# Earth Mover's Distance (EMD) between two Gaussians

Is there a closed-form formula for (or some kind of bound on) the EMD between $x_1\sim N(\mu_1, \Sigma_1)$ and $x_2 \sim N(\mu_2, \Sigma_2)$?

$$\DeclareMathOperator\EMD{\mathrm{EMD}} \DeclareMathOperator\E{\mathbb{E}} \DeclareMathOperator\Var{Var} \DeclareMathOperator\N{\mathcal{N}} \DeclareMathOperator\tr{\mathrm{tr}} \newcommand\R{\mathbb R}$$The earth mover's distance can be written as $$\EMD(P, Q) = \inf \E \lVert X - Y \rVert$$, where the infimum is taken over all joint distributions of $$X$$ and $$Y$$ with marginals $$X \sim P$$, $$Y \sim Q$$. This is also known as the first Wasserstein distance, which is $$W_p = \inf \left( \E \lVert X - Y \rVert^p \right)^{1/p}$$ with the same infimum.

Let $$X \sim P = \N(\mu_x, \Sigma_x)$$, $$Y \sim Q = \N(\mu_y, \Sigma_y)$$.

Lower bound: By Jensen's inequality, since norms are convex, $$\E \lVert X - Y \rVert \ge \lVert \E (X - Y) \rVert = \lVert \mu_x - \mu_y \rVert,$$ so the EMD is always at least the distance between the means (for any distributions).

Upper bound based on $$W_2$$: Again by Jensen's inequality, $$\left( \E \lVert X - Y \rVert \right)^2 \le \E \lVert X - Y \rVert^2$$. Thus $$W_1 \le W_2$$. But Dowson and Landau (1982) establish that $$W_2(P, Q)^2 = \lVert \mu_x - \mu_y \rVert^2 + \tr\left( \Sigma_x + \Sigma_y - 2 (\Sigma_x \Sigma_y)^{1/2} \right) ,$$ giving an upper bound on $$\EMD = W_1$$.

A tighter upper bound: Consider the coupling \begin{align} X &\sim \N(\mu_x, \Sigma_x) \\ Y &= \mu_y + \underbrace{\Sigma_x^{-\frac12} \left( \Sigma_x^\frac12 \Sigma_y \Sigma_x^\frac12 \right)^\frac12 \Sigma_x^{-\frac12}}_A (X - \mu_x) .\end{align} This is the map derived by Knott and Smith (1984), On the optimal mapping of distributions, Journal of Optimization Theory and Applications, 43 (1) pp 39-49 as the optimal mapping for $$W_2$$; see also this blog post. Note that $$A = A^T$$ and \begin{align} \E Y &= \mu_y + A (\E X - \mu_x) = \mu_y \\ \Var Y &= A \Sigma_x A^T \\&= \Sigma_x^{-\frac12} \left( \Sigma_x^\frac12 \Sigma_y \Sigma_x^\frac12 \right)^\frac12 \Sigma_x^{-\frac12} \Sigma_x \Sigma_x^{-\frac12} \left( \Sigma_x^\frac12 \Sigma_y \Sigma_x^\frac12 \right)^\frac12 \Sigma_x^{-\frac12} \\&= \Sigma_x^{-\frac12} \left( \Sigma_x^\frac12 \Sigma_y \Sigma_x^\frac12 \right) \Sigma_x^{-\frac12} \\&= \Sigma_y ,\end{align} so the coupling is valid.

The distance $$\lVert X - Y \rVert$$ is then $$\lVert D \rVert$$, where now \begin{align} D &= X - Y \\&= X - \mu_y - A (X - \mu_x) \\&= (I - A) X - \mu_y + A \mu_x ,\end{align} which is normal with \begin{align} \E D &= \mu_x - \mu_y \\ \Var D &= (I - A) \Sigma_x (I - A)^T \\&= \Sigma_x + A \Sigma_x A - A \Sigma_x - \Sigma_x A \\&= \Sigma_x + \Sigma_y - \Sigma_x^{-\frac12} \left( \Sigma_x^\frac12 \Sigma_y \Sigma_x^\frac12 \right)^\frac12 \Sigma_x^{\frac12} - \Sigma_x^{\frac12} \left( \Sigma_x^\frac12 \Sigma_y \Sigma_x^\frac12 \right)^\frac12 \Sigma_x^{-\frac12} .\end{align}

Thus an upper bound for $$W_1(P, Q)$$ is $$\E \lVert D \rVert$$. Unfortunately, a closed form for this expectation is surprisingly unpleasant to write down for general multivariate normals: see this question, as well as this one.

If the variance of $$D$$ ends up being spherical (e.g. if $$\Sigma_x = \sigma_x^2 I$$, $$\Sigma_y = \sigma_y^2 I$$, then the variance of $$D$$ becomes $$(\sigma_x - \sigma_y)^2 I$$), the former question gives the answer in terms of a generalized Laguerre polynomial.

In general, we have a simple upper bound for $$\E \lVert D \rVert$$ based on Jensen's inequality, derived e.g. in that first question: \begin{align} \left( \E \lVert D \rVert \right)^2 &\le \E \lVert D \rVert^2 \\&= \lVert \mu_x - \mu_y \rVert^2 + \tr\left( \Sigma_x + \Sigma_y - A \Sigma_x - \Sigma_x A \right) \\&= \lVert \mu_x - \mu_y \rVert^2 + \tr\left( \Sigma_x \right) + \tr\left( \Sigma_y \right) - 2 \tr\left( \Sigma_x^{-\frac12} \left(\Sigma_x^\frac12 \Sigma_y \Sigma_x^\frac12 \right)^\frac12 \Sigma_x^{\frac12} \right) \\&= \lVert \mu_x - \mu_y \rVert^2 + \tr\left( \Sigma_x \right) + \tr\left( \Sigma_y \right) - 2 \tr\left( \left(\Sigma_x^\frac12 \Sigma_y \Sigma_x^\frac12 \right)^\frac12 \right) \\&= W_2(P, Q)^2 .\end{align} The equality at the end is because the matrices $$\Sigma_x \Sigma_y$$ and $$\Sigma_x^\frac12 \Sigma_y \Sigma_x^\frac12 = \Sigma_x^{-\frac12} (\Sigma_x \Sigma_y) \Sigma_x^{\frac12}$$ are similar, so they have the same eigenvalues, and thus their square roots have the same trace.

This inequality is strict as long as $$\lVert D \rVert$$ isn't degenerate, which is most cases when $$\Sigma_x \ne \Sigma_y$$.

A conjecture: Maybe this closer upper bound, $$\E \lVert D \rVert$$, is tight. Then again, I had a different upper bound here for a long time that I conjectured to be tight that was in fact looser than the $$W_2$$ one, so maybe you shouldn't trust this conjecture too much. :)