g-RIPS Sendai 2022
Mitsubishi-A Group

Final Report

Authors:
Tomoya Akamatsu1
Gabriel Gress2
Katelynn Huneycutt3
Seiya Omura 4

Mentors:
Dr. Shunsuke Kano+
Dr. Masashi Yamazaki*

1

Osaka University
University of New Mexico
3
Ohio State University
4
Nagoya University+
Academic Mentor, Tohoku University,
RACMaS*
Industrial Mentor, MITSUBISHIElectric Corp
2

August 17, 2022

Page 1 of 34

g-RIPS Sendai 2022, Mitsubishi-A

Contents
1 Introduction

2

2 Problem Summary and Statement2

3 Background
3.1 Geometric Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Hidden Markov Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Application to Map-Matching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
4
5

4 Mathematical formulation for our Approach

5

5 Wasserstein Distance Method
5.1 A brief introduction of (discrete) optimal transport theory . . . . . . . . . . . . . . . . . . . .
5.2 The case in which input data is only the trajectory coordinates and timestamps . . . . . . . .
5.3 Incorporating speed and direction information . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 Our motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.2 Our strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.3 The definition of weights including speed and direction information . . . . . . . . . . .
5.3.4 Briefly explanation of route determination in Figure 5.13 using Wasserstein method .
5.4 Future problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7
8
10
10
10
11
12
14

6 Methods by Physical Analogy
6.1 “Electric” Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1.1 Details of method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1.2 Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 “Harmonic oscillator” method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.1 Details of method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.2 Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15
15
15
16
17
18
19

7 Implementation
7.1 Preprocessing Routes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.1 Restricting to a Sub-graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.2 Obtaining Candidate Routes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.3 Candidate Route Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 Wasserstein Distance Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3 Electrical Method and Harmonic Oscillator Implementation . . . . . . . . . . . . . . . . . . .
7.4 Evaluation (Error) Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.5 Preliminary Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.6 Computational Complexity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.7 Metric-based Algorithm Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.8 Data Fusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.9 Other Utilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20
21
21
22
22
22
23
24
25
25
26
26
28

8 Future Work (Implementation)
8.1 IMU Inclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.2 Incomplete Map Matching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.3 Extension to Higher Dimensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28
29
29
29

9 Conclusion

29

Appendix A Geometric background of Wasserstein method
A.1 Ricci curvature of graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.2 “Ricci curvature” between the trajectory and route . . . . . . . . . . . . . . . . . . . . . . . .

29
29
31

1

Page 2 of 34

g-RIPS Sendai 2022, Mitsubishi-A

1

Introduction

A fundamental question when dealing with geospatial information is, given GPS trajectory data and aroad map, how can one determine which route on a map this trajectory corresponds to? This problem iscalled the map-matching problem. Data-driven methods, such as hidden Markov models, are currently themost prominent approach to the map-matching problem; however, they are inflexible changes in the data. In
contrast, a mathematical model will have the benefit of being more adaptable to changes in the data, such as
the extension from two to three-dimensional data or the inclusion of speed and angle data. Our project aims
to develop new solutions to the map-matching problem, favoring mathematical formulations, rather than adata-driven models, that are stable under small perturbations in the road map and trajectory data. We
propose three novel methods: Wasserstein distance method, electrical force method, and harmonic oscillator
method. Each of these methods has an associated loss (objective) function whose minimum is the “matched”
map. The loss function of each of the three methods employs the mathematical object for which they are
named. We implement the Wasserstein distance method and electrical force method in Python and evaluate
their performance using both theoretical and numerical techniques.

2

Problem Summary and Statement

Here we will describe the mathematical formulation of our problem. We will adopt the problem statementfrom definitions 2.1-2.4 in [CXHZ] with slight modifications. The geospatial data points will be modeled bya trajectory, see Definition 2.1, and map data will be encapsulated by a road network Definition 2.2.
Let us fix N ∈ N, N ≥ 2, (but almost everywhere we consider the case N = 2).
Definition 2.1 (Trajectory). A trajectory T r is a sequence p = (p1 , p2 , . . . , pn ) of points in RN equippedwith• a sequence, t(p) = (t1 , . . . , tn ), of positive numbers satisfying t1 < t2 < · · · < tn , called the timestampof p,• a sequence, spd(p) = (spd1 , . . . , spdn ), of positive numbers called the speed of p (optional),• a sequence, u(p) = (u1 , . . . , un ), of unit vectors in RN , called the direction of p (optional).
Definition 2.2 (Road Network). A road network (also known as map) is a directed graph G = (V, E)
consists of the set V (resp. E) of vertices (resp. edges) with an embedding ϕ : |G| → RN of the geometricrealization |G| of G. We will identify G and the image ϕ(|G|) by ϕ as long as there is no confusion.
Definition 2.3 (Local Road Network). A local road network is a directed connected subgraph of G = (V, E).
Definition 2.4 (Route). A route r on a road network G = (V, E) is a sequence of connected edges(e1 , e2 , . . . , en ) ⊂ E, i.e. the head of ei coincides with the tail of ei+1 for each i = 1, 2, . . . , n − 1. Let
R denote the set of all routes.
Definition 2.5 (Candidate Routes). For the local road network graph as H of the road network G, we define
CRH = CR := {routes on a local road network graph H},
Definition 2.6 (Map-Matching). Given a road network G = (V, E) and a trajectory T r, the map-matching,
MRG (T r), is the route that is the argument of the minimum of some function L : CR → R+ , called the loss
function.
In the map-matching problem, the trajectory and road network are given as input data. Preprocessing
can narrow the scope of the problem to compare candidate routes on the local road network. The main
goal of this project is to find a suitable loss function such that MRG (T r) is the route that is the minimal
distance from the actual route traveled by a vehicle or pedestrian taken concerning a chosen metric (see § 7.4
for details on chosen metrics). A flow chart of our approach to the map-matching problem can be seen in
Figure 2.1.

2

Page 3 of 34

g-RIPS Sendai 2022, Mitsubishi-A

Road network and trajectory

Pre-processing

Candidate routes

Minimize loss function

Optimal route
Figure 2.1: Steps of the map-matching process3

BackgroundIn this section, we will review three existing map-matching methods, the point-to-point method, thepoint-to-curve method, and the Hidden Markov Model method, and evaluate their strengths and weaknesses.
Surveys of these and many other existing map-matching algorithms are available in [CXHZ, QON].
3.1

Geometric Model

Descriptions of several geometric methods including point-to-point, point-to-curve, and curve-to-curve
methods can be found in [BK]. The simplest geometric map-matching algorithm is the point-to-point method(for this method only, consider a route to instead be a sequence of vertices in the road network). Suppose
there are n trajectory points, pi with 1 ≤ i ≤ n and m vertices vj with 1 ≤ j ≤ m. The procedure for the
point-to-point method is given in [BK] to be
1. For each trajectory point, pi , compute the distance from pi to v for each v ∈ V ,2. Find v ∗ ∈ V such that the euclidean distance between pi and v ∗ is minimal.
3. Let the route found be the sequence of the v ∗ ’s found for each pi , 1 ≤ i ≤ n.
Similarly, the procedure for the point-to-curve method is
1. Project each trajectory point onto each edge in the road network,
2. Compute the distance between each trajectory point and all of its projections,
3. Determine which edge minimizes this distance for each trajectory point
4. Let the route found be the sequence of the edges found in step 3.
Geometric methods are easily implemented and have low time complexity [BK]. Techniques such as
creating Voronoi diagrams can be used to further decrease computational time. For example, one can partitiona subset of R2 based on which points are closer to a given vertex on the graph than all others. Thispartition could be computed once and used to map-match different trajectories on the same road network.
See Figure 3.2 for an example of a Voronoi diagram of a road network. The black lines and circles in the
figure correspond to the road network while the blue lines partition the space into regions whose points areclosest to the given vertex.

3

Page 4 of 34

g-RIPS Sendai 2022, Mitsubishi-A

However, these geometric methods, can be sensitive to measurement errors and are highly dependent on
the network structure [BK]. Bernstein and Kornhauser presented the example in Figure 3.3 to demonstratethis. Suppose the sequence of red points labeled x1 , x2 , x3 is our trajectory and the lines and circles represent
the vertices and edges of the underlying road network. Clearly, the trajectory appears to be closest to thepath from vertex v4 to v5 , but the closest node to trajectory point x2 is v2 .

v1
x1

Figure 3.2: Voronoi Example
3.2

v2
x2

v4

v3
x3

v5

Figure 3.3: Geometric Method Error
Hidden Markov Model

The use of Hidden Markov Models (HMM) is one of the most popular approaches to the map-matchingproblem [CXHZ]. An open source example of a map-matching algorithm using HMM is GraphHopper [GH].
y1

y2

···

yN

x1

x2

···

xN

Figure 3.4: Hidden Markov Model (HMM)
A Markov chain is a probabilistic model for sequential events subject to the condition that the probabilityof a given event depends on the on the previous event alone, i.e. it is a sequence of random variable z1 , . . . , znsatisfyingp(zn |z1 , . . . , zn−1 ) = p(zn |zn−1 ).
A Hidden Markov Model assumes that observations, x1 , x2 , . . . , xn are generated by a Markov chain of
unobserved states y1 , . . . , yn , seen in Figure 3.4. The joint probability of the observed and unobserved statesisp(x1 , x2 , . . . , xn , y1 , . . . , yn ) = p(y1 )p(x1 |y1 )nYp(yi |yi−1 )p(xi |yi ).

i=2

The probability p(xi |yi ) is called the emission probability, p(yi |yi−1 ) is the transition probability, and p(y1 )is the initial distribution. If there are a finite number of states each xi and yi can take on, then we can form
emission, transition, and initial distribution matrices. Each probability is a parameter in our model. Oncethese parameters are established, one of several existing algorithms can be used to compute the probabilities
of some sequences of yi values occurring given the sequence of observed xi ’s. This method can be used tofind the sequence of events with the highest probability given a sequence of observations.
4Page 5 of 34

g-RIPS Sendai 2022, Mitsubishi-A

Example 3.1. A simple application of an HMM can be used to determine how busy a teacher is, given
the observed lecture quality. Take the observed states x1 , . . . , xn ∈ {good, bad} to be the quality of thelecture on days 1, . . . , n, and the unobserved states to be y1 , . . . , yn ∈ {busy, not busy}. In this example, theemission matrix, A, and transition matrix, B, are given below based on the quantities decided in Figure 3.5:




.7 .3
.2 .8
A=
,
B=
.
.6 .4
.9 .1

.6
.7

Busy

Not Busy

.4

.3
.2

.8

.9

Good

.1

Bad

Figure 3.5: Transition and Emission probabilities for Example 3.1.

3.2.1

Application to Map-Matching

From the survey from Chao [CXHZ], HMM models for map-matching have the following setup. The
observed quantities in the HMM for the map-matching problem are the sequence of geospatial data points
and the hidden variables are the possible edges on the road network. The transition probabilities describe
the likelihood of the route containing some edge, given the previous edge, but the emission probabilities
describe the probability of observing a trajectory point given some edge in the road network. The choice of
these parameters differentiates existing HMM models for map-matching. After the parameters are chosen,
one can compute the route in the road network with the highest probability given the observed trajectory.
Unfortunately, this approach to map-matching is limited, because a two-dimensional hidden Markov model
for map-matching cannot be extended to include three-dimensional model or be changed to include speed
information without training on the new data altogether.

4

Mathematical formulation for our Approach

We introduce the two methods each employing one of the following mathematical tool and physicalconcept:• Wasserstein Distance,• Coulomb forces and harmonic oscillator.
We will give our situation settings used in Section 5 and Section 6, with mathematical formulations.
First, we explain how to handle the GPS error. Preprocessing can reduce the effect of GPS error, but
it cannot eliminate it. Although the average value of the GPS error can be obtained, the value of the error
at each trajectory point is often not known. However, for simplicity, we impose and discuss here a rather
strong setting in which the GPS error at each point is obtained. The error Err is used to mean that the true
location is likely to exist within radius Err of the location where the trajectory was observed. For simplicity,
we handle the error Err in the sense that there is always a true position within a radius Err. In addition,in areas with low signal reception, such as skyscrapers, GPS data is less reliable, so speed and direction

5

Page 6 of 34

g-RIPS Sendai 2022, Mitsubishi-A

information is obtained from the IMU1 on the mover2 . In this case, the GPS data has a large error rate,
but the accuracy of the speed and direction information remains high. We make the following assumptions
taking these into account.
Assumption 4.1 (The GPS error).
• Set the GPS error as Err : p → R≥0 . Then, suppose that the spherically-symmetric probability measure
γp is given such that supp γp = B p; Err(p) := v ∈ RN dRN (v, p) ≤ Err(p) . We assume that p istruly located at v with (conditional) probability γp (v) (see Figure 4.6).
• Suppose that there are NO errors for both the speed and direction information.
Remark 4.2.
• On spherically-symmetric probability measure: For example, consider the Gaussian measureparameterized so that the total measure is closeenough to 1 in the bounded domain B p; Err(p) . One
can then cut off the measure by B p; Err(p) and normalize it so that the total measure is 1.
• Validity of the error: Notice that for any p ∈ p, there exists e ∈ E such that B p; Err(p) ∩ e ̸= ∅.
supp γp
Err(p)
p

vp1

vp2

Figure 4.6: supp γp which represents the true region of existence of p. The mover exists at vpi with probability
γp (pi ) for i = 1, 2. then γp (p) > γp (vp1 ) > γp (vp2 ) holds.

If we try to handle GPS errors more rigorously, we face the following problem.
Problem 4.3 (Renormalization of (γp )p∈p ). Although in Remark 4.2, we defined supp γp as a ball in RN
for each p ∈ p, this is unreasonable, given that the true location of the mover cannot exist outside of theroad (see Figure 4.7). Hence, supp γp should be restricted to supp γp ∩ E (⊂ RN ) and renormalized sothat the total probability of being there is 1. However, since we do not consider the width of the road in
our setting, a high-level mathematical theory such as geometric measure theory would be needed to justify
such a renormalization. We, therefore, used a simplified setup, as in Remark 4.2, in this project. This
renormalization will be an important procedure when considering the width of the road (see Figure 4.8).
supp γp

supp γp

e1
p

e1
p

e3

e3

e2

e2

Figure 4.7: The situation where supp γp andedges e1 , e2 , e3 are intersecting. Since thewidth of the edges is not considered in our
setup, the true position of p is the wholesupp γp .

Figure 4.8: The renormalized supp γp by restricting it to where the mover can be locatedin supp γp . This renormalization will be necessary especially when considering the width of
the road.

1 Inertial Measurement Unit
2 see https://www.onosokki.co.jp/HP-WK/products/keisoku/vehicle/lc8_principle.htm,
https://www.digikey.jp/ja/articles/use-inertial-measurement-units-to-enable-precision-agriculture (both articles
are in Japanese).

6Page 7 of 34

g-RIPS Sendai 2022, Mitsubishi-A

It is unreasonable to use RN -distance dRN in quantifying how far away a trajectory point is from the edge
due to GPS error. To address this, we introduce a distance taking into account the GPS error between a
trajectory point and an edge.
Definition 4.4 (The “distance” with error between p ∈ p and e ∈ E). We define the “distance” with error
dErr between p ∈ p (with errors) and e ∈ E (without errors) as
Z
dErr (p, e) := dRN (vp , e) dγp (vp ).
vp ∈B p;Err(p)Remark 4.5. Note that this dErr is a quantity that describes “how far is it p ∈ p to e ∈ E can be considered
to be” and not used as a “distance function”.
Observation 4.6. If p ∈ e, i.e. dRN (p, e) = 0, then dErr (p, e) = 0. It is only in this case that dErr (p, e) = 0.

5

Wasserstein Distance Method

In this section, we explain the method using Wasserstein distance which comes from optimal transport
theory. We will briefly describe discrete optimal transport theory (see also [FG, Sa, Vi], etc.).
5.1

A brief introduction of (discrete) optimal transport theory

Consider transporting sand from a sand pile at x1 , . . . , xn to a hole at y1 , . . . , ym . Note that n, m ∈ N
are independent. Suppose that each sand pile x1 , . . . , xn has µ(x1 ), . . . , µ(xn ) mass of sand, respectively, andeach hole y1 , . . . , ym can contain ν(y1 ), . . . , ν(ym ) mass of sand, respectively. Moreover, we assume that thecost of transporting from a sand pile of xi to a hole of yj is linearly dependent on their distance d(xi , yj ): it
costs d(xi , yj )π(xi , yj ) to transport sand of mass π(xi , yj ) from the sand pile at xi to the hole at yj . Since
the sum of the mass of sand transported from each xi equals the sum of the mass of sand placed in each holeyj , the following holds:nXµ(xi ) =i=1mX

ν(yj ).

(5.1)

j=1

For simplicity, we normalize both sides of (5.1) to be 1. Then, µ, ν can be regarded as probability measures,
respectively. Optimal transport problem is the problem such minimizing total cost for transporting. In otherwords, we consider the following the minimizing problem.
minimizen XmXd(xi , yj )π(xi , yj )i=1 j=1

subject toπ(xi , yj ) ≥ 0mXµ(xi ) =π(xi , yj )for i = 1, 2, . . . , n, j = 1, 2, . . . , m,
for i = 1, 2, . . . , n,

(5.2)

j=1

µ(yj ) =nXπ(xi , yj )for j = 1, 2, . . . , m.

i=1

We call the map π that satisfies these conditions the optimal transport plan or the coupling of µ, ν. Π(µ, ν)denotes a set of all couplings of µ, ν.
Definition 5.1 ((L1 -)Wasserstein distance). For probability measures µ, ν with supp µ = {x1 , . . . , xn } andsupp ν = {y1 , . . . , ym }, we define
W1 (µ, ν) :=infπ∈Π(µ,ν)n X
mXd(xi , yj )π(xi , yj ).

(5.3)

i=1 j=1

Under the appropriate conditions, W1 is a distance function on the probability measure space. We call this(L1 -)Wasserstein distance between probability measures µ and ν.
7

Page 8 of 34

g-RIPS Sendai 2022, Mitsubishi-A

A probability measure space often contains geometrical information about its underlying space. Therefore,
it is useful to analyze the probability measure space using Wasserstein distance to understand the geometricalstructure of its underlying space.
Remark 5.2. Although Π(µ, ν) is an infinite set since the mass of sand transported can be continuously
varied, this is known to be a bounded closed, i.e. compact set. Hence, the right-hand side of (5.3) attains
the minimum value.
5.2

The case in which input data is only the trajectory coordinates and timestamps

First, we introduce probability measures that are determined only from the coordinate information of the
trajectory points.
Definition 5.3 (Probability measures associated with the trajectory and candidate route). Let G = (V, E)
be the local road network graph and m ∈ N be large enough. We divide the candidate route A ∈ CRG into
m + 1 equal parts and let a1 , . . . , am denote the m threshold points, then set V (A, m) := {a1 , a2 , . . . , am }.
Then, we define probability measures µp and νA,m = νA associated with the trajectory p and A ∈ CRG ,respectively, as follows:X11X
δp , νA,m = νA :=
δa .
µp :=
n p∈p
m
a∈V (A,m)Remark 5.4.
• Independence of n and m: The number of trajectory points n and the number of divisions m are
independent. On the other hand, for each candidate route A ∈ CR to be compared, the number ofdivisions m is the same: for any A, B ∈ CR, # supp νA = # supp νB .
• Size of m: Although the threshold points a1 , . . . , am is placed in the image of a discrete approximationof the route A, we still do not know exactly what the minimum magnitude of m should be (see also
Problem 5.18).
• On the definition of µp : Considering the GPS error, γp should be adopted as the definition of µp
instead of δp at each p ∈ p, but for simplicity, we defined µp in this way (see also § 5.4).
Assumption 5.5. In the following, unless otherwise noted, we proceed with the discussion in this section
assuming that the dividing point of the candidate route is m.
Examine µp and ν, i.e. νA,m = νA and νB,m = νB in the following example as the simplest case where is
it is difficult to determine the true route.
Example 5.6 (Square model). Consider the case where the edges of the local road network graph (⊂ R2 )form a square, and suppose that the trajectory points p1 , . . . , pn is observed near the diagonal as shown in
Figure 5.9. Although the trajectory points are drawn in a zigzag pattern in Figure 5.9, our approach canbe used for anything near the diagonal line. We assume that we know from the information of timestamps
that the trajectory began near v1 and ended near v4 . It is difficult to determine from this trajectory along
whether (v1 , v2 ) → (v2 , v4 ) (we call this route A) or (v1 , v3 ) → (v3 , v4 ) (we call this route B) is the true route
traveled. Therefore, we propose the Wasserstein method, which is the method of route determination basedon the optimal transport theory described in § 5.1.
We will find µp and νA,m and νB,m in this situation, respectively. From the definition, µp can be written
directly as
n

1X
δp .
µp =
n i=1 iHence, we put weights of 1/n on each point p ∈ p as shown in Figure 5.10. To find νA,m and νB,m , we firstdivide the route A and B into m + 1 equal parts. These are obtained by putting weight 1/m on each of them threshold points in this way, respectively (see also Figure 5.11 and Figure 5.12):m

νA,m = νA =

m

1 X
δa ,
m j=1 jνB,m = νB =

8

1 X
δb .
m j=1 jPage 9 of 34

g-RIPS Sendai 2022, Mitsubishi-A

v1

v1

v2

v2

p1
p2

1/n

p3

1/n

1/n

Route B

Route A
pn−1
pn

v3

1/n
1/n

v3

v4

v4

Figure 5.10: A probability measure µp on the
trajectory points. The weight 1/n is placed on
the trajectory points p1 , . . . , pn respectively.

Figure 5.9: The case where the trajectory
points p1 , . . . , pn is observed near the diagonalof a square road network.

v1
v1

v2

v2

a1 a2

b1
b2

Route B
Route A
am−1
am

v3
v3

bm−1 bmv4

v4

Figure 5.11: A probability measure νA,m associated with the route A. The edges (v1 , v2 )and (v2 , v4 ) forming the route A are divided
into m + 1 equal parts, each with 1/m weighton its threshold points a1 , . . . , am .

Figure 5.12: A probability measure νB,m associated with the route B. The edges (v1 , v3 )and (v3 , v4 ) forming the route B are divided
into m + 1 equal parts, each with 1/m weighton its threshold points b1 , . . . , bm .

Remark 5.7. If m is odd and can be expressed as m = 2ℓ + 1, then the threshold points aℓ and bℓ coincidewith the nodes v2 and v3 , respectively.
If input data is only the trajectory coordinates and timestamps, then the true route is determined byWasserstein distance between µp and ν.
Remark 5.8. Although φm (A) take finite values for any A ∈ CR and any m ∈ N, they might not converge
as m → ∞. Although we would conjecture that they converge as m approaches infinity in an even or odd
state, respectively, their values may be different, i.e. their values may oscillate as m → ∞.
Remark 5.9. Even if all trajectory points are on the candidate route A ∈ CR, W1 (µp , νA,m ) ̸= 0 holds ingeneral. Therefore, we should apply this method only when it is non-trivial which route was passed such as
the trajectory points are clustered near the diagonal.
Remark 5.10. Note that the trajectory p is a finite set, and a route is a curve, that is, a set of infinite points.
The point-to-curve method projects from the trajectory point to the edges that compose the candidate route,so it is a “point-to-(point on curve)” iteration. In contrast, our method transports weights from trajectory
p to a candidate route A ∈ CR, in other words, it is a “trajectory-to-route” process.

9Page 10 of 34

g-RIPS Sendai 2022, Mitsubishi-A

v1

v2
p1 •
p2 •
p3 •
•
p4

Route B

•
p5

v3

•
p6

Route Av4

Figure 5.13: A local road network graph G. The trajectory is close to the route A, while the unit vector issimilar to the shape of the route B. We should judge the route B in such a case.

5.3

Incorporating speed and direction information

In this subsubsection, we propose a method that takes speed and direction information into account.
5.3.1

Our motivation

Consider an extreme example such as Figure 5.13. Suppose that each pi is located between edges v1 v3and v2 v4 for i = 1, 2, 3, and each pj is located between edges v1 v2 and v3 v4 for j = 4, 5, 6. Consider µp , νA
and νB in this case, and determine the route from only location information using the method described in§ 5.2. From the situation,
W1 (µp , νA ) > W1 (µp , νB )(5.4)

holds and We would have to select the route A. In this case, however, the GPS data is probably noisy due toskyscrapers and other factors, making the location data unreliable. Also, since direction should have a small
margin of error even in such cases, we would like to determine the route from the direction information as
much as possible. Thus, in this case, the true route would be the route B.
5.3.2

Our strategy

The first idea is to use speed and direction information to modify the transport method, i.e., W1 distance.
In this case, however, we would modify the optimal transport plan, not the loss function for W1 distance.
Although the modification itself may not be difficult, the modified W1 distance that is constructed as a result
is likely to be very hard to handle, in particular, the triangle inequality may no longer hold. Then, we modify
the probability measures µp and ν using speed and direction information. The following is a brief descriptionof the flow of the method.
1. Place a weight on each threshold point a ∈ V (A, m) that takes speed and direction information into
account. This weight is separated into two parts: one wA which is uniformly placed on each route andanother wp,A which is affected by speed and direction.
2. According to the weights, the local road network graph is updated to a weighted graph. This weighted
graph is a graph that is defined for each route A and we denote it as GA .
10

Page 11 of 34

g-RIPS Sendai 2022, Mitsubishi-A

3. Perturb µp and νA on the weighted graph GA by ε, where 0 < ε ≪ 1. We denote the probabilityεεmeasures obtained in this way as µεp and νA, respectively. We similarly define νB
for B ∈ CR.
4. Then, comparing the valuesεε
W1 (µεp , νA)
W1 (µεp , νB)
and,
W1 (µp , νA )W1 (µp , νB )determine the route with the smallest value to be the true route. In other words, we adopt the valueε
W1 (µεp , νA)
W1 (µp , νA )(5.5)

as the quantity for route determination.
εRemark 5.11 (The quantity to be used to determine route). Note that since 0 < ε ≪ 1, W1 (µεp , νA)εεand W1 (µp , νB ) are not yet valid quantities for determining route because of the strong influence of location
information. There is also a technique to make this ε larger, but in that case, the question is how to proportion
it to the location information. Therefore, we adopted quantities normalized by W1 (µp , νA ) for 0 < ε ≪ 1.
This quantity (5.5) is very useful in Riemannian geometry on graphs. See Appendix A for more details.

5.3.3

The definition of weights including speed and direction information

We will put additional weighted edges in this vertex set V . Our method quantifies the distance between
p and each route separately. Therefore, it is necessary to consider two types of weighted edges, since each ofthem is evaluated from a different point of view.
The first type is a “natural” weight in the sense that it decides which routes to consider.
Definition 5.12 (Natural weights for each route). For each trajectory point p ∈ p and each candidate routeA ∈ CR, we define wp,A , wp,B : RN → {0, 1} as(
1 (v ∈ V (A)),
wp,A (v) :=0 (otherwise),
V (A) is the threshold points set of A ∈ CR.
Next, we introduce a weight including speed and direction information. This weight is a quantity that
depends on the input data at each trajectory point p ∈ p, and is commonly defined independently of whichroute is considered.
We will make some preparations to define this weight.
Definition 5.13 (Weights depending on speed and direction information).
• For each p ∈ p and e ∈ E, we define
Cp,e := spd(p) ·⟨up , e⟩,
exp dErr (p, e)

where e is the unit vector parallel to e for each e ∈ E. Then, we define Cp :=Pe∈E Cp,e .

• With the above preparations, we define wp : RN → R≥0 for each p ∈ p as follows:
wp (v) :=

X

Cp,e
1
.
·
C
#
e
∩
V (A)
p,e
e∈E

P
e∈E;
e∋v

Observation 5.14 (Weighted local road network graph). The weight introduced in Definition 5.12 andDefinition 5.13 modifies the local road network G = (V, E) to the weighted graph G = (V, E, w). Noticethat the weights of each eA ∈ E(A) are all zero, so these edges vanished. Moreover, note that the11

Page 12 of 34

g-RIPS Sendai 2022, Mitsubishi-A

vertex set VA of GA determined by the weights wp,A and wp is VA = p ∪ V (A). In order to consider the
perturbation (of the probability measure µp associated with the trajectory) at G = (V, E, w), let us now
calculate the weighted degree dp,A,m = dp,A of each p ∈ p, which can be considered in #CR types dependingon the different weights wp,A for each A ∈ CR, but the values are the same:


X
X
X X
C
1
 = m + 1.
P p,edp,A :=
wp,A (v) + wp (v) = ·1 +e∈E Cp,e # e ∩ V (A)v∈Ve∈E;v∈V (A)v∈V (A)e∋v

This right-hand side does not depend on A ∈ CR. For any A ∈ CR, we have dp,A = m + 1.
Finally, we get to the definition of ε-perturbations µεp,· and ν·ε .
Definition 5.15 (The ε-perturbations of the probability measures associated with the trajectory and routes).
For each ε ∈ [0, 1], p ∈ p, candidate route A ∈ CR and a ∈ V (A, m), we define as follows:(
1−εµεp,A (v) :=w(x)+w (x)ε · p,A dp,A p1 − εενa (v) := ε · n10(v = p),(v ∈ V (A, m)),(v = a),(v ∈ p),(otherwise),

Next, we define the ε-perturbations of the probability measures associated with the trajectory and routes asfollows, respectively:µεp,A :=

1X εµ ,n p∈p p,A

ε
:=
νA

1
m

X

νaε .

a∈V (A,m)εεcan be simply expressed by the definition asand νA
Remark 5.16. νA

(
ε
νA
(v) =

1−ε
m
ε
n

(v ∈ V (A, m)),(v ∈ p).

00= νB hold.
= νA and νB
Observation 5.17. It is clear from the definition that µ0p,A = µp = µ0p,B , νA

5.3.4

Briefly explanation of route determination in Figure 5.13 using Wasserstein method

Consider Figure 5.13 again and use the Wasserstein method to determine the route.
First, we find each C. For simplicity, suppose that spd(p), Err(p) are constant at each p ∈ p (spd(p) =: S),→−−→
vec(p(i)) ⊥ −v−
1 v3 for i = 1, 2, 3 and vec(p(j)) ⊥ v1 v2 for j = 4, 5, 6. Then, we obtain for i = 1, 2, 3
S⟨upi , e⟩= , Cpi ,v1 v2 = Cpi ,v3 v4 = 0,
exp dErr (p, e)
exp dErr (p, e)where exp dErr (p, e) := exp dErr (p, v1 v3 ) = exp dErr (p, v2 v4 ) . We also obtain for j = 4, 5, 6
Cpi ,v1 v3 = Cpi ,v2 v4 = spd(pi ) ·Cpj ,v1 v3 = Cpj ,v2 v4 = 0,

Cpj ,v1 v2 = Cpj ,v3 v4 = spd(pj ) ·⟨upj , e⟩S=,
exp dErr (p, e)
exp dErr (p, e)where exp dErr (p, e) := exp dErr (p, v1 v2 ) = exp dErr (p, v3 v4 ) . Then, we have for i = 1, 2, 3 and j = 4, 5, 6
Cpi ,v1 v3
1
Cp ,v vCpi ,v1 v2
Cp ,v v
= = i 2 4,=0= i 3 4,
Cpi
2
Cpi
Cpi
Cpi
Cpj ,v1 v3
Cpj ,v2 v4
Cpj ,v1 v2
Cpj ,v3 v4
1
=0=
,
= =
.
Cpj
Cpj
Cpj
2
Cpj
Suppose that we divide the routes A and B into m = 2ℓ + 1 equal parts. Then, the weight wp (·) is distributed
as follows:12Page 13 of 34

g-RIPS Sendai 2022, Mitsubishi-A(
wp1 (v) = wp2 (v) = wp3 (v) =0(

wp4 (v) = wp5 (v) = wp6 (v) =1
1
2 · ℓ+1

1
1
2 · ℓ+10(v = aℓ+1 , . . . , a2ℓ+1 , b1 , . . . , bℓ+1 ),(otherwise),(v = a1 , . . . , aℓ+1 , bℓ+1 , . . . , b2ℓ+1 ),(otherwise).

Then, it holds that µεp,A = µεp,B for any ε ∈ [0, 1]. The supports of probability measures µp , νA , νB , µεp,A (=εεµεp,B ), νA
and νBare illustrated in Figure 5.14, Figure 5.15, Figure 5.16,Figure 5.17, Figure 5.18 and Figure 5.19, respectively. Let α be α := 1 − ε.

Figure 5.14: supp µp . • =1/6.

Figure 5.15: supp νA . • =1/m.

Figure 5.16: supp νB . • =1/m.

Figure 5.17: supp µεp,A • =α/6, • = ε/4(ℓ + 1) and •= ε/2(ℓ + 1).

εFigure 5.18: supp νA. •=ε/6 and • = α/m.

εFigure 5.19: supp νB. •=ε/6 and • = α/m.

Now note that the value (5.5) can be decomposed as follows:ε
W1 (µεp,A , νA)A(ε) + δ(ε)
A(ε)
δ(ε)
=
=
+
,
W1 (µp , νA )W1 (µp , νA )W1 (µp , νA ) W1 (µp , νA )ε
W1 (µεp,B , νB)
B(ε) + δ(ε)
B(ε)
δ(ε)
=
=
+
,
W1 (µp , νB )W1 (µp , νB )W1 (µp , νB ) W1 (µp , νB )

(5.6)
(5.7)

where A(ε) and B(ε) are the transport costs from the trajectory to each route, respectively, and δ(ε) isthe transport cost from one route to the other. Note that δ(ε) is a common quantity in both transports.
Moreover, we can see that
A(ε)
B(ε)
≒
W1 (µp , νA )W1 (µp , νB )13

Page 14 of 34

g-RIPS Sendai 2022, Mitsubishi-A

is valid when ε is close enough to 0. Then, the magnitude relationship between (5.6) and (5.7) is determinedby the magnitude relationship between W1 (µp , νA ) and W1 (µp , νB ).
Hence, since (5.4) holds, we obtainεε
W1 (µεp,A , νA
W1 (µεp,B , νB))>.
W1 (µp , νA )W1 (µp , νB )(5.8)

Therefore, the Wasserstein method outputs the route B as the true route in a situation like Figure 5.13.
5.4

Future problems

We list the problems that should be addressed to develop the Wasserstein method in the future in thefollowing aspects:• the definition of a probability measure associated with the trajectory (Definition 5.3),• how to distribute weights by speed and direction information (Definition 5.13),• the property to be shown for implementation:– the case where determination is based only on location information (§ 5.2),– the case where determination is based on speed and direction information as well (§ 5.3).
The definition of a probability measure associated with the trajectory (Definition 5.3)
Inspired by Definition A.6, we defined the probability measure µp associated with the trajectory as
Definition 5.3, but it may be useful to make the definition correspond to (γp )p∈p , which represents the rangeof GPS error see also § A.1. This is natural considering the original meaning of “coarse Ricci curvature”,
which was introduced in [Ol], as defined on metric measure spaces. However, the concrete calculations would
be very complicated beyond our setup.
How to distribute weights by speed and direction information (Definition 5.13)
After all, the distribution of weights in Definition 5.13 is the most important point in the map-matchingproblem. The weight in Definition 5.13 was simply assigned equally to all the edges of the local road network,
but I think it should be more sophisticated. For example, the weight value to be placed on the thresholdpoint of each edge can be different depending on the location of each p ∈ p and the spd(p).
The property to be shown for implementation
This is a common problem for both § 5.2 and § 5.3 in the sense that the goal is to derive the property
corresponding to Theorem A.5 in each setting. However, the properties to be shown in each situation are
somewhat different and will be explained separately.
The case where determination is based only on location information (§ 5.2)
For each candidate route A ∈ CR, as m is increased, supp νA,m approaches to the route A respectively.
Our idea is that we can determine which route is closer to the trajectory points p by computing W1 (µp , νA,m )for a sufficiently large m ∈ N. Then, we need to show the following.
Problem 5.18. For any A, B ∈ CR, there exists a me ∈ N such that W1 (µp , νA,m ) − W1 (µp , νB,m ) is amonotone function for any m ∈ N such that m ≥ m.
e
We can determine which route the trajectory points is closer to by checking the sign of W1 (µp , νA,m ) −W1 (µp , νB,m ) for a sufficiently large m if we can prove Problem 5.18. Hence, if W1 (µp , νA,m ) < W1 (µp , νB,m )holds, then we can conclude that the trajectory points is closer the route A, which is the true route.
Remark 5.19. Although W1 (µp , νA,m ) take finite values for any A ∈ CR and any m ∈ N, they might notconverge as m → ∞. Although we would conjecture that they converge as m approaches infinity in an evenor odd state, respectively, their values may be different, i.e. their values may oscillate as m → ∞.
14

Page 15 of 34

g-RIPS Sendai 2022, Mitsubishi-A

The case where determination is based on speed and direction information as well (§ 5.3)
We could not show the piecewise linearity like Theorem A.5 in LLY–Ricci curvature yet within the periodof this project (see also Remark A.12):
Problem 5.20. Is the function ε 7→ κ(ε; p, A) piecewise linear on [0, 1]?
This is due to the complexity of our setting, and I conjecture κ(ε; p, A) will almost certainly be piecewise
linear on ε ∈ [0, 1] for each candidate route A ∈ CR. However, the number of segments will depend on n andm, not 3. Therefore, we will calculate κ(ε; p, A)/ε by assigning a value sufficiently close to 0 to ε. It would
have been much more useful in terms of implementation of this first linear part ε ∈ [0, ε̃(n, m)] could havebeen derived.

6

Methods by Physical Analogy

We introduce the “electric method” and the “harmonic oscillator method” as a way to define our new
loss function.
6.1“Electric” Method

First, we propose a method to measure the closeness of trajectory points and routes by the interactionof Coulomb forces. In this method, we consider not only trajectory points, but also the entire polyline
connecting the trajectory points, and compare it with the entirety of each route. That is, we consider a“trajectory segments-to-route method.”
6.1.1

Details of method

We describe the specific procedure of this “electric” method. The abstract of the step is as follows:
1. Connecting trajectory points with segments and parameterizing them and each candidate route.
2. Giving the polyline and routes opposing electrical charge.
3. Calculating electric forces between the polyline and routes.
4. Choosing “the closest route” as the true route.
First, we connect trajectory points with segments and parametrize the candidate route. We assume that
this polyline can move but the road map is fixed.
Definition 6.1 (polyline). We define the polyline lp : [0, 1] → RN associated with a trajectory p = (pi )0≤i≤n
by


i+1
i
≤t≤
.
(6.1)
lp (t) := (1 − t)pi + tpi+1nnDefinition 6.2 (parameterized route). Let r = (e1 , . . . , em ) be a route. We define the parametrised edgeei : [0, 1] → RN by parametrising |ei | by arclength parameter and the parametrised route r(t) : [0, m] → RN
by union of the parametrised route ei .
Remark 6.3. We assume that route is piecewise smooth and regular, thus arclength parameter exists oneach interval.
Next, we give the candidate routes and our trajectory polyline opposing electrical charges. For simplicity,
we assume that electric charge density is uniform.
Definition 6.4 (electric charge). We define electric charge density as uniform probability density function
on the polyline ρpolyline : [0, 1] → R and on the candidate route ρroute : [0, 1] → R.

15

Page 16 of 34

g-RIPS Sendai 2022, Mitsubishi-A

Then, all of the candidate routes exert an electric force on our polyline. The magnitude of Coulomb forces
is represented as follows:
Z
[0,m]

˙
(l(t) − r(s)) · l(t)ρpolyline (t)ρroute (s) dtds.
2[0,1] dErr (l(t), r(s))

Z

(6.2)

We choose the route which exerts the most force on the polyline as the true route.

Route A
p1

p2
p4

p3

Figure 6.21: Connecting traRoute B
Figure 6.20: Given trajectory
points and candidate routes.

jectory points and giving opposing electric charges.
The
red colored area have positive
charge and the blue colored area
have negative charge.

6.1.2

Figure 6.22:

The polyline
moving to “closer” route by
Coulomb forces.

Observations• Strength:
By using this method, we can consider both the whole polyline and the entirety of each route all at
once, because the electric forces are determined by the whole body of each.
• Problems: We have two problems with this “electric” method;– The theories about static electrodynamics and rigid body dynamics do not tell us how to treat
speed and direction data.
– Coulomb forces are proportional to square inverse of distance, so when the polyline intersects withthe route they diverge. We can avoid divergence by adding a constant to the denominator. But,
if a constant is small, the effects of intersection points are still dominant. Thus it is difficult tochoose the correct route which has fewer intersection points than other wrong routes like the below
figure.
p2 p4 p6

• • •p8
p1 p3 p5 p7• • • •

p9
p10
p11

• • •
•
•
•
•

• • • •••••

Also, due to the nature of floating-point arithmetic, it is difficult for the computer to distinguishroutes that share intersection points like the below figure.

16Page 17 of 34

g-RIPS Sendai 2022, Mitsubishi-A

p2 p4 p6

p7
p1 p3 p5

p8
p9
p10

• • •

• • •

• • • ••••• • • ••••

As we see in above figures, consequently, we must take the appropriate constant. But because it
is dependent on the number of intersection points, we cannot choose it arbitrarily. This problem
was found in the implementation step.
• Strategy to solve this problem:
Our option to avoid the above issue is to instead invert the problem: that is, we replace maximizing
inverse squared distance maxA∈CR r12 to minimizing squared distance minA∈CR r2 . The latter term issimilar to the potential of harmonic oscillators. We consider the analogy of this model.
6.2“Harmonic oscillator” method

Based on the consideration of the “electric” method, we developed a policy to minimize the square of the
distance, and this equation inspired us to consider an analogy with the analytical dynamical description ofharmonic oscillators.
Harmonic oscillators are a system in which a material point is subjected to a force proportional to its
distance from some fixed point. An example is the motion of a mass connected to a wall by a spring.

k

k

x(t)

m

v(t)

If mass m has velocity v(t) and displacement x(t) from the natural length of the spring at a time t, thissystem has momentum 12 mv(t)2 and potential 12 kx(t)2 where k is a spring constant. Then the Hamiltonian
H(t) and the energy E of this system are defined by

H(t) :=

1
1
mv(t)2 + kx(t)2 ,
2
Z2
E :=

L(t) dt.

We formulate our map-matching problem as the following setup: we take the true root as the location where the potential of the system is minimal. Then, assume that the trajectory point has a minutedisplacement from the point of minimum potential due to error.
We now assume that our assumptions about error place more emphasis on speed and orientation thanon location information. Therefore, we formulate the distance to the edge to be small. This operation isconsistent with the picture that the harmonic oscillator approximates the system at small displacements fromthe point where the potential is smallest.

17Page 18 of 34

g-RIPS Sendai 2022, Mitsubishi-A

6.2.1

Details of method

We describe the specific procedure of this “harmonic oscillator” method. The abstract of the step is asfollows:
1. Connecting trajectory points by segments
2. Giving speed, direction, and error information to points on polyline connecting trajectory points.
3. Defining a score of edges included in the candidate route and choosing the edge which minimizes thescore.
4. Connecting points on the polyline and anchors on the “closest” edge by “spring”.
5. Defining “momentum” and “potential” of each point on the polyline.
6. Calculating “energy” of “Hamiltonian” for the polyline and each candidate route.
Definition 6.5. Let p = {p0 , . . . , pn } be trajectory points
1. We define polyline lp : [0, 1] → RN by
lp (t) := (1 − t)pi + tpi+1

i+1
i
≤t≤
,
n
n


0≤i≤n−1 .

(6.3)

2. We define the speed, direction, and error of each point on the polyline by piecewise linear interpolation:
ii+1• (speed) spdlp (t) := (1 − t)spdi + tspdi+10≤i≤n−1 ;n ≤t≤ n ,i+1i
0≤i≤n−1 ;• (direction) ulp (t) := (1 − t)ui + tui+1
n ≤t≤ n ,
i+1i
• (error) Errlp (t) := (1 − t)Erri + tErri+10≤i≤n−1 .
n ≤t≤ n ,Definition 6.6 (Score of edges included in candidate route). Let p be a trajectory and lp denotes the polyline
formed by p Let p be a trajectory, lp : [0, 1] → RN be the polyline formed by p and A = {e ∈ E} ∈ CR bea candidate route.
1. We define the score Slp : [0, 1] × A → R by
Slp (t, e) := ⟨spdlp (t) ulp (t), ⃗e⟩RN · dRN (lp (t), e),(6.4)

where ⃗e ∈ RN is the unit vector induced naturally from timestamp information of trajectory.
2. We define the anchor map Anclp ,A : [0, 1] → A; t 7→ Anclp ,A (t), so that S(t, Anclp ,A (t)) = maxe∈A Slp (t, e).
Definition 6.7. Let p be a trajectory and lp denotes the polyline formed by p. We define the “mass” and“spring constant” of this system then use them to define the “momentum” and the “potential ”.
1. We define the mass mlp : [0, 1] → R by
mlp (t) :=1.
Errlp (t) + 1

2. We define the momentum Mlp : [0, 1] × CR → R by*
Mlp (t, A) := mlp (t)spdlp (t)ulp (t)log(1 + spdlp (t))+2−−−−−−−→, Anclp ,A (t)

(6.5)

RN

−−−−−−−→
where Anclp ,A (t) ∈ RN is the unit direction vector of Anclp ,A (t) induced naturally from timestamp
information of trajectory.
18

Page 19 of 34

g-RIPS Sendai 2022, Mitsubishi-A

3. We define the spring constant klp : [0, 1] → R byklp (t) := exp −Errlp (t) .

(6.6)

4. We define the potential Plp : [0, 1] × CR → R by2
Plp (t, A) := klp (t) · dErrlp (t) lp (t), Anclp ,A (t)
.

(6.7)

Remark 6.8.
1. We have now decided to focus on speed and direction rather than positional informationin our assumptions, and have formulated the distance to be small. Therefore, the spring constants
involved in the potential, of which distance is a major component, are allowed to decay exponentiallywith respect to the error. On the other hand, for mass, it is similarly affected by errors, but the choiceof the edges connecting the springs is based on the score determined to direction, and since the direction
is assumed to have no error, the correction is assumed to decay about a first order polynomial. 1 isadded to eliminate singularity.
2. Velocities vary widely among pedestrians, bicyclists, and motor vehicles, and momentum can be very
different for them, which can lead to significant differences in results. This is undesirable in light of thevalidity of the method. The addition of 1 is to eliminate singularity.
Definition 6.9 (loss function). Let T r = p be a trajectory. We define a loss function ET r : CR → R byZ
ET r (A) :=Mlp (t, A) + Plp (t, A) dt.
(6.8)
t∈[0,1]

When several candidate routes are given, we calculate the energy of each route on the trajectory points
and choose the route that has the least energy as the true route. For example, trajectory points and candidate
routes are given as below figure, we calculate ET r (A) and ET r (B). Then if ET r (B) < ET r (A) holds, wechoose route B as the true route.
Route A
p1

p2
p3

p4

Route B

6.2.2

Observations• Strength:– There is no divergence problem in contrast to the “electric” method.
– We can take into account information about speed, direction, and error naturally.
• Weakness:
The edges connecting the springs were formulated so that only one edge with the smallest score isselected and the springs are connected to that edge, but this method is inadequate to reflect the shapeof the trajectory points and the shape of the entire route. Consider the case where a route and trajectory
points are given such that a similar shape is repeated.

19

Page 20 of 34

g-RIPS Sendai 2022, Mitsubishi-A

Figure 6.23: Given trajectory points and candidate routes.

Figure 6.24: The candidate route that seems to be the true route.

Figure 6.25: The candidate route that is chosen to be the true route by “harmonic oscillator” method inwhich the only one route is chosen to be connected to trajectory points

In this case, the true route should be the one that has the same number of peaks. However, the methodof selecting one edge and connecting it to the spring would select a route with only one peak.
To improve this, we want to connect all edges with springs and give them appropriate weights, so thatthe shape can be more reflective of the shape.

7

Implementation

After formulating the proposed mathematical methods into robust map-matching algorithms, we implemented them in Python to evaluate their performance numerically. The code written for the project can befound at: https://github.com/gjgress/G-RIPS-2022-Mitsubishi-A. We will show the resulting solutionof our algorithms on two cornerstone examples:
1. the Square road network with two candidate routes to decide between and a trajectory in a neighborhood of the diagonal Figure 7.26
2. the Sendai map example with road network from OpenStreetMap and trajectory collected from walkaround Sendai. Figure 7.27
The square example has a total of 6 trajectory points and four edges in its road network. The left andbottom edges (in blue) make up route 1, and the top and right edges (in yellow) make up route 2. The
Sendai example has 645 trajectory points (in blue) and 191 edges (in white).
The dataset we used for evaluation of the electrical and harmonic oscillator methods was the Dataset
for testing and training of map-matching algorithms [KCMMN]. This dataset was a good fit for the project
because it included ground truths, despite only possessing GPS data. We were unable to apply the Wasserstein20

Page 21 of 34

g-RIPS Sendai 2022, Mitsubishi-A

Figure 7.26: Square Example
Figure 7.27: Sendai Example
method to this data set due to issues discussed in § 7.2. Although we pre-processed the BDD100K open
dataset provided by Berkeley for evaluation, it does not contain ground truths, so we were unable to utilize
it for this project. [YCWXCLMD].
7.1

Preprocessing Routes

As the number of edges in a directed graph increases, the number of routes between two nodes grows at a
rate similar to the factorial (and can be infinite if loops are allowed); moreover, the determination of all suchroutes is NP-hard. As a result, it is not sufficient to simply determine a suitable loss function and evaluate all
possible routes. Finally, our loss functions are often expressed as integrals, so even once we have our routes,
we need to interpolate each route so that the sum sufficiently approximates the integral. In practice, we doall this by restricting to a sub-graph, using an algorithm to obtain a restricted set of candidate routes, andthen interpolating the resulting routes. We refer to these steps as preprocessing.
Of course, the preprocessing process need only be done once for a given dataset. That is, once we
preprocess our graph with our trajectory, we can utilize the processed information for all of our algorithms.
Furthermore, all the preprocessed information can readily be saved to disk, so for larger datasets, this need
only be done once.
7.1.1

Restricting to a Sub-graph

The method we use to restrict to a sub-graph is relatively simple. We choose a k ∈ Z+ , and for eachpoint pi in the trajectory, we find the set of k-nearest edges with their respective nodes (or all the possibleedges, in the rare case that there are fewer than k edges), denoted Epki . Our restricted sub-graph G̃ is thengiven by[G̃ =Epki .
(7.1)
pi

Unlike k-nearest nodes, k-nearest neighbors are more computationally difficult. A k-d tree is not possible
because we have to consider all points around an edge. As a result, we are forced to rely on geometric tools
to find this set. We do this by choosing an r > 0 and creating a square of size 2r centered at each pi , andfinding the geometric intersection with the road network. If at a given point the square does not intersect atleast k edges, we double the value of r and repeat the process, until we have found enough edges or we havedoubled the value ten times. This method was inspired by the implementation of FMM [YG].
This method is quite inefficient. In the future, we hope to explore using mathematical geometric tools tobetter find our sub-graph.
21

Page 22 of 34

g-RIPS Sendai 2022, Mitsubishi-A

7.1.2

Obtaining Candidate Routes

Once we have obtained our sub-graph G̃, we utilize a variant on Dijkstra‘s algorithm to generate a listof short paths between our source and target. First we identify the source and target node on G̃ by finding
the closest neighbor to p1 and pn . Then we follow the standard procedure of Dijkstra’s algorithm by visiting
the nodes connected to the route with lowest cost. Unlike Dijkstra’s algorithm, we keep each route obtained,even if we have found a shorter route to get to the same intermediate node. We continue Dijkstra’s algorithm
until we have either enumerated every route or have obtained the number of routes desired. We do returnthe cost of each route as well, but we do not use this information in our algorithms– the shortest route is notnecessarily the best candidate route!
This method is a reasonably fast way to take a breadth-first search approach, but nevertheless is still quite
slow. Furthermore, while Dijkstra’s algorithm can be modified to allow multiple source and target nodes, wechoose not to apply this as it would require us to generate many more candidate routes to ensure the true
route is within our set.
There are some edge cases where this algorithm will fail. For one, if the number of nearest edges is
not sufficiently high, the sub-graph will not be able to find any candidate route. In this situation one
has to increase the number of candidate edges and start the algorithm from scratch– which can be very
computationally expensive. The algorithm is also not designed to handle routes where the same edge istraversed more than once; these cases are exceedingly rare, and so we did not find it necessary to modify thealgorithm for this purpose.
This algorithm also is difficult to run in parallel. Conceptually, it should be possible to search multiple
routes at once, but this requires a more robust task assignment process. We did not have the time to explore
this option, and more importantly, it may not be worth it if better candidate route generation methods areavailable.
Despite these shortcomings, the candidate routes obtained are usually high quality. Unlike some datadriven implementations, we also know that when we find our optimal route, it will be a valid route.
7.1.3

Candidate Route Interpolation

The interpolation process is rather simple. Because our objects are Shapely geometries, we can interpolate
each edge using Shapely’s built-in interpolation. We choose constants nr , nt ∈ Z+ and interpolate each edgeof the candidate route into nr segments and each edge of the trajectory into nt segments.
One concern is how to determine the constants nr , nt a priori. We found that choosing nr ≈ 100 andnt ≈ 10 produced the best results; however, in the future we wish to implement a more dynamic approach
which interpolates based on the relative density of the points.
Typically with simpler loss functions such as the electric method or harmonic oscillator method, increasing the interpolation points does not significantly increase computation time. However, in the case of theWasserstein method, it had a big impact on computation time. Thus the optimal choice of nr and nt depend
greatly on the method.
A shortcoming of this method is that it interpolates every edge equitably. This means that very smalledges will have a high density of points, while a long edge will be more sparse. This may give an unfair
weighting to small edges. In the future we wish to explore a dynamic choice of interpolation points based on
the lengths of an edge.
7.2

Wasserstein Distance Implementation

Let {pi } be the set of trajectory points and discretize A ∈ CRG into m + 1 equal parts to obtain the set ofthreshold points V (A, m) := {a1 , a2 , . . . , am }. Then, from § 5.1, we know can find the Wasserstein distance22

Page 23 of 34

g-RIPS Sendai 2022, Mitsubishi-A

by solving the following linear programn X
mX

minimized(pi , aj )π(pi , aj )i=1 j=1

subject toπ(pi , aj ) ≥ 0mXµ(pi ) =π(pi , aj )for i = 1, 2, . . . , n, j = 1, 2, . . . , m,
for i = 1, 2, . . . , n,

(7.2)

j=1

ν(aj ) =nXπ(pi , aj )for j = 1, 2, . . . , m.

i=1

After selecting n, m, we solve this linear program in python using the linprog command from Scipy’s
optimization package. The constraint matrix constructed will be size (n + m) × nm, but only 2mn elements
will be nonzero. Thus, the constraint matrix can be saved as a sparse matrix. For the computation, we
remove one of the constraints, because the system is overdetermined as written. In contrast, it is likely thatno d(pi , aj ) will be zero. Therefore, for each candidate route, we must compute nm distances.
For the simple example of a trajectory contained in a square road network, see Figure 7.28, the loss
computed for each route is the Wasserstein distance. The chosen route in this instance is Route 1. We cansee that the loss for Route 1 is smaller, and is the route we would have chosen to be correct through visual
analysis.

Figure 7.28: Wasserstein solution to the Square Example
In contrast, when we apply Wasserstein distance to the example of the Sendai map, the route that minimizes Wasserstein distance is longer and more complex than the true route. In figure Figure 7.29, we see theroute with minimal Wasserstein distance. Not only does it stray from the trajectory substantially, but it also
contains a loop. At this time, we are unsure if this is an issue with the method itself or some implementation
error. One obstacle to determining this is the computational time. Computing the Wasserstein distances on
the set of candidate routes for the Sendai map alone can take several hours.
7.3

Electrical Method and Harmonic Oscillator Implementation

In comparison, the electric method and harmonic oscillator are simpler to compute. Let Tr = {pi } be theset of trajectory points and A = {qi } be the set of points along a candidate route. For each pi we determinepi
the k-nearest neighbors from the set {qi }; that is, for each pi we obtain a subset {qj }j∈1,...,k⊂ {qi }. Then

23

Page 24 of 34

g-RIPS Sendai 2022, Mitsubishi-A

Figure 7.29: Wasserstein solution to the Sendai Example
we compute the sums:
ETr (A) =kXXpi

ActTr (A) =

−1
d(p
,
qj )2 + ε
i
j=1

k
XX

d(pi , qj )2(Electrical Method)(7.3)(Harmonic Oscillator Method)(7.4)

pi j=1

where 0 < ε ≪ 1 is a small constant chosen to prevent divergence.
In particular, for the electric method, we choose k = # {qi } (that is, we compare to all the points along
the candidate route), and for the harmonic oscillator method, we choose k = 1 (that is, we only compare
to the nearest neighbor). On one hand, choosing large k means you take into account more of the polyline;
however, points far away from the polyline may have a large influence on the loss. Typically this does notchange the relative losses between two candidate routes, but it does make the losses numerically much closer–
and in extreme cases, we may lose this relative information due to floating point errors. Conversely, for small
k we are only considering the part of the candidate route most relevant to a given point; however, outliers in
the trajectory may cause the loss of a candidate route to be lower than expected. Therefore, it is likely thatthe optimal value for k lies somewhere in between. Of course, it is difficult to determine a priori what k isappropriate for a given case.
7.4

Evaluation (Error) Method

How do we measure the accuracy of our prediction? We use the following formula presented by Newson
and Krumm [NK] to calculate the error:
d− + d+Err =d0where d0 is the length of the correct route, d− is the length the prediction erroneously subtracted fromthe ground truth, and d+ is the length the prediction erroneously added outside the ground truth. See
Figure 7.30.
24

Page 25 of 34

g-RIPS Sendai 2022, Mitsubishi-A

d0
Ground Truth
Prediction

d−

d+
+
Err = d−d+d0d0 = length of ground truthd− = length of prediction route erroneously subtracted
d+ = length of prediction route erroneously added
Figure 7.30: Error Formula by Newson and Krumm

Figure 7.31: FMM evaluated on the annotated dataset
7.5

Preliminary Results

When tested on the dataset, we found that FMM had an error of 14.7% on average 7.31.
Unfortunately, due to processing power and time constraints, we were unable to evaluate our algorithms
on the entirety of the dataset. In particular, the preprocessing stage was too slow on the much larger networks
to allow for testing. As a result we can only provide preliminary results for two of the methods 7.32 7.33.
We anticipate that the geometric methods will perform more accurately than FMM on the dataset, but
at the cost of (significantly) greater computation time.
7.6

Computational Complexity

We tested the computational complexity of the preprocessing stage and our methods. In particular, thepreprocessing stage scales linearly with the number of candidate routes to find, but scales exponentially with
the number of candidate edges provided 7.34. Both the Electric Method and the Harmonic Oscillator Method
seem to scale linearly O(n) with the number of candidate routes 7.357.36. When the number of distancecalculations were increased (i.e. more nodes in the candidate route and trajectory), runtime increased, but notsignificantly– certainly slower than linear growth. Unfortunately, the Wasserstein method grew too rapidly
when distance calculations were increased, and so a similar figure could not be produced for the method. We
suspect that it scales as a polynomial of degree > 1 with respect to distance calculations.

25

Page 26 of 34

g-RIPS Sendai 2022, Mitsubishi-A

Figure 7.32: Electric Method evaluated onannotated dataset

Figure 7.33: Harmonic Oscillator Method
evaluated on the annotated dataset

The metric algorithm framework can parallelize over the number of candidate routes, which significantly
improves the computation complexity. It is unlikely that improvements can be made within the loss functions
to parallelize further, as Electric Method and Harmonic Oscillator Method utilize array operations, andWasserstein utilizes linear programming.
While candidate route generation has limited room for parallelization, one can easily parallelize the processacross the entirety of the dataset. The simulators can also be easily run in parallel on the entirety of thedataset.
7.7

Metric-based Algorithm Framework

Because the basic formulation of our proposed map matching algorithms are all purely metric-based, we
realized it was more efficient to write the algorithm in a modular capacity. This led to the creation of thegeneralized metric mm algorithm class.
To create a metric based map matching algorithm simulator, one first creates a loss function in one of
two ways:• Providing simple interior and exterior functions which operate in stages on the distance arrays• Providing a singular loss function which takes as arguments the candidate route distance array and the
trajectory distance array.
Often it is simpler to provide the latter, but the former may be convenient in the case that one wishes tomake small modifications to the loss function procedure in a systematic capacity.
If preprocessing has already been performed by another simulator, the interpolated candidate route nodes
and trajectory can be provided directly.
The class then provides a preprocessing method if preprocessing has not already been performed. It’s
only required argument is a GeoDataFrame consisting of LineStrings of the trajectory (this can be generated
from a sequence of points using a utility function provided in another module). If candidate routes are notprovided, it applies the variant on Dijsktra’s algorithm described in 7.1.2.
To run the simulator, one simply has to call the run method. This will return the best candidate route(and its loss) or the loss of all candidate routes depending on parameters passed.
7.8

Data Fusion

In addition to map matching, the Jupyter notebooks providing a framework for implementing and testingalgorithms on driving datasets. One issue that pervades this field of research is the high variety in file formatting, and in particular, these formats are not well-suited for containing IMU data. Our notebooks provides
a rudimentary fusion method to align asynchronous data and incorporates IMU data into a GeoDataFrame

26

Page 27 of 34

g-RIPS Sendai 2022, Mitsubishi-A

Figure 7.34: Preprocessing Computation Runtimes

Figure 7.35: Electric Method Computation Runtimes

27

Page 28 of 34

g-RIPS Sendai 2022, Mitsubishi-A

Figure 7.36: Harmonic Oscillator Method Computation Runtimes(which can be exported to GPX, GeoJSON, KML, etc.) in a sensible manner. In particular, this method
can sample any discrete-time data such as speed, accelerometer, and gyroscope, and merge it with the GPSdata sequence.
7.9

Other Utilities

In the process of creating this framework, several utility functions and classes have been provided. Within
mm utils.py we provide functions such as:• Basic plotting methods• An evaluation method to calculate error• A method to create trajectory edges (LineStrings) from a sequence of coordinates (Points)• A method which provides a road network from OSM (as a MultiDiGraph or GeoDataFrame) from the
trajectory data only
• A method to find k nearest points using k-d trees• A method to find k nearest edges using the method described in 7.1.1

8

Future Work (Implementation)

We briefly mention again improvements discussed in the previous section:• Modify k-nearest edge search with more efficient methods• Modify or replace Dijsktra’s algorithm for candidate route generation• Investigate and improve upon Wasserstein method

28

Page 29 of 34

g-RIPS Sendai 2022, Mitsubishi-A

8.1

IMU Inclusion

While IMU approaches were proposed theoretically, we did not have the time to include IMU information
into our algorithm implementations. Currently the generalized metric-based map matching framework does
not pass IMU data to the loss functions, but it is preserved in all the intermediate stages, so it should be
simple to modify it so that the relevant data is passed along.
8.2

Incomplete Map Matching

One area of interest with map matching problems is incomplete maps. Our assumption is that our road
network includes every possible path in reality. Of course, some older roads or even new roads may not
be incorporated into our road network information. One direction of interest is then to explore if we candetermine when roads are missing from our network based on the loss values along the trajectory.
8.3

Extension to Higher Dimensions

Including the z-axis into our implementation faces one major hurdle: OpenStreetMaps does not include
elevation. Because this information is not included within our network, we cannot test our algorithms in
3-space. One can circumvent this by merging elevation data from an external source. However, highways
and roads are often not incorporated into these data sets, and so it is approximate at best. Instead, it isprobably best to source network information from more detailed data sets, such as proprietary sources, but
results can vary greatly depending on the region.
Fortunately, the Python framework developed is quite modular, so if one does have elevation data intheir network and test data, it is easy to implement within our notebooks. Hopefully, with time, data setssuch as OpenStreetMaps will gain detail, and this line of inquiry becomes more accessible to the scientificcommunity.

9

Conclusion

We developed 3 new loss functions for solving the map-matching problem: Wasserstein, Electrical, andHarmonic Oscillator. Within each theoretical formulation, not only can we determine a route from distanceinformation, but also by accounting for speed and direction information. For the Wasserstein distancemethod, this achieved through a perturbation of the probability measures on both the trajectory data andthe threshold points on the candidate routes. For the harmonic oscillator, the speed and direction areincorporated through the momentum term of the Hamiltonian. Each method was also implemented in
Python, but using only the distance information. The results for Wasserstein distance suggests either theformulation or the implementation requires more investigation. However, the electric and harmonic oscillators
results are promising. Although their computational time is greater than that of FMM, in the experimentsthat we have run so far, the accuracy seems much improved.

Appendix A

Geometric background of Wasserstein method

This section provides the geometric background of the Wasserstein method in Section 5. Wassersteinmethod is based on the concept of “Ricci curvature of graph”, which was introduced by [Ol, LLY]. This
concept is a metric to measure the strength of cohesion between two vertices and has attracted attention asa new tool for graph analysis and has already been applied to real problems ([JL, NLGGS, NLLG], etc.). We
modified Ricci curvature of graphs to quantify the “strength of cohesion” between the trajectory and eachroute.
A.1

Ricci curvature of graphs

In this section, we briefly describe the Ricci curvature of graphs. We considerPa weighted graph. Denote
the weight of edge e as we . In this case, the node degree dx of vertex x is dx = y∼x wxy . We call a graph29

Page 30 of 34

g-RIPS Sendai 2022, Mitsubishi-A

with a weight of 1 on all edges a combinatorial graph. First, we introduce a transition probability measure
on a node. In this section, we use x, y as symbols that denote nodes.
Definition A.1 (Random walk with idleness parameter ε ∈ [0,1]). For each node x ∈ V and ε ∈ [0, 1], we
define a random walk µεx as follows:1−ε(y = x),wxyεµx (y) := ε · dx (y ∼ x),
0(otherwise).
The Ricci curvature with idleness parameter ε is defined as follows.
Definition A.2 (ε-Ollivier–Ricci curvature; [Ol, LLY]). Let ε ∈ [0, 1]. The ε-Ollivier–Ricci curvatureκ(ε; x, y) between two nodes x, y ∈ V is defined as
W1 µεx , µεyκ(ε; x, y) := 1 −.
d(x, y)When ε = 0, κ(0; x, y) = 0 for any nodes x, y. Then, in [LLY], they defined a curvature as the (right)
limit limε↓0 κ(ε; x, y)/ε instead of simply assigning ε = 0 (Definition A.6). To confirm the existence of thislimit, the following two lemmas were shown.
Lemma A.3 ([LLY, Lemma 2.2]). For any ε ∈ [0, 1] and x, y ∈ V , we have|κ(ε; x, y)| ≤2ε.
d(x, y)Lemma A.4 ([LLY, Lemma 2.1]). For two vertices x, y, κ(ε; x, y) is concave in ε ∈ [0, 1].
The shape of κ(ε; x, y) in ε ∈ [0, 1] was later examined in more detail.
Theorem A.5 ([BCLMP, Theorem 3.4], [CK, Theorem 3.2]). For any connected, locally finite and simple
graph G = (V, E) and any of its nodes x, y, the function φ : ε 7→ κ(ε; x, y) is concave and piecewise linear on[0, 1]. Moreover, the number of its partitions is at most 3. In particular, φ(ε) := κ(ε; x, y)/ε is constant for εsufficiently close to 0.
In [BCLMP], the case of x ∼ y was shown, then in [CK] established the general case by using the
Kantorovich duality of W1 distance.
From Lemma A.3 and Lemma A.4, we can see the existence of the (right) limit limε↓0 κ(ε; x, y)/ε.
Definition A.6 (LLY–Ricci curvature; [LLY]). We define the LLY–Ricci curvature κ(x, y) between two
nodes x, y ∈ V as
κ(x, y) := limε↓0κ(ε; x, y).
ε

Although LLY–Ricci curvature κ(x, y) is a limit value, it can be obtained by linear programming ([CKLLS])
thanks to Theorem A.5. Moreover, the Graph Curvature Calculator3 has been developed to calculate LLY–
Ricci curvature of each edge of a graph by inputting node and edge information.
3 https://www.mas.ncl.ac.uk/graph-curvature/ On this page, you can select from the tabs in the lower left corner tocalculate various types of curvatures. The tab “Lin–Lu–Yau Curvature” allows you to calculate the curvatures used in thisstudy. Remark that these are calculated on unweighted graphs, and the following discussion considers weighted graphs.

30

Page 31 of 34

g-RIPS Sendai 2022, Mitsubishi-A

A.2“Ricci curvature” between the trajectory and routeAccording to Definition A.6, we introduce a “curvature” between the trajectory and each route. It is thisquantity that determines which routes the trajectory is closer to.
The denominator d(x, y) used in the Definition A.6 (i.e. Definition A.2) was the distance between thetwo vertices x, y of the target. In this setting, the target to be measured is not two vertices but two “setsof vertices”. Therefore, it is necessary to first consider what the quantity corresponding to d(x, y) should
be in this setting. In conclusion, we adopt W1 (µp , νA ) (with respect to the route A) as described in themidterm presentation. This is because W1 (µp , νA ) was quantifying the distance between the trajectory androute using only the location information of the trajectory. We next need to consider the amount of thenumerator of Definition A.2. In the idea of Definition A.6, they perturbed the Dirac measures4 at x, y byε along the graph structure. In the setting of Wasserstein method (§ 5.2), the vertices of the local road
network are only the trajectory and the divided points {a1 , . . . , am , b1 , . . . , bm } of each route A and B, andthus it is a problem of how to introduce the edges. Note that we define the transition probabilities forroutes as combinatorial graphs. Then, we consider introducing weighted edges using the input speed
and direction at each trajectory point and we set up and defined them as described in § 5.3.
Now that ε-perturbation measures have been defined, we can define the “curvature” between the trajectoryand routes as in Definition A.2 and Definition A.6. In the following, we will write about route A only, since
the definition is the same for both routes A and B.
Definition A.7 (ε-Ricci curvature between the trajectory and route). Let ε ∈ [0, 1]. The ε-Ricci curvatureκ(ε; p, A) between the trajectory p and route A is defined as:ε
W1 (µεp,A , νA).
W1 (µp , νA )

κ(ε; p, A) := 1 −(A.1)

We next want to prove the corresponding properties for Lemma A.3 and Theorem A.5. However, we prove
only the concavity of the function ε 7→ κ(ε; p, A) corresponding to Lemma A.4 (Proposition A.9) because thepiecewise linearity of h like Theorem A.5 is not yet clearly known (see also § 5.4).
Proposition A.8 (ε-boundedness of κ(ε; p, A)). For any ε ∈ [0, 1], we have|κ(ε; p, A)| ≤2ε.
W1 (µp , νA )Proof. By the triangle inequality of W1 , we obtain
ε
0
0
ε
0
W1 (µεp,A , νA
) ≤ W1 (µεp,A , µ0p ) + W1 (µ0p,A , νA
) + W1 (νA, νA) = W1 (µ0p,A , νA
) + 2ε = W1 (µp , νA ) + 2ε,
ε
0
0
ε
0
W1 (µεp,A , νA
) ≥ W1 (µ0p,A , νA
) − W1 (µ0p,A , µεp,A ) − W1 (νA, νA) = W1 (µ0p,A , νA
) − 2ε = W1 (µp , νA ) − 2ε.

This implies the following:−ε
W1 (µεp,A , νA)2ε2ε
≤ κ(ε; p, A) := 1 −≤.
W1 (µp , νA )W1 (µp , νA )W1 (µp , νA )Proposition A.9 (Concavity of κ(ε; p, A)). The function h : [0, 1] ∋ ε 7→ κ(ε; p, A) ∈ R is concave.
Proof. Let ε1 , ε2 , ε3 be 0 ≤ ε1 < ε2 < ε3 ≤ 1 and t := (ε3 − ε2 )/(ε3 − ε1 ). Then, ε2 = tε1 + (1 − t)ε3 holds.
We show thatκ(ε2 ; p, A) ≥ tκ(ε1 ; p, A) + (1 − t)κ(ε3 ; p, A)(A.2)

holds. First, we show the following.
4 Here, notice that d(x, y) can be transformed to d(x, y) = W (δ , δ ) and µ0 = δ , µ0 = δ . This means that the fraction
x
y
1 x y
x
y
Definition A.2: W1 (µεx , µεy )/W1 (δx , δy ) measures the fundamental probability measures δx , δy in the denominator and the εperturbations µεx , µεy of them in the numerator, with W1 between them, respectively.

31

Page 32 of 34

g-RIPS Sendai 2022, Mitsubishi-Aεεjand νAj , with respect to j = 1, 3. Then, π2 :=Claim Let πj be the optimal coupling between µp,Aε2
ε2
tπ1 + (1 − t)π3 is a coupling between µp,A and νA .
Proof of Claim
By the definition of π2 , we obtain
X
X
X
ε1
ε3
π2 (v1 , v2 ) = t
π1 (v1 , v2 ) + (1 − t)
π3 (v1 , v2 ) = tνA(v2 ) + (1 − t)νA
(v2 )

v1 ∈V

v1 ∈V

v1 ∈V

1 X ε1
1 X ε3=t·νa (v2 ) + (1 − t) ·νa (v2 ),mma∈V (A)a∈V (A)X
X
X
13π2 (v1 , v2 ) = t
π1 (v1 , v2 ) + (1 − t)
π3 (v1 , v2 ) = tµεp,A(v1 ) + (1 − t)µεp,A
(v1 )
v2 ∈V

v2 ∈V

(A.3)

v2 ∈V

1 X ε1
1 X ε3=t·µp,A (v1 ) + (1 − t) ·µ (v1 ).
n p∈p
n p∈p p,A

(A.4)

ε2
2
It is sufficient to check that the right hand side of (A.3) and (A.4) coincide with νA(v2 ) and µεp,A(v1 ),
respectively.
(A.3) ( i ) In the case of v2 ∈ V (A): it holds that

(A.3) = t ·1
1
1
1 X ε2ε2(v2 ).
(1 − ε1 ) + (1 − t) · (1 − ε3 ) = (1 − ε2 ) =νa (v2 ) = νAmmmma∈V (A)(ii) In the case of v2 ∈ p: it holds that






1
1
1
1
1
1
1 X ε2ε2(A.3) = t ·ε1 ·· m + (1 − t) ·
ε3 ·
·m=
ε2 ·
·m=
(v2 ).
νa (v2 ) = νA
m
n
m
n
m
n
m
a∈V (A)ε2(v2 ).
(iii) In the case of v2 ∈ V (B): it holds that (A.3) = 0 = νA(A.4) (iv) In the case of v1 ∈ p: it holds that

(A.4) = t ·1
1
1
1 X ε22(1 − ε1 ) + (1 − t) · (1 − ε3 ) = (1 − ε2 ) =µ (v1 ) = µεp,A(v1 ).
nnnn p∈p p,A(v) In the case of v1 ∈ V (A): it holds that!
!
!
1X
1 + wp (v1 )1X
1 + wp (v1 )1X
1 + wp (v1 )(A.4) = t ·ε1 ·+ (1 − t) ·
ε3 ·
=
ε2 ·
n p∈p
m+1
n p∈p
m+1
n p∈p
m+1
X
1
2
µε2 (v2 ) = µεp,A(v2 ).
=
n p∈p p,A(vi) In the case of v1 ∈ V (B): it holds that!
!
!
1X
wp (v1 )1X
wp (v1 )1X
wp (v1 )(A.4) = t ·ε1 ·+ (1 − t) ·
ε3 ·
=
ε2 ·
n p∈p
m+1
n p∈p
m+1
n p∈p
m+1
1 X ε22=µ (v2 ) = µεp,A(v2 ).
n p∈p p,A
This concludes the proof of Claim. ■From Claim, we obtain
X
X
X
2
W1 (µεp,A, νA ) ≤π2 (v1 , v2 )d(v1 , v2 ) = t
π1 (v1 , v2 )d(v1 , v2 ) + (1 − t)
π3 (v1 , v2 )d(v1 , v2 )v1 ,v2 ∈Vv1 ,v2 ∈V
ε1
ε3= tW1 (µp,A , νA ) + (1 − t)W1 (µp,A , νA ).

32

v1 ,v2 ∈VPage 33 of 34

g-RIPS Sendai 2022, Mitsubishi-A

This yields the following:κ(ε2 ; p, A) := 1 −ε2
2W1 (µεp,A, νA)

W1 (µp , νA )

≥t 1−

ε1
1
W1 (µεp,A
, νA
)

!

W1 (µp , νA )+ (1 − t) 1 −ε33W1 (µεp,A
, νA
)

!

W1 (µp , νA )= tκ(ε1 ; p, A) + (1 − t)κ(ε3 ; p, A).
This proves (A.2).
Definition A.10 (Ricci curvature between the trajectory and route). The existence of the value of thelimit limε↓0 κ(ε; p, A)/ε is guaranteed by Proposition A.8 and Proposition A.9. We define this value as Ricci
curvature κ(p, A) between the trajectory and route A:
κ(p, A) := limε↓0κ(ε; p, A).
εObservation A.11. We determined the route with the smallest (5.5) to be the output as the true route.
Note that if (5.8) holds, then we obtainκ(ε; p, A) := 1 −εεW1 (µεp,A , νA)
W1 (µεp,A , νA)
<1−=: κ(ε; p, B).
W1 (µp , νA )W1 (µp , νA )Therefore, the Wasserstein method is a method that outputs the route with the largest Ricci curvature withthe trajectory as the true route.
Remark A.12. As mentioned in § 5.4, we have not shown the piecewise linearity like Theorem A.5 in
LLY–Ricci curvature yet.
Remark A.13. This value may require modifications in the way the weights (Definition 5.13) are assigned,
as we have not yet calculated the concrete examples.

References[BCLMP]

D. P. Bourne, D. Cushing, S. Liu, F. Münch and N. Peyerimhoff, Ollivier-Ricci idleness functions ofgraphs, SIAM J. Discrete Math. 32(2) (2018), 1408–1424.

[BK]

D. Bernstein and A. Kornhauser, An introduction to map-matching for personal navigation assistants,(1996).

[CK]

D. Cushing and S. Kamtue, Long-scale Ollivier Ricci curvature of graphs, Anal. Geom. Metr. Spaces.
7(1) (2019), 22–44.

[CKLLS]

D. Cushing, R. Kangaslampi, V. Lipiäinen, S. Liu and G. W. Stagg, The Graph Curvature Calculator
and the curvatures of cubic graphs, Exp. Math. (2019), 13pp.
https://doi.org/10.1080/10586458.2019.1660740[CXHZ]

P. Chao, Y. Xu, W. Hua and X. Zhou, A survey on map-matching algorithms, Springer, Cham.
(2020).

[FG]

A. Figalli and F. Glaudo, An invitation to optimal transport, Wasserstein distances, and gradient
flows, EMS Textbooks in Mathematics. EMS Press (2021).

[GH]

GraphHopper Github Repository, https://github.com/graphhopper/graphhopper.

[JL]

J. Jost and S. Liu, Ollivier’s Ricci Curvature, Local Clustering and Curvature-Dimension Inequalities
on Graphs, Discrete Comput. Geom. 51(2) (2014), 300–322.

[Jo]

J. Jost, Riemannian geometry and geometric analysis, 7th edn. Springer, Berlin (2017).

[KCMMN]

M. Kubic̆ka, A. Cela, P. Moulin, H. Mountier and S. I. Niculescu, Dataset for testing and training of
map-matching algorithms, In 2015 IEEE Intelligent Vehicles Symposium (IV), 1088–1093 (2015).

[LLY]

Y. Lin, L. Lu and S.-T. Yau, Ricci curvature of graphs, Tohoku Math. J. (2) 63(4) (2011), 605–627.

[LSG]

X. Li, Y. Shi and I. Gutman, Graph energy, Springer, New York. (2012).

33

Page 34 of 34

g-RIPS Sendai 2022, Mitsubishi-A[NK]

P. Newson and J. Krumm, Hidden Markov map matching through noise and sparseness, In Proceedings of the 17th ACM SIGSPATIAL international conference on advances in geographic information
systems, 336–343 (2009).

[NLGGS]

C.-C. Ni, Y.-Y. Lin, J. Gao, X. D. Gu and E. Saucan, Ricci Curvature of the Internet Topology, In:
IEEE Conference on Computer Communications. (2015), 2758–2766.

[NLLG]

C.-C. Ni, Y.-Y. Lin, F. Luo and J. Gao, Community Detection on Networks with Ricci Flow, Sci Rep
9, 9984 (2019), 12pp.
https://doi.org/10.1038/s41598-019-46380-9[Ol]

Y. Ollivier, Ricci curvature of Markov chains on metric spaces, J. Funct. Anal. 256(3) (2009), 810–
864.

[QON]

M. A. Quddus, W. Y. Ochieng and R. B. Noland. Current map-matching algorithms for transportapplications: State-of-the art and future research directions, Transportation research part c: Emerging
technologies, 15(5), 312–328 (2007).

[QOZN]M. A. Quddus, W. Y. Ochieng, L. Zhao and R. B. Noland, A general map-matching algorithm fortransport telematics applications, GPS Solutions 7, 157–167 (2003).

[Sa]

F. Santambrogio, Optimal transport for applied mathematicians. Calculus of variations, PDEs, andmodeling, Progress in Nonlinear Differential Equations and their Applications, Birkhäuser/Springer,
Cham. (2015).

[Vi]

C. Villani, Optimal Transport. Old and New, Springer-Verlag. (2008).

[YG]

C. Yang and G. Gidófalvi, Fast map matching, an algorithm for integrating a hidden Markov modelwith precomputation, International Journal of Geographical Information Science. Taylor & Francis,
32(3), 547–570 (2018).

[YCWXCLMD] F. Yu, H. Chen, X. Wang, W. Xian, Y. Chen, F. Liu, V. Madhavan and T. Darrell, BDD100K:
A Diverse Driving Dataset for Heterogeneous Multitask Learning, In Proceedings of the IEEE/CVF
conference on computer vision and pattern recognition, 2636–2645 (2020).

34