# Efficient algorithm for Wasserstein-1 distance in graph setting

I first asked this question on math.stackexchange, but I think this question is high-level enough that is better suited here.

I'm looking for an efficient algorithm to calculate the Wasserstein-1 distance in the following setting:

Let $G = (V, E)$ be a (locally) finite, undirected graph, with associated weight function $w : E \to \mathbb{R}_{\geq 0}$. From the weights we derive the geodesic distance $d : V \times V \to \mathbb{R}_{\geq 0}$ as either $d(x, y) = 1 / w(x, y)$ when $(x, y) \in E$ or the shortest path over existing edges between $x$ and $y$ otherwise.

Also associated to the graph $G$ is a simple random walk, where in each vertex $x \in V$ the SRW has probability \begin{align*} \mathbb{P}_x(y) = \begin{cases} \frac{w(x, y)}{\sum_{z \sim x}w(x, z)} & \text{ if } y \sim x \\ 0 & \text{otherwise} \end{cases} \end{align*} to jump to vertex $y \in V$.

The Wasserstein-1 distance $\mathcal{W}_1$ between probability measures $\mathbb{P}_x$ and $\mathbb{P}_y$ is defined as: \begin{align*} \mathcal{W}_1(\mathbb{P}_x, \mathbb{P}_y) = \inf_{\gamma \in \prod(\mathbb{P}_x, \mathbb{P}_y)} \left\{\sum_{(a, b) \in V \times V} d(a, b) \gamma(a, b)\right\} \end{align*} where $\prod(\mathbb{P}_x, \mathbb{P}_y)$ denotes the space of couplings between $\mathbb{P}_x$ and $\mathbb{P}_y$.

Essentially, this is the problem of modifying one mass configuration to the other, over a cost function $d$. Thus, $\gamma(x, y)$ represents the amount of mass that is sent from vertex $x$ to vertex $y$, and $d(x, y) \gamma(x, y)$ is the cost of transporting this mass.

Since this is a discrete setting, calculating $\mathcal{W}_1$ isn't too difficult, but I want to calculate $\mathcal{W}_1$ for every edge in the graph, and for larger graphs this quickly becomes very time consuming.

My first implementation was to simply use linear programming, since the coupling properties can be used as constraints and the coupling $\gamma$ can be interpreted as a collection of variables.

The problem here is that the constraints require an enormous amount of memory: take 200 vertices, then in the worst case scenario you need 40,000 variables, and get a 40,000 by 400 matrix of constraints at least.

Another method is to model a minimal cost flow problem. I've made a greedy implementation of this problem, but the problem here is that with sorting all the paths the algorithm isn't very efficient: only $O(n^3)$.

I realise that my attempts have not been very clever, so I'd be very grateful if someone can help me with a better algorithm or point to a C++ or R package with a good implementation!

• My guess is that in the worst case scenario an exact computation will be very costly. But since the mesures you consider are supported on the neighborhood of a vertex and you restrict to adjacent vertices, for sparse enough graphs things are simpler. The measures are supported on at most $d$ vertices ($d$ the max degree); to compute distances between a neighbor of $x$ and a neighbor of $y$ you have an obvious candidate path (through $x,y$) to avoid testing paths above some length (if the weigths do not have roughly the correct order distance computation might be costly). – Benoît Kloeckner Sep 9 '17 at 15:54
• Also, if by any chance this has to do with Ollivier-Ricci curvature, it often yields better results to consider lazy random walks. – Benoît Kloeckner Sep 9 '17 at 15:56
• Actually, it does have to do with Ollivier-Ricci curvature! Could you explain why lazy random walks would give better results? – JustSomeGuy Sep 9 '17 at 16:14
• With a simple walk, the cube has curvature 0 while a lazy walk has positive curvature, implying many nice propertied. This is all well explained in Ollivier's paper. – Benoît Kloeckner Sep 10 '17 at 6:54
• Ah right, I have actually seen another paper where they assume a lazy walk as well. I'll have to look into it, thanks for the tip! As for your algorithm suggestion, I'm not sure I understand the implication, isn't this a consideration for the calculation of the distances? – JustSomeGuy Sep 11 '17 at 15:30

See this PLOS paper for algorithm and extensive survey. (Gottschlich and Schumacher, 2014)

• The funny thing is that I already tried the R package supplied by the authors when it appeared a few months ago, but I missed the 'shortlist' option with the algorithms they provided. Having read this paper I tried that option, and the results are astounding! – JustSomeGuy Sep 21 '17 at 10:43
• Just for comparison's sake: calculating one Wasserstein-1 distance for the averaged stochastic block model with 2 blocks of 250 vertices each with parameters $p_{in} = 0.12$ and $p_{out} = 0.1$ takes about 0.13 seconds with the shortlist method, while it takes about 3.22 and 3.27 second with my LP and flow implementations. – JustSomeGuy Sep 21 '17 at 10:46
• For a single realization of the SBM with the same parameters, it takes about 25 second to calculate all the Wasserstein-1 distances for the shortlist method, while it takes 51 and 117 seconds for the LP and flow implementations. The difference grows even larger for bigger graphs, while it decreases for smaller graphs (due to my parallel implementation I suspect) – JustSomeGuy Sep 21 '17 at 10:48
• So thank you @IgorRivin, you've just made my life a whole lot easier! – JustSomeGuy Sep 21 '17 at 10:49
• @Thomas glad it helped! – Igor Rivin Sep 21 '17 at 14:01

I could be wrong, but this seems to be addressed in this 2013 paper.

• You are. First, that paper deals with finding barycenters (i.e., minimization of weighted sums of distances to a given collection of measures) rather than with finding the distance per se. Second, it deals with 2-distance rather than 1-distance the OP is asking about. – R W Sep 9 '17 at 15:59
• @RW The second part of the paper deals with the regularization of the Wasserstein distance to make it more computationally tractable. – Igor Rivin Sep 9 '17 at 16:48
• @IgorRivin I'll make sure to look over the second part, thanks for the tip! – JustSomeGuy Sep 21 '17 at 10:40

To begin with, it is wrong, both historically and conceptually, to call the transportation metric Wasserstein 1-distance. Historically, because it was much earlier introduced by Kantorovich and Rubinshtein. Conceptually, because the presence of duality makes 1-distance quite different from other Wasserstein $p$-distances. More precisely, the distance between two measures $\alpha$ and $\beta$ coincides with $$\sup_h \bigl[ \langle h,\alpha \rangle - \langle h, \beta \rangle \bigr] \;,$$ where the $\sup$ is taken over all Lipschitz functions on the state space. Moreover, a coupling $\Lambda$ of the measures $\alpha$ and $\beta$ realizes the transportation distance if and only if there exists a Lipschitz function $h$ such that $h(x)-h(y)=d(x,y)$ for any $(x,y)\in\mathop{\rm supp}\Lambda$. The latter property is quite useful in concrete situations.

• This should be a comment I guess? – Federico Poloni Sep 9 '17 at 17:50
• Why? Isn't it clear enough that I suggest to look at the dual problem? – R W Sep 9 '17 at 21:35
• It seems like there is a lot of distance between this suggestion and a practical algorithm... – Federico Poloni Sep 10 '17 at 6:08
• Yeah, I kind of agree with @FedericoPoloni. In my field the focus is on Wasserstein metrics (of any order), and in order 1 calculating the metric is equivalent to solving a transportation problem. I'm aware of the duality, but I don't see how this can be utilized to derive an efficient algorithm? – JustSomeGuy Sep 11 '17 at 15:34
• My point is that order 1 metric has properties absent in other orders, and therefore it might be sensible to use these properties rather than ignoring them. – R W Sep 11 '17 at 15:52