# Simplify Wasserstein distance between Gaussians with binary cost function

Let $$\mu_1$$ and $$\mu_2$$ be 1D gaussian distributions with means $$m_1$$ and $$m_2$$ respectively and common variance $$\sigma$$. Let $$\Omega$$ be a closed subset of $$\mathbb R^2$$, and consider the cost function $$c_\Omega:\mathbb R \times \mathbb R \rightarrow \{0, 1\}$$ defined by $$c_\Omega(x',x) = 1$$ if $$(x', x) \in \Omega; c(x',x) = 0$$ otherwise. Further suppose that $$\Omega$$ is symmetric in the sense that $$(x',x) \in \Omega \iff (x,x') \in \Omega$$.

Now, consider the Wasserstein distance: $$W_c(\mu_1,\mu_2) := \min_{\gamma \in \Pi(\mu_1,\mu_2)}\int_{(x',x) \in \mathbb R^2}c_\Omega(x',x)d\gamma(x',x) = \min_{\gamma \in \Pi(\mu_1,\mu_2)}\int_{(x',x) \in \Omega} d\gamma(x',x),$$ where $$\Pi(\mu_1,\mu_2)$$ is the set of all couplings of $$\mu_1$$ and $$\mu_2$$.

# Question

• Is there an analytic formula for $$W_c(\mu_1,\mu_2)$$ (perhaps involving Gaussian integrals) ?

• Same question for the special case $$\Omega = \{(x',x) \in \mathbb R^2 \mid |x'-x| \ge \alpha\}$$, for some $$\alpha > 0$$.

• I wrote an answer to you second question below. Regarding the first one, I doubt there is an answer with no further assumptions on $\Omega$. At the very least, one should assume that $\Omega = \{(x_1, x_2) : |x_1 - x_2| \geqslant \phi(x_1 + x_2)\}$ for some $1$-Lipschitz $\phi$. – Mateusz Kwaśnicki Jul 30 '19 at 21:09

This is an answer to the second question, for the particular choice of $$\Omega$$. In fact, in order to make the notation slightly simpler, let us assume that $$\Omega = \{(x_1, x_2) : |x_1 - x_2| > \alpha\}$$, with a strict inequality. By a simple approximation argument one easily sees that this change does not influence the answer.

Obviously, if $$|m_1 - m_2| \leqslant \alpha$$, then $$W_c(\mu_1, \mu_2) = 0$$: simply take $$\gamma$$ to be the distribution of $$(X + m_1, X + m_2)$$, where $$X$$ has Gaussian distribution with mean $$0$$ and variance $$\sigma^2$$.

Suppose that $$|m_1 - m_2| > \alpha$$, and assume with no loss of generality that $$m_2 > m_1$$. Let $$f_m(x)$$ be the probability density function of a Gaussian distribution with mean $$m$$ and variance $$\sigma^2$$. Define $$\gamma$$ in the following way (see the edit at the end of this answer for some intuition): \begin{aligned} \gamma(dx_1, dx_2) & = \Bigl(\min\{f_{m_1}(x_1), f_{m_2 - \alpha}(x_1)\} \delta_{x_1 + \alpha}(dx_2) \\ & \qquad + \max\{0, f_{m_1}(x_1) - f_{m_2 - \alpha}(x_1)\} \delta_{m_1 + m_2 - x_1}(dx_2)\Bigr) dx_1 . \end{aligned} This is a distribution concentrated on the line $$x_2 = x_1 + \alpha$$ and the half-line $$\begin{cases} x_2 = m_1 + m_2 - x_1, \\ x_2 \geqslant x_1 + \alpha. \end{cases}$$ It is immediate to see that $$\gamma(dx_1, \mathbb{R}) = f_{m_1}(x_1) dx_1$$, and it takes a while to find that \begin{aligned} & \gamma(dx_1, dx_2) \\ & = \Bigl(\min\{f_{m_1}(x_1), f_{m_2 - \alpha}(x_1)\} \delta_{x_2 - \alpha}(dx_1) \\ & \qquad + \max\{0, f_{m_1}(x_1) - f_{m_2 - \alpha}(x_1)\} \delta_{m_1 + m_2 - x_2}(dx_1)\Bigr) dx_2 \\ & = \Bigl(\min\{f_{m_1}(x_2 - \alpha), f_{m_2 - \alpha}(x_2 - \alpha)\} \delta_{x_2 - \alpha}(dx_1) \\ & \qquad + \max\{0, f_{m_1}(m_1 + m_2 - x_2) - f_{m_2 - \alpha}(m_1 + m_2 - x_2)\} \delta_{m_1 + m_2 - x_2}(dx_1)\Bigr) dx_2 \\ & = \Bigl(\min\{f_{m_1 + \alpha}(x_2), f_{m_2}(x_2)\} \delta_{x_2 - \alpha}(dx_1) \\ & \qquad + \max\{0, f_{m_2}(x_2) - f_{m_1 + \alpha}(x_2)\} \delta_{m_1 + m_2 - \alpha - x_2}(dx_1)\Bigr) dx_2\end{aligned} so that $$\gamma(\mathbb{R}, dx_2) = f_{m_2}(x_2) dx_2$$. It follows that \begin{aligned} W_c(\mu_1, \mu_2) & \leqslant \int_{-\infty}^\infty \max\{0, f_{m_1}(x_1) - f_{m_2 - \alpha}(x_1)\} dx_1 \\ & = \int_{-\infty}^{(m_1 + m_2 - \alpha) / 2} (f_{m_1}(x_1) - f_{m_2 - \alpha}(x_1)) dx_1 \\ & = \mathbb{P}(X < \tfrac{m_2 - m_1 - \alpha}{2}) - \mathbb{P}(X < -\tfrac{m_2 - m_1 - \alpha}{2}) \\ & = \mathbb{P}(|X| < \tfrac{m_2 - m_1 - \alpha}{2}) , \end{aligned} where $$X$$ has Gaussian distribution with mean $$0$$ and variance $$\sigma^2$$.

We claim that the above bound is optimal. Indeed, suppose that $$(X_1, X_2)$$ is a vector with marginals $$f_{m_1}(x_1) dx_1$$ and $$f_{m_2}(x_2) dx_2$$. Then \begin{aligned} \mathbb{P}(|X_2 - X_1| > \alpha) & \geqslant \mathbb{P}(X_2 - X_1 > \alpha) \\ & \geqslant \mathbb{P}(X_2 > \tfrac{m_1 + m_2 + \alpha}{2}, \, X_1 < \tfrac{m_1 + m_2 - \alpha}{2}) \\ & \geqslant \mathbb{P}(X_2 > \tfrac{m_1 + m_2 + \alpha}{2}) - \mathbb{P}(X_1 \geqslant \tfrac{m_1 + m_2 - \alpha}{2}) \\ & = \mathbb{P}(X > \tfrac{m_1 - m_2 + \alpha}{2}) - \mathbb{P}(X \geqslant \tfrac{-m_1 + m_2 - \alpha}{2}) \\ & = \mathbb{P}(|X| < \tfrac{m_2 - m_1 - \alpha}{2}) , \end{aligned} with $$X$$ as above. Thus, $$W_c(\mu_1, \mu_2) \geqslant \mathbb{P}(|X| < \tfrac{m_2 - m_1 - \alpha}{2}) .$$

This gives the desired expression for $$W_c(\mu_1, \mu_2)$$: $$W_c(\mu_1, \mu_2) = \mathbb{P}(|X| < \tfrac{m_2 - m_1 - \alpha}{2}) = 2 \Phi(\tfrac{m_2 - m_1 - \alpha}{2 \sigma}) - 1 ,$$ where $$\Phi$$ is the cumulative distribution function of the standard Gaussian distribution.

Edit: The idea behind the coupling $$\gamma$$ constructed above is in fact quite simple. We like to have $$X_1 = X_2 - \alpha$$ with as high probability as possible. The "probability" that $$X_1 = x_1$$ must equal to $$f_{m_1}(x_1)$$, while the "probability" that $$X_2 - \alpha = x_1$$ must be equal to $$f_{m_2 - \alpha}(x_1)$$. Thus, we set $$X_1 = X_2 - \alpha = x_1$$ with "probability" $$\min\{f_{m_1}(x_1), f_{m_2 - \alpha}(x_1)\}$$.

This is now extended to a full coupling in a completely arbitrary way. It is quite easy to see that this is indeed possible (any sub-probability distribution $$\gamma_0$$ with marginals less than $$f_{m_1}(x_1) dx_1$$ and $$f_{m_2}(x_2) dx_2$$ can be extended to a full coupling of these measures), and in fact one easily finds an explicit expression by exploiting symmetry of Gaussians.

• Thanks. Please could you explain the intuition behind the construction of your coupling $\gamma$ ? I really looks mysterious. – dohmatob Jul 31 '19 at 8:56
• @dohmatob: I have just edited in some intuitions. The idea is quite simple, I think I rather failed to describe it clearly then. – Mateusz Kwaśnicki Jul 31 '19 at 9:09
• Thanks :). In the meanwhile, I'd come up with a somewhat synthetic solution to a generalization of the second part of my original question. See post below. Would like to know what you think. N.B.: It's possible my answer has a loophole due to an unjustified step (b). – dohmatob Jul 31 '19 at 10:42

I agree with the comment by Mateusz that a simple expression is unlikely to exist for general sets $$\Omega=:C$$. For such sets, the best result is apparently as follows: $$W_c(\mu,\nu)=\sup\{\mu(A)-\nu(A^C)\colon A\subseteq\mathbb R, A\text{ closed}\},$$ where $$C$$ is a nonempty open subset of $$\mathbb R^2$$, $$c=1_C$$, $$A^C:=\{y\in\mathbb R\colon\;\exists x\in A\ (x,y)\notin C\},$$ and $$\mu$$ and $$\nu$$ are any (not necessarily Gaussian) probability measures over $$\mathbb R$$; see e.g. Theorem 1.27, p. 44.

• Thanks. I had already managed to obtain this general expression, from the "Strassen formula" you referred me to on another question :) – dohmatob Aug 2 '19 at 8:24

# Alternative olution to generalized version of second part

For simplicity of notation, let $$\varepsilon := \alpha/2$$. I'll only consider the second part of the problem, but generalized so that the base space is a Banach space $$\mathcal X$$ with norm $$\|.\|$$ (in the statement of my problem, I considered $$\mathcal X = \mathbb R$$, equiped with the absolute-value norm. Also, $$\Omega := \{(x',x) \in \mathcal X^2 \mid \|x'-x\| > 2\varepsilon\}$$.

Lemma. We have the identity $$\begin{eqnarray} c_\Omega(x',x) = \inf_{\|z\| \le \varepsilon}\delta(x'-z,x+z),\;\forall x',x \in \mathcal X. \end{eqnarray}$$

Proof. We show that $$\|x'-x\| > 2\varepsilon$$ iff $$x'-z \ne x + z$$ for all $$z \in \mathbb B_\varepsilon(0)$$. Indeed, if $$x'-z = x + z$$ for some $$z \in \mathcal B_\varepsilon(0)$$, then $$\|x'-x\| = 2\|z\| \le 2\epsilon$$.

Conversely, if $$\|x'-x\| \le 2\varepsilon$$, we may take $$z=(x'-x)/2 \in \mathbb B_\varepsilon(0)$$ and note that $$x'-x = x + z = (x' + x) / 2$$. $$\quad\quad\quad\quad\Box$$

Given a distribution $$P$$ on $$\mathcal X$$ and a point $$z \in \mathcal X$$, let $$P+z$$ be the translation of $$P$$ by $$z$$. Alternatively, if $$X$$ is a random variable with distribution $$P$$, then $$P+z$$ corresponds to the distribution of the r.v $$X + z$$.

Theorem. We have the inequality $$\begin{eqnarray} c_\Omega(P_1,P_2) \le \inf_{\|z\| \le \varepsilon}\text{TV}(P_1 - z, P_2 + z) \end{eqnarray}$$

Proof. One computes $$\begin{split} c_\Omega(P_1,P_2) &:= \inf_{\gamma \in \Pi(P_1,P_2)}\int_{\mathcal X^2} c_\Omega(x',x)d\gamma(x',x) \overset{(a)}{=} \inf_{\gamma \in \Pi(P_1,P_2)}\int_{\mathcal X^2}\inf_{\|z\| \le \varepsilon}\delta(x'-z,x+z)d\gamma(x',x)\\ &\overset{(b)}{\le}\inf_{\gamma \in \Pi(P_1,P_2)}\inf_{\|z\| \le \varepsilon}\int_{\mathcal X^2}\delta(x'-z,x+z)d\gamma(x',x)\\ &\overset{(c)}{=}\inf_{\|z\| \le \varepsilon}\inf_{\gamma \in \Pi(P_1,P_2)}\int_{\mathcal X^2}\delta(x'-z,x+z)d\gamma(x',x)\\ &\overset{(d)}{=}\inf_{\|z\| \le \varepsilon}\inf_{\gamma \in \Pi(P_1-z,P_2+z)}\int_{\mathcal X^2}\delta(x',x)d\gamma(x',x) \overset{(e)}{=} \inf_{\|z\| \le \varepsilon}\text{TV}(P_1-z,P_2+z), \end{split}$$ where

• (a) is by the above Lemma
• (b) inequality due to interchange of integral and inf
• (c) is nothing special
• (d) is a simple change-of-variable formula
• (e) is simply the definition of total-variation distance. $$\quad\quad\quad\quad\Box$$

## Special case of multivariate Gaussians

The following is a direct corollary.

Corollary. In $$\mathbb R^p$$ with $$\Omega := \{(x',x) \in \mathbb R^p \mid \|x'-x\| > 2\varepsilon\}$$, we have $$\begin{eqnarray} c_\Omega(\mathcal N(\mu_1,\Sigma),\mathcal N(\mu_2,\Sigma)) = 2\Phi(\Delta(\varepsilon))-1, \end{eqnarray}$$ where $$\Delta(\varepsilon):=\inf_{\|z\| \le \varepsilon}\|z-\mu\|_{\Sigma^{-1}}$$, $$\mu := \mu_1-\mu_2$$, and $$\|z-\mu\|_{\Sigma^{-1}} := \sqrt{(z-\mu)^T\Sigma^{-1}(z-\mu)}$$.

Proof. Follows from the above Theorem and the well-known fact that fact that $$\text{TV}(\mathcal N(\mu_1,\Sigma),\mathcal N(\mu_2,\Sigma)) = \Phi(\|\mu\|_{\Sigma^{-1}})-\Phi(-\|\mu\|_{\Sigma^{-1}}) = 2\Phi(\|\mu\|_{\Sigma^{-1}})-1. \quad\quad\Box$$

## Special case of 1d Gaussians

In one dimension ($$p=1$$), $$\Sigma = \sigma^2$$, the common variance, and we have the closed-form formula $$\Delta(\varepsilon) = \sigma^{-1}(\mu-\varepsilon)_+$$, which is the result obtained by user Mateusz.

• Your Theorem is false without further assumptions. Not only "$\sup$" should be "$\inf$" (both in Lemma and Theorem; I assume that $\delta(x,y) = 0$ if $x = y$ and $\delta(x,y) = 1$ otherwise). Consider $\alpha = 1$, $P_1 = \delta_0$ and $P_2 = \tfrac{1}{2} \delta_{-1} + \tfrac{1}{2} \delta_1$. The unique coupling $\gamma = P_1 \times P_2$ gives $W_c(P_1, P_2) = 0$, but $\operatorname{TV}(P_1 - z, P_2 + z)$ is always at least $1$. – Mateusz Kwaśnicki Jul 31 '19 at 19:56
• That said, Theorem is likely true when $P_1$ and $P_2$ are translations of a given isotropic unimodal distribution. The upper bound for $W_c(P_1, P_2)$ can be constructed in the same way as in my answer (and in fact the inequality $W_c(P_1, P_2) \geqslant \inf \operatorname{TV}(P_1 - z, P_2 + z)$ is always true). For the lower bound, certain modifications are necessary, but if I am not mistaken, the idea is the same. – Mateusz Kwaśnicki Jul 31 '19 at 20:02
• Thanks. The $sup$'s were a typo / thinko in the Lemma which propagated all through. Fixed. With the typo out of the way, it seems the bug in the theorem proof is my step (d). Are you saying this step is justified only when the $P_k$ are translations os unimodal isotropic distributions ? Also, why would we need isotropy here ? – dohmatob Jul 31 '19 at 20:21
• I only meant that the theorem should be fine if $P_1$ and $P_2$ are translations of an isotropic unimodal distribution. This is not a necessary condition, the result may well be true under more general conditions.. – Mateusz Kwaśnicki Jul 31 '19 at 20:47