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)$?
 A: $\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. :)
