Question about Wasserstein metric Let $\mu$ and $\nu$ be two probability measures on $\mathbb R^n$ with finite first moment. Denote by $d:=W_1(\mu,\nu)$, where $W_1(\cdot,\cdot)$ stands for the Wasserstein distance of order $1$. 
My question is the following: Let $X$ be a random variable defined on some probability space (rich enough) with law $\mu$, could we find a measurable function $f:\mathbb R^n\times \mathbb R^n\to\mathbb R^n$ and a random variable $G$ independent of $X$ s.t.
$$Y:=f(X,G)~\sim~\nu~~~~~~ \mbox{ and }~~~~~~ \mathbb E[|X-Y|]~\le ~2d~?$$
Thought 1: Let $d_0:=\rho(\mu,\nu)$, where $\rho(\cdot,\cdot)$ denotes the Prokhorov distance. Then we have a measurable function $f_0:\mathbb R^n\times \mathbb R^n\to\mathbb R^n$ and a random variable $G_0$ independent of $X$ s.t.
$$Y_0:=f_0(X,G)~\sim~\nu~~~~~~ \mbox{ and }~~~~~~ \mathbb P[\{|X-Y_0|\ge2d_0\}]~\le ~2d_0.$$
The above construction is from the paper On a representation of random variables by Skorokhod, but I can't find this paper. 
Thought 2: Let $\pi(dx,dy)$ be the optimal transport plan, i.e. $\pi(A\times\mathbb R^n)=\mu(A)$ and $\pi(\mathbb R^n\times A)=\nu(A)$ for all measurable $A\subset\mathbb R^n$. Disintegration w.r.t. the first coordinate $x$, one has $\pi(dx,dy)=\mu(dx)\otimes \lambda_x(dy)$, where $(\lambda_x)_{x\in\mathbb R^n}$ denotes the r.c.p.d. (regular conditional probability distribution). But I've no idea how to recover the function $f$ using $\lambda_x$.
Any answer, help or comment is highly appreciated. Thanks a lot!
 A: Yes, you can even ensure $\mathbb{E}(|X-Y|)=d$ and all you need for $G$ is that it has an atomless law, $\rho$ say.
Take the disintegration $(\lambda_x)$ of the optimal coupling $\Pi$ with respect to $\mu$ (note your formula has a $dx$ where one should read a $dy$). What you want is that $f(x,\cdot)$ sends the law $\rho$ of $G$ to $\lambda_x$. As soon as $\rho$ has no atom, this can be done. Then you only have to check that you can collect the $f(x,\cdot)$ into a measurable map $f$ (you need to be a bit careful in the construction for this to hold; in dimension $1$ using distribution functions should help, and in higher dimension you can use the fact that all standard space are isomorphic; alternatively you can possibly take advantage of the extra bit of margin you took in the question).
A: $\newcommand{\R}{\mathbb R}
\newcommand{\B}{\mathcal B}
\newcommand{\la}{\lambda}
\newcommand{\Si}{\Sigma}
\renewcommand{\c}{\circ}
\newcommand{\tr}{\operatorname{tr}}$
The desired function $f$ and random variable (r.v.) $G$ can be built recursively, by induction, using the increasing rearrangement/inverse transformation method: If $F$ is any cumulative distribution function (cdf), 
\begin{equation}
 F^{-1}(u):=\inf\{x\in\mathbb R\colon F(x)\ge u\}
\end{equation}
for $u\in(0,1)$, and $U\sim\mathcal U(0,1)$ (a r.v. uniformly distributed on the interval $(0,1)$), then the cdf of the r.v. $F^{-1}(U)$ is $F$. 
Indeed, for $j=0,\dots,n$, let $\pi_j$ be the push-forward image of the probability measure $\pi$ under the projection of $\R^n\times\R^n$ onto $\R^n\times\R^j$, so that $\pi_j(A\times B_j)=\pi(A\times B_j\times\R^{n-j})$ for $A$ in the Borel sigma-algebra $\B(\R^n)$ and $B_j$ in $\B(\R^j)$; naturally, $\R^0=\{0\}$, and we identify $\R^j\times\R^{n-j}$ with $\R^n$. 
Similarly, let $\nu_j$ be the push-forward image of the probability measure $\nu$ under the projection of $\R^n$ onto $\R^j$, so that $\nu_j(B_j)=\nu( B_j\times\R^{n-j})$ for $B_j$ in $\B(\R^j)$. 
Then $\pi_0 =\mu$, $\pi_n=\pi$, $\nu_0$ is the only probability measure on $\B(\R^0)=\B(\{0\})$, and $\nu_n=\nu$. 
Write $Y=(Y_1,\dots,Y_n)$ and let $Y_{1;j}:=(Y_1,\dots,Y_j)$, with $Y_{1;0}:=0$; the r.v.'s $Y_1,\dots,Y_n$ are to be constructed. 
Accordingly, for $y=(y_1,\dots,y_n)\in\R^n$ let $y_{1;j}:=(y_1,\dots,y_j)$, with $y_{1;0}:=0$. For $j=0,\dots,n$, let $G_j:=(U_1,\dots,U_j)$, where $U_1,\dots,U_n$ are independent $\mathcal U(0,1)$ r.v.'s. In particular, $G_0=0$. 
For $j=1,\dots,n$, we are going to construct, by induction, a function $f_j\colon\R^n\times(0,1)^j\to\R^j$ such that for $Y_{1;j}:=f_j(X,G_j)$ the distribution of $(X,Y_{1;j})$ is $\pi_j$ and hence the distribution of $Y_{1;j}$ is $\nu_j$. 
To complete the basis of induction, let $f_0(x,0):=0$ for all $x\in\R^n$, so that $Y_0=f_0(X,G_0)$. 
Take now any $j=1,\dots,n$. Let 
$$\R^n\times\R^{j-1}\times\B(\R)\ni(x,y_{1;j-1},C)\longmapsto
\la_{x,y_{1;j-1}}(C)$$ 
be a regular version of the conditional distribution of $Y_j$ given $(X,Y_{1;j-1})$ assuming that the joint distribution of $(X,Y_{1;j})$ is $\pi_j$, so that 
$$\pi_j(dx\times dy_{1;j-1}\times dy_j)=\pi_{j-1}(dx\times dy_{1;j-1})\la_{x,y_{1;j-1}}(dy_j).$$ 
For each $(x,y_{1;j-1})\in\R^n\times\R^{j-1}$, let $F_{x,y_{1;j-1}}$ be the cdf of the probability measure $\la_{x,y_{1;j-1}}$ on $\B(\R)$, and define the function $f_j\colon\R^n\times(0,1)^j\to\R^j$ by the formula 
\begin{equation}
 f_j(x,u_{1;j}):=\big(y_{1;j-1},F^{-1}_{x,y_{1;j-1}}(u_j)\big)\quad\text{with}\quad 
 y_{1;j-1}=f_{j-1}(x,u_{1;j-1})  
\end{equation}
for $(x,u_{1;j})\in\R^n\times(0,1)^j$, where the notation $u_{1;j}$ is of course quite similar to $y_{1;j}$. 
Let now $Y_{1;j}:=f_j(X,U_{1;j})=f_j(X,G_j)$, which is in agreement with the definition $Y_{1;j-1}:=f_{j-1}(X,U_{1;j-1})=f_{j-1}(X,G_{j-1})$ at the previous step of the induction process. 
Then the distribution of $(X,Y_{1;j})$ is $\pi_j$ and hence the distribution of $Y_{1;j}$ is $\nu_j$. 
In particular, the distribution of $(X,Y)=(X,Y_{1;n})$ is $\pi_n=\pi$ and hence the distribution of $Y=Y_{1;n}$ is $\nu$. Moreover, $Y=Y_{1;n}=f_n(X,G_n)$, as desired. 
A: I've a solution but it's not perfectly satisfying. Assume that 
$$V~~~:=~~~\int |x|^pd\mu(x)~+~\int |x|^pd\nu(x)~~~<~~~+\infty$$
for some fixed $p>1$. It follows from Thought 1 that, there exists $f_0$ and $G$ s.t. 
$$Y_0~:=~f_0(X,G)~\sim~\nu~~~~\mbox{ and }~~~~\mathbb P[|X-Y_0|\ge 2d_0]~\le~2d_0.$$
It follows that
\begin{eqnarray}
\mathbb E[|X-Y_0|]\quad&=&\quad\mathbb E[|X-Y_0|{\bf 1}_{\{|X-Y_0|\ge 2d_0\}}]~+~\mathbb E[|X-Y_0|{\bf 1}_{\{|X-Y_0|< 2d_0\}}] \\
&\le&\quad 2d_0 ~+~ \big(E[|X-Y_0|^p]\big)^{1/p}\cdot \mathbb P[|X-Y_0|\ge 2d_0]^{1/q} \\
&\le&\quad 2d_0 ~+~ \left[\big(E[|X|^p]\big)^{1/p}~+~\big(E[|Y_0|^p]\big)^{1/p}\right]\cdot (2d_0)^{1/q} \\
&\le&\quad 2d_0 ~+~ 2V^{1/p}\cdot (2d_0)^{1/q},
\end{eqnarray}
where $q>1$ denotes the conjugate number of $p$, i.e. $1/p+1/q=1$. Let $\alpha(x):=2x ~+~ 2V^{1/p}\cdot (2x)^{1/q}$, then one has 
$$\mathbb E[|X-Y_0|]\quad\le\quad \alpha(d_0)\quad\le \quad \alpha(d),$$
as $d_0\le d$.
My question is that could we remove the assumption $V<+\infty$.
