Truncating Shortest Path Search for Efficient Map-Matching

1 Introduction

The Global Positioning System (GPS) is heavily used in devices such as cell phones and car navigation systems, and a huge amount of GPS traces, or sequences of GPS points, are being recorded every day. For example, OpenStreetMap provides public GPS traces 1 , and some city governments such as Dublin 2 and Rio de Janeiro 3 provide GPS traces of buses to public. Moreover, as Internet of Things emerges, more devices are to come with GPS sensors and be connected to the Internet. This gives us the opportunities to leverage the knowledge that can be extracted from such a huge amount of GPS traces. For example, one might analyze the traffic conditions such as congestion from the GPS traces.

However, the GPS traces often involve noise and are sparse (interval between consecutive points are long and irregular), and it is nontrivial even to extract such basic information as the route that the GPS device has traversed. This task of determining the route corresponding to a given GPS trace on a

1 http://www.openstreetmap.org/traces

2 http://www.dublinked.ie/

3 http://data.rio/

Takayuki Osogami and Rudy Raymond

IBM Research - Tokyo 19-21 Nihonbashi Hakozaki-cho, Chuo-ku, Tokyo, 103-8510, Japan ¶osogami,rudyhar♦@jp.ibm.com

road network is known as map-matching and has been studied extensively in the literature. Figure 1 shows examples of GPS traces from OpenStreetMap in Rio de Janeiro, Brazil. Here, we can see that the lines connecting the GPS traces (red line segments) are not always on the road in the map. The left panel in Figure 2 illustrates the discrepancy between a road network (black line segments) and a GPS trace (red points), where blue dashed line segments connect consecutive points in the GPS trace. The right panel in Figure 2 illustrates the sequence of road segments (green line segments) that actually corresponds to the GPS trace in the left. Observe that, connecting road segments closest to the GPS points, we would obtain an unrealistic winding route, which is drastically different from the actual route shown in the right panel.

Recently, some popular algorithms for map-matching are based on Hidden Markov Models (HMMs) [Newson and Krumm, 2009; Osogami and Raymond, 2013]. Here, a hidden state corresponds to (a point on) a road segment, and an observed variable corresponds to a GPS point. Map-matching with HMMs is known to achieve high accuracy and be robust against the noise and sparsity in the GPS traces. Unfortunately, these algorithms are computationally expensive and not suitable for dealing with a large amount of GPS traces.

We observe, from preliminary experiments, that the bottleneck of existing algorithms for map-matching with HMMs is in the calculation of transition probabilities. Specifically, the transition probability from a hidden state to another is calculated based on the length of the shortest path between the two states on a road network. Notice that the length of a path in the model can take into account various measures such as the number of turns and the expected travel time in addition to the length of traversed roads. If we consider k hidden states for each GPS point, the existing algorithms invoke shortest path search for k2 times at every step of the Viterbi algorithm. Our goal is to speed up the map-matching with HMMs in order to enable highly accurate map-matching with HMMs for a huge amount of GPS traces. To this end, we have two contributions.

First, we propose a technique to reduce the computational time needed for the shortest path search in HMM-based mapmatching. Our technique is to run one-to-many shortest path search on a reversed network, where the direction of every directed arc is reversed, and truncate this shortest path search early without exploring all the shortest paths. We will see that this truncation produces no loss of accuracy.

Our second contributions are in numerical experiments, which show that the proposed approach can drastically reduce computational time without losing accuracy. Specifically, our proposed approach reduces the computational time by a factor of at least 5.4 for 24 GPS traces (1.18 million GPS points) in the T-drive data set [Yuan et al., 2010; 2011].

1.1 RelatedWork

Existing methods of map-matching can be grouped into four categories [Quddus et al., 2007; Zheng, 2015]: geometrical, topological, probabilistic, and advanced methods that incorporate the advantages from multiple categories, for example, an HMM-based approach [Lamb and Thi´ebaux, 1999; Hummel, 2006; Newson and Krumm, 2009; Osogami and Raymond, 2013]. The geometrical approach focuses only on the shapes of the routes. The topological approach considers the topological aspects such as the connectivity of the road network. The probabilistic approach defines the confidence region around the GPS points and looks for the best possible route. An HMM-based approach takes advantage of a probabilistic model of measurement error, similar to Kalman filters [Kim et al., 2000] and Particle filters [Pink and Hummel, 2008]. It thus has the advantages of probabilistic methods, which have been successfully applied for recovering lanes of roads [Sehestedt et al., 2007] and for localizing autonomous vehicles [Thrun et al., 2001]. Unlike purely probabilistic methods, an HMM-based approach also takes into account the topological information of the road network to achieve high accuracy, similar to topological methods [Brakatsoulas et al., 2005]. HMM-based approaches can also run in real-time [Goh et al., 2012] with guaranteed accuracy using an online decoding algorithm [Wang and Zimmermann, 2014], unlike geometrical methods [White et al., 2000; Yanagisawa, 2010], which often suffer from large computational cost.

Our work builds upon the recent work on an HMM-based method [Osogami and Raymond, 2013], which extends the method by Newson and Krumm [2009] by taking into account drivers' preferences about routes. Although the running time of the HMM-based methods is usually not an issue for a single GPS trace, it becomes unacceptable for many GPS traces. The bottleneck is in the computational cost for calculating the shortest paths between many hidden states.

Despite the long history of research on the shortest path problem, only the standard algorithm of one-to-one shortest path has been used in the existing HMM-based mapmatching. We use one-to-many shortest path search and investigate when we can truncate the search particularly in the context of HMM-based map-matching. One might further improve the efficiency of HMM-based map-matching by exploiting sophisticated techniques of shortest path search, including Floyd-Warshall [Floyd, 1962; Warshall, 1962] for multi-source and multi-destination as well as advanced heuristics such as A* [Hart et al., 1968].

Another approach to accelerate shortest path search is preprocessing [Goldberg and Harrelson, 2005; Geisberger et al., 2008; Sturtevant et al., 2009]. Preprocessing produces the auxiliary data to accelerate the shortest path search. Many of preprocessing approaches rely on a particular fixed metric of the shortest path, or can only deal with small changes to the metric. Those preprocessing approaches need to recalculate the auxiliary data whenever we change the metric. To avoid recalculation, one can apply an advanced preprocessing approach by Delling et al. [2015] that essentially produces auxiliary data supporting multiple metrics. This preprocess partitions the road network into connected cells, generates an overlay graph connecting the cells, and runs shortest path search on the overlay graph.

However, we do not adopt the preprocessing approaches in our map-matching because of the following reasons. Firstly, we use the weighted sum of the travel distance and the turn cost, where the weight can vary drastically for each GPS trace, depending on who traveled along that trace. Secondly, consecutive GPS points on a trajectory of map matching that we consider are likely to be located closely because their intervals are usually short (up to several hundreds seconds). In this case, because the source and destination points are often close, the preprocessing by Delling et al. [2015] is likely to perform a shortest path search in the same cell, or at most between a pair of neighboring cells. Thus, the shortest path search with preprocessing is essentially the same as a shortest path search on the original road network. Our proposed method exploits techniques for fast shortest path computations in the HMM context similar to that of Felzenszwalb et al. [2004]. Our approach can also accelerate HMMs for highvolume data processing beyond map matching applications, as long as the transition probabilities are determined by some distance metrics on a directed graph whose nodes represent hidden states.

2 Map-matching Algorithm

We first introduce the HMM for map-matching briefly and then explain the details of the novel technique to reduce the computational time for calculating transition probabilities. We suppose that a set of GPS traces is given as input. The points in each GPS trace are ordered according to the time they are observed, but the exact time is not necessarily recorded. Our goal is to output a set of sequences of road segments corresponding to the GPS traces.

2.1 Map-matching with a HiddenMarkovModel

An HMM consists of the observable variables, the hidden states, the initial state probability, the emission probability, and the transition probability. For map-matching, an observable variable is the GPS point, i.e., the longitude and latitude. In the following, we describe the other components of the HMM for map-matching, except the initial probability, which is given analogously to the emission probability. The approach described in this section closely follows Osogami and Raymond [2013].

A hidden state corresponds to a road segment of a road network. Here, we define the location of a hidden state by the midpoint of the corresponding road segment where we assume that road segments are straight line segments. We adopt a simpler form of a hidden state than the model of Osogami and Raymond [2013], where an arbitrary point on a road segment can become a hidden state. A GPS point, zt at time t, is associated with k hidden states, where k is a parameter to be fixed. The k hidden states correspond to the k nearest road segments to zt.

The emission probability of a GPS point follows from a model of measurement error with the GPS. Let z be the location of the GPS point and s be the location of the hidden state, namely the midpoint of the road segment corresponding to the hidden state, and σ be the standard deviation of the measurement error. Then the emission probability follows the normal distribution whose probability density function is given by

where k · k denotes a distance measure such as the great circle distance or the Euclidean distance.

The transition probability from a hidden state s to another s0 is given by the probability density function of an exponential distribution. Let d⇤ (s↪ s0 ) be the minimum value of a weighted travel cost from s to s0 , namely, the shortest path from s to s0 . Following Osogami and Raymond [2013], we use the weighted sum of the length of traversed road segments and the cost of turns for d⇤ (s↪ s0 ). To calculate the travel cost, we construct a network G = (V↪E) whose nodes V correspond to road segments, and edges E connect adjacent road segments. Each edge is weighted by the corresponding travel cost, as is defined in the following. Let v and v0 be adjacent road segments (v0 follows v) and dv and dv0 be the length of the road segments. The cost u of turn is defined as follows.

where ✓v↪v0 is the angle between the road segments v and v0 . We add an edge e 2 E connecting v and v0 on the network G and define the cost ce of e by

where wturn is a parameter. Note that (2) contains costs with multiple metrics. Osogami and Raymond [2013] used Inverse Reinforcement Learning [Ziebart et al., 2008] to estimate the weight wturn. The optimal path that minimizes the weighted cost between a source node and a destination node can be found by applying a shortest path algorithm such as the Dijkstra algorithm [Dijkstra, 1959].

The transition probability from s to s0 is given, with a parameter β, by

where ks-s0 k denotes the distance between the midpoints of the road segments that correspond to the hidden states s and s0 . Note that the following transition probability is used by Osogami and Raymond [2013]:

We omit the absolute value, because the shortest path distance d⇤ (s↪ s0 ) is longer than the distance ks - s0 k in most cases, and our condition for truncating shortest path search, to be explained later, is suitable for the form (3).

Having all HMM elements, we compute the most likely sequence of hidden states by the Viterbi algorithm [Viterbi, 1967]. Let z1↪ ▷ ▷ ▷ ↪ zn be the sequence of GPS points and St = ¶st↪1↪ ▷ ▷ ▷ ↪ st↪k♦ be the set of hidden states at time t, where k is the parameter denoting the number of the hidden states. The Viterbi algorithm recursively computes the maximum log likelihood LL⇤ (s) for each hidden state, s:

2.2 Truncating Shortest Path Search

We observe, from preliminary experiments, that the bottleneck of the existing algorithms for map-matching with HMMs is in the calculation of transition probabilities whose values depend on d⇤ (s↪ s0 ), the shortest paths between two consecutive hidden states (road segments). Computation of these shortest paths is time consuming for long and irregular interval of hidden states. The existing algorithms invoke shortest path search for k2 times, one for each pair of hidden states, at every step of the Viterbi algorithm.

Instead of many runs of one-to-one shortest path search, we consolidate k2 runs of one-to-one shortest path search from all s 2 St-1 to all s0 2 St into k runs of one-to-many shortest path search from all s0 2 St to all s 2 St-1 on a reversed network, where the travel cost from a node s0 to another s is the travel cost from s to s0 on the original network. The duplicated one-to-one traversal is avoided by the one-to-many traversal, and the number of the shortest path searches at every step in Viterbi algorithm is reduced from O(k2 ) to O(k). We call it baseline algorithm and compare it with our technique. Note the reversed network is necessary for our technique.

We propose a technique to truncate the one-to-many shortest path search in map-matching without loss of accuracy of the Viterbi algorithm. Consider a case where there is a pair of a source node and a destination node that are not reachable or very far. The shortest path search takes long time in such a case, because it traverses all the nodes on the road network. Note that the transition probability (3) for a long distant pair of nodes is so small that we can safely truncate its shortest path search. For example, see Figure 3, where there are two highways, (s1↪ s2↪ s3) and (s4↪ s5), and three GPS points, (zt-2↪ zt-1↪ zt). The first highway lies above the second one, but they are not connected on the map under consideration. Let s1 be a hidden state of zt-2, and s2↪ s4↪ s5 be the hidden states of zt-1. In the computation of the transition probability, the shortest path search would attempt to compute d⇤ (s1↪ s4) and d⇤ (s1↪ s5) by traversing all nodes on the road network, even though the s4 and s5 are not (directly) reachable from s1, and the corresponding transition probabilities are negligible.

In order to truncate such redundant shortest path search, we carefully study the log likelihood (5) of the Viterbi algorithm. At one stage of (5), the Viterbi algorithm has found LL⇤ (st-1↪j) for all j and seeks to find LL⇤ (st↪i) for an i. To do so, we compute the shortest path from each of st-1↪j for 1  j  k to st↪i. This many-to-one shortest paths can be computed by finding one-to-many shortest paths on the reversed network.

We then show how to exactly, i.e., without losing accuracy, truncate the search for one-to-many shortest paths. This exact truncation is made possible by exploiting the property of the Viterbi algorithm on the reversed network. Let

for s 2 St-1. Plugging (6) into (5), we obtain

We have already computed LL⇤ (st-1↪j) for every j in the t - 1-st iteration of the Viterbi algorithm. In the beginning of the t-th iteration of the Viterbi algorithm, we compute and store the value of Gt↪i(st-1↪j) for all j. In (7), a j cannot become the maximizer if there exists an ` such that

Thus, we can stop searching for the shortest path from st↪i to st-1↪j on the reversed network if the following inequality holds for an ` that has been searched in previous steps:

The procedure for the exact truncation to compute (7) is as follows. It invokes shortest path search from st↪i to St-1 = ¶st-1↪i ♣ 1  i  k♦ on the reversed network by the Dijkstra algorithm. Starting from the source node st↪i, the algorithm extends the search gradually and determines the shortest path distance. Let W be a subset of all nodes V whose shortest path distances from st↪i have been determined and, let F ⌘ St-1 \ W. Every time the algorithm determines the shortest path distance to a new node u, which is not necessarily in St-1, it updates W and F and checks the following inequality (9) to see if we can truncate the shortest path search:

The inequality (8) holds if inequality (9) holds because the undetermined nodes St-1 ∩ F are further than u. Note that ` = argmax¶Gt↪i(s) - d⇤ (s↪ st↪i) ♣ s 2 F♦ in (8). If the inequality (9) holds, the procedure truncates the Dijkstra algorithm and computes (7) using ` as the maximizer.

3 Experimental Results

We conduct two types of experiments. The first is to evaluate the quality of the results of our map-matching, while the second is to compare the computational time against the baseline algorithm. We implemented our algorithm in Python as a single threaded program and run it on PyPy runtime version 2.6.1 on a PC with 3.3 GHz Intel Xeon E5-2643 CPU.

3.1 Comparison of Accuracy

We apply our map-matching algorithm to a benchmark set that is used and made publicly available by Newson and Krumm [2009]. It contains GPS points for a route of approximately 80 km in Seattle and the ground truth of the traversed route, or simply, the true route. We compare the accuracy of the outcomes by our algorithm against that by Osogami and Raymond [2013] 4 . We change the sampling interval of GPS trace data from 30 seconds to 600 seconds. Throughout we fix σ = 10▷0 meters and β = 0▷1 meters in (1) and (3) and we use k = 5 hidden states at every time step in accordance with Osogami and Raymond [2013].

In order to measure the accuracy of results, we use the route mismatched fraction (error) introduced by Newson and Krumm [2009]. Roughly speaking, the error of a route is defined by the ratio of the total length of the erroneously added or subtracted road segments and the length of the correct route. We also calculate the precision and recall of the resulting route. Let R be a set of road segments of a resulting route and T be a set of road segments of the true route. The precision and recall of the route are defined by ♣R\ T♣◁♣R♣↪ ♣R\ T♣◁♣T♣, respectively.

We divide the data set into 12 data sets in accordance with Osogami and Raymond [2013] and use the wturn inferred by Osogami and Raymond [2013], which ranges from 117.1 to 2,257. We plot the average values of the metrics of the error, precision and recall in Figure 4 among the 12 experiments, where the horizontal axis denotes the sampling interval. The lines labeled with 'Ours' and 'OR' represent the average values of the corresponding metrics by our algorithm and Osogami and Raymond [2013], respectively.

Overall, the accuracy of our algorithm is comparable to that of Osogami and Raymond [2013]. This first experiment confirms the soundness of the implementation of our algorithm. Recall that the only essential difference between our algorithm and that of Osogami and Raymond [2013] is in the way transition probabilities (or shortest paths) are calculated, and our truncation of shortest path search does not lose accuracy.

We also apply our map-matching without the truncation of shortest path search and observe that the outputs with and without the truncation are equivalent. It demonstrates that the truncation does not lose accuracy of the Viterbi algorithm as we discussed in Section 2.2 .

3.2 Comparison of Speed

We now apply our algorithm and the baseline algorithm to a big data set to compare their computational time. Recall that the baseline performs one-to-many shortest path search, without truncation, to calculate the transition probabilities.

We use the first 24 GPS traces from the T-drive data set [Yuan et al., 2010; 2011], which contains a large amount of GPS trace data in Beijing but not the true routes. In our experiment, we remove the GPS points with the following properties. First, we delete the GPS points with the exact same coordinates as the previous time step. Second, we eliminate the GPS points which are located too far (more than 100 km) from ones at the previous time step. Third, we delete the GPS points whose estimated speed from ones at the previous time step is too fast (more than 200 km/h). The total number of GPS points in the resulting data is approximately 1.18 million. We extract the road network data of Beijing from OpenStreetMap, which contains approximately 301,000 nodes and 461,000 edges.

4 We also performed the comparison against the results of Newson and Krumm [2009] and obtained the results similar to those reported in Osogami and Raymond [2013].

We apply our algorithm and the baseline algorithm to the data set and compare the computation time. Note that we do not compare the precision as in the previous subsection, because the data set does not include the ground truth data. We fix σ = 10▷0 meters, β = 0▷1 meters and wturn = 1000▷0 in (1), (3) and (2), respectively.

Table 1 shows the total computational time of the baseline algorithm and our algorithm with k = 5 and 10. We set the time limit for each GPS trace by 10,000 seconds. There are 4 GPS traces (among 24 traces) in which the baseline with k = 10 exceeded the time limit. Our algorithm finished the computation within the time limit for all 24 GPS traces in the two cases. We observe that our algorithm is 5.4 times faster than the baseline with k = 5 and at least 6.4 times faster with k = 10. See Figure 5 for the computation time for each GPS trace. Our algorithm runs faster than the baseline by a factor of 1.8 to 20.5 (or even faster), depending on the trace and k.

The advantages of the proposed algorithm over the baseline can increase with a larger value of k. While k = 5 appears to be sufficient for the data set of Newson and Krumm [2009], a larger k can improve the accuracy for other data sets. Our algorithm indeed finds different routes with a larger k for the T-drive data set, although the changes in the accuracy cannot be evaluated without true routes.

We additionally test an approach to reuse the intermediate data structure of the shortest path search such as the shortest path tree and the priority queue of nodes. Because we deal with geospatial temporal GPS data sets, we expect the locality property to hold: namely, recently accessed hidden states and its neighbors are likely to be accessed soon. Therefore, storing the intermediate data structure in a Least Recently Used (LRU) cache, which is fast and easy to implement, we can avoid duplicate traversal in the shortest path search. "baseline + LRU" and "ours + LRU" denote the computation time of the baseline algorithm with an LRU cache and our algorithm with an LRU cache, respectively. We observe that the LRU cache accelerated our algorithm and the baseline algorithm by a factor of 3.7 and 5.2, respectively, with k = 5. Our algorithm with the cache runs 3.9 times faster than the baseline algorithm with the cache with k = 5.

Finally, we compare the outputs of "baseline + LRU" and "ours" with k = 5 and 10 in Figure 5 and observe that they are equivalent for all traces. The difference between the algorithms is the LRU and the truncation of shortest path search, and it is clear that the LRU does not affect outputs of mapmatching. Thus, the observation confirms that the truncation does not lose accuracy of the Viterbi algorithm. Note that we use "baseline + LRU" instead of "baseline" because "baseline" exceeded the time limit 10,000 seconds for trace ID 2237, 3105, 3557, and 6275 with k = 10.

4 Conclusion

We propose a new technique to speed up HMM-based mapmatching algorithms. Our technique computes the transition probabilities in the HMM by applying one-to-many shortest path search on the reversed network and truncate the shortest path search early without exploring all the shortest paths by taking into account the log likelihood of the Viterbi algorithm. It was shown effective to significantly reduce the computational time without loss of accuracy. We observe that the proposed method runs 5.4 times faster than the baseline algorithm through the computational experiments with k = 5 and even faster with larger k.

We used the proposed method in the framework of Osogami and Raymond [2013], but the proposed method can also be applied to a larger set of HMM-based map-matching algorithms. For example, Osogami and Raymond [2013] define the cost of the path only from the length of the traversed road segments and the number of turns, but one could also take into account other distance measures such as the travel time, type of turns (as some drivers might prefer less number of left turns than right ones), distance to major landmarks, and so on. This is important in modeling the real-world driving, especially in the context of imitation learning, where the utility function of the drivers to be recovered often depends on many features or distance measures. In future work, we plan to evaluate the proposed model with such rich features against large data sets that come with the ground truth data, such as, bus trajectories that are usually fixed.

Acknowledgments

A part of this research was supported by CREST, JST.