  We investigate the problem of pairwise multi-marginal optimal transport, that
is, given a collection of probability distributions $\{P_\alpha\}$ on a Polish
space $\mathcal{X}$, to find a coupling $\{X_\alpha\}$, $X_\alpha\sim
P_\alpha$, such that $\mathbf{E}[c(X_\alpha,X_\beta)]\le r\inf_{X\sim
P_\alpha,Y\sim P_\beta}\mathbf{E}[c(X,Y)]$ for all $\alpha,\beta$, where $c$ is
a cost function and $r\ge1$. In other words, every pair $(X_\alpha,X_\beta)$
has an expected cost at most a factor of $r$ from its lowest possible value.
This can be regarded as a locality sensitive hash function for probability
distributions, and has applications such as robust and distributed computation
of transport plans. It can also be considered as a bi-Lipschitz embedding of
the collection of probability distributions into the space of random variables
taking values on $\mathcal{X}$. For $c(x,y)=\Vert x-y\Vert_2^q$ on
$\mathbb{R}^n$, where $q>0$, we show that a finite $r$ is attainable if and
only if either $n=1$ or $0<q<1$. As $n\to\infty$, the growth rate of the
smallest possible $r$ is exactly $\Theta(n^{q/2})$ if $0<q<1$. Hence, the
metric space of probability distributions on $\mathbb{R}^n$ with finite $q$-th
absolute moments, $0<q<1$, with the earth mover's distance (or 1-Wasserstein
distance) with respect to the snowflake metric $c(x,y)=\Vert x-y\Vert_2^q$, is
bi-Lipschitz embeddable into $L_1$ with distortion $O(n^{q/2})$. If we consider
$c(x,y)=\Vert x-y\Vert_2$ (i.e., $q=1$) on the grid $[0..s]^n$ instead of
$\mathbb{R}^n$, then $r=O(\sqrt{n}\log s)$ is attainable, which implies the
embeddability of the space of probability distributions on $[0..s]^n$ into
$L_1$ with distortion $O(\sqrt{n}\log s)$, and improves upon the $O(n\log s)$
result by Indyk and Thaper. The case of the discrete metric cost
$c(x,y)=\mathbf{1}\{x\neq y\}$ and more general metric and ultrametric costs
are also investigated.
