Deep Representation Learning for Trajectory
Similarity Computation
Xiucheng Li† , Kaiqi Zhao† , Gao Cong† , Christian S. Jensen‡ , Wei Wei◦†

Nanyang Technological Univeristy, ‡ Aalborg University, ◦ Huazhong University of Science and Technology{xli055@e., kzhao002@e., gaocong@}ntu.edu.sg, csj@cs.aau.dk, weiw@hust.edu.cn

Abstract—Trajectory similarity computation is fundamental
functionality with many applications such as animal migration
pattern studies and vehicle trajectory mining to identify popularroutes and similar drivers. While a trajectory is a continuous curve in some spatial domain, e.g., 2D Euclidean space,trajectories are often represented by point sequences. Existing
approaches that compute similarity based on point matching
suffer from the problem that they treat two different point
sequences differently even when the sequences represent thesame trajectory. This is particularly a problem when the point
sequences are non-uniform, have low sampling rates, and havenoisy points. We propose the first deep learning approach tolearning representations of trajectories that is robust to low
data quality, thus supporting accurate and efficient trajectory
similarity computation and search. Experiments show that our
method is capable of higher accuracy and is at least one order ofmagnitude faster than the state-of-the-art methods for k-nearest
trajectory search.

I. I NTRODUCTION
With the proliferation of GPS-enabled devices, trajectorydata is being generated at an unprecedented speed. A trajectoryis typically represented as a sequence of discrete locations,or sample points, which describes the underlying route of amoving object over time. Quantifying the similarity between
two trajectories is a fundamental research problem and isa foundation for many trajectory based applications, suchas tracking migration patterns of animals [1], mining hotroutes in cities [2], trajectory clustering [3] and moving group
discovery [4], [5]. The importance of measuring trajectory
similarity has also been recognized by researchers, and anumber of classic methods have been proposed, such as
dynamic time wrapping (DTW) [6], longest common subsequence (LCSS) [7], edit distance with real penalty (ERP) [8],and edit distance on real sequences (EDR) [9].
These existing methods usually try to form a pairwise
matching the sample points of two trajectories and identify the best alignment using dynamic programming. Morespecifically, these pairwise point-matching methods implicitly
partition the space into cells according to a threshold  andmatch two points if they fall into the same cell. Dynamic
programming is then used to find an alignment that minimizes
a matching cost.
The pairwise point-matching methods for computing trajectory similarity often suffer in three scenarios: the first is when
the sampling rates of trajectories are non-uniform. Sampling
rates often vary across devices due to different device settings,
battery constraints, and communication failures. Even for thesame device, the sampling rates may vary. For example, a taxi
driver may alter the default device sampling rate to reduce
the power consumption [10]. This may result in sampling
points alternating between sparse and dense episodes. Thisposes challenges to the existing methods—if two trajectories
represent the same underlying route, but are generated at
different sampling rates, it is difficult for these methods toidentify them as similar trajectories. The second scenario
where it is challenging to align the sample points of similar
trajectories is when the sampling rate is low. For example, thesampling rates for trajectories generated from online “checkins” (e.g., Foursquare, Facebook), geo-tagged tweets, geotagged photo albums, and call detail records are low andinherently non-uniform [11]. Third, the performance of thesemethods may be degraded when the sample points are noisy.
Such noise may occur in GPS points when moving in urban
canyons.
Moreover, the existing methods rely on local point matching
and identify an optimal alignment using dynamic programming, which leads to quadratic computational complexity
O(n2 ), where n is the mean length of the trajectories. We
argue that a good trajectory similarity measure should be bothaccurate and scalable to large datasets.
We illustrate the problems caused by non-uniform andlow sampling rates with two examples. First, consider thetwo trajectories Ta = [a1 , a2 , a3 ] and Tb = [b1 , . . . , b6 ] in
Figure 1a that are sampled from the same underlying routewith different sampling rate. Assuming that the cell threshold is 1, a pairwise point-matching method such as EDR will
match (a1 , b1 ) and (a3 , b6 ) while the remaining points will be
unmatched, which yields a cost (edit distance) of 5. Although
the two trajectories represent the same underlying route,
they differ when compared using the pairwise point-matching
methods. Second, to see the challenge of low sampling rates,consider the trajectory in Figure 1b consists of sample points[a1 , . . . , a5 ]. Due to the low sampling rate in the middle partof the trajectory, the pairwise point-matching methods may
falsely match the two trajectories because four point pairs arematched.
The essential problem underlying the two examples is how
to infer the true route from the observed trajectory points in thepresence of non-uniform sampling rates, low sampling rates,
and noise. As shown in Figure 1b, it is hard to tell whichroute, e.g., RA , RB is represented by the set of sample points[a1 , . . . , a5 ]. Recent studies [10], [12] reveal that the transitingFig. 1: Examples of the challenge in quantifying trajectory
similarity.
patterns between certain locations are often highly skewed,
i.e., some routes are more likely to be traveled than others,
and these transition patterns, which are accumulated in thespatiotemporal databases [10], [12], [13], hold the potential tohelp quantify trajectory similarity. To harness such transition
patterns, Su et al. [13] proposed an Anchor Points based
Method (APM). APM first learns the transfer relationships
among a fixed set of spatial objects (such as Points ofInterest), called anchor points, from dense (i.e., with highsampling rate) historical trajectories by using Hidden Markov
Models (HMM). Then sparse trajectories are calibrated tothe anchor points, such that the existing pairwise pointmatching methods, e.g., DTW, EDR, LCSS, can be employed
to compute similarities of trajectories more accurately. APM
requires the availability of a large amount of POIs and suffersfrom the limitations inherent in HMM like requiring explicit
dependency assumptions to make inferences tractable [14].
Moreover, even after being trained on a historical dataset,
APM still cannot reduce the O(n2 ) time complexity of the
pairwise point-matching methods.
In this paper, we propose a novel approach, called t2vec(trajectory to vector), to inferring and representing the underlying route information of a trajectory based on deep
representation learning. The learned representation is designed
to be robust to non-uniform and low sampling rates, andnoisy sample points for trajectory similarity computation. Thisis achieved by taking advantage of the archived historicaltrajectory data and a new deep learning framework. With thelearned representation, it only takes a linear time O(n + |v|)(|v| is the length of vector v) to compute the similarity
between two trajectories, while all the existing approaches
take O(n2 ) time. To the best of our knowledge, this is thefirst deep learning based solution for computing the similarityof trajectories.
To learn trajectory representations, it is natural to consider
the use of Recurrent Neural Networks (RNNs), which areable to embed a sequence into a vector. However, simply
applying RNNs to embed trajectories is impractical. First,
the representation obtained using RNNs is unable to reveal
the most likely true route of a trajectory when uncertainty
arises due to low sampling rates or noise. Second, the existing
loss functions used to train RNNs fail to consider spatial
proximity, which is inherent in spatial data. Thus, they can-not guide the model to learn consistent representations for
trajectories generated by the same route. To overcome thefirst challenge, we propose a sequence-to-sequence (seq2seq)
based model to maximize the probability of recovering thetrue route of trajectory. To contend with the second challenge,
we design a spatial proximity aware loss function and acell pretraining algorithm that encourage the model to learn
consistent representations for trajectories generated from thesame route. We also propose an approximate loss function
using Noise Contrastive Estimation [15] to boost the training
speed. Overall, the paper makes the following contributions:• We propose a seq2seq-based model to learn trajectory
representations, for the fundamental research problem oftrajectory similarity computation. The trajectory similarity based on the learned representations is robust to nonuniform, low sampling rates and noisy sample points.
Our solution computes the similarity of two trajectories
in linear time.
• For the purpose of learning consistent representations, wedevelop a new spatial proximity aware loss function anda cell representation learning approach that incorporate
the spatial proximity into the deep learning model. To
further speed up training, we propose an approximate loss
function based on Noise Contrastive Estimation.
• We conduct extensive experiments on two real-world
trajectory datasets that offer evidence that the proposed
method is capable of outperforming the existing trajectory
similarity measure techniques in terms of both accuracy
and efficiency.
The rest of the paper is organized as follows. In Section II,
we discuss the related work. The problem definition andpreliminaries are given in Section III. Section IV presents thedetails of our method. The experimental results are presented
in Section V. Finally, we summarize the paper and discuss
future research directions in Section VI.
II. R ELATED W ORK
We briefly review the related work on trajectory similarity
computation and deep representation learning.
A. Trajectory similarity computation
Computing the similarity between two trajectories is fundamental functionality in many spatiotemporal data analysis
tasks. Not surprisingly, the problem of accurately and efficiently measuring the similarity of trajectories has been studied
extensively [6]–[9]. DTW [6] was a first attempt at tackling
the local time shift issue for computing trajectory similarity.
ERP [8], EDR [9], DISSIM [16], and the model-driven
approach MA [17] were developed to further improve theability of capturing the spatial semantics in trajectories. Wanget al. [18] studied the effectiveness of these similarity methods
according to their robustness to noise, varying sampling rates,and shifting. All of these methods focus on identifying the
optimal alignment based on sample point matching, and thus
they are inherently sensitive to variation in the sampling rates.
To solve this issue, APM [13] and EDwP [11] are proposed.

As discussed in the introduction, APM solves this issue bylearning transition patterns of anchor points from historicaltrajectories. To compute the similarity of two trajectories,
EDwP computes the cheapest set of replacement and insertion
operations using linear interpolation to make them identical.
Our solution is very different from APM and EDwP in that
we aim to learn a vector that represents a trajectory andto then compute similarity using the new representation. In
experiments, we compare with EDwP and not APM for two
reasons: i) The implementation of APM partly requires anabundance of POIs which are not available in our datasets; ii)
the more recent EDwP has been reported to perform betterin similarity analysis of trajectories with non-uniform and low
sampling rates. Our work is also related to the inference ofhidden routes from partial observations. Zheng et al. [10] firststudied the problem, and Banerjee et al. [12] further exploredit using Bayesian posterior inference to estimate the top-k
most likely routes. Our work differs from these two in that
our ultimate goal is to learn representations of the trajectories
rather than solely inferring the most possible routes.
Moreover, all the aforementioned existing measures for
trajectory similarity are based on the dynamic programming
technique to identify the optimal alignment which leads toO(n2 ) computation complexity. Given the complexity, it will
be computationally expensive if we want to apply thesesimilarity measures to cluster a large trajectory database. In
contrast, our method has a linear time complexity O(n + |v|)
to measure the similarity of two trajectories, which is able tosupport analysis on big trajectory data, such as clustering trajectories. We can also offer near-instantaneous response times
that support interactive use, while the competition cannot.

B. Representation learning
Learning representations for specific tasks has been alongstanding open problem in machine learning. Recently,
inspired by the success of word2vec [19], the idea of learning
general representation has been extended to paragraphs [20],networks [21], [22], etc. To capture the sequential order information emerging in the sequence processing tasks, Recurrent
Neural Networks (RNNs) based encoder-decoder models havebeen developed, such as sequence to sequence learning [23]–[25], and skip-thought vectors [26]. Our model is based on
the general sequence encoder-decoder framework. However,these existing sequence encoder-decoder models were initially
proposed for natural language processing to deal with discrete
tokens (i.e., words, punctuations). Our scenarios is different inthat the tokens inherently share the spatial proximity relation.
Our model therefore differs from the above sequence encoderdecoders in two ways: i) we design a spatial proximity aware
loss function and a cell representation pretraining approach toincorporate the spatial proximity into the deep representation
model, and ii) we also propose an approximate loss function
based on Noise Contrastive Estimation to accelerate the training.

III. D EFINITIONS AND P RELIMINARIESIn this section, we present definitions and preliminaries
essential to understand the problem addressed and the sequence encoder-decoder model used in our solution. For easeof reference, frequently used notation is given in Table I.
TABLE I: Frequently used notation.
SymbolT (or Ta )
R
x
xt
x1:t
|T | (or |x|)v
htV
r1 (r2 )Definition
Trajectory
Underlying route
A sequence of tokens
Token at position t
A sequence of tokens from position 1 to t
Length of T (or x)
Embedded vector
Hidden state (vector)
Vocabulary
Dropping rate (distorting rate)

A. Definitions
We next define the notions of underlying route and trajectory.
Definition 1. (UNDERLYING ROUTE) An underlying routeof a moving object is a continuous spatial curve (e.g., in thelongitude-latitude domain), indicating the exact path taken by
the object.
The underlying route is only a theoretical concept as location acquisition techniques do not record moving locations
continuously.
Definition 2. (TRAJECTORY) A trajectory T is a sequence ofsample points from the underlying route of a moving object.
In practice, an underlying route can be represented by enormous trajectories, depending on the specifics of the moving
objects and the sampling strategies used. Each generated trajectory can be considered as a representative of an underlying
route. In the rest of this paper, a trajectory is also referred toas trip, depending on the context.
Problem statement. Given a collection of historical trajectories, we aim to learn a representation v ∈ Rn (n is thedimension of a Euclidean space) for each trajectory T suchthat the representation can reflect the underlying route of the
trajectory for computing trajectory similarity. The similarity of
two trajectories based on the learned representations must be
robust to non-uniform, low sampling rates and noisy sample
points.
Our proposed method for solving the problem is based on
deep representation learning techniques. Specifically, we adapt
the sequence encoder-decoder framework for the first time tocompute trajectory similarity (the motivation for adapting thatparticular framework is explained in Section IV-A).
B. Preliminaries of sequence encoder-decoders
We briefly present the sequence encoder-decoder frame|x||y|work. Consider two sequences x = hxt it=1 and y = hyt it=1

where each xt and yt denotes token (e.g., a word or punctuation mark in a natural language sentence) and |x| and|y| represent lengths. We next illustrate how to build theconditional probability P(y|x) in the framework.
The sequence encoder-decoder model has two main
components—an encoder and a decoder, as depicted in Figure 2. The encoder is used to encode sequence x into a fixeddimensional vector v that preserves the sequential informationin x, and then the decoder decodes out sequence y conditionedon the encoded representation v. Since Recurrent Neural
Networks (RNNs) [27], [28] accept input in the form of realvalued vectors, a token embedding layer is added to embedthe discrete token in a vector. The token embedding layer isform of a feed forward neural network [29].

yt−1 is further encoded into ht along with y1 , y2 , . . . , yt−2 .
Finally, P(yt |y1:t−1 , v) can be modeled as follows.
P(yt = u|y1:t−1 , v) = P(yt = u|ht ) = P

exp(Wu> ht )>v∈V exp(Wv ht )Here, W > is the projection matrix that projects ht from thehidden state space into the vocabulary space, Wu> denotes its
u-th row, and V is the vocabulary.
IV. P ROPOSED M ETHOD
We first present the underlying motivation for the proposed
method and the challenges to be addressed when developing
it. Then we describe how to handle non-uniform and low
sampling rates, as well as noisy sample points in our model.
In Section IV-C, we discuss the details of the proposed spatial
proximity aware loss function and cell representation learning
approach in order to encourage the sequence encoder-decoder
to learn consistent representations for trajectories. We discuss
the time complexity of using our model to compute trajectory
similarity in Section IV-D.
A. Motivations and challenges

Fig. 2: Sequence encoder-decoder model. The model reads
the sequence x and outputs the sequence y, where EOS isa special token indicating the end of a sequence and v is therepresentation of x. The hidden state ht captures the sequential
information in [x, y1 , y2 , . . . , yt−1 ].
To see how the model builds the probability P(y|x), we firstrewrite it by the chain rule, i.e.,P(y1 , . . . , y|y| |x1 , . . . , x|x| ) = P(y1 |x)|y|Y

P(yt |y1:t−1 , x),t=2where y1:t−1 represents y1 , y2 , . . . , yt−1 . The encoder reads
and encodes sequence x into the fixed-dimensional vector v.
Since v encodes the sequential information in x, we haveP(yt |y1:t−1 , x) = P(yt |y1:t−1 , v)
The decoder builds the probability P(yt |y1:t−1 , v) at every
position t by squashing v and y1:t−1 into the hidden state
ht , which is simply a forward computation. More specifically,
ht is computed from the output of the previous position ht−1
and token yt−1 , as follows:(f (v, EOS)
t=1
ht =
(1)
f (ht−1 , yt−1 ) t ≥ 2,
where f (·, ·) indicates the RNNs forward computation and
EOS is the special token that signals the end of a sequence,
which is necessary in order to support variable-length sequences [24].
Note that we recursively compute ht = f (ht−1 , yt−1 ): ht−1
encodes the information of v along with y1 , y2 , . . . , yt−2 , andRecent work [10], [12] reveals that the transition patterns
between road network locations are often highly skewed.
This implies that some routes have higher probability of
being traveled than others. Rich movement patterns have beenarchived in spatiotemporal databases (we have 20 different
data sources at Daisy1 and different data sources have different
sampling rates), which then afford us the opportunity to learn
representations for trajectories that overcome the shortcomingsof the point-matching methods.
Recall that we intend to learn a representation v ∈ Rn for a
trajectory T that is robust to non-uniform, low sampling rates
and noisy sample points when using it to compute trajectory
similarity. To learn a representation for sequential data, it isvery natural to consider RNNs [27], [28] since these havebeen shown successful at handling sequences in natural language processing, including generating sequences [30], neuralmachine translation [24], [25] and paraphrases detection [26].
However, by simply applying RNNs we will not be able toaccomplish our goal, which is due to two difficulties: First,
the trajectory representation obtained using RNNs is unable
to represent the most likely underlying route that generates
the trajectory when uncertainty arises due to low sampling
rates and noise. Second the loss functions employed in natural
language processing do not consider the spatial proximity
information inherent in the spatial data, so simply adapting
existing loss functions to our scenario will fail to guide RNNs
to learn consistent representations for trajectories that sharethe same underlying route.
To tackle the first difficulty, the desired model must be able
to maximize the conditional probability P(R|T ), i.e., finding
the most likely underlying route R for the given trajectory
T . However, the underlying route R is often not available in
1 http://www.daisy.aau.dk/

practice. To circumvent this, we exploit two observations: i)both a non-uniform, relatively low sampling rate trajectory,
denoted by Ta , and a relatively high sampling rate trajectory,
denoted by Tb , are paraphrases of their underlying route,and ii) a relatively high sampling rate trajectory Tb is closer
to their true underlying route R than is Ta , and it has
lower uncertainty. These observations cause us to replace
the objective of maximizing P(R|Ta ) into the objective of
maximizing P(Tb |Ta ) and to build the model using a sequence
encoder-decoder framework. The encoder embeds Ta into
its representation v, and the decoder will try to recover its
counterpart Tb with relatively high sampling rate conditioned
on v by optimizing its parameters. When the model is trainedusing the real-world trajectories, the transition patterns hidden
in historical data will be learned by the model. To overcome
the second difficulty, we propose a new spatial proximity
aware loss function and cell (token) representation pre-training
method to incorporate spatial proximity into the model. To
accelerate the training, we develop an approximate spatial
proximity aware loss function based on Noise Contrastive
Estimation [15].
B. Handling varying sampling rates and noise
Based on the above analysis, given a collection of sampling(i)
rate trajectories {Tb }N
i=1 (where N is the cardinality of thecollection), we create a collection of pairs (Ta , Tb ), where Tb isan original trajectory and Ta is obtained by randomly dropping
sample points from Tb with dropping rate r1 . By doing so,
each down-sampled Ta is also non-uniformly sampled andthus represents a real-life trajectory with non-uniform and low
sampling rate. The start and end points of Tb are preservedin Ta to avoid changing the underlying route of the downsampled trajectory. To illustrate, consider generating the subtrajectories for Tb in Figure 1b, we randomly drop points inb2:5 , i.e., all generated sub-trajectories will start with b1 andend with b6 . After the generating procedure, we maximize
the joint probability of all (Ta , Tb ) pairs with the sequence
encoder-decoder model:

influence of noisy sample points, we only keep the cells which
are hit by more than δ sample points. These cells are referred
to as hot cells and form the final vocabulary V (in the rest of
paper, we will interchangeably use token and cell to refer toan element V ). Sample points are represented by their nearest
hot cell. To make the learned representations more robust to
the noisy data, we further distort each downsampled Ta basedon a distorting rate r2 to create the distorted variants, i.e., werandomly sample a fraction of the points (size indicated by r2 )that are then distorted. Point (px , py ) is distorted by adding a
Gaussian noise with a radius 30 (meters) as follows,
px = px + 30 · dx ,dx ∼ Gaussian(0, 1)

py = py + 30 · dy ,dy ∼ Gaussian(0, 1)(3)

We can optimize the same objective as shown in Equation 2
where Ta is both downsampled and distorted.
C. Learning consistent representations
The original sequence encoder-decoder does not model the
spatial correlation between cells, which is important in orderto learn consistent representations for trajectories drawn fromthe same route. To address this, we propose a novel spatial
proximity aware loss function (in Section IV-C1) and a newcell representation pretraining approach that takes into accountspatial proximity (in Section IV-C2) To further improve thetraining, an approximate loss function based on Noise Contrastive Estimation [15] is also proposed (in Section IV-C1).
1) Spatial proximity aware loss function: To train a sequence encoder-decoder, we need a loss function to characterize the optimization objective. This is important, as differencesin the loss function would encourage the model to learn
different representations [32]. When the sequence encoderdecoder is employed in natural language processing, e.g.,in neural machine translation [23], [24], [33], Negative Log
Likelihood (NLL) loss is chosen to minimize the negative log
likelihood function for tokens in the target sentence as follows,
Y
L1 = − logP(yt |y1:t−1 , x)(4)tmaximizeNY(i)P(Tb |Ta(i) )

(2)

i=1

In the sequence encoder-decoder model, the inputs should
be sequences of discrete tokens. Therefore, we need to finda way to map the continuous coordinates (i.e., longitude,
latitude) into discrete tokens (analogous to words in natural
language). We adopt a simple strategy that is used commonlyin spatial data analytics, i.e., we partition the space into cellsof equal size [31] and treat each cell as a token. All sample
points falling into the same cell are then mapped to the same
token.
The above helps mainly to overcome the problems of nonuniform and low sampling rates. However, the realistic trajectories also may have noisy sample points. For example, whena GPS receiver is in an urban canyon and satellite visibility
is poor, inaccurate locations may result. To eliminate the

However, simply adopting this loss function is problematicfor the spatiotemporal data. Recall that our original purpose
is to maximize P(R|Ta ). Due to the unavailability of theunderlying route R, we use the original trajectory Tb torepresent R and instead maximize P(Tb |Ta ). In practice, even
original trajectories may not cover their underlying route R
well. For example, Figure 3 has two trajectories Tb and Tb0
generated from the same route R (after transforming thecoordinates to cells, their corresponding sequences are y and
y 0 , respectively). The sample points of the two trajectories
interleave the cells in our space partitioning. Let Ta and Ta0
denote the sub-trajectories of Tb and Tb0 , respectively. Ideally,
the representations learned for Ta and Ta0 should be similar, as
they are both generated from route R. The NLL loss function
in Equation 4 differentiates Tb and Tb0 as two identical target
trajectories, so that it cannot discover the similarity between
Ta and Ta0 .

Fig. 3: Tb and Tb0 are two trajectories generated from anunderlying route R. After transforming the coordinates tocells, their corresponding sequences are y, y 0 respectively. The
sample points in the two trajectories interleave on the route.

The reason is that the loss function in Equation 4 penalizes
the output cells with equal weight. Intuitively, the output cells
that are closer to the target should be more acceptable thanthose that are far way. For example, if the decoded target cell
is y3 (in Figure 3), the loss function penalizes the outputs y30
and y1 equally. This is not a good penalty strategy. Since y30
is spatially closer to y3 , it is more acceptable for the decoder
to output y30 rather than to output y1 .
The intuition behind the proposed spatial proximity aware
loss function is that we assign a weight for each cell when wetry to decode a target cell yt from the decoder. The weightof cell u ∈ V is inversely proportional to its spatial distanceto the target cell yt , so the closer the cell is to yt the larger
weight we will assign to it. The spatial proximity aware loss
is given as follows.
|y|L2 = −

XX

wuyt log P

t=1 u∈V

where

exp(Wu> ht ),>v∈V exp(Wv ht )(5)

exp (−||u − yt ||2 /θ)v∈V exp(−||v − yt ||2 /θ)is the spatial proximity weight for cell u when decoding target
yt , and ||u − yt ||2 denotes the Euclidean distance between
the centroid coordinates of the cells. Here θ > 0 is a spatial
distance scale parameter. A small θ penalizes far away cells
heavily, and when θ → 0, the loss function will be reduced to
the NLL loss function in Equation 4. The exponential kernel
function is chosen as it decays fast at the tail, which would
encourage the model to learn to output cells near the target
cell yt .
Although the spatial proximity aware loss function in Equation 5 helps us learn consistent representations for trajectories
generated from the same routes, it requires us to sum over theentire vocabulary twice every time we decode a target yt :
!
X
P
>
>
(6)
Wu ht −
exp(Wv ht )u∈V wuytv∈V{zlog probability

t=1 u∈NK (yt )

v∈N O

(7)
where
exp (−||u − yt ||2 /θ)v∈NK (yt ) exp(−||v − yt ||2 /θ)wuyt = Pwuyt = P|

are close to target cell yt ; ii) it is not necessary to calculate
the exact value of the log probability in Equation 6, as longas we can encourage the decoder to assign the probability
to the cells that are close to the target cell. Based on the firstobservation, we can use just the K nearest cells of yt , denoted
as NK (yt ), instead of using the whole vocabulary in the firstsum in Equation 6. Based on the second observation, we
can use Noise Contrastive Estimation (NCE) [15] to computethe log probability in Equation 6. NCE was developed by
Gutmann et al. [15] to differentiate data from noise by training
a logistic regression. We can use it to approximately maximize
the log probability of cells in NK (yt ) by randomly samplinga small set of cells from V −NK (yt ) as noise data, denoted asO(yt ). In our experiments, we find that 500 randomly sampled
noise cells can give a very good approximation, and thus
the time complexity is reduced from O(|y| × |V |) to O(|y|).
In summary, our approximate spatial proximity aware loss isgiven as follows.
!
|y|
X
X
X
exp(Wv> ht ) ,L3 = −
wuyt Wu> ht −}Thus, the cost of decoding a trajectory y is O(|y|×|V |). Whenthe vocabulary size |V | is large, it will be expensive to train
the model.
Approximate spatial proximity aware loss function. To
reduce the training cost, we design an approximate spatial
proximity aware loss function based on the following two
observations: i) most of wuyt are very small except cells thatN O = NK (yt ) ∪ O(yt )2) Pre-training cell representations: To further guaranteethat the trajectories generated by the same route have close
representations in the latent space, we propose a cell representation learning algorithm to pre-train the cells in theembedding layer of the model. The intuition is that the encoder
squeezes a sequence of cells covered by the trajectory to get
the trajectory representation v, and thus the representations of
two trajectories along the same route will be close in theirlatent space if we can learn similar representations for cells
that are spatially close. For example, if each yi has a similarcell representation as that of yi0 in Figure 3, the trajectory
representations of Ta and Ta0 will be close in the latent spacesince they are encoded by the same encoder.
Two straightforward representations exist for the cells, theone-hot representation [32] and the centroid coordinates of
the cells (GPS coordinates). However, both representations
have limitations. The one-hot representation loses the spatial
distance relation of the cells as all the cells are treatedindependently. As a result, the proposed model may takemore training time to discover spatial relations in the cell
embedding layer. It would help accelerate the training if
the input cell representations provide the prior knowledgeof spatial proximity. Next, the centroid coordinates of thecells naturally encode the spatial proximity for the cells but
restrict the representations in a two-dimensional space, whichmake it difficult for the loss function to further optimize therepresentations in their parameter space.
Based on the above analysis, we propose to feed the
distributed cell representations, which capture the cell spatial

proximity relation, to the embedding layer of the model. We
achieve this by borrowing the key idea from skip-grams [19].
The intuition behind skip-grams is that words with similar
meanings tend to appear together in the same contexts, and if
we use the representation of a word to predict its surroundingwords then we can embed the words into a Euclidean space ina way that captures the semantic distances between the words.
Towards this end, we create the context for a given cell u ∈ V
by randomly sampling its neighbor u0 ∈ NK (u) (we also only
consider its K-nearest neighbors) according to the following
cell sampling distribution:exp(−||u0 − u||2 /θ)v∈NK (u) exp(−||v − u||2 /θ)P(u0 ) = P(8)The cell sampling distribution is similar to the spatial proximity weight in Equation 5. Note that their θ values do not have tobe equal. For each cell u ∈ V , the cell sampling distribution
tends to sample the cells that are spatially close to it as its
context. In this fashion, we are able to create the context foreach cell and to learn the cell representation efficiently with
the negative sampling algorithm [34] by maximizing the log
probability of observing the neighboring cells in its context,
C(u), given cell u:Xmaximizelog P(C(u)|g(u))(9)
u∈VHere, g(u) denotes the representation of cell u. The learned
cell representations will be used to initialize the embedding
layer in the model, but we do not fix their values. Thus
they can still be further optimized by the loss function in
Equation 7. The learning algorithm is shown in Algorithm 1.
Algorithm 1: CellLearning (CL)
Input: The dimension of the learned representations d, context window
size lOutput: The learned cell representations g(u)1 for u ∈ V do
2
C(u) ← ∅ ;3while |C(u)| < l do
4
u0 ∼ P(u0 ) according to Equation 8;
5
C(u) ← C(u) ∪ u0 ;6
7

Optimizing the loss function in Equation 9;
return g(u) for u ∈ V ;D. Complexity of similarity computation
Our model can be trained completely unsupervised with the
SGD (Stochastic Gradient Descent) algorithm. Given a trainedmodel, it only requires O(n) time (as shown in Equation 1,
where f (·, ·) indicates the encoder-RNN and h0 is a zero
vector) to embed a trajectory into a vector which is fast andcan be done using GPUs. Then we can use the Euclidean
distance of the two vectors to measure the similarity of two
trajectories, with a time complexity of O(|v|). Therefore thetime complexity of measuring the similarity between two
trajectories is O(n + |v|).

TABLE II: Dataset statistics.
DatasetPorto
Harbin#Points
74,269,739
184,809,109#Trips
1,233,766
1,527,348

Mean length
60
121

V. E XPERIMENTS
We study the effectiveness and scalability of our proposed
method on two real-world taxi datasets. The experimental
setup and parameter settings are presented in Sections V-A
and V-B respectively. Then we evaluate the accuracy of different methods using most similar search, cross-similarity, and knn queries in Sections V-C1 to V-C3, respectively. Scalability
is covered in Section V-D. The proposed loss functions andcell learning approach are evaluated in Section V-E. We end by
evaluating the impact of the cell size, the hidden state size ofthe encoder, and the training data size on the learned trajectory
representations in Sections V-F and V-G.
A. Experimental setup
Dataset. The experiments are conducted on two real-worldtaxi datasets. The first dataset2 is collected in the city of Porto,
Portugal over 19 months and contains 1.7 million trajectories.
Each taxi reports its location at 15 second intervals. We remove
trajectories with length less than 30, which yields 1.2 million
trajectories. The second dataset contains trajectories collectedfrom 13,000 taxis over 8 months in Harbin, China. We selecttrajectories with length at least 30 and time gaps between consecutive sample points being less than 20 second. This yields
1.5 million trajectories. We partition both sets into training
data and testing data based on the starting timestamp of the
trajectories. For both sets, the first 0.8 million trajectories areused for training, and the remaining trajectories are used fortesting. Statistics of the two sets are shown in Table II.
To create the training trajectory pairs as described in Section IV-A, we perform two kinds of transformations, downsampling and distortion. For each trajectory Tb we first downsample it with a dropping rate r1 varied in [0, 0.2, 0.4, 0.6]
to create its 4 sub-trajectories Ta . And we further distort eachdown-sampled Ta based on a distorting rate r2 (as described in
Equation 3) varied in [0, 0.2, 0.4, 0.6]. As a result, 16 training
pairs (Ta , Tb ) are created for each original trajectory Tb .
Benchmarking Methods: We compare t2vec with threeother methods for measuring the trajectory similarity, namely
EDR [9], LCSS [7], and EDwP [11]. LCSS and EDR are
two of the most widely adopted trajectory similarity measuresin spatiotemporal data analyses. EDwP is the state-of-the-art
method for measuring similarity of non-uniform and low sampling rate trajectories. We do not include DTW as it has been
demonstrated to be consistently inferior to EDR in trajectory
similarity computation [11]. Moreover, we compare with thevanilla RNN (vRNN) [35] and the common set representation(CMS). The vanilla RNN serves as an embedding baseline
method, and the common set representation is used to measure
2 http://www.geolink.pt/ecmlpkdd2015-challengethe similarity of two trajectories based on their common set
after they have been mapped to cells. We discuss the reasonsfor comparing with the two baselines in Section V-C1.
Evaluation Platform: Our method3 is implemented in Julia [36] and PyTorch, and trained using a Tesla K40 GPU.
The baseline methods are written in Java4 . The platform runs
the Ubuntu 14.04 operating system with an Intel Xeon E51620 CPU.

TABLE III: Mean rank versus the database size using the Portoand Harbin datasets.
DB size
EDR
LCSS
CMS
vRNN
EDwP
t2vec

20k
25.73
31.95
62.18
32.73
6.78
2.30

Porto
40k
60k
50.70
76.07
59.20
95.85
112.84173.3461.24
100.20
11.48
16.08
3.45
4.73

80k
104.01130.40231.55135.22
23.02
6.35

100k
130.98150.67291.26163.1028.90
7.67

DB size
EDR
LCSS
CMS
vRNN
EDwP
t2vec

20k
30.37
35.49
97.41
34.30
12.80
5.10

Harbin
40k
60k
57.90
85.72
63.20
105.46141.04209.3765.24
103.05
20.64
29.10
7.50
9.62

80k
118.02
137.20271.45140.25
35.20
12.51

100k
149.01160.67316.81
162.1045.30
15.70

B. Parameter settings and training details
Cell size: The default cell size in the experiments is 100
meters. After removing the cells hit by less than 50 points(i.e., δ = 50), we get 18,866 hot cells for the Porto dataset
and 22,171 hot cells for the Harbin dataset.
RNN units: In our model, GRU [35] with 3 layers is chosen
as the computational unit because it has been shown to be as
good as LSTM [37] in sequence modeling tasks, while it ismuch more efficient to compute [35].
Gradient clipping: Although RNNs tend not suffer from
gradient vanishing problem, they may have exploding gradients [38]. Hence, we clip the gradients by enforcing amaximum gradient norm constraint [30], which is set to 5
in our experiments.
Terminating condition: We randomly select 10,000 trajectories as a validation dataset from the test dataset (the selectedtrajectories are removed from the test dataset). The training
is terminated if the loss in the validation dataset does notdecrease in 20,000 successive iterations.
In addition, both the hidden layer size in GRU and thedimension of the learned cell representation d are set to 256,
the context window size l in the cell learning algorithm is set
to 10. For simplicity, θ in Equations 5 and 8 is fixed at 100(meters). Parameter K and the size of O(yt ) in Section IV-C
are set to 20 and 500, respectively. We adopt Adam stochastic
gradient descent [39] with an initial learning rate of 0.001 to
train the model. We evaluate the training time in Sections V-E
and V-F.
To set the parameter  of the baseline methods EDR and
LCSS, we adopt the strategies described in the studies [7], [9]proposing the two methods; the parameters of vRNN are set
to be the same as our encoder-RNN except that it is trained by
predicting the next cell based on the cells that it has alreadyseen.
C. Performance evaluation
The lack of ground-truth dataset makes it a challenging
problem to evaluate the accuracy of trajectory similarity. Tworecent studies [11], [13] propose to evaluate the accuracy of
methods for computing trajectory similarity using the selfsimilarity and cross-similarity comparisons and assessmentsof the precision of finding the k-nearest neighbors. Currently,
this is the best evaluation methodology, and we adopt thismethodology in the experiments. In addition, we also design
a new experiment, called most similar search (which can beconsidered as a sort of self-similarity measure) to evaluate theeffectiveness of different methods, in Section V-C1.
One of the most important tasks in trajectory analysis issimilar trajectory search. To overcome the lack of groundtruth, we design three experiments to evaluate the performanceusing different methods for this task.
We randomly choose 10,000 trajectories from the test
dataset, denoted as Q, and then we choose another m (aparameter to be evaluated) trajectories, denoted by P . For
each trajectory Tb ∈ Q, we create two sub-trajectories fromit by alternately taking points from it, denoted as Ta and Ta0(see Figure 4), and we use them to construct two datasets
0
DQ = {Ta } and DQ= {Ta0 }. We perform the sametransformation for the trajectories in P to get DP and DP0 .
Then for each query Ta ∈ DQ , we retrieve its top-k most0
∪ DP0 and calculate thesimilar trajectories from database DQ
rank of Ta0 . Ideally Ta0 is ranked at the top since it is generated
from the same original trajectory as Ta . The reason for using00
DQ∪ DP0 as the database instead of DQ∪ P is that the query
trajectory will have similar mean length as the trajectories inthe database5 . Moreover, to evaluate whether RNN encodes
two sequences of cells into two similar vectors simply because
the two sequences have the same starting or ending cells, orjust because they have sufficient numbers of common cells,
we include another two baselines, vRNN and CMS. If theaforementioned reason is true, vRNN and CMS should also
give good performance in the task.

Fig. 4: Creating two sub-trajectories Ta and Ta0 from trajectory
Tb by alternately taking points from it.
1) Most similar trajectory search: Experiment 1 We firststudy the performance of the different methods when we
increase m, the size of P , from 20,000 to 100,000. Table III

3 https://github.com/boathit/research-papers/tree/master/t2vec
4 The authors of EDwP give us access to their compiled jar file.

5 Similar results were found using database D 0 ∪ P .
Q

TABLE IV: Mean rank versus the down-sampling rate r1 using
the Porto and Harbin datasets.
r1
EDR
LCSS
CMS
vRNN
EDwP
t2vec

0.2
160.03168.02296.56173.45
29.10
7.88

Porto
0.3
0.4
208.01235.60173.45
187.60317.70430.00179.58190.24
30.50
31.64
8.00
9.48

0.5
285.10188.40387.90200.13
39.67
12.70

0.6
340.68192.20446.50210.20
61.72
15.99

r1
EDR
LCSS
CMS
vRNN
EDwP
t2vec

0.2
183.02178.80330.70176.45
47.32
15.92

Harbin0.3
0.4
231.01265.50193.50195.30376.20450.04184.83
191.42
49.80
51.39
17.21
19.87

0.5
316.30208.40460.90203.13
57.91
21.74

0.6
380.86210.30476.50250.2081.81
30.95

TABLE V: Mean rank versus the distorting rate r2 using thePorto and Harbin datasets.
r2
EDR
LCSS
CMS
vRNN
EDwP
t2vec

0.2
132.40210.30296.16212.16
30.10
9.10

Porto
0.3
0.4
133.10135.60215.70214.60317.27337.31220
217.30
30.16
32.63
9.20
9.52

0.5
134.90215.05327.90220.61
31.23
9.49

0.6
139.10228.03346.05235.70
33.53
10.80

r2
EDR
LCSS
CMS
vRNN
EDwP
t2vec

0.2
142.44230.32306.62222.41
46.41
16.43

Harbin0.3
0.4
143.69145.67235.93244.65327.87329.31
227.20236.37
47.97
49.32
17.28
17.52

0.5
146.90245.15339.34250.62
50.68
18.51

0.6
152.11251.60349.51245.75
51.10
21.08

shows the mean rank of the 10,000 queries in DQ using
different methods when using the Porto and Harbin datasets.
As the size of P grows, the performance of all methods
degrades. CMS performs the worst, as it ignores the sequential
information in the trajectories. vRNN and LCSS demonstrate
similar performance, and this is reasonable because bothmethods can be considered as enhanced version of CMS that
preserve the order of the points in the sequence. EDwP performs the best among all baseline methods. t2vec outperforms
the other methods significantly, and even when the size of
database P reaches 100,000, it gives a low mean rank for thequeries.
Experiment 2 Next, we study the impact of down-sampling on0
the methods with a fixed database size |DQ∪DP0 | = 100, 000.
We vary the dropping rate r1 from 0.2 to 0.6 and down0∪DP0 based on thesample the trajectories in both DQ and DQ
dropping rate. Table IV depicts the mean rank for the queries in
DQ of the different methods using Porto and Harbin datasets.
EDR degrades quickly when we increase the dropping rate r1 ,while LCSS and vRNN are not very sensitive to variationsin r1 . EDwP shows relatively consistent performance when

TABLE VI: Mean cross-distance deviation for varying dropping rate r1 and distorting rate r2 .
r1
t2vec
EDwP
EDR
r2
t2vec
EDwP
EDR

0.1
0.057
0.059
0.130
0.1
0.010
0.010
0.012

0.2
0.010
0.010
0.190
0.2
0.013
0.018
0.019

0.4
0.016
0.024
0.380
0.4
0.018
0.031
0.033

0.6
0.025
0.039
0.580
0.6
0.021
0.038
0.039

the down-sampling rate r1 varies between 0.2 and 0.5, butwhen raising r1 to 0.6, its mean rank increases markedly. Thisimplies that the linear interpolation may not be able to handle alow sampling rate effectively if the dropping rate is large. The
reason is that the assumption made by the linear interpolationthat the object would move along a straight line between two
consecutive sample points is no longer true. t2vec consistently
outperforms the other methods by a large margin.
Experiment 3 Finally, we evaluate the effect of noise on re0∪ DP0 | = 100, 000 and distort the points
sults. We still fix |DQ0∪ DP0 withof trajectories in both query DQ and database DQ
distorting rate r2 , as described in Equation 3. The results areshown in Table V using the Porto and Harbin datasets. Unlikewith down-sampling, we observe that all methods are not very
sensitive to point distortion. Even if the distorting rate is set to
0.6, no method shows obvious performance degradation. We
observe that t2vec achieves the best performance.
In subsequent experiments, we observes similar results on
the two datasets, and we only report the results when using
the Porto dataset due to the space limitation.
2) Cross-similarity comparison: A good similarity measure
should be able to not only recognize the trajectory variantsof the same underlying route (self-similarity), but should
also preserve the distance between two different trajectories,
regardless of the sampling strategy. We adopt an evaluation
criterion from the literature [13], [18], namely cross distancedeviation, which is calculated as follows:|d(Ta (r), Ta0 (r)) − d(Tb , Tb0 )|,d(Tb , Tb0 )where Tb and Tb0 represent two distinct original rate trajectories, and Ta (r) and Ta0 (r) are their variants obtainedby randomly dropping (or distorting) sample points with the
dropping (or distorting) rate r. We randomly select 10,000
trajectory pairs (Tb , Tb0 ) from the test dataset to calculate
their mean cross distance deviation. A small cross distancedeviation indicates that the evaluated distance is much close
to the ground truth. The mean cross distance deviation of thedifferent methods is described in Table VI, where we vary
the dropping rate r1 and the distorting rate r2 . We notice thatt2vec outperforms the other two methods in terms of crossdistance deviation for different dropping and distorting rates.
3) k-nn queries: The similarity measure is a fundamental
operation that can be used in many applications, e.g., similarity search, clustering, and classification. In this experiment,

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 5: (a)-(c) k-nn results when varying the dropping rate for k = 20, 30, 40. (d)-(f) k-nn results when varying the distorting
rate for k = 20, 30, 40.
we evaluate the performance of different similarity searchmethods. To contend with the lack of ground-truth, we follow the methodology used in previous work [11], [13]. We
first randomly choose 1000 trajectories as query and 10,000
trajectories as the target database from the test dataset. We
apply each method to find the k-nearest-neighbors (k-nn) of
each query trajectory from the target database as its groundtruth. Next, we transform queries and database trajectories
by randomly dropping (resp. distorting) certain sample points
according to the dropping (resp. distorting) rate r1 (resp. r2 ).
Finally, for each transformed query, we find its k-nn from thetarget database using each method and then compare the resultwith the corresponding ground-truth. The rationale behind
this methodology is that a robust distance measure should
adapt to non-uniform and relatively low sampling rates (resp.
distortion) and yield results close to those for relatively highsampling rate (resp. non-distorted) counterparts.
Figure 5 shows the precision (the proportion of true knn trajectories) of the different similarity methods when the
dropping rate (as shown in Figure 5a-5c) and the distorting
rate (as shown in Figure 5d-5f) is varied for k = 20, 30, 40.
The precision of all methods decreases when the dropping
rate or distorting rate increases. EDR and LCSS show similar
performance, but the precision of EDR drops rapidly when
the dropping rate reaches 0.6. EDwP surpasses them by a fairmagnitude, and t2vec consistently performs the best.
D. Scalability
The time complexity of LCSS and EDR for determining
the similarity of two trajectories Ta , Tb is O(|Ta | × |Tb |).
EDwP has complexity O((|Ta | + |Tb |)2 ). These methods rely
on intricate pruning techniques [9], [11] to answer k-nn queries

(a)

(b)

Fig. 6: (a) k-nn query efficiency versus database size using thePorto dataset. (b) k-nn query efficiency versus database size
using the Harbin dataset.
on large datasets. As discussed in Section IV-D, once ourmodel has been trained offline, it takes linear time to encodea trajectory into a vector v, and the encoding process can also
be done offline. As an example, we can encode the 1.7 million
trajectories in the Porto dataset into vectors within 30 minutes
using one Tesla K40 GPU. After encoding the trajectories into
vectors offline, its online complexity is O(|v|). The linear
complexity makes t2vec scale well on large datasets. Notethat model training is done offline; we will study the training
time in Section V-E.
Although it is obvious that our method is much more
efficient in terms of the complexity analysis, we conduct anexperiment to compare the efficiency of t2vec with those of
EDR and EDwP empirically. Figures 6a and 6b show thatthe query time grows with the size of the target databasefor the k-nn query (k = 50) using the Porto and Harbin
dataset, respectively. Although both EDR and EDwP employ

TABLE VII: Mean rank and training time (hours) for themodel equipped with loss functions L1 , L2 , L3 , L3 +CL using
the Porto dataset.
Loss
MR@r1 = 0.4
MR@r1 = 0.5
MR@r1 = 0.6
Time

L1
46.56
55.72
68.49
26

L2
21.34
27.30
32.01
120

L3
9.70
13.50
16.52
22

L3 + CL
9.48
12.70
15.99
14

TABLE VIII: The impact of the cell size on the model using
the Porto dataset.
Cell size#Cells
MR@r1 = 0.5
MR@r1 = 0.6
MR@r2 = 0.5
MR@r2 = 0.6
Time25
60,004216.23234.18291.57302.9137

50
35,335
15.21
19.21
9.49
10.87
25

100
18,866
12.70
15.99
9.49
10.80
14

150
12,425
12.70
16.03
9.51
11.03
8

carefully designed pruning and indexing techniques, t2vec is
at least one order of magnitude faster than both methods.
t2vec offers near-instantaneous response times that support
interactive use and analysis on big trajectory data, such as
trajectory clustering, while the competition cannot. A response
in less that 200 ms is perceived as instantaneous.
E. Evaluation on the loss function
In this experiment we evaluate the effectiveness of the
proposed loss function and the cell representation learning (CLin Algorithm 1) approach on most similar trajectory search byusing the same setting in Section V-C1. The database size is0∪DP0 | = 100, 000. Table VII shows the mean rank
fixed at |DQ(MR) w.r.t dropping rates r1 = 0.4, 0.5, 0.6 and the training
time of different loss functions using the Porto dataset. The L2
loss, is very expensive to compute, and since the model does
not converge even after training for over 5 days (120 hours),
we terminate the training process before it converges. L3 loss
is capable of improving the mean rank significantly whencompared to L1 . The cell representation learning approach
further improves the mean rank and reduces the training timeby 1/3.
F. Effect of the cell size and the hidden layer size
Intuitively, a small cell size provides a higher resolution of
the underlying space, but it also generates more cells (tokens),
which leads to higher training complexity since the model
complexity is linear in the number of tokens [40]. We evaluate
the influence of the cell granularity on the method performancein answering most similar search with r1 = 0.5, r2 = 0.5. As
TABLE IX: The impact of the hidden layer size on the modelusing the Porto dataset.
|v|
MR@r1 = 0.5
MR@r1 = 0.6
MR@r2 = 0.5
MR@r2 = 0.6

64400.01
431.11
390.27397.22

128
50.21
63.71
48.36
50.26

256
12.70
15.99
9.49
10.80

484
10.24
16.70
8.01
9.27

512
11.26
17.42
9.09
10.05

Fig. 7: The impact of training dataset size on the model using
the Porto dataset.
shown in Table VIII, a cell size of 100 gives the best mean
rank. The smallest cell size (25) performs the worst, which
is probably because it has the highest model complexity andthus is much more difficult to train.
Another important parameter that determines the quality of
the learned representation v is the dimension of the hidden
layer in the encoder. A high dimension of the hidden layer
is typically much more expressive, but requires more training
data to avoid overfitting. Table IX summarizes the impact of
the hidden layer size on the most similar trajectory search withthe same settings as in the above experiment. It is obviousthat increasing |v| from 64 to 256 significantly enhances thequality of the learned representations, while the performancedrops when we increase it further.
G. Effect of the training data size
In this experiment, we evaluate the effect of training data
size on the accuracy of the trajectory similarity search. The
setting is similar to the one in Section V-E, with exceptionthat we vary the training data size and fix the dropping rate
r1 at 0.6.
Figure 7 shows the effect of the training data size on themost similar search on the Porto dataset. The mean rank drops
rapidly as we increase the training data size from 200,000 to
600,000, and the decrease slows down when we continue toenlarge the training data size. When we increase the training
data size from 200,000 to 600,000, more transition patterns
can be learned for representing the trajectories, and thus themodel quality is enhanced significantly. However, when we
further increase the training data size, the marginal benefit ofusing larger training data is less pronounced.
VI. C ONCLUSIONS AND F UTURE W ORK
To the best of our knowledge, we present the first solution
for using learning for creating representations of trajectories
that are well suited for use in trajectory similarity computation.
Most existing methods for trajectory similarity computation
measure similarity by means of pairwise point-matching,which may not capture the underlying route information of
the trajectories, thus being sensitive to non-uniform and low
sampling rates, and noise sample points. Moreover, all ofthem have a quadratic time complexity. Motivated by theseobservations, we propose a seq2seq-based method to learn
representations for trajectories that enable accurate and efficient trajectory similarity search that is robust to sampling ratevariations and noisy sample points. The method is evaluatedempirically with favorable results in terms of both accuracy and efficiency. It consistently outperforms the baseline
methods by a large margin in the most similarity tasks. The
method is at least one order of magnitude faster than theother methods, thus it can support big trajectory analysis andinteractive use while the competition cannot.
This work sheds light on several new research directions:
1) Employing the learned representations to explore more
downstream tasks, e.g., trajectory clustering and popularroutes search. 2) Extending the proposed method to moregeneral time series data beyond trajectories. 3) Developing
indexing techniques like Locality-Sensitive Hashing [41] tofurther speed up the proposed method.
Acknowledgments. This work is supported in part by the
Rapid-Rich Object Search (ROSE) Lab at Nanyang Technological University. The ROSE Lab is supported by the National
Research Foundation, Prime Minister’s Office, Singapore, under its IDM Futures Funding Initiative, administered by theInteractive and Digital Media Programme Office. This work
was also supported by the MOE Tier-2 grant MOE2016-T21-137, MOE Tier-1 grant RG31/17, NSFC under the grant
61772537, and a grant from Microsoft. C. S. Jensen was
supported by the DiCyPS project and by a grant from the Obel
Family Foundation. W. Wei was supported by NSFC under the
grant 61602197, NSF under the grant 2016CFB192.
R EFERENCES[1] Z. Li, J. Han, M. Ji, L. A. Tang, Y. Yu, B. Ding, J. Lee, and
R. Kays, “Movemine: Mining moving object data for discovery ofanimal movement patterns,” ACM TIST, vol. 2, no. 4, pp. 37:1–37:32,
2011.
[2] Z. Chen, H. T. Shen, and X. Zhou, “Discovering popular routes fromtrajectories,” in ICDE, 2011, pp. 900–911.
[3] C.-C. Hung, W.-C. Peng, and W.-C. Lee, “Clustering and aggregatingclues of trajectories for mining trajectory patterns and routes,” VLDBJ,
vol. 24, no. 2, pp. 169–192, 2015.
[4] H. Jeung, M. L. Yiu, X. Zhou, C. S. Jensen, and H. T. Shen, “Discoveryof convoys in trajectory databases,” PVLDB, vol. 1, no. 1, pp. 1068–
1080, 2008.
[5] X. Li, V. Ceikute, C. S. Jensen, and K.-L. Tan, “Effective online group
discovery in trajectory databases,” IEEE TKDE, vol. 25, no. 12, pp.
2752–2766, 2013.
[6] B.-K. Yi, H. Jagadish, and C. Faloutsos, “Efficient retrieval of similar
time sequences under time warping,” in ICDE, 1998, pp. 201–208.
[7] M. Vlachos, G. Kollios, and D. Gunopulos, “Discovering similar multidimensional trajectories,” in ICDE, 2002, pp. 673–684.
[8] L. Chen and R. Ng, “On the marriage of lp-norms and edit distance,”in PVLDB, 2004, pp. 792–803.
[9] L. Chen, M. T. Özsu, and V. Oria, “Robust and fast similarity search
for moving object trajectories,” in SIGMOD, 2005, pp. 491–502.
[10] K. Zheng, Y. Zheng, X. Xie, and X. Zhou, “Reducing uncertainty of
low-sampling-rate trajectories,” in ICDE, 2012, pp. 1144–1155.
[11] S. Ranu, P. Deepak, A. D. Telang, P. Deshpande, and S. Raghavan,“Indexing and matching trajectories under inconsistent sampling rates,”in ICDE, 2015, pp. 999–1010.
[12] P. Banerjee, S. Ranu, and S. Raghavan, “Inferring uncertain trajectories
from partial observations,” in ICDM, 2014, pp. 30–39.

[13] H. Su, K. Zheng, H. Wang, J. Huang, and X. Zhou, “Calibrating
trajectory data for similarity-based analysis,” in SIGMOD, 2013, pp.
833–844.
[14] A. Graves, S. Fernández, F. Gomez, and J. Schmidhuber, “Connectionist temporal classification: labelling unsegmented sequence data withrecurrent neural networks,” in ICML, 2006, pp. 369–376.
[15] M. Gutmann and A. Hyvärinen, “Noise-contrastive estimation: A new
estimation principle for unnormalized statistical models,” in AISTATS,
2010, pp. 297–304.
[16] E. Frentzos, K. Gratsias, and Y. Theodoridis, “Index-based most similar
trajectory search,” in ICDE, 2007, pp. 816–825.
[17] S. Sankararaman, P. K. Agarwal, T. Mølhave, J. Pan, and A. P.
Boedihardjo, “Model-driven matching and segmentation of trajectories,”in SIGSPATIAL, 2013, pp. 234–243.
[18] H. Wang, H. Su, K. Zheng, S. Sadiq, and X. Zhou, “An effectivenessstudy on trajectory similarity measures,” in Australian Database Conference, 2013, pp. 13–22.
[19] T. Mikolov, K. Chen, G. Corrado, and J. Dean, “Efficient estimation ofword representations in vector space,” arXiv preprint arXiv:1301.3781,
2013.
[20] Q. V. Le and T. Mikolov, “Distributed representations of sentences anddocuments.” in ICML, 2014, pp. 1188–1196.
[21] B. Perozzi, R. Al-Rfou, and S. Skiena, “Deepwalk: Online learning ofsocial representations,” in SIGKDD, 2014, pp. 701–710.
[22] J. Tang, M. Qu, M. Wang, M. Zhang, J. Yan, and Q. Mei, “Line: Largescale information network embedding,” in WWW, 2015, pp. 1067–1077.
[23] K. Cho, B. Van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares,
H. Schwenk, and Y. Bengio, “Learning phrase representations using
rnn encoder-decoder for statistical machine translation,” arXiv preprint
arXiv:1406.1078, 2014.
[24] I. Sutskever, O. Vinyals, and Q. V. Le, “Sequence to sequence learningwith neural networks,” in NIPS, 2014, pp. 3104–3112.
[25] O. Vinyals, Ł. Kaiser, T. Koo, S. Petrov, I. Sutskever, and G. Hinton,“Grammar as a foreign language,” in NIPS, 2015, pp. 2773–2781.
[26] R. Kiros, Y. Zhu, R. R. Salakhutdinov, R. Zemel, R. Urtasun, A. Torralba, and S. Fidler, “Skip-thought vectors,” in NIPS, 2015, pp. 3294–
3302.
[27] D. Williams and G. Hinton, “Learning representations by backpropagating errors,” Nature, vol. 323, no. 6088, pp. 533–538, 1986.
[28] P. J. Werbos, “Backpropagation through time: what it does and how to
do it,” Proceedings of the IEEE, vol. 78, no. 10, pp. 1550–1560, 1990.
[29] Y. Bengio, R. Ducharme, P. Vincent, and C. Jauvin, “A neural probabilistic language model,” JMLR, vol. 3, pp. 1137–1155, 2003.
[30] A. Graves, “Generating sequences with recurrent neural networks,” arXiv
preprint arXiv:1308.0850, 2013.
[31] R. H. Güting and M. Schneider, “Realm-based spatial data types: therose algebra,” VLDBJ, vol. 4, no. 2, pp. 243–286, 1995.
[32] Y. Bengio, A. Courville, and P. Vincent, “Representation learning: A
review and new perspectives,” TPAMI, vol. 35, no. 8, pp. 1798–1828,
2013.
[33] D. Bahdanau, K. Cho, and Y. Bengio, “Neural machine translation by
jointly learning to align and translate,” arXiv preprint arXiv:1409.0473,
2014.
[34] T. Mikolov, I. Sutskever, K. Chen, G. S. Corrado, and J. Dean,“Distributed representations of words and phrases and their compositionality,” in NIPS, 2013, pp. 3111–3119.
[35] J. Chung, C. Gulcehre, K. Cho, and Y. Bengio, “Empirical evaluation of
gated recurrent neural networks on sequence modeling,” arXiv preprint
arXiv:1412.3555, 2014.
[36] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh
approach to numerical computing,” SIAM review, vol. 59, no. 1, pp.
65–98, 2017.
[37] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neuralcomputation, vol. 9, no. 8, pp. 1735–1780, 1997.
[38] R. Pascanu, T. Mikolov, and Y. Bengio, “On the difficulty of training
recurrent neural networks,” in ICML (3), 2013, pp. 1310–1318.
[39] D. Kingma and J. Ba, “Adam: A method for stochastic optimization,”
arXiv preprint arXiv:1412.6980, 2014.
[40] S. Jean, K. Cho, R. Memisevic, and Y. Bengio, “On using verylarge target vocabulary for neural machine translation,” arXiv preprint
arXiv:1412.2007, 2014.
[41] A. Andoni and P. Indyk, “Near-optimal hashing algorithms for approximate nearest neighbor in high dimensions,” in FOCS, 2006, pp. 459–468.