# How do you compute numerically the Earth mover's distance (EMD)?

I was trying to compute numerically (write a program) that calculated the EMD for two probability distribution $p_X$ and $q_X$. However, I had a hard time finding an outline of how to exactly compute such a distance.

I was wondering if there existed a closed form equation for EMD or if there existed an outline of an algorithm to compute it.

Just like there is very compact equation for say, the KL-divergence $D(p_X||q_X) = \sum_{x \in X}p_X(x) log\frac{p_X(x)}{q_X(x)}$

Is there one such equation for EMD or an algorithm for EMD such that it can be easily numerically computed?

What I have found so far is the following, that $EMD(p_X,q_X)$ between distributions $p_X$ and $q_X$ is:

$$EMD(P,Q) = \frac{\sum^m_{i=1}\sum^n_{j=1} f_{ij}d_{ij}}{\sum^m_{i=1}\sum^b_{j=1} f_{ij}}$$

However, to find $d_{ij}$ one needs to define some distance measure I guess and to find $f_{ij}$ one needs to solve the transportation linear programming defined in:

http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/RUBNER/emd.htm#RUBNER98A

Isn't there a solution to the transportation problem such that there is an easy way of evaluating EMD?

Also, I found a wikipedia section to compute it but, the pseudo-code was to ambiguous for me to understand how to actually compute the EMD. If someone understands that section better and can explain it, it would awesome!

http://en.wikipedia.org/wiki/Earth_mover's_distance

or you can just see the pseudo-code right here:

The Pseudo-code/wikipedia section:

If the domain D is discrete, the EMD can be computed by solving an instance transportation problem, which can be solved by the so-called Hungarian algorithm. In particular, if D is a one-dimensional array of "bins" the EMD can be efficiently computed by scanning the array and keeping track of how much dirt needs to be transported between consecutive bins. For example:

$EMD_0 = 0 \\$

$EMD_{i+1} = ( A_i + EMD_i ) - B_i \\$

$TotalDistance = \sum | EMD_i | \\$

I guess what I don't understand is what $A_i$ and $B_i$ are and when the algorithm even stops running. Its just to unclear to me what its doing and I guess I don't understand EMD well enough to derive it myself (If I could I would!)

As an extension to my question, is it a problem if the sample space for the probability distributions is infinite?

I probably won't accept answer that are not as general as possible, but an answer clarifying what wikipedia article is trying to say, would definitively get an up vote.

• It is not entirely clear to me whether you are looking for a solution for general metric spaces, or for the special case of subsets of $\mathbb{Z}$ with the absolute difference metric. The algorithm from Wikipedia applies only to the latter case. With some work it can be extended to suitable subsets of $\mathbb{R}$, but not much beyond that. – Niels J. Diepeveen Mar 20 '14 at 17:25
• I am not 100% sure if I understand everything you said but the only thing I want to know is, given two probability distributions functions for two different distributions, how do I compute the EMD distance between them. Does that help? (sorry for my ignorance on the topic if my clarification didn't help, I will keep trying until it makes sense :) ) – Charlie Parker Mar 20 '14 at 19:55
• The question is: probability distributions on what space? What is the geometry of the problem? For arbitrary distance functions I am fairly sure that you will just have to solve the corresponding transportation problem. The algorithm given in Wikipedia only works for points that are all in a straight line. – Niels J. Diepeveen Mar 20 '14 at 20:33
• An answer specifying how to calculate it for different cases is an acceptable answer, if it makes sense of course. – Charlie Parker Mar 22 '14 at 9:47

The Wikipedia algorithm is correct, but a little terse. A simple example should clear it up if we stick to actual earth moving. Histogram $A$ is three piles of dirt (earth) in a row: $A_0$, $A_1$, and $A_2$, and histogram $B$ is three piles of dirt $B_0$, $B_1$, and $B_2$ in a row.

Suppose the number of shovelfuls in each dirt pile is as follows.

$A_0 = 3$, $A_1 = 2$, $A_2 = 1$

$B_0 = 1$, $B_1 = 2$, $B_2 = 3$

So to make $A$ look like $B$, we really need to take two shovelfuls of dirt from $A_0$ and put them in $A_2$; i.e., move two shovelfuls two units of distance: that's four units of work. The algorithm has us move the two shovelfuls from $A_0$ to $A_1$, and then two from $A_1$ to $A_2$, using $EMD_{i+1} = (A_i + EMD_i) - B_i$.

$$EMD_0 = 0 \\ EMD_1 = (3 + 0) - 1 = 2 \\ EMD_2 = (2 + 2) - 2 = 2 \\ EMD_3 = (1 + 2) - 3 = 0 \\ \Rightarrow EMD = 2 + 2 + 0 = 4$$

• Can this algorithm be used for 1D real numbers that sum up to 1? – Asterisk Jan 19 '18 at 8:59

In the discrete algorithm you quote, $A_i$ seems to be the amount "of earth" at position $i$ in the first distribution and $B_i$ the amount "of earth" at position $i$ in the second distribution, with $i$ taking values in the non-negative integers. The algorithm as stated needs just one pass, and so long as the two total amounts "of earth" are the same and the ranges of the distributions are bounded, the algorithm will terminate.

For two probability distributions on the real numbers, where the total amounts "of probability" are each $1$ by definition, suppose their cumulative distribution functions are $F_a(x)$ and $F_b(x)$. Then the Earth Mover Distance becomes $$\int_{x=-\infty}^{\infty} |F_a(x) - F_b(x)|\, dx .$$

If the two distributions have densities $f_a(x)$ and $f_b(x)$ this is equivalent to $$\int_{x=-\infty}^{\infty} \left|\int_{y=-\infty}^x(f_a(y) - f_b(y))\,dy\,\right|\, dx.$$

As Niels Diepeveen says in the comments, this changes with distributions on other metric spaces, but I get the impression you are not asking about that.

In general, yes you need to defined a "cost" function $c(x,y)$ and then do some calculation. The continuous case ($x\in\mathbb{R}, y\in\mathbb{R}$) is not so straightforward except in certain cases where the formulae have already been worked out. Generally continuous case requires some numeric routines to build out the flows.

For the finite, discrete case, as mentioned in the post, you can setup a linear program (LP) to solve the Earth Moving Distance (EMD) problem. If $c_{ij}=c(x_i, y_j)$ denotes the cost and $A$ your matrix of constraints and $\vec{b}$ your vector of the marginal "earths" (the stuff you move from $X \to Y$). Then you setup the LP where your variables optimizing over are $p_{ij}$ and the LP is written:

$minimize_{p}\{\vec{c}^T\vec{p} | A\vec{p} = \vec{b} \}$

For a simple case of 2 probability measures on discrete marginal spaces I've worked this at in a post over on cross-validated (second half of the post is on EMD). The neat thing about this solution is that I show how the EMD is related to the Total Variation (TV) distance of probability measures. The TV distance requires a specific cost function $c(x_i, y_i) =\mathbb{1}(y_i\neq x_i)$. The solution to the EMD and the TV distance are the same in this case. On that linked post there is also an example set of R code to do the calculation. You can set this up on most scientific/mathematical computing languages that have support for LP solutions.