Distributed Subtrajectory Join on Massive Datasets

TL;DR

This paper presents three distributed algorithms for subtrajectory join processing on massive datasets using MapReduce, significantly improving performance over existing methods in large-scale trajectory data analysis.

Contribution

It introduces three novel distributed solutions, including indexing, for efficient subtrajectory join processing on big datasets, addressing a key challenge in mobility data analytics.

Findings

DTJi is up to 16x faster than DTJb

DTJi is 10x faster than DTJr

DTJi outperforms the state-of-the-art algorithm

Abstract

Joining trajectory datasets is a significant operation in mobility data analytics and the cornerstone of various methods that aim to extract knowledge out of them. In the era of Big Data, the production of mobility data has become massive and, consequently, performing such an operation in a centralized way is not feasible. In this paper, we address the problem of Distributed Subtrajectory Join processing by utilizing the MapReduce programming model. Compared to traditional trajectory join queries, this problem is even more challenging since the goal is to retrieve all the "maximal" portions of trajectories that are "similar". We propose three solutions: (i) a well-designed basic solution, coined DTJb, (ii) a solution that uses a preprocessing step that repartitions the data, labeled DTJr, and (iii) a solution that, additionally, employs an indexing scheme, named DTJi. In our experimental study, we utilize a 56GB dataset of real trajectories from the maritime domain, which, to the best of our knowledge, is the largest real dataset used for experimentation in the literature of trajectory data management. The results show that DTJi performs up to 16x faster compared with DTJb, 10x faster than DTJr and 3x faster than the closest related state of the art algorithm.

Figures

Tables

Table 1. Table 1: Parameters and default values (in bold)

Parameter	Values
${ϵ}_t$ (in minutes)	10, 15, 20, 25, 30
${ϵ}_{sp}$ (%)	10%, 20%, 30%, 40%, 50%
$\delta t$ (in minutes)	10, 15, 20, 25, 30
max # of points per cell (%)	1%, 2%, 3%, 4%, 5%
# of Nodes	12, 24, 36, 48

Full text

K****eywords (Sub)Trajectory Join, Distributed Join Processing, MapReduce

1 Introduction

During the recent years, the proliferation of GPS enabled devices has led to the production of enormous amounts of mobility data. This “explosion” of data generation has posed new challenges in the world of mobility data management. One of these challenges is the so-called trajectory join problem, which aims to find all pairs of “similar” (i.e. nearby in space-time) trajectories in a dataset [5, 7, 9, 30]. An even more interesting and challenging problem is the subtrajectory join query [3], where, for each pair of trajectories, we want to retrieve all the “portions” of trajectories that are “similar”. However, the subtrajectory join is a processing-intensive operation. Centralized algorithms do not scale with the size of today’s trajectory data, thus parallel and distributed algorithms are necessary in order to provide efficient processing of subtrajectory join, an issue largely overlooked in the related research.

Several modern applications that manage trajectory data could benefit from such an operation. For instance, in the urban traffic domain, carpooling is becoming increasingly popular. More concretely, consider a mobile application which tries to match users that can share a ride based on their past movements. Here, given a set of trajectories we want to find all the pairs of users that can share a ride for a portion of their everyday routes without significantly deviating (spatially and temporally) from their daily routine (i.e. retrieve all pairs of maximal subtrajectories that move close in space and time). Another interesting scenario concerns the identification of suspicious movement by a governmental security agency. For instance, given a set of trajectories that depict the movement of suspicious individuals, we would like to retrieve all the pairs of moving objects that move “close” to each other for more than a threshold (moving together for small periods of time could be considered as coincidental) as candidates for illegal activity. Moreover, such a query is in fact the building block for a number of operations than aim to identify mobility patterns, such as co-movement patterns (e.g. flocks [15], convoys [17], swarms [18]). An even more challenging problem is that of subtrajectory clustering  [22, 1]. An interesting application scenario of subtrajectory clustering is network discovery, where given a set of trajectories (e.g from the maritime or the aviation domain) we want to identify the underlying network of movement by grouping subtrajectories that move “close” to each other and use cluster representatives/medoids as network edges. One of the main goals of subtrajectory clustering is to segment trajectories to subtrajectories. Finally, trajectory segmentation techniques [20, 22], can directly benefit from the subtrajectory join query since their input, for each trajectory, is the number of objects that were located close to it at any given time. However, the bottleneck in all these applications is the underlying processing cost of the join operation, which calls for parallel and distributed solutions that scale beyond the limitations of a single machine.

Inspired by the above application scenarios, the problem that we address in this paper is as follows: given two sets of trajectories (or a single set and its mirror in the case of self-join), identify all pairs of maximal “portions” of trajectories (or else, subtrajectories) that move close in time and space w.r.t a spatial threshold $\epsilon_{sp}$ and a temporal tolerance $\epsilon_{t}$, for at least some time duration $\delta t$. To illustrate this informal definition, as depicted in Figure 1(a), given two trajectories $r$ and $s$, the pair of their maximal matching “portions” is ($\{r_{4},r_{5},r_{6},r_{7},r_{8}\},\{s_{3},s_{4},s_{5},s_{6},s_{7}\}$). Each point of a trajectory defines a spatio-temporal ’neighborhood’ area around it, a cylinder of radius $\epsilon_{sp}$ and height $\epsilon_{t}$. In order for a pair of subtrajectories to be considered “matching”, each point of a subtrajectory must have at least one point of the other subtrajectory in its “neighborhood”, thus making the result symmetrical. A pair of matching subtrajectories is maximal if there exists no superset of either subtrajectories that can replace them and the pair still qualifies as a “matching” pair.

There have been some efforts to tackle variations of this problem in a centralized way [3, 5, 7]. However, these solutions discover pairs of entire trajectories and cannot identify matching sub-trajectories. In [4], all pairs of “matching” (w.r.t a spatial threshold) subtrajectories of exactly $\delta t$ duration are retrieved in contrast with the problem addressed in this paper, where the goal is to identify maximally “matching” subtrajectories, which is vital for exploiting the output in subsequent steps, e.g. the mining operations mentioned above. Moreover, applying these centralized solutions to a parallel and distributed environment is not straightforward and is often impossible if radical changes to the methods/algorithms do not take place, since there are several non-trivial issues that arise. For instance, how to partition the data in such a way so that each partition can be processed independently and be of even size.

In a recent effort, in  [25] the authors try to tackle the problem of trajectory similarity join in spatial networks in parallel. The solution proposed in  [25] handles each trajectory separately and all the data have to be replicated for each trajectory and, consequently, to each node. Due to this fact, such a solution cannot scale to terabytes of data, thus making it inapplicable to Big Data. Furthermore, such an approach assumes that the underlying network is known in advance, hence it cannot support datasets of moving objects that move freely in space (e.g. from the maritime or the aviation domain). As a result, a scenario where the goal is to identify the underlying network cannot be supported. Finally, the output of  [25] is pairs of trajectories and not subtrajectories, which is significantly different than the problem addressed in this paper. More recently, in  [26] the authors try to tackle the problem of trajectory similarity join. Specifically, given two sets of trajectories, a similarity function (e.g. DTW) and a similarity threshold, they aim to identify all pairs of trajectories that exceed this similarity threshold. Again, the problem addressed in [26] is to retrieve pairs of trajectories in contrast with the problem that we try to tackle in this paper, which is to retrieve all pairs of “maximally matching” subtrajectories.

It is straightforward to claim that an integral part of any algorithm that tries to address the subtrajectory join query is to identify all pairs of points that move “close enough” in time and space w.r.t a spatial threshold $\epsilon_{sp}$ and a temporal tolerance $\epsilon_{t}$, e.g. $r_{4}$ and $s_{3}$ in Figure 1(a). In that sense, another line of research that is closely related to our problem is that of MapReduce-based spatial [33, 2, 11] and multidimensional joins [28, 19, 13], where the goal is to identify such points. A generic solution which could form the basis for any MapReduce-based spatial (or spatiotemporal) join algorithm is presented in [33], where the input data are partitioned into small, disjoint tiles at Map stage and get joined at the Reduce stage by performing a plane sweep algorithm along with a duplication avoidance technique. However, all of the above approaches try to solve a problem that is significantly different from ours since our problem is not to join spatial or multidimensional objects but identify all pairs of “maximally matching” subtrajectories.

In this paper, we provide efficient solutions for the Distributed Subtrajectory Join processing problem, as it is formally defined in Section 2. To the best of our knowledge, this problem has not been addressed in the literature yet. Our main contributions are the following:

• •

We formally define the problem of Distributed Subtrajectory Join processing, investigate its main properties, and discuss its main challenges.

• •

We present a well-designed algorithm, called $\mathit{DTJb}$, that solves the problem of Distributed Subtrajectory Join processing by employing two MapReduce phases.

• •

We propose an improvement of $\mathit{DTJb}$, termed $\mathit{DTJr}$, which is equipped with a repartitioning mechanism that achieves load balancing and collocation of temporally adjacent data.

• •

To boost the performance of query processing even further, we introduce $\mathit{DTJi}$, which extends $\mathit{DTJr}$ by exploiting an indexing scheme that speeds up the computation of the join.

• •

We compare with an appropriately modified state of the art MapReduce spatial join algorithm and show that our solution performs several times better.

• •

We study the performance of the proposed algorithms by using, to the best of our knowledge, the largest real trajectory dataset (56GB) used before in the relevant literature, thus demonstrating the scalability of our algorithms.

The rest of the paper is organized as follows. In Section 2 we introduce the problem. In Section 3, we present $\mathit{DTJb}$. Subsequently, in Section 4 we propose $\mathit{DTJr}$ that utilizes a preprocessing step. In Section 5, we introduce $\mathit{DTJi}$ that boosts the performance of the join processing. In Section 6, we provide our experimental study. Finally, we present an overview of the related work in Section 7, and we conclude the paper in Section 8.

2 Problem Statement

Given a set $R$ of trajectories moving in the xy-plane, a trajectory $r\in R$ is a sequence of timestamped locations $\{r_{1},\dots,r_{N}\}$. Each $r_{i}=(x_{i},y_{i},t_{i})$ represents the $i$-th sampled point, $i\in{1,\dots,N}$ of trajectory $r$, where $N$ denotes the length of $r$ (i.e. the number of points it consists of). The pair $(x_{i},y_{i})$ and $t_{i}$ denote the 2D location in the xy-plane and the time coordinate of point $r_{i}$ respectively. A subtrajectory $r_{i,j}$ is a subsequence $\{r_{i},\dots,r_{j}\}$ of $r$ which represents the movement of the object between $t_{i}$ and $t_{j}$ where $i<j$.

Given a pair $(r,s)$ of trajectories (the same holds for subtrajectories) with $r\in R$ and $s\in S$, the common lifespan $w_{r,s}$ is defined as the time interval $[max(r_{1}.t,s_{1}.t),min(r_{N}.t,s_{M}.t)]$, where $r_{1}$ ($s_{1}$) is the first sample of $r$ ($s$, respectively) and $r_{N}$ ($s_{M}$) is the last sample of $r$ ($s$, respectively). The duration of the common lifespan $w_{r,s}$ is $\Delta w_{r,s}$ = $min(r_{N}.t,s_{M}.t)$ - $max(r_{1}.t,s_{1}.t)$

Further, let $DistS(r_{i},s_{j})$ denote the spatial distance between two points $r_{i}$, $s_{j}$, which is defined as the Euclidean distance in this paper, even though other distance functions are also applicable. Also, let $DistT(r_{i},s_{j})$ denote the temporal distance, defined as $|r_{i}.t-s_{j}.t|$.

Definition 1**.**

(Matching subtrajectories) Given a spatial threshold $\epsilon_{sp}$, a temporal tolerance $\epsilon_{t}$ and a time duration $\delta t$, a “match” between a pair of subtrajectories $(r^{\prime},s^{\prime})$ occurs iff $\Delta w_{r^{\prime},s^{\prime}}\geq\delta t-2\epsilon_{t}$, and $\forall r^{\prime}_{i}\in r^{\prime}$ there exists at least one $s^{\prime}_{j}\in s^{\prime}$ so that $DistS(r^{\prime}_{i},s^{\prime}_{j})\leq\epsilon_{sp}$ and $DistT(r^{\prime}_{i},s^{\prime}_{j})\leq\epsilon_{t}$, and $\forall s^{\prime}_{j}$ there exists at least one $r^{\prime}_{i}$ so that $DistS(s^{\prime}_{j},r^{\prime}_{i})\leq\epsilon_{sp}$ and $DistT(s^{\prime}_{j},r^{\prime}_{i})\leq\epsilon_{t}$.

Definition 2**.**

(Maximally matching subtrajectories) Given a pair of “matching” subtrajectories $(r^{\prime},s^{\prime})$ which belong to trajectories $r,s$ respectively, this pair is considered a “maximal match” iff $\nexists$ superset $r^{\prime\prime}$ of $r^{\prime}$ or $s^{\prime\prime}$ of $s^{\prime}$ where the pair $(r^{\prime\prime},s^{\prime})$ or $(r^{\prime},s^{\prime\prime})$ or $(r^{\prime\prime},s^{\prime\prime})$ would be “matching”.

At this point, we should clarify that two trajectories may have more than one “maximal matches” (i.e. pairs of subtrajectories). Having provided the above background definitions, we can define the subtrajectory join query between two sets of trajectories.

Definition 3**.**

(Subtrajectory join) Given two sets of trajectories $R$ and $S$, a spatial threshold $\epsilon_{sp}$, a temporal tolerance $\epsilon_{t}$ and a time duration $\delta t$, the subtrajectory join query searches for all pairs $(r^{\prime},s^{\prime})$, $r^{\prime}\in r\in R$ and $s^{\prime}\in s\in S$, which are “maximally matching” subtrajectories.

2.1 A Closer Look at the Subtrajectory Join Problem

An integral part of any algorithm addressing the subtrajectory join query, as defined in Defn. 3 above, is to identify all pairs of joining points ($r_{i}$,$s_{j}$), where $r_{i}\in r$ and $s_{j}\in s$, which satisfy the following property: $DistS(r_{i},s_{j})\leq\epsilon_{sp}$ and $DistT(r_{i},s_{j})\leq\epsilon_{t}$. In fact, the set of joining points is the outcome of the inner join $R\bowtie S$, where the evaluated join predicates are the ones mentioned above. However, as it will be explained next, these pairs of points do not suffice to return the correct query result.

Let $\mathcal{A}$ denote the class of correct algorithms for the subtrajectory join problem. A naive algorithm $A\in\mathcal{A}$ would require the Cartesian product $R\times S$ to produce the correct result. We claim that $R\times S$ can be represented by three sets of points, the set of joining points ($JP$), the breaking points ($BP$) and the non-joining points ($NJP$). Formally, $R\times S=JP\cup BP\cup NJP$. The definitions of these sets follow, and the discussion is aided by Figure 1(b), which is a variation of Figure 1(a) in order to emphasize the distinction between $JPs$, $BPs$ and $NJPs$.

The first set of points that needs to be identified, besides $JP$, contains all points $r_{i}\in r\forall r\in R$ that do not “match” with any other point in $S$. We call such points as breaking points.

Definition 4**.**

(Breaking points) A point $r_{i}\in r\in R$ is a breaking point iff it is not a joining point with any other point $s_{j}\in S$:

$\nexists s_{j}\in S$: $DistS(r_{i},s_{j})\leq\epsilon_{sp}$ $\land$ $DistT(r_{i},s_{j})\leq\epsilon_{t}$.

As it will be shown later, the lack of information about $BPs$ can make an algorithm $A\in\mathcal{A}$ to falsely identify a pair of subtrajectories as “matching”. The set of $BP$ along with the set of $JP$ is actually the outcome of the full outer join of $R$ and $S$. Figure 1(b) depicts the case where $r_{1}$ is a breaking point of $r$ ($r_{2}$, $r_{3}$ and $r_{10}$ are also breaking points), since it does not “match” with any other point of any trajectory. Obviously, breaking points are never reported as part of the answer set and the portion of $r$ that could possibly contribute to the result is subtrajectory $r_{4,9}$.

The last set of points that is necessary to be identified consists of the pairs of points that do not “match”, coined non-joining points, since some of them might indicate the start or the end of “maximally matching” subtrajectories.

Definition 5**.**

(Non-joining points) A point $r_{i}\in r\in R$ is a non-joining point w.r.t. $s_{j}\in s\in S$ iff: (a) $r_{i}$ and $s_{j}$ are not breaking points, and (b) $r_{i}$ is not a joining point with $s_{j}\in S$:

$r_{i},s_{j}\notin BP$ $\land$ ($DistS(r_{i},s_{j})>\epsilon_{sp}$ $\lor$ $DistT(r_{i},s_{j})>\epsilon_{t}$).*

This case is illustrated in Figure 1(b), where $r_{5}$ is a non-joining point w.r.t. $s_{5}$.

Actually, if we remove condition (a) from Definition 5 then it is obvious that a breaking point $r_{i}$ is a special case of non-joining point where $r_{i}\in R$ is a non-joining point with every other point $\in S$. However, we differentiate breaking points from non-joining points so as to reduce the amount of information that needs to be kept, i.e. instead of keeping multiple non-joining points we only keep one breaking point. In the section that follows, we investigate the theoretical properties of an efficient algorithm in class $\mathcal{A}$.

2.2 Properties of Subtrajectory Join

In this section, we provide the theoretical properties for designing efficient algorithms for the subtrajectory join problem. The properties shown below essentially determine which pairs of points from the sets $BP$ and $NJP$ are necessary for a correct algorithm in class $\mathcal{A}$.

Lemma 1**.**

The set of breaking points is necessary in order to produce the correct result set for the Subtrajectory Join problem.

This result indicates that breaking points cannot be ignored by an algorithm, without compromising the correctness of the result. The remaining question is whether all non-joining points are also necessary. In the following, we define a subset of non-joining points points $sNJP\subseteq NJP$, and show that this subset is actually necessary.

Definition 6**.**

(Subset $sNJP$ of non-joining points) A non-joining point $s_{j}\in S$ w.r.t. a point $r_{i}\in R$ belongs to $sNJP$, if at least one of its adjacent points $s_{j-1}$ or $s_{j+1}$ is a joining point with any point $r_{p}\in r$, with $p\neq i$ and $\nexists$ a point $r_{q}$, with $q\neq i$, such that $DistT(r_{q},s_{j})\leq DistT(r_{i},s_{j})$.

Returning to the example of Figure 1(b), $s_{5}$ does not “match” with any point in $r$, even though all points of $r_{4,8}$ “match” with a point in $s_{3,7}$. Again, failure to identify pairs of points such as $(s_{5},r_{6})$ would result in erroneously identifying larger “matching” subtrajectories.

Lemma 2**.**

The set $sNJP$ of pairs of non-joining points is necessary in order to produce the correct result set for the Subtrajectory Join problem.

In summary, our main finding is that a typical join algorithm that identifies only the set of $JP$ is not enough in order to address the subtrajectory join problem. Additionally to the set of $JP$, an algorithm needs to identify both the set of $BP$ and the subset $sNJP$ during the join processing, in order to ensure correctness.

2.3 Distributed Subtrajectory Join

Given two sets $R$ and $S$ of trajectories, the typical approach for parallel join processing consists of two main phases: (a) data repartitioning, in order to create pairs of partitions $R_{i}\subset R$ and $S_{j}\subset S$, such that part of the join can be processed using only $R_{i}$ and $S_{j}$, and (b) join processing, where a join algorithm is performed on partitions $R_{i}$ and $S_{j}$.

Problem 1**.**

(Distributed Subtrajectory Join) Given two distributed sets of trajectories, $R=\cup R_{i}$ and $S=\cup S_{j}$, compute the subtrajectory join (Defn. 3) in a parallel manner.

In this setting, the main challenges are the following: (a) ensure that the created partitions are sufficient to produce parts of the total join without additional data, (b) generate even-sized partitions in order to balance the load fairly to multiple nodes, (c) handle the problem of potential duplicate existence in the join results, which may arise due to the way partitions are created, and (d) process the actual join on the partitions in an efficient way. The first challenge sets the foundations for parallel processing, as it identifies pairs of partitions that can be processed together, without any additional data, and produce a subset of the final join result. The second challenge is about load balancing and determines the efficiency of parallel processing, which is not straightforward, since processing uneven work units in parallel may lead to sub-optimal performance (as the slowest task will determine the query execution time). The third challenge, labeled duplicate avoidance, is about avoiding to generate duplicate results which typically occurs in parallel join processing. Finally, the fourth challenge, labeled efficient join, refers to the efficiency of the (centralized) algorithm used to join two partitions.

Clearly, solving the above problem is quite challenging in a distributed setting, as multiple challenges need to be addressed at the same time. In the following sections, we present a well designed solution solution to the Distributed Subtrajectory Join problem along with two improved versions, following the popular MapReduce paradigm.

3 The Basic Subtrajectory Join Algorithm

3.1 Preliminaries

One of the prevalent technologies for dealing with Big Data and offline analytics, is the MapReduce programming paradigm [8] and its open-source implementation Hadoop [27]. A lot of efforts have been made as far as it concerns join processing through this technology and a survey on limitations of MapReduce/Hadoop, also related to join processing, is conducted in [10]. In more detail, Hadoop is a distributed system created in order to process large volumes of data which are usually stored in the Hadoop Distributed File System (HDFS). In more detail, when running a MapReduce (MR) job, each Mapper processes (in parallel) an input split, which is a logical representation of data. An input split typically consists of a block of data (the default block size is 128MB) but it can be adjusted according to the users’ needs by implementing a custom FileInputFormat along with the corresponding FileSplitter and RecordReader. Subsequently, for each record of the split the “map” function is applied. The output of the Map phase is sorted and grouped by the “key” and written to the local disk. Successively, the data is partitioned to Reducers based on a partitioning strategy (also known as shuffling), and each Reducer receives a partition (group) of data and applies the “reduce” function to the specific group. Finally, the output of the Reduce phase is written to HDFS.

3.2 The DTJb Algorithm

Our first algorithm, named $\mathit{DTJb}$, consists of three phases: (a) the Partitioning phase, where input data is read and partitioned, (b) the Join phase, where the sets $JP$, $BP$ and $sNJP$ are identified in each partition, and (c) the Refine phase, where these sets are grouped by trajectory and sorted by time in order to identify all the pairs of “maximally matching” subtrajectories111For the sake of simplicity, from now on, we are going to consider the case of self-join. The transition to the problem of joining two relations is straightforward..

3.2.1 Partitioning Phase

The first challenge is how to partition the input data in order to satisfy the requirement for parallel processing. Partitioning the data into $N$ disjoint temporal partitions $R=\cup_{i=1}^{N}part_{i}$, where $R$ is the set of trajectories, cannot guarantee the correctness during parallel processing, due to the temporal tolerance parameter $\epsilon_{t}$. Hence, we define a partitioning where each $part_{i}$ is expanded by $\epsilon_{t}$, thus expanded partitions can be processed independently in parallel. Let $expPart_{i}$ denote such an expanded partition. Processing each $expPart_{i}$ individually guarantees correctness, but at the cost of having duplicates due to the point replication in temporally overlapping partitions. To address this duplication avoidance challenge, we supplement each point with a flag partFlag that indicates whether this point belongs to the original partition (i.e. not expanded by $\epsilon_{t}$) or not.

Lemma 3**.**

An expanded partition $expPart_{i}$ is sufficient in order to produce the sets of $JP$ and $BP$ for $part_{i}$

Unfortunately, an expanded partition $expPart_{i}$ is not sufficient in order to produce the set of $sNJP$ since, according to Defn. 6, for each pair $(r_{j},s_{k})$ that belongs to $NJP$ we need to examine $r_{j-1}$ and $r_{j+1}$, which may span to other partitions. However, the set of $sNJP$ can be identified at the Refine phase, where all the pairs concerning a trajectory are grouped together.

In more detail, we choose to partition the data into uniform temporal partitions, where for each pair of partitions $(part_{i},part_{j})$, with $i\neq j$ and $i,j\in[1,N]$, it holds that $DistT(t_{e}^{part_{i}},t_{s}^{part_{i}})=DistT(t_{e}^{part_{j}},t_{s}^{part_{j}})$. Typically, the duration of a partition is larger than the maximum interval between two consecutive points of any trajectory. As illustrated in Figure 2, in the Map phase we access each data point and assign it to the expanded partition with which it intersects, essentially applying a temporal range partitioning. Then, the data is grouped by expanded partition, sorted by time and fed to the Reduce phase, where the Join procedure takes place.

3.2.2 Join Phase

Figure 2 shows that each Reducer task takes as input an expanded partition and performs the Join operation. At this point, the duplication avoidance technique is applied, by employing the aforementioned flag and emitting only pairs where at least one point belongs to the original partition. The input of this phase is a set of tuples of the form $\left\langle t,x,y,trajID,partFlag\right\rangle$. The output of this MR job is a set of (a) $JP$, (b) $BP$ and (c) candidate $sNJP$.

In more detail, we apply a plane sweep technique in order to perform the Join, by sweeping the temporal dimension. We choose to employ such a technique due to the fact that is much more efficient than a nested-loop join approach, since our data already arrive sorted by the temporal dimension, as illustrated in Figure 2. A typical plane sweep algorithm would emit only the set of $JP$, which is not enough in our case. For this reason, we devised and implemented a modified plane sweep technique, named TRJPlaneSweep, depicted in Figure 3, which also reports the sets of $BP$ and candidate $sNJP$.

Algorithm 1 presents how the Join processing is performed. Each accessed point is inserted to an array $D$, which contains points sorted in increasing time. After point insertion, (Algorithm 2 is invoked for the currently accessed point (say $D[i]$) if $D[i]$ belongs to the original partition. TRJPlaneSweep examines the previously accessed points for the previous $\epsilon_{t}$ window (line 4). The role of this function is threefold. First, it identifies joining points with $D[i]$, e.g., point $D[j]$, and emits them in the form $((D[i],D[j]),True)$ (lines 5-10). Depending on the outcome of the duplicate avoidance technique, pairs $((D[j],D[i]),True)$ are also output. Second, it discovers points that belong to the candidate $sNJP$ set by examining whether the previous trajectory point (getPrevTrPoint)) of $D[j]$ (and $D[i]$), say $D[k]$, is a $NJP$ (FindMatch) with each point $\in$ $D[i].trajID$ ($D[j].trajID$, respectively) (lines 11-17). In case such points are identified, they are output with a different flag $((D[i],D[k]),False)$ to differentiate them from $JP$. Third, it discovers the points that belong to $BP$. In more detail, in lines 18–19, a breaking point $D[i]$ is added to the breaking points set $BP$ and in lines 7 and 10 is removed if a point has a “match”. The remaining points in $BP$ are reported as breaking points, using the following form: $((D[i],null),True)$ (Algorithm 1 lines 7–8).

By examining only the previous point of a $JP$ in a trajectory, we might not examine a possible temporary adjacent point that might lie after the last $JP$ of a trajectory in each partition. For this reason, we post-process the last $JP$s in order to check for candidate $sNJP$s by invoking the TreatLastTrPoints function (Algorithm 1 line 6).

Example 1**.**

As illustrated in Figure 3(a), we suppose that the current point inserted into $D$ is $q_{2}$. In Figure 3(b), assuming that $DistS(q_{2},r_{2})\leq\epsilon_{sp}$, we get a “match” and pair $((q_{2},r_{2}),True)$ is reported (the symmetric pairs are omitted for simplicity). Subsequently, we need to find the previous point of $r$ and in order to achieve this we should traverse our data backwards until we find it, as presented in Figure 3(c). When we find $r_{1}$, we need to check whether it is a $NJP$ for each point $\in q$, as illustrated in Figure 3(c). If there exists a point $\in q$ that “matches”, in our case $q_{1}$, nothing is reported and we proceed to examine whether $q_{2}$ and $p_{2}$ are $JP$s. If $DistS(q_{2},p_{2})\leq\epsilon_{sp}$ then we output the pair $((q_{2},p_{2}),True)$, as shown in Figure 3(d). Subsequently, we need to find $p_{1}$ and check whether it is $NJP$ for each point $\in q$. As depicted in Figure 3(e) there is no “match” between $p_{1}$ and any of the points of $q$. For this reason, we report the pair $((q_{2},p_{1}),False)$. The same procedure is continued to the next point inserted to memory as delineated in Figure 3(f) until there are no more points inserted.

The complexity of the Join procedure is $O(|D|\cdot a\cdot((1-b)\cdot|D|+b\cdot|D|\cdot(2\cdot(L+2\cdot a\cdot|D|))))$, where $|D|$ is the number of points, $a$ is the selectivity of $\epsilon_{t}$ and $b$ is the selectivity of $\epsilon_{sp}$. $L$ is the number of points that have to be traversed in order to find the previous point of a specific trajectory. It is obvious that when $a$ tends to reach 1 the complexity tends to reach $O(|D|^{2})$. In the worst case, the complexity can be analogous to $O(|D|^{3})$, when both $a$ and $b$ tend to reach 1. However, for a typical analysis task $\epsilon_{t}$ and $\epsilon_{sp}$ are much smaller than the dataset duration and the dataset diameter respectively. Roughly, we can say that the complexity is $O(a\cdot b\cdot|D|^{2})$.

3.2.3 Refine Phase

The output of the Join phase is actually pairs of points. From now on, let us refer to the left point of such a pair as reference point and the trajectory that it belongs to, reference trajectory. The Refine phase consists of a second MR job that reads the output of the Join step and groups points by the reference trajectory. Each Reduce task receives all pairs of points belonging to a specific trajectory, sorted first by the reference point’s time and the by the non-reference trajectory ID. Figure 4 shows an example where the output pairs of points from the Join step are grouped, sorted and fed as input to three Reduce tasks (for trajectories $p$, $q$, and $r$ respectively). The general idea here is to scan the set of $JP$ in a sliding window fashion so as to identify “maximally” matching subtrajectories and at the same time “consult” the sets of $BP$ and $sNJP$ in order to avoid false identifications, as described in Section 2.2.

Hence, each Reducer accesses all the pairs of a reference trajectory (say $p$) sorted by time, i.e., $\{p_{1},p_{2},\dots,p_{n}\}$. Algorithm 3 describes the pseudo-code of the Refine phase which aims to identify all the “maximally matching” pairs of subtrajectories of $p$ with other subtrajectories of any trajectory $x$ ($x\neq p$). For each accessed pair (($p_{i},x_{j}$),flag), the algorithm assigns it in one of the two structures that it maintains: the MatchList and the FalseList. All $JP$ and $BP$ will be kept in the MatchList, whereas the candidate $sNJP$ is kept in the FalseList (lines 10–13). Again, this is more clearly depicted in the example of Figure 4. Also, notice that for each reference point in the MatchList, we maintain a list of points sorted by trajectory ID.

Lemma 4**.**

The set of candidate $sNJP$ is sufficient so as to identify the set of $sNJP$ at the Refine phase.

The algorithm proceeds as follows: as soon as all pairs of points of a specific reference point $p_{i}$ have been accessed, it initiates processing on the MatchList. The processing takes place only if the first and last point of $p$ in MatchList have temporal distance greater than or equal to $\delta t$ (line 5). The processing essentially identifies points of other trajectories that join with points of $p$ in the whole temporal window. This is performed by intersecting the lists in MatchList and excluding points existing in the FalseList (line 6). List intersection is efficiently performed in linear time to the length of the lists, since the lists are sorted by trajectory ID. Figure 5 depicts the result of this processing as resultT.

Subsequently, the points in resultT are processed as follows. We start from the first point and take into consideration all points with temporal distance at most $\delta t-2\epsilon_{t}$ from the first point. From this set of points, we derive the subtrajectories that “match” for the entire $\delta t-2\epsilon_{t}$ window, and insert them in resultF (line 7). The temporary results of the resultF structure are added to the final result structure resultFinal, if not already contained in it (line 8). Then, a new set of points is considered, of temporal distance at most $\delta t-2\epsilon_{t}$ from the second point of resultT and the process is repeated, similarly to a sliding a window of duration $\delta t-2\epsilon_{t}$ on resultT. In the end, the first entry of the MatchList ($p_{1},\{q_{1},r_{1},s_{1}\}$) is removed (line 9), as all potential results containing $p_{1}$ have already been produced. The algorithm terminates when the entire trajectory is traversed, the resultFinal is returned and each element of this list is emitted.

Example 2**.**

Figure 5 presents a working example of the Refine algorithm given the specific MatchList and FalseList of trajectory $p$.Assuming that $DistT(p_{1}.t,p_{7}.t)\geq\delta t$, we intersect all the lists contained in the specific window of the MatchList and we pass the result to resultT. In this way, the list of the last entry of resultT will contain only the points that belong to the subtrajectories that move “close” enough with $p$ for the whole $\delta t$ window. During list intersection, we take into account the FalseList structure in order to deal with points that belong to $sNJP$. Specifically, even though for each $p_{i}$, with $i\in$ [1,7] $\exists$ a “match” with $q$, however $q_{6}$ has no “match” with $p$, as depicted in the FalseList. For this reason, $q$ should be excluded from resultT after $p_{5}$. Then, a sliding $\delta t-2\epsilon_{t}$ window is created that traverses resultT, and for each such window we intersect all lists and the result is stored in resultF. For the first $\delta t-2\epsilon_{t}$ window, as depicted in Figure 5, subtrajectories $r_{1,5}$ and $q_{1,5}$ are identified. The reason for this is to discover the subtrajectories that move “close” enough, with $p$ for the whole $\delta t-2\epsilon_{t}$ window. Subsequently, before proceeding to the next $\delta t$ window, the contents of resultF are inserted to the final result, if not already contained.

The complexity of the Refine procedure is $O(T\cdot SW\cdot dt\cdot l)$, where $T$ is the average number of points in a trajectory, $SW$ is the number of points contained in the $\delta t$ window, $dt$ is the number of points contained in the $\delta t-2\epsilon_{t}$ window and $l$ is the size of the list. The complexity, here, clearly depends on the average number of points per trajectory, the $\epsilon_{t}$ and $\delta t$ parameter, and the number of pairs emitted by the join phase which in turn depends on $\epsilon_{t}$ and $\epsilon_{sp}$.

4 Subtrajectory Join with Repartitioning

Even though the $\mathit{DTJb}$ algorithm provides a correct solution to the Distributed Subtrajectory Join problem, it has some limitations. In particular, it does not address the load balancing challenge, since it does not handle the case of temporally skewed data. Also, due to the two chained MR jobs, the intermediate output of the first job is written to HDFS and must be read again by the second job, which imposes a significant overhead as its size is comparable and can be even bigger than the original dataset.

Motivated by these limitations, we propose an improved two-step algorithm ($\mathit{DTJr}$), which consists of the repartitioning and the query step. Each step is implemented as a MR job. However, the repartitioning step is considered a preprocessing step, since it is performed once and is independent of the actual parameters of our problem, namely $\epsilon_{sp}$, $\epsilon_{t}$, and $\delta t$.

4.1 Repartitioning

The aim of the repartitioning step is to split the input dataset in $M$ equi-sized, temporally-sorted partitions (files), which are going to be used as input for the join algorithm. This is essential for two reasons: (a) it will provide the basis for load balancing, by addressing the issue of temporal skewness in the input data, and (b) it will result in temporal collocation of data, thus drastically reducing processing and network communication costs.

The repartitioning step is performed by means of a MR job as follows. We sample the input data, using Hadoop’s InputSampler, and construct an equi-depth histogram on the temporal dimension. The histogram contains $M$ equi-sized bins, i.e. the numbers of points in any two bins are equal, where the borders of each bin correspond to a temporal interval $[t_{i},t_{j})$.

The equi-depth histogram is exploited by the Map phase in order to assign each incoming data object in the corresponding histogram bin, based on the value of its temporal dimension. Each “map” function outputs each data object using as key a value $[1,M]$ that corresponds to the bin that the object belongs to. During shuffling, all data objects that belong to a specific bin are going to be sorted in time and will be collected by a single “reduce” function (thus having $M$ “reduce” functions). As a result, each “reduce” function writes an output file to HDFS that contains all data objects in a specific temporal interval $[t_{i},t_{j})$ sorted by increasing time. A graphical view of the MR job is provided in Figure 6(a).

A subtle issue is how to determine the number $M$ of bins (and, consequently, output files). A small value of $M$, smaller than the number of nodes in the cluster, would be opposed to the collocation property because data would have to be transferred through the network. On the other hand, a large value of $M$ would result to many small files, smaller than the HDFS block size, and would lead to inefficient use of resources as well as increasing the management cost of these HDFS files. A good compromise is to have files of equal size to the HDFS block. Hence, the number of files can be calculated as $M=\lceil\frac{InputTotalSize}{hdfsblocksize}\rceil$. Collocation can be further improved by placing temporally adjacent files to the same nodes. This can be achieved by grouping together $k$ consecutive files, where $k=\lceil\frac{M}{N}\rceil$, with $N$ being the number of nodes, and using their group id as the partition key.

4.2 The DTJr Algorithm

In order to minimize the I/O cost, the MR job that implements the proposed algorithm performs the Join procedure in the Map phase, and the Refine in the Reduce phase. To achieve this, we need to provide to a Map task as input, a data partition that contains all necessary data in order to perform part of the Join procedure independently from other Map tasks. Thus, an HDFS block produced by the repartitioning phase is expanded with additional points that exist at time (+/-)$\epsilon_{t}$, and this is the process of InputSplits creation. In this way, points are duplicated to other HDFS blocks, which means that the same point may be output by two different Map tasks. To avoid this pitfall, a different duplicate avoidance mechanism is introduced which practically determines that a point is going to be output only by a single Map task; the Map task processing the HDFS block where the point belongs to.

As already mentioned, each data partition (InputSplit) that is fed to a Map task should contain all the data needed to perform the join of points for the specific partition, i.e. data for the period $[t_{s}^{part}-\epsilon_{t},t_{s}^{part}+\epsilon_{t}]$. However, an output file produced by the repartitioning step is not sufficient due to the temporal tolerance $\epsilon_{t}$, thus we need to augment these output files with extra data points, so that they form independent data partitions. At technical level, we devised and implemented a new FileInputFormat called BloatFileInputFormat, along with the corresponding FileSplitter and RecordReader, which selectively combines different files in order to create splits that carry all the necessary data points. Furthermore, during the creation of input splits we augment (as metadata) each split with the starting and ending time of the original partition of each split, termed $t_{s}^{base}$ and $t_{e}^{base}$. The utility is to provide us with a simple way to perform duplicate avoidance at the Join phase.

Figure 6(b) shows that each Map task takes as input a split and performs the join at the level of point for a specific data partition. The input of this phase is a set of tuples of the form $\left\langle t,x,y,trajID\right\rangle$ sorted in ascending time $t$ order. Since the data are already sorted w.r.t. the temporal dimension, we can apply the Join procedure, presented in Section 3.2.2. The output of the Map phase will be the $JP$, $BP$ and $sNJP$ sets. Finally, the Refine procedure presented in Section 3.2.3 can be performed at the Reduce phase.

5 Index-based Subtrajectory Join with Repartitioning

The Join step of the previous algorithms is common and operates on the array $D$ that contains temporally sorted points. However, it can be improved in two ways. First, by employing spatial filtering in order to avoid attempting to join points that are far away. Second, by having an index structure that given a point $p_{i}$ can efficiently locate the (temporally) previous point $p_{i-1}$ of $p$. Motivated by these observations, we devised and implemented an indexing scheme in order to speed up the processing of the join.

5.1 Indexing Scheme

As illustrated in Figure 7, this scheme consists of 3 levels. We already covered the first level in Section 4.1, where the initial data are partitioned to equi-sized temporal partitions (Section 4). At the second level, we partition the space. In order to have load balanced partitions we utilize the spatial partitioning provided by QuadTrees. More specifically, an “empty” QuadTree is created once, by sampling the original data, as in  [11], and is written to HDFS. It is important to mention here that the QuadTree contains only the spatial partitions and not the actual points. Then, when a new query is posed, the QuadTree is loaded into Hadoop’s distributed cache in order to be accessible by all the nodes. Moreover, at the same level, we employ two indexes. The first index is a spatial index (SpI) which enables pruning of points based on their spatial distance, thus decreasing significantly the number of points that need to be examined within the $\epsilon_{t}$ window. The second index is an index that keeps track of the representation of each individual trajectory within the temporally sorted structure $D$ (TrI), thus providing an efficient way to access the previous trajectory point. The two indexes are created gradually, as the data are read from HDFS. Finally, at the third level, we have the temporally sorted data that correspond to the specific temporal partition.

5.1.1 Spatial Index (SpI)

The spatial index, called SpI, utilizes a given space partitioning, in our case QuadTrees. For each spatial partition of the QuadTree, SpI keeps a temporally sorted array where each entry is the position of a point that is contained in the given partition expanded by $\epsilon_{sp}$. SpI is implemented as a HashMap with key the partition id and value the sorted array. Thus, a partition can be accessed in O(1), while a point in a partition can be accessed in $O(logP_{i})$, where $P_{i}$ here is the number of points in the corresponding sorted array. The construction of SpI has $O(|{D}|\cdot h)$ complexity, where $|{D}|$ is the number of points in the specific temporal partition and $h$ is the height of the QuadTree, since for each point we need to traverse the QuadTree in order to find out in which expanded partition it is contained. Note that each point is enriched with the id of its original (i.e. not expanded) spatial partition, thus consisting of $\left\langle trajID,x,y,t,PartitionID\right\rangle$.

5.1.2 Trajectory Index (TrI)

The TrI index keeps track of each individual trajectory within $D$. TrI is also implemented as a HashMap with key the trajectory id. For each trajectory, the value is a temporally sorted array, where each entry corresponds to a point of a trajectory, and the value of the entry is an integer indicating the point’s position in $D$. Thus, a trajectory point can be efficiently accessed in $O(logT)$, where $T$ is the number of points of a trajectory. To exemplify, the first element of the array holds the position of the first point of the trajectory inside $D$ and so on. The construction of this index has $O(T)$ time complexity since the data is already sorted in time.

5.2 The DTJi Algorithm

Having these two indexes at hand we can utilize them in order to perform the join operation in an efficient way. Algorithm 4, presents the index-enhanced plane sweep procedure. Initially, the QuadTree is loaded into memory from the distributed cache (line 3) and then, each accessed point is inserted not only to an array $D$, which contains points sorted in increasing time, but also to the SpI and TrI indexes. Finally, the TRJPlaneSweepI() algorithm is invoked for each accessed point (lines 4–7).

Algorithm 5, presents the TRJPlaneSweepI() algorithm. Here, given a point $p_{i}\in p$, instead of scanning the whole $\epsilon_{t}$ window before it, in order to find “matches”, we perform a search in SpI and get only the points that belong to the same partition as $p_{i}$ by invoking the getCandidatePoint() method (line 4). The partition id is retrieved in $O(1)$ and then binary search is performed in the temporally sorted list of points in order to find the position of $p_{i}$ inside it. Having that, we can get the previous element, which will be the previous point in time that lies within the same partition, and check if the temporal and spatial constraint are satisfied. If they are satisfied, we have a “match”, we proceed to the previous element of SpI and so on and so forth. Assuming that we have a “match” with $q_{j}$ that belongs to trajectory $q$ we need to find the previous point of $q$. This is achieved by invoking getPrevTrPointI, which performs a search in TrI in order to retrieve in O(1) the entry of $q$ (lines 11, 14). Then, by performing binary search in the temporally sorted list, we can find the position of $q_{j}$ and can easily get $q_{j-1}$. Having that, we need to find if it “matches” with any point that belongs to $p$. Here, instead of scanning the whole $2\epsilon_{t}$ window of $q_{j-1}$ in order to check for “matches” with $p$, we perform a search in TrI in order to get the points of $p$ that exist “close” to the time of $q_{j-1}$ (lines 12, 15). Then, if the spatial and temporal constraints are satisfied we have a “match” and the FindMatchI() method returns True. Otherwise, the whole procedure continues, until the temporal constraint is not satisfied anymore.

The complexity of the index-based solution is $O(|D|\cdot h\cdot(log_{2}P_{i}\cdot a\cdot P_{i}((1-b)\cdot P_{i}+b\cdot P_{i}\cdot(2\cdot(log_{2}T+(log_{2}T+a\cdot T)))))))$, with $|D|$ being the number of points, $h$ the height of the QuadTree, $a$ and $b$ the selectivity of $\epsilon_{t}$ and $\epsilon_{sp}$ respectively. $P_{i}$ is the number of points within the $i$-th partition expanded by $\epsilon_{sp}$, where $P_{i}\ll|D|$, and $T$ is the number points per trajectory. In the worst case, where $a$ and $b$ tend to 1, the complexity can reach $O(|D|\cdot log_{2}P_{i}\cdot P_{i}^{2})$. However, again this only occurs for values of $\epsilon_{t}$ and $\epsilon_{sp}$ that are comparable to the dataset’s duration and diameter respectively. Roughly speaking, the complexity drops to $O(|D|\cdot(log_{2}P_{i}\cdot a\cdot b\cdot P_{i}^{2}))$, which clearly shows the benefit attained when employing the proposed indexing scheme.

6 Experimental Study

In this section, we provide our experimental study on the comparative performance of the three variations of our solution, namely (1) $\mathit{DTJb}$ that uses two MR jobs (Section 3), (2) $\mathit{DTJr}$ that employs repartitioning and a single job to perform the join (Section 4), and (3) $\mathit{DTJi}$ that additionally uses the SpI and TrI indexes for more efficient join processing (Section 5). Furthermore, we compare our solution with the work presented in [33].

The experiments were conducted in a 49 node Hadoop 2.7.2 cluster, provided by * okeanos*222https://okeanos.grnet.gr/home/, an IAAS service for the Greek Research and Academic Community. The master node consists of 8 CPU cores, 8 GB of RAM and 60 GB of HDD while each slave node is comprised of 4 CPU cores, 4 GB of RAM and 60 GB of HDD. Our configuration enables each slave node to launch 4 containers, thus resulting that at a given time the cluster can run up to 192 jobs (Map or Reduce). The real dataset employed for our experiments is IMIS333IMIS dataset has been kindly provided by IMIS Hellas for research and educational purposes. It is available for downloading at http://chorochronos.datastories.org, which consists of 699,031 trajectories of ships moving in the Eastern Mediterranean for a period of 3 years. This dataset contains approximately 1.5 billion records, 56GB in total size.

Our experimental methodology is as follows: Initially, we verify the scalability of our algorithms by varying (a) the dataset size, and (b) the number of cluster nodes (Section 6.1). Then, we examine the benefits of the repartitioning step as well as the associated cost (Section 6.2). Successively, we compare our solution with the work presented in [33] (Section 6.3). Subsequently, we perform a sensitivity analysis in order to evaluate the effect of different parameters to our algorithms (Section 6.4). Finally, we perform a set of experiments so as to examine the creation time and the size of the proposed indexes w.r.t. to varying the number of spatial partitions and $\epsilon_{sp}$ (Section 6.5).

Table 1 shows the experimental setting, where we vary the following parameters: $\epsilon_{t}$, $\epsilon_{sp}$, $\delta t$, the maximum number of points per cell, and the number of cluster nodes, which are the main parameters affecting the performance of our algorithms. In more detail, the values of $\epsilon_{sp}$ were calculated as a percentage of the diameter of the smallest cell produced by the QuadTree, whereas the maximum number of points per cell is calculated as a percentage over the total population.

6.1 Scalability

Initially, we vary the size of our dataset and measure the execution time of our algorithms. To study the effect of dataset size, we created 4 portions (20%, 40%, 60%, 80%) of the original dataset. As the dataset size increases and the number of nodes remains the same, it is expected that the execution time will increase. In order to measure this, for each portion $D_{i}$ of the dataset with $i\in[1,5]$, we calculate $SlowDown=\frac{T_{D_{i}}}{T_{D_{1}}}$, where $T_{D_{1}}$ is the execution time of the first portion (i.e. 20%) and $T_{D_{i}}$ the execution time of the current one. As shown in Figure 8(a), as the size of the dataset increases, the $\mathit{DTJr}$ linear and the $\mathit{DTJi}$ appears to have linear behaviour, with $\mathit{DTJi}$ presenting better scalability. On the other hand, $\mathit{DTJb}$ appears to have a somehow “abnormal” behaviour. This can be justified if we study Figure 9, which presents the standard deviation of the different portions of the dataset. In fact, we can observe that $\mathit{DTJb}$ in Figure 8(a) is affected by how imbalanced is the partitioning in each portion of the dataset, as depicted in $\mathit{DTJb}$ Figure 9. To investigate further the performance of the different algorithms, we measure separately the execution time of the Join and Refine phases for all algorithms. Concerning the Join phase, as illustrated in Figure 8(b), $\mathit{DTJi}$ outperforms $\mathit{DTJb}$ by 16$\times$ and the $\mathit{DTJr}$ by almost one order of magnitude. Regarding the Refine phase, as depicted in Figure 8(b), $\mathit{DTJr}$ and $\mathit{DTJi}$, perform exactly the same, as anticipated, since they use an identical algorithm. Instead, $\mathit{DTJb}$ performs worse due to the fact that the Refine phase is implemented as a second MR job, which means that the output of the Join phase, which is typically several times larger than the input data, needs to be read from HDFS and get sorted, grouped and shuffled to the Reduce tasks.

Subsequently, we keep the size of the dataset fixed (at 100%) and vary the number of nodes. As the number of nodes increase and the dataset size remains the same, it is expected that the execution time will decrease. In order to measure this, for each portion $N_{i}$ of the dataset with $i\in[1,5]$, we calculate $SpeedUp=\frac{T_{N_{i}}}{T_{N_{1}}}$, where $T_{N_{1}}$ is the execution time when using the minimum number of nodes (i.e. 12) and $T_{N_{i}}$ the execution time of the current one. In this experiment, as illustrated in Figure 8(c) and (d), we observe that all three approaches present linear scaling, with $\mathit{DTJi}$ demonstrating slightly better scalability. The reason why the behaviour is different here, is that in this experiment the dataset that was employed was fixed (100%). This means that, despite the fact that we vary the number of nodes, the effect of the different algorithms over the data is the same. On the contrary, when we increase the amount of data, as already shown in Figure 8(b), we can see than the performance of $\mathit{DTJb}$ is affected by the skewness of the different portion of the dataset that were used. As depicted in Figure 9, we can observe that the standard deviation of $\mathit{DTJb}$ affect significantly $\mathit{DTJb}$ in Figure 8(a) and (b).

6.2 Repartitioning and Load Balancing

In this set of experiments, we evaluate the cost of the repartitioning step employed by $\mathit{DTJr}$ and $\mathit{DTJi}$. In Figure 9(a), we compare $\mathit{DTJb}$ (which does not use this step) against $\mathit{DTJr}$ and $\mathit{DTJi}$ after including in the latter two algorithms the time needed for repartitioning. The result shows that even for a single query both algorithms outperform $\mathit{DTJb}$. Obviously, for multiple queries with different query parameters ($\epsilon_{t}$, $\epsilon_{sp}$, $\delta t$), the gain is multiplied, since the repartitioning cost needs to be paid only once, before processing the first query. This experiment justifies the use of the repartitioning step, while demonstrating its low overhead in the case of a single query, which in the case of multiple queries becomes negligible.

In order to quantify whether the work allocation of the Join is balanced to the different parallel tasks, we compare the input size of the Join phase of $\mathit{DTJb}$, against the Join phase of $\mathit{DTJr}$. In Figure 9(b), we report the standard deviation of the size of input data for the various tasks. Smaller values of the standard deviation, indicate that the different tasks are assigned with similar-sized input data, thus the load is more fairly balanced. $\mathit{DTJr}$ demonstrates significantly lower standard deviation, approximately one order of magnitude, than $\mathit{DTJb}$. This also partly justifies the overall better performance of $\mathit{DTJr}$ illustrated in Figure 8(b) and Figure 8(d).

6.3 Comparative Evaluation

As already mentioned, the problem of Distributed Subtrajectory Join has not been addressed yet in the literature and it is not straightforward (if and) how state of the art solutions to trajectory similarity search and trajectory join can be adapted to solve the problem. However, if we utilize only a specific instance of our problem, when $\delta{t}=0$, then we only need to identify the set of $JP$ during the Join phase. Based on this observation, we select to compare with the work presented in  [33], called $SJMR$, a state of the art MapReduce-based spatial join algorithm, which is able to identify efficiently the set of $JP$ that will be passed to the Refine procedure and produce the desired result. The reason why $SJMR$ was chosen is that it is a generic solution which could form the basis for any distributed spatial join algorithm and thus required the minimum amount of modifications so as to match with our problem specification.

More specifically, $SJMR$ repartitions the data at the Map phase and Joins them at the Reduce phase by performing a plane sweep join. For the sake of comparison, we modified $SJMR$ by injecting time as a third dimension and introducing parameters $\epsilon_{sp}$ and $\epsilon_{t}$. In more detail. at the Map phase the spatiotemporal space is divided to tiles using a fine grained grid. Then, each data point is expanded by $\epsilon_{sp}$ and $\epsilon_{t}$ and is assigned to the tiles with which it intersects. Subsequently, the tiles are mapped to partitions using the method described in  [33]. At the Reduce phase, the points are grouped by partition and sorted by one of the dimensions (we chose the temporal dimension so as to be aligned with our solution). Finally, we sweep through the time dimension and report the set of $JP$.

So, in this set of experiments we compare $\mathit{DTJi}$-Join, which outperforms $\mathit{DTJb}$-Join and $\mathit{DTJr}$-Join, with $SJMR$. In more detail, we vary the size of our dataset and measure the execution time of the two algorithms. The results, as illustrated in Figure 10 show that $\mathit{DTJi}$-Join not only performs significantly better than $SJMR$ but more importantly, the gain of $\mathit{DTJi}$-Join over $SJMR$ increases for larger data sets. The reason for this behaviour lies mainly due to the utilization of the indexing structure of $\mathit{DTJi}$ ( [33] uses no indexes) and the fact that $\mathit{DTJi}$-Join is a Map-only job where the repartitioning cost is “paid” only once (as a preprocessing step), unlike $SJMR$, where this cost is “paid” every time at the Map phase, as explained earlier.

6.4 Sensitivity Analysis

In the following experiment, we perform a sensitivity analysis of algorithms $\mathit{DTJr}$ and $\mathit{DTJi}$. We exclude $\mathit{DTJb}$ from this set of experiments, as it consistently performs significantly worse than the other two algorithms.

Initially we vary the value of $\epsilon_{t}$ while retaining fixed the values of the other parameters. As shown in Figure 11(a), the execution time of both algorithms, as expected, increases with $\epsilon_{t}$. In more detail, the Join phase of $\mathit{DTJr}$ is more sensitive to the fluctuation $\epsilon_{t}$ than the Join phase of $\mathit{DTJi}$, due to the fact that the latter is utilizing the SpI index which, for a given $\epsilon_{t}$, performs spatial filtering instead of scanning the entire space in order to find “matching” pairs of points. What is more interesting is that as $\epsilon_{t}$ increases the difference between the two approaches increases, which means that for higher values of $\epsilon_{t}$ the difference, in terms of execution time, will be higher than one order of magnitude. As far as it concerns the Refine step, both approaches present the same increasing behavior when $\epsilon_{t}$ increases, since both of them employ the same algorithm, due to the fact that the higher the $\epsilon_{t}$, the larger the sliding window that is created.

Then, we set different values to $\epsilon_{sp}$ while keeping the values of the other parameters fixed. As illustrated in Figure 11(b), $\epsilon_{sp}$ affects directly the Join and indirectly the Refine phase of both approaches. More specifically, the Join phase of $\mathit{DTJr}$ is slightly affected by setting different values to $\epsilon_{sp}$ due to the fact that, for a given $\epsilon_{t}$, $\mathit{DTJr}$ will search the whole space in order to find “matches”. Hence, $\epsilon_{sp}$ will only affect the number of “matches”. On the other hand, $\mathit{DTJi}$ does not search the whole space but utilizes the SpI index which, consequently, makes it more sensitive to $\epsilon_{sp}$. The only case where the Join phase of $\mathit{DTJr}$ performs the same as $\mathit{DTJi}$ is when $\epsilon_{sp}$ spans the whole dataset space. Regarding the Refine step, as expected, both approaches perform the same and the higher the $\epsilon_{sp}$, the higher the execution time. The reason for this behaviour is that as $\epsilon_{sp}$ increases, the product of the Join phase increases.

Finally, we vary the values of $\delta{t}$ while keeping the values of the other parameters fixed. As presented in Figure 11(c), this parameter affects only the Refine phase, as anticipated. More specifically, the higher the $\delta{t}$ the slightly higher the execution time of both approaches. This takes place due to the fact that as $\delta{t}$ increases, the sliding window gets larger.

6.5 Indexing

In order to measure the effect of having different number of spatial partitions in spatial index size and spatial index construction time, we perform a final set of experiments. More specifically, we vary the maximum number of points per cell parameter of the QuadTree, and we measure the index creation time and the index size. As illustrated in Figure 12(a), the SpI construction time increases as the maximum number of points per cell decrease, while the TrI construction time is, as expected, not affected by that. This occurs due to the fact that as the maximum number of points per cell decreases, the number of spatial partitions increases. Furthermore, as depicted in Figure 12(a), the fewer the maximum number of points per cell the smaller the execution time of the Join algorithm. It is worth mentioning that the overall index construction time as a percentage over the execution time of the Join algorithm varies only between 4% and 11%.

As far as the size of the indexes is concerned, Figure 12(b) illustrates how the size is affected when varying the maximum number of points per cell. As expected, the TrI index is not affected, whereas the SpI index slightly increases its size as the number of partitions increase. At this point, we should mention that compared to the size of $D$, the percentage of the total size of the indexing scheme over $D$ varies only between 24% and 28%.

Another parameter that can affect the SpI index is $\epsilon_{sp}$ due to the fact that each spatial partition is enlarged by $\epsilon_{sp}$. As depicted in Figure 12(c), the size of SpI increases as $\epsilon_{sp}$ increases.

7 Related Work

The work performed in this paper is closely related to three topics in the spatial and spatio-temporal database management literature. These are, (a) centralized trajectory joins, (b) MapReduce spatial and multidimensional joins, and (c) distributed trajectory joins. In the following paragraphs we outline the related work and try to highlight their differences from our work.

Centralized Trajectory Joins: In [3] the effort focuses in identifying pairs of trajectories that move close enough, with respect to a spatial threshold, during a user specified temporal window, unlike our approach where we identify “matches” of different duration (at least $\delta t$) during the whole lifespan of the datasets. Furthermore, in [3], no temporal tolerance is considered which can lead in missing pairs of trajectories that move with some temporal displacement. In  [5], the authors try to solve the same problem in a streaming environment. In [4] the authors extend their work by not binding to the temporal dimension the interval in which two trajectories should move “together”. Hence, all pairs of “matching” (w.r.t a spatial threshold) subtrajectories of exactly $\delta t$ duration will be returned. This definition, although more general, it still suffers from the rest of the aforementioned problems. Moreover, in these approaches there is an assumption made, that all trajectories have the same number of points which are synchronized. However, an assumption like that is not realistic in real life scenarios and supposes a preprocessing step that can be prohibitive when dealing with Big Data. A slightly different definition is provided in [7] where the goal is to identify all pairs of moving objects that, for some time intervals, move closer than a given spatial threshold. Here, the duration of the “matches” is not fixed. The most significant difference with our approach is that, here, they identify only the joining points, which leads in a trajectory join definition which is asymmetric. Moreover, time relaxation is not considered, so the distance between two objects refers to their distance at the same time point $t$. Furthermore, the minimum duration of the “matches” cannot be limited which can lead to pairs with very small duration that might not be useful for some applications. Finally, the solution provided is focused to an instantiation of the problem, called Window Trajectory Distance Join, which limits the problem to a user-specified time window.

A similar but different problem is the one of trajectory similarity join, where the goal is to retrieve all pairs of trajectories that exceed a given similarity threshold as in  [30] and  [9]. However, both of them return as a result pairs of trajectories and not subtrajectories, thus they cannot support some of the scenarios presented in Section 1. An approach very similar to ours is presented in  [6], where, given a pair of trajectories they try to perform partial matching, finding the most similar subtrajectories between these two trajectories. Different variations of the problem are presented, where the duration of the “match” is specified beforehand or not. Nevertheless, the problem in [6] is not a join operation and temporal tolerance is not considered. To sum up, all of the above approaches are centralized and applying them to a parallel and distributed environment is non-trivial.

MapReduce Spatial & Multidimensional Joins: A special class of joins which is very relative to our problem is that of spatial join. There have been several efforts to tackle this issue using the MapReduce framework. In particular, [33], which is based on the traditional PPBSM algorithm [21], partitions the input data into small evenly disjoint tiles at Map stage and joins them at Reduce stage by further partitioning the data into strips and performing a plane sweeping algorithm along with a duplication avoidance technique. Other works that try to deal with the problem of spatial join are  [2] and [11]. In more detail,  [2] first partitions the data by focusing on recursively breaking high density tiles into smaller ones. Objects that exceed the borders of a partition are replicated and a post-processing step is employed in order to eliminate duplicate results. The join process takes place at the Reduce phase by utilizing the R*-Tree indexes that are created and loaded in-memory at query time for each of the relations. Another system that copes with the problem of spatial join is [11]. This approach first re-partitions the data by taking into account load balancing and spatial locality. In more detail the data are sampled and an index (Grid File, R-Tree or R+-Tree) is created which will set the boundaries of each partition. During, the re-partitioning phase global and local indexes are created and stored in HDFS. The join takes place in the Map phase by utilizing the local indexes. Concerning borderline objects, a duplicate avoidance method is applied. An approach that enhances [11] with the functionality of identifying closest pairs of points is presented in [14]. Multidimensional similarity join is also related to our work. In  [29, 28] the problem of distance range join is studied, which is probably the most common case of similarity join. In this approach the data is iteratively partitioned similarly to the Quickjoin algorithm [16], which results in having multiple MR jobs in order to get the final results. The problem of high dimensional similarity joins on massive datasets using MapReduce is tackled in [19]. In [23] the problem of $\epsilon$-distance similarity self-join on vector data is tackled by employing in the Map phase a fixed size grid with cell width $\epsilon$ and assign the data to the corresponding cells. In order to compute the $\epsilon$-neighborhood of each cell only the adjacent cells are needed. In order to reduce the replication of data they avoid taking into account all the adjacent cells. Despite this, the algorithm used to perform the join in the Reduce phase is still a Nested Loop Join. [13] is an extension of [23] for medium- to high-dimensional spaces where the full $d$-dimensional space is broken down to $k$ dimension groups, the join is performed in each group and then the results are merged. Unlike [23], they try to cope with the skewness of data by starting with a very fine grid and merging cells until a balanced grid is created. Similarly to [23] the join process is still a Nested Loop Join. However, all of the above approaches try to solve a problem that is significantly different from ours since our problem is not to join spatial or multidimensional objects but to join sequences of spatial points. In more detail, the output of the above approaches would be the set of $JP$. However, as already stated, in order to get the desired result we also need to the sets of $BP$ and $SNJP$, which cannot be acquired by post-processing the result of these approaches.

Distributed Trajectory Joins: The approaches mainly focus to the $k$-nn join or the $k$ most similar trajectories join. More specifically,  [32] and  [31] address the problem where given a reference trajectory and an integer $k$ they want to discover the $k$ most similar trajectories to the reference trajectory. The problem of $k$-nn join by using the MapReduce framework is tackled in  [12]. More specifically, given two sets of trajectories $R$ and $M$, an integer $k$ and a time interval $[t_{s},t_{e}]$, the algorithms proposed there return the k nearest neighbors from $R$ for each object in $M$ during this interval. The above approaches, address an entirely different problem than the one presented in this paper. Recently, the algorithms proposed in [24, 25] find all pairs of network-constrained trajectories that exceed a similarity threshold in a parallel manner. However, the parallelization proposed there handles each trajectory separately by assuming that all data need to be replicated for each trajectory, which makes such a solution inapplicable to the Big Data setting. Finally, these approaches (a) assume that the underlying network is known, which is not something trivial in some domains (e.g. maritime or aviation) and (b) work at the entire trajectories and cannot identify matching sub-trajectories.

8 Conclusions

In this paper, we introduced the Distributed Subtrajectory Join query, an important operation in the spatiotemporal data management domain, where very large datasets of moving object trajectories are processed for analytic purposes. To address this problem in an efficient manner following the MapReduce programming model, we initially provided a well-designed basic solution which is used as a baseline in order to propose two efficient improvements, called $\mathit{DTJr}$ and $\mathit{DTJi}$ which can boost the performance by up to 16$\times$ and 10$\times$, respectively. Our experimental study was performed on a very large real dataset of trajectories from the maritime domain,consisting of 56 GB of data (or 1.5 billion time-stamped locations). As for future work, the solution to the Distributed Subtrajectory Join problem that we provided can be utilized as the basis for efficiently identifying various mobility patterns (e.g. group behaviours) and discovering clusters of moving objects over massively distributed data.

9 Acknowledgements

This work was partially supported by projects datACRON (grant agreement No 687591), Track&Know (grant agreement No 780754) and MASTER (Marie Sklowdoska-Curie agreement N. 777695), which have received funding from the EU Horizon 2020 R&I Programme. We would like to thank the Greek Research and Technology Network (GRNET) for providing access to the Okeanos cloud infrastructure and IMIS Hellas for donating the IMIS dataset.

Appendix A Appendix

Proof of Lemma 1

Proof.

It suffices to construct an instance of the problem where an algorithm $A$ that operates solely on joining points and is unaware of breaking points would produce erroneous results, thus $A\notin\mathcal{A}$. Let $r^{\prime}=r_{1,n}$ and $s^{\prime}=s_{1,(n+1)}$ denote two subtrajectories, such that there exist $n$ pairs of joining points ($r_{i}$,$s_{j}$) and $\Delta w_{r^{\prime},s^{\prime}}=\delta t$. However, let us assume that there exists a point $s_{k}\in s^{\prime}$ which is a breaking point. Based on the problem definition (Defn. 1), if algorithm $A$ was unaware of breaking points, it would falsely identify $r^{\prime}$ and $s^{\prime}$ as “matching” subtrajectories. ∎

Proof of Lemma 2

Proof.

The proof is similar to the proof of Lemma 1, only using a non-joining point instead of a breaking point in the constructed instance of the problem. ∎

Proof of Lemma 3

Proof.

Joining points: By contradiction. Let us assume that an expanded partition $expPart_{i}$ is not sufficient to produce the set of $JP$ in $part_{i}$. Then, there must exist a pair of joining points $r_{j}\in r$ and $s_{k}\in s$, such that $r_{j}$ belongs to partition $part_{i}$, whereas $s_{k}$ does not belong to $expPart_{i}$. Based on the definition of expanded partitions, it follows that $DistT(r_{j},s_{k})>\epsilon_{t}$. Thus, $r_{j}$ and $s_{k}$ cannot be joining points, which is a contradiction. Breaking points: By contradiction. Let us assume that an expanded partition $expPart_{i}$ is not sufficient to produce the set of $BP$ for trajectories in $part_{i}$. Then, there must exist a point $r_{*}$ (that belongs to $part_{i}$) of trajectory $r$, such that $r_{*}$ is a breaking point. To identify if $r_{*}$ is a breaking point, we need to examine whether there exists a point $s_{j}\in s$ of any other trajectory $s$ with $DistT(r_{*},s_{j})\leq\epsilon_{t}$ (Defn. 4). However, based on the definition of expanded partitions, such a point $s_{j}$ must belong to $expPart_{i}$, which contradicts with the assumption that $expPart_{i}$ is not sufficient. ∎

Proof of Lemma 4

Proof.

Consider two pairs of $JP$, $(r_{i},s_{j})$ and $(r_{i+1},s_{j+n})$, where $r_{i},r_{i+1}\in r$ and $s_{j},s_{j+n}\in s$. Moreover, let us assume that for each point $s_{m}$ where $m\in[j+1,j+n-1]$, $s_{m}$ is a non-joining point for each point of $r$. Then, points $r_{i}$, $r_{i+1}$, $s_{j}$, $s_{m}$ and $s_{j+n}$ will span at most to two consecutive temporal partitions: $part_{k}$ and $part_{k+1}$. This means that either $s_{j},s_{m}\in part_{k}$ or $s_{m},s_{j+n}\in part_{k+1}$. In both cases $s_{m}$ will be recognized as a non-joining point for each point of $r$, by the procedure that generates the set of $sNJP$(Defn. 6). ∎