On the Computation of Distances between 2D-Histograms by Minimum - - PowerPoint PPT Presentation

1/25 Introduction Formulations Results Conclusions References

On the Computation of Distances between 2D-Histograms by Minimum Cost Flows

Stefano Gualandi Federico Bassetti and Marco Veneroni

Università di Pavia, Dipartimento di Matematica

email: stefano.gualandi@unipv.it twitter: @famo2spaghi

2/25 Introduction Formulations Results Conclusions References

Let X = {x1, . . . , xn} and Y = {y1, . . . , ym} be two discrete spaces. Let µ = (µ(x1), . . . , µ(xn)) probability vector on X ν = (ν(y1), . . . , ν(ym)) probability vector on Y c : X × Y → R+ a cost function.

2/25 Introduction Formulations Results Conclusions References

Let X = {x1, . . . , xn} and Y = {y1, . . . , ym} be two discrete spaces. Let µ = (µ(x1), . . . , µ(xn)) probability vector on X ν = (ν(y1), . . . , ν(ym)) probability vector on Y c : X × Y → R+ a cost function. Definition 1 (Kantorovich-Rubinshtein Functional (San15; LS17)) The Kantorovich-Rubinshtein functional in the discrete setting is the following LP problem (a special case of the Hitchcock Problem) Wc(µ, ν) = min

• x∈X
• y∈Y

c(x, y)π(x, y) s.t.

• y∈Y

π(x, y) = µ(x) ∀x ∈ X

• x∈X

π(x, y) = ν(y) ∀y ∈ Y π(x, y) ≥ 0.

3/25 Introduction Formulations Results Conclusions References

Definition 2 (Wasserstein distance (San15)) When X = Y and c(x, y) = dp(x, y), where d is a ground distance on X, we define the Wasserstein distance of order p as: Wp(µ, ν) := Wdp(µ, ν)min(1/p,1) which is a distance on the simplex of probability vectors on X.

3/25 Introduction Formulations Results Conclusions References

Definition 2 (Wasserstein distance (San15)) When X = Y and c(x, y) = dp(x, y), where d is a ground distance on X, we define the Wasserstein distance of order p as: Wp(µ, ν) := Wdp(µ, ν)min(1/p,1) which is a distance on the simplex of probability vectors on X. OUR CONTRIBUTION: For Wasserstein distances of order p = 1 with the following ground distances d1(x, y) = ||x − y||1 d2(x, y) = ||x − y||2 d∞(x, y) = ||x − y||∞

3/25 Introduction Formulations Results Conclusions References

Definition 2 (Wasserstein distance (San15)) When X = Y and c(x, y) = dp(x, y), where d is a ground distance on X, we define the Wasserstein distance of order p as: Wp(µ, ν) := Wdp(µ, ν)min(1/p,1) which is a distance on the simplex of probability vectors on X. OUR CONTRIBUTION: For Wasserstein distances of order p = 1 with the following ground distances d1(x, y) = ||x − y||1 → Exact method (easy) d2(x, y) = ||x − y||2 d∞(x, y) = ||x − y||∞

3/25 Introduction Formulations Results Conclusions References

Definition 2 (Wasserstein distance (San15)) When X = Y and c(x, y) = dp(x, y), where d is a ground distance on X, we define the Wasserstein distance of order p as: Wp(µ, ν) := Wdp(µ, ν)min(1/p,1) which is a distance on the simplex of probability vectors on X. OUR CONTRIBUTION: For Wasserstein distances of order p = 1 with the following ground distances d1(x, y) = ||x − y||1 → Exact method (easy) d2(x, y) = ||x − y||2 d∞(x, y) = ||x − y||∞ → Exact method (easy)

3/25 Introduction Formulations Results Conclusions References

Definition 2 (Wasserstein distance (San15)) When X = Y and c(x, y) = dp(x, y), where d is a ground distance on X, we define the Wasserstein distance of order p as: Wp(µ, ν) := Wdp(µ, ν)min(1/p,1) which is a distance on the simplex of probability vectors on X. OUR CONTRIBUTION: For Wasserstein distances of order p = 1 with the following ground distances d1(x, y) = ||x − y||1 → Exact method (easy) d2(x, y) = ||x − y||2 → Exact and approximation methods (tricky) d∞(x, y) = ||x − y||∞ → Exact method (easy)

3/25 Introduction Formulations Results Conclusions References

Definition 2 (Wasserstein distance (San15)) When X = Y and c(x, y) = dp(x, y), where d is a ground distance on X, we define the Wasserstein distance of order p as: Wp(µ, ν) := Wdp(µ, ν)min(1/p,1) which is a distance on the simplex of probability vectors on X. OUR CONTRIBUTION: For Wasserstein distances of order p = 1 with the following ground distances d1(x, y) = ||x − y||1 → Exact method (easy) d2(x, y) = ||x − y||2 → Exact and approximation methods (tricky) d∞(x, y) = ||x − y||∞ → Exact method (easy) All the methods rely on solving an Uncapacitated Min Cost Flow problem (but on different networks)

4/25 Introduction Formulations Results Conclusions References

2D Histograms

2D Histograms 2D histograms can be seen as discrete measures on a finite set of points in R2. To represent 2D histograms with N × N equally spaced bins, we take X = LN := {i = (i1, i2) : i1 = 0, . . . , N − 1, i2 = 0, . . . , N − 1} One can think of each point (i1, i2) as the center of a bin.

5/25 Introduction Formulations Results Conclusions References

Computing Wasserstein Distance of order 1

Computing W1 distances between 2D histograms with n = N2 bins reduces to an Uncapacitated Min Cost Flow problem on a bipartite graph with 2n nodes and n2 arcs, and can be solved in O(n3 log n) time (Orl93; GTT89). Related to Earth Mover Distance (RTG00; PW09; Cut13)

5/25 Introduction Formulations Results Conclusions References

Computing Wasserstein Distance of order 1

Computing W1 distances between 2D histograms with n = N2 bins reduces to an Uncapacitated Min Cost Flow problem on a bipartite graph with 2n nodes and n2 arcs, and can be solved in O(n3 log n) time (Orl93; GTT89). Related to Earth Mover Distance (RTG00; PW09; Cut13)

5/25 Introduction Formulations Results Conclusions References

Computing Wasserstein Distance of order 1

Computing W1 distances between 2D histograms with n = N2 bins reduces to an Uncapacitated Min Cost Flow problem on a bipartite graph with 2n nodes and n2 arcs, and can be solved in O(n3 log n) time (Orl93; GTT89). Related to Earth Mover Distance (RTG00; PW09; Cut13) ground distance: d2(i, j) =

• (i1 − j1)2 + (i2 − j2)2

6/25 Introduction Formulations Results Conclusions References

W1 distances: Recomputation from the Literature (SSG17)

Best results with the Network Simplex of the Lemon Graph Library v1.3.1 Similar results with CPLEX 12.7 and Gurobi 7.0 32x32 vertices: 2 048 arcs: 1 048 576 runtime: 0.4 s

6/25 Introduction Formulations Results Conclusions References

W1 distances: Recomputation from the Literature (SSG17)

Best results with the Network Simplex of the Lemon Graph Library v1.3.1 Similar results with CPLEX 12.7 and Gurobi 7.0 32x32 64x64 vertices: 2 048 8 192 arcs: 1 048 576 16 777 216 runtime: 0.4 s 10.9 s

6/25 Introduction Formulations Results Conclusions References

W1 distances: Recomputation from the Literature (SSG17)

Best results with the Network Simplex of the Lemon Graph Library v1.3.1 Similar results with CPLEX 12.7 and Gurobi 7.0 32x32 64x64 128x128 vertices: 2 048 8 192 32 768 arcs: 1 048 576 16 777 216 268 435 456 runtime: 0.4 s 10.9 s

• ut-of-memory

6/25 Introduction Formulations Results Conclusions References

W1 distances: Recomputation from the Literature (SSG17)

Best results with the Network Simplex of the Lemon Graph Library v1.3.1 Similar results with CPLEX 12.7 and Gurobi 7.0 32x32 64x64 128x128 256x256 vertices: 2 048 8 192 32 768 131 072 arcs: 1 048 576 16 777 216 268 435 456 4 294 967 296 runtime: 0.4 s 10.9 s

• ut-of-memory
• ut-of-memory

6/25 Introduction Formulations Results Conclusions References

W1 distances: Recomputation from the Literature (SSG17)

Best results with the Network Simplex of the Lemon Graph Library v1.3.1 Similar results with CPLEX 12.7 and Gurobi 7.0 32x32 64x64 128x128 256x256 512x512 vertices: 2 048 8 192 32 768 131 072 524 288 arcs: 1 048 576 16 777 216 268 435 456 4 294 967 296 68 719 476 736 runtime: 0.4 s 10.9 s

• ut-of-memory
• ut-of-memory
• ut-of-memory

7/25 Introduction Formulations Results Conclusions References

Min Cost Flow Formulation on Kn

For computing Wdh distances between µ, ν ∈ X, we define the complete flow network Kn = (V , A): Nodes: V = LN

7/25 Introduction Formulations Results Conclusions References

Min Cost Flow Formulation on Kn

For computing Wdh distances between µ, ν ∈ X, we define the complete flow network Kn = (V , A): Nodes: V = LN Arcs: A = {(i, j) | ∀i, j ∈ V , i = j}

7/25 Introduction Formulations Results Conclusions References

Min Cost Flow Formulation on Kn

For computing Wdh distances between µ, ν ∈ X, we define the complete flow network Kn = (V , A): Nodes: V = LN Arcs: A = {(i, j) | ∀i, j ∈ V , i = j} Arc costs: cij = dij = ||i − j||h , ∀(i, j) ∈ A Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V

7/25 Introduction Formulations Results Conclusions References

Min Cost Flow Formulation on Kn

For computing Wdh distances between µ, ν ∈ X, we define the complete flow network Kn = (V , A): Nodes: V = LN Arcs: A = {(i, j) | ∀i, j ∈ V , i = j} Arc costs: cij = dij = ||i − j||h , ∀(i, j) ∈ A Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V Let fij be the flow sent from bin i to bin j, we denote by FKn,c(b) = min

• (i,j)∈A

cijfij, s.t.

• ij∈A

fij −

• ji∈A

fji = bi, fij ≥ 0

8/25 Introduction Formulations Results Conclusions References

Proposition 2.1 Let µ and ν be two probability vectors on X. (a) If c(x, x) = 0 for every x ∈ X, and b = µ − ν, then Wc(µ, ν) ≥ FKn,c(b). (b) If, in addition, c(x, y) ≤ c(x, z) + c(z, y) for every x, y, z in X, then Wc(µ, ν) = FKn,c(b).

8/25 Introduction Formulations Results Conclusions References

Proposition 2.1 Let µ and ν be two probability vectors on X. (a) If c(x, x) = 0 for every x ∈ X, and b = µ − ν, then Wc(µ, ν) ≥ FKn,c(b). (b) If, in addition, c(x, y) ≤ c(x, z) + c(z, y) for every x, y, z in X, then Wc(µ, ν) = FKn,c(b). For any directed graph G ⊂ Kn one has that FG,c(b) ≥ FKn,c(b), and EG(b) := FG,c(b) − FKn,c(b) FKn,c(b) ≥ 0. Ideally, we look for a subgraph G ⊂ Kn for which EG(b) = 0.

9/25 Introduction Formulations Results Conclusions References

d1 Ground Distance

We define the flow network G0 = (V , A0) as follows: Nodes: V = LN V0 = {(1, 0), (0, 1), (−1, 0), (0, −1)}

9/25 Introduction Formulations Results Conclusions References

d1 Ground Distance

We define the flow network G0 = (V , A0) as follows: Nodes: V = LN V0 = {(1, 0), (0, 1), (−1, 0), (0, −1)} Arcs: A0 = {(i, i + j) : i ∈ LN, j ∈ V0, i + j ∈ LN} Arc costs: cij = d1(i, j) = |i1 − j1| + |i2 + j2|= 1

9/25 Introduction Formulations Results Conclusions References

d1 Ground Distance

We define the flow network G0 = (V , A0) as follows: Nodes: V = LN V0 = {(1, 0), (0, 1), (−1, 0), (0, −1)} Arcs: A0 = {(i, i + j) : i ∈ LN, j ∈ V0, i + j ∈ LN} Arc costs: cij = d1(i, j) = |i1 − j1| + |i2 + j2|= 1

9/25 Introduction Formulations Results Conclusions References

d1 Ground Distance

We define the flow network G0 = (V , A0) as follows: Nodes: V = LN V0 = {(1, 0), (0, 1), (−1, 0), (0, −1)} Arcs: A0 = {(i, i + j) : i ∈ LN, j ∈ V0, i + j ∈ LN} Arc costs: cij = d1(i, j) = |i1 − j1| + |i2 + j2|= 1 Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V Let FG0,d1(b) be the min cost flow value.

10/25 Introduction Formulations Results Conclusions References

d1 Ground Distance

We define the flow network G0 = (V , A0) as follows: Nodes: V = LN V0 = {(1, 0), (0, 1), (−1, 0), (0, −1)} Arcs: A0 = {(i, i + j) : i ∈ LN, j ∈ V0, i + j ∈ LN} Arc costs: cij = d1(i, j) = |i1 − j1| + |i2 + j2|= 1 Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V Let FG0,d1(b) be the min cost flow value. Proposition 2.2 EG0,d1 = 0 and hence FG0,d1(b) = FK,d1(b). In particular, Wd1(µ, ν) = FG0,d1(µ − ν) for every couple of probability measures µ and ν on LN.

10/25 Introduction Formulations Results Conclusions References

d1 Ground Distance

We define the flow network G0 = (V , A0) as follows: Nodes: V = LN V0 = {(1, 0), (0, 1), (−1, 0), (0, −1)} Arcs: A0 = {(i, i + j) : i ∈ LN, j ∈ V0, i + j ∈ LN} Arc costs: cij = d1(i, j) = |i1 − j1| + |i2 + j2|= 1 Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V Let FG0,d1(b) be the min cost flow value. NOTE: |A0| = O(4n) Proposition 2.2 EG0,d1 = 0 and hence FG0,d1(b) = FK,d1(b). In particular, Wd1(µ, ν) = FG0,d1(µ − ν) for every couple of probability measures µ and ν on LN.

11/25 Introduction Formulations Results Conclusions References

d∞ Ground Distance

We define the flow network G1 = (V , A1) as follows: Nodes: V = LN V1 = V0 ∪ {(1, 1), (−1, 1), (−1, −1), (1, −1)}

11/25 Introduction Formulations Results Conclusions References

d∞ Ground Distance

We define the flow network G1 = (V , A1) as follows: Nodes: V = LN V1 = V0 ∪ {(1, 1), (−1, 1), (−1, −1), (1, −1)} Arcs: A1 = {(i, i + j) : i ∈ LN, j ∈ V1, i + j ∈ LN} Arc costs: cij = d∞(i, j) = max{|i1 − j1|, |i2 − j2|}= 1

11/25 Introduction Formulations Results Conclusions References

d∞ Ground Distance

We define the flow network G1 = (V , A1) as follows: Nodes: V = LN V1 = V0 ∪ {(1, 1), (−1, 1), (−1, −1), (1, −1)} Arcs: A1 = {(i, i + j) : i ∈ LN, j ∈ V1, i + j ∈ LN} Arc costs: cij = d∞(i, j) = max{|i1 − j1|, |i2 − j2|}= 1

11/25 Introduction Formulations Results Conclusions References

d∞ Ground Distance

We define the flow network G1 = (V , A1) as follows: Nodes: V = LN V1 = V0 ∪ {(1, 1), (−1, 1), (−1, −1), (1, −1)} Arcs: A1 = {(i, i + j) : i ∈ LN, j ∈ V1, i + j ∈ LN} Arc costs: cij = d∞(i, j) = max{|i1 − j1|, |i2 − j2|}= 1 Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V Let FG1,d∞(b) be the min cost flow value.

12/25 Introduction Formulations Results Conclusions References

d∞ Ground Distance

We define the flow network G1 = (V , A1) defined as follows: Nodes: V = LN V1 = V0 ∪ {(1, 1), (−1, 1), (−1, −1), (1, −1)} Arcs: A1 = {(i, i + j) : i ∈ LN, j ∈ V1, i + j ∈ LN} Arc costs: cij = d∞(i, j) = max{|i1 − j1|, |i2 − j2|}= 1 Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V Let FG1,d∞(b) be the min cost flow value. Proposition 2.3 EG1,d∞ = 0 and hence FG1,d∞(b) = FK,d∞(b). In particular, Wd∞(µ, ν) = FG1,d∞(µ − ν) for every couple of probability measures µ and ν on LN.

12/25 Introduction Formulations Results Conclusions References

d∞ Ground Distance

We define the flow network G1 = (V , A1) defined as follows: Nodes: V = LN V1 = V0 ∪ {(1, 1), (−1, 1), (−1, −1), (1, −1)} Arcs: A1 = {(i, i + j) : i ∈ LN, j ∈ V1, i + j ∈ LN} Arc costs: cij = d∞(i, j) = max{|i1 − j1|, |i2 − j2|}= 1 Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V Let FG1,d∞(b) be the min cost flow value. NOTE: |A1| = O(8n) Proposition 2.3 EG1,d∞ = 0 and hence FG1,d∞(b) = FK,d∞(b). In particular, Wd∞(µ, ν) = FG1,d∞(µ − ν) for every couple of probability measures µ and ν on LN.

13/25 Introduction Formulations Results Conclusions References

d2 Ground Distance: Exact Method

We define the flow network GL = (V , AL) as follows: Nodes: V = LN VL = V0 ∪

(i1, i2) ∈ Z0 × Z0 : |i1| ≤ L, |i2| ≤ L,

and |i1| and |i2| are coprime

13/25 Introduction Formulations Results Conclusions References

d2 Ground Distance: Exact Method

We define the flow network GL = (V , AL) as follows: Nodes: V = LN VL = V0 ∪

(i1, i2) ∈ Z0 × Z0 : |i1| ≤ L, |i2| ≤ L,

and |i1| and |i2| are coprime

• Arcs:

AL = {(i, i + j) : i ∈ LN, j ∈ VL, i + j ∈ LN} Arc costs: cij = d2(i, j) =

• (i1 − j1)2 + (i2 − j2)2

13/25 Introduction Formulations Results Conclusions References

d2 Ground Distance: Exact Method

We define the flow network GL = (V , AL) as follows: Nodes: V = LN VL = V0 ∪

(i1, i2) ∈ Z0 × Z0 : |i1| ≤ L, |i2| ≤ L,

and |i1| and |i2| are coprime

• Arcs:

AL = {(i, i + j) : i ∈ LN, j ∈ VL, i + j ∈ LN} Arc costs: cij = d2(i, j) =

• (i1 − j1)2 + (i2 − j2)2

Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V Let FGL,d∞(b) be the min cost flow

14/25 Introduction Formulations Results Conclusions References

d2 Ground Distance: Exact Method

We define the flow network GL = (V , AL) as follows: Nodes: V = LN VL = V0 ∪

(i1, i2) ∈ Z0 × Z0 : |i1| ≤ L, |i2| ≤ L,

and |i1| and |i2| are coprime

• Arcs:

AL = {(i, i + j) : i ∈ LN, j ∈ VL, i + j ∈ LN} Arc costs: cij = d2(i, j) =

• (i1 − j1)2 + (i2 − j2)2

Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V Let FGL,d∞(b) be the min cost flow. Proposition 2.4 For L = N − 1, we have that EGN−1,d2 = 0. In particular, Wd2(µ, ν) = FGN−1,d2(µ − ν) for every couple of probability measures µ and ν on LN.

14/25 Introduction Formulations Results Conclusions References

d2 Ground Distance: Exact Method

We define the flow network GL = (V , AL) as follows: Nodes: V = LN VL = V0 ∪

(i1, i2) ∈ Z0 × Z0 : |i1| ≤ L, |i2| ≤ L,

and |i1| and |i2| are coprime

• Arcs:

AL = {(i, i + j) : i ∈ LN, j ∈ VL, i + j ∈ LN} Arc costs: cij = d2(i, j) =

• (i1 − j1)2 + (i2 − j2)2

Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V Let FGL,d∞(b) be the min cost flow. BAD NEWS: |AN−1| = O(n2) Proposition 2.4 For L = N − 1, we have that EGN−1,d2 = 0. In particular, Wd2(µ, ν) = FGN−1,d2(µ − ν) for every couple of probability measures µ and ν on LN.

15/25 Introduction Formulations Results Conclusions References

d2 Ground Distance: Approximate Method

We define the flow network GL = (V , AL) as follows: Nodes: V = LN VL = V0 ∪

(i1, i2) ∈ Z0 × Z0 : |i1| ≤ L, |i2| ≤ L,

and |i1| and |i2| are coprime

• Arcs:

AL = {(i, i + j) : i ∈ LN, j ∈ VL, i + j ∈ LN} Arc costs: cij = d2(i, j) =

• (i1 − j1)2 + (i2 − j2)2

Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V Let FGL,d∞(b) be the min cost flow. Proposition 2.5 For L < N − 1, we have that EGN−1,d2≥0. In particular, Wd2(µ, ν)≤FGN−1,d2(µ − ν) for every couple of probability measures µ and ν on LN.

15/25 Introduction Formulations Results Conclusions References

d2 Ground Distance: Approximate Method

We define the flow network GL = (V , AL) as follows: Nodes: V = LN VL = V0 ∪

(i1, i2) ∈ Z0 × Z0 : |i1| ≤ L, |i2| ≤ L,

and |i1| and |i2| are coprime

• Arcs:

AL = {(i, i + j) : i ∈ LN, j ∈ VL, i + j ∈ LN} Arc costs: cij = d2(i, j) =

• (i1 − j1)2 + (i2 − j2)2

Flow balance: bi = µ(xi) − ν(yi), ∀i ∈ V Let FGL,d∞(b) be the min cost flow. NOTE : |AL| = O(kLn) Proposition 2.5 For L < N − 1, we have that EGN−1,d2≥0. In particular, Wd2(µ, ν)≤FGN−1,d2(µ − ν) for every couple of probability measures µ and ν on LN.

16/25 Introduction Formulations Results Conclusions References

d2 Ground Distance: Approximations

Theorem 3 EGL,d2 ≤ 1 − g(1/L), with g(1/L) :=

• 1

2 + L 2 √ 1 + L2 for every 1 ≤ L ≤ N − 1. The 2nd order power expansion of g is g(1/L) = 1 − 1 8L2 + o

1

L3

• for L → +∞.

16/25 Introduction Formulations Results Conclusions References

d2 Ground Distance: Approximations

Theorem 3 EGL,d2 ≤ 1 − g(1/L), with g(1/L) :=

• 1

2 + L 2 √ 1 + L2 for every 1 ≤ L ≤ N − 1. The 2nd order power expansion of g is g(1/L) = 1 − 1 8L2 + o

1

L3

• for L → +∞.

Indeed, the worst case error we get on Wd2(µ, ν) depends only on L:

L 1 2 3 4 5 10 15 1 − g(1/L) 7.61% 2.68% 1.29% 0.75% 0.49% 0.12% 0.06% |AL| O(8n) O(16n) O(32n) O(48n) O(80n) O(256n) O(576n)

17/25 Introduction Formulations Results Conclusions References

Definition of set VL for L ∈ {1, 2, 3, 4, 5}

L=1

17/25 Introduction Formulations Results Conclusions References

Definition of set VL for L ∈ {1, 2, 3, 4, 5}

L=2

17/25 Introduction Formulations Results Conclusions References

Definition of set VL for L ∈ {1, 2, 3, 4, 5}

L=3

17/25 Introduction Formulations Results Conclusions References

Definition of set VL for L ∈ {1, 2, 3, 4, 5}

L=4

17/25 Introduction Formulations Results Conclusions References

Definition of set VL for L ∈ {1, 2, 3, 4, 5}

L=5

18/25 Introduction Formulations Results Conclusions References

Implementation Details

Tests implemented in C++ using the Network Simplex for Minimum Cost Flow implemented in the Lemon Graph library v1.3.1 DOTmark instances (SSG17): 10 classes of 2D histograms (gray scale images), 10 instances for classes Distances computed for each pair of 2D histograms of the same class Test on a laptop, CPU i7 6700HQ, 16Gb RAM, single thread

19/25 Introduction Formulations Results Conclusions References

EXACT methods

32x32 vertices: 1 024 d1: 0.006 s arcs: 3 968 d∞: 0.015 s arcs: 7 812 d2: 0.29 s arcs: 638 692

19/25 Introduction Formulations Results Conclusions References

EXACT methods

32x32 64x64 vertices: 1 024 4 096 d1: 0.006 s 0.037 s arcs: 3 968 16 128 d∞: 0.015 s 0.048 s arcs: 7 812 32 004 d2: 0.29 s 7.63 s arcs: 638 692 10 205 236

19/25 Introduction Formulations Results Conclusions References

EXACT methods

32x32 64x64 128x128 vertices: 1 024 4 096 16 384 d1: 0.006 s 0.037 s 0.38 s arcs: 3 968 16 128 65 024 d∞: 0.015 s 0.048 s 0.32 s arcs: 7 812 32 004 129 540 d2: 0.29 s 7.63 s

• ut-of-memory

arcs: 638 692 10 205 236

19/25 Introduction Formulations Results Conclusions References

EXACT methods

32x32 64x64 128x128 256x256 vertices: 1 024 4 096 16 384 65 536 d1: 0.006 s 0.037 s 0.38 s 6.4 s arcs: 3 968 16 128 65 024 261 120 d∞: 0.015 s 0.048 s 0.32 s 3.8 s arcs: 7 812 32 004 129 540 521 220 d2: 0.29 s 7.63 s

• ut-of-memory
• ut-of-memory

arcs: 638 692 10 205 236

19/25 Introduction Formulations Results Conclusions References

EXACT methods

32x32 64x64 128x128 256x256 512x512 vertices: 1 024 4 096 16 384 65 536 262 144 d1: 0.006 s 0.037 s 0.38 s 6.4 s 156.6 s arcs: 3 968 16 128 65 024 261 120 1 046 528 d∞: 0.015 s 0.048 s 0.32 s 3.8 s 79.5 s arcs: 7 812 32 004 129 540 521 220 2 091 012 d2: 0.29 s 7.63 s

• ut-of-memory
• ut-of-memory
• ut-of-memory

arcs: 638 692 10 205 236

20/25 Introduction Formulations Results Conclusions References

d2 ground distance: APPROXIMATE method

32x32 vertices: 1 024 L = 1: 0.006 s arcs: 3 968 L = 4: 0.027 s arcs: 42 844 L = 5: 0.039 s arcs: 68 332 L = 10: 0.091 s arcs: 185 468 L = 15: 0.174 s arcs: 341 396

20/25 Introduction Formulations Results Conclusions References

d2 ground distance: APPROXIMATE method

32x32 64x64 vertices: 1 024 4 096 L = 1: 0.006 s 0.037 s arcs: 3 968 16 128 L = 4: 0.027 s 0.15 s arcs: 42 844 182 772 L = 5: 0.039 s 0.23 s arcs: 68 332 299 884 L = 10: 0.091 s 0.55s arcs: 185 468 888 572 L = 15: 0.174 s 1.11 s arcs: 341 396 1 828 628

20/25 Introduction Formulations Results Conclusions References

d2 ground distance: APPROXIMATE method

32x32 64x64 128x128 vertices: 1 024 4 096 16 384 L = 1: 0.006 s 0.037 s 0.4 s arcs: 3 968 16 128 65 024 L = 4: 0.027 s 0.15 s 1.2 s arcs: 42 844 182 772 760 540 L = 5: 0.039 s 0.23 s 1.8 s arcs: 68 332 299 884 1 254 508 L = 10: 0.091 s 0.55s 3.9 s arcs: 185 468 888 572 3 867 644 L = 15: 0.174 s 1.11 s 7.4 s arcs: 341 396 1 828 628 8 342 036

20/25 Introduction Formulations Results Conclusions References

d2 ground distance: APPROXIMATE method

32x32 64x64 128x128 256x256 vertices: 1 024 4 096 16 384 65 536 L = 1: 0.006 s 0.037 s 0.4 s 6.4 s arcs: 3 968 16 128 65 024 261 120 L = 4: 0.027 s 0.15 s 1.2 s 14.5 s arcs: 42 844 182 772 760 540 3 093 724 L = 5: 0.039 s 0.23 s 1.8 s 19.3 s arcs: 68 332 299 884 1 254 508 5 129 836 L = 10: 0.091 s 0.55s 3.9 s 44.5 s arcs: 185 468 888 572 3 867 644 16 117 244 L = 15: 0.174 s 1.11 s 7.4 s 72.4 s arcs: 341 396 1 828 628 8 342 036 35 524 628

20/25 Introduction Formulations Results Conclusions References

d2 ground distance: APPROXIMATE method

32x32 64x64 128x128 256x256 512x512 vertices: 1 024 4 096 16 384 65 536 262 144 L = 1: 0.006 s 0.037 s 0.4 s 6.4 s 157 s arcs: 3 968 16 128 65 024 261 120 1 046 528 L = 4: 0.027 s 0.15 s 1.2 s 14.5 s 178 s arcs: 42 844 182 772 760 540 3 093 724 12 478 684 L = 5: 0.039 s 0.23 s 1.8 s 19.3 s 240 s arcs: 68 332 299 884 1 254 508 5 129 836 20 744 812 L = 10: 0.091 s 0.55s 3.9 s 44.5 s 561 s arcs: 185 468 888 572 3 867 644 16 117 244 65 782 268 L = 15: 0.174 s 1.11 s 7.4 s 72.4 s 1 021 s arcs: 341 396 1 828 628 8 342 036 35 524 628 146 512 916

21/25 Introduction Formulations Results Conclusions References

Comparison with the bound given by Theorem 4

Results for histograms with N=32

|AL| runtime (avg) EGL (practice) EGL (theory) L = 2 15252 0.006 s 1.295% 2.680% L = 3 29404 0.010 s 0.428% 1.290% L = 5 68332 0.018 s 0.098% 0.490% L = 10 185468 0.045 s 0.011% 0.120% L = N − 1 638692 0.139 s

Results for histograms with N=64

|AL| runtime (avg) EGL (practice) EGL (theory) L = 2 63252 0.053 s 1.309% 2.680% L = 3 124252 0.086 s 0.466% 1.290% L = 5 299884 0.152 s 0.116% 0.490% L = 10 888572 0.341 s 0.015% 0.120% L = N − 1 10205236 3.131 s Each row gives the average over 450 istances

22/25 Introduction Formulations Results Conclusions References

Conclusions

Contribution We are able to compute Wasserstein distances of order 1 in short computation time and we scale up to 512 × 512 2D histograms Future work Implementation of uncapacitated min cost flow algorithms that exploits the network structure Generalization to histograms with d > 2 Different approaches for Wasserstein distances of order p > 1

23/25 Introduction Formulations Results Conclusions References

THANKS!! 1

1Draft paper ready, to be submitted

24/25 Introduction Formulations Results Conclusions References

[Cut13] Marco Cuturi, Sinkhorn distances: Lightspeed computation of optimal transport, Advances in neural information processing systems, 2013,

• pp. 2292–2300.

[GTT89] Andrew V Goldberg, Éva Tardos, and Robert Tarjan, Network flow algorithm, Tech. report, Cornell University Operations Research and Industrial Engineering, 1989. [LS17] Bruno Levy and Erica Schwindt, Notions of optimal transport theory and how to implement them on a computer, arXiv preprint arXiv:1710.02634 (2017). [Orl93] James B Orlin, A faster strongly polynomial minimum cost flow algorithm, Operations research 41 (1993), no. 2, 338–350. [PW09] Ofir Pele and Michael Werman, Fast and robust earth mover’s distances, Computer vision, 2009 IEEE 12th international conference on, IEEE, 2009,

• pp. 460–467.

[RTG00] Yossi Rubner, Carlo Tomasi, and Leonidas J Guibas, The earth mover’s distance as a metric for image retrieval, International journal of computer vision 40 (2000), no. 2, 99–121. [San15] Filippo Santambrogio, Optimal transport for applied mathematicians, Springer, 2015. [SSG17] Jörn Schrieber, Dominic Schuhmacher, and Carsten Gottschlich, Dotmark–A benchmark for discrete optimal transport, IEEE Access 5 (2017), 271–282.

25/25 Introduction Formulations Results Conclusions References

Complexity

Proposition 5.1 (Theorem 5.3.7 in (GTT89)) Given a network G with n nodes and m arcs, Orlin’s algorithm solves the transshipment problem in O(n log (min{n, D})S(n, m)), where D = maxi{|bi|} and S(n, m) is the complexity of finding a shortest path in G. Using Dijkstra algorithm S(n, m) = O(n log n + m).