We investigate the problem of pairwise multi-marginal optimal transport, that is, given a collection of probability distributions $\{P_α\}$ on a Polish space $\mathcal{X}$, to find a coupling $\{X_α\}$, $X_α\sim P_α$, such that $\mathbf{E}[c(X_α,X_β)]\le r\inf_{X\sim P_α,Y\sim P_β}\mathbf{E}[c(X,Y)]$ for all $α,β$, where $c$ is a cost function and $r\ge1$. In other words, every pair $(X_α,X_β)$ 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 $Θ(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.