Efficient algorithm for Wasserstein-1 distance in graph setting

Question

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!

Answer 1

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

Answer 2

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

Answer 3

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.