Proof of Couplings on Euclidean Space: Product Coupling, Swap, Finiteness of the Cost, Push-Forward Couplings, Modifying One Marginal, Quantisation, Gluing over a Finitely Supported Measure, and the Lipschitz Bound

Each result cited is universally quantified over the data in its own statement. "The integral theorem" refers to Linearity and Monotonicity of the Lebesgue Integral, and the Integrability Criterion (claim 1: additivity, homogeneity with a factor in $[0,\infty)$, and monotonicity of the nonnegative integral; claim 2: linearity, $|\int f|\le\int|f|$ and monotonicity for integrable functions), and "change of variables" to the formula of Probability Measures on Euclidean Space and Random Vectors: Standing Notation (iv). Let $\varphi:\mathbb{R}^{d+d}\to\mathbb{R}$ be the Borel map $z\mapsto\lVert\mathrm{pr}_{1}(z)-\mathrm{pr}_{2}(z)\rVert^{2}$ and $\varphi_{0}$ the Borel map $z\mapsto\lVert\mathrm{pr}_{1}(z)-\mathrm{pr}_{2}(z)\rVert$, both of Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (iv), so that $I(\pi)=\int\varphi\,d\pi$ for every $\pi\in\mathcal{P}(\mathbb{R}^{d+d})$ by Couplings of Two Probability Measures on Euclidean Space and Their Quadratic Cost (ii), each such $\pi$ being a coupling of its own push-forwards by Couplings of Two Probability Measures on Euclidean Space and Their Quadratic Cost (i). Two facts are used repeatedly. (i) For Borel $G:\mathbb{R}^{m}\to\mathbb{R}^{l}$ and $H:\mathbb{R}^{l}\to\mathbb{R}^{k}$ and $\lambda\in\mathcal{P}(\mathbb{R}^{m})$, one has $(H\circ G)_{\#}\lambda=H_{\#}(G_{\#}\lambda)$, since both sides assign to a Borel $B$ the value $\lambda(G^{-1}(H^{-1}(B)))$. (ii) For $u,v:\mathbb{R}^{m}\to\mathbb{R}^{d}$ Borel, $\mathrm{pr}_{1}\circ(u,v)=u$ and $\mathrm{pr}_{2}\circ(u,v)=v$ by Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (i), and $\varphi((u,v)(x))=\lVert u(x)-v(x)\rVert^{2}$, $\varphi_{0}((u,v)(x))=\lVert u(x)-v(x)\rVert$. Finally, if $\lambda$ is a probability measure on a measurable space, a nonnegative measurable real function with finite integral is integrable with the same integral, by Integrable Functions and Their Lebesgue Integral (its negative part is $0$); this identifies the two readings of such integrals below.

Claim 1. By Probability Measures on Euclidean Space and Random Vectors: Standing Notation (v) and Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (iii) (with $q=p=d$), $\mu\boxtimes\nu\in\mathcal{P}(\mathbb{R}^{d+d})$ and its push-forwards by $\mathrm{pr}_{1}$ and $\mathrm{pr}_{2}$ are $\mu$ and $\nu$; so $\mu\boxtimes\nu\in\Pi(\mu,\nu)$ by Couplings of Two Probability Measures on Euclidean Space and Their Quadratic Cost (i).

Claim 2. The swap $\sigma$ is Borel by Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (ii), so $\sigma_{\#}\pi\in\mathcal{P}(\mathbb{R}^{d+d})$. By Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (i), $\mathrm{pr}_{1}\circ\sigma=\mathrm{pr}_{2}$ and $\mathrm{pr}_{2}\circ\sigma=\mathrm{pr}_{1}$, and $\sigma(\sigma(z))=\iota(\mathrm{pr}_{2}(\sigma(z)),\mathrm{pr}_{1}(\sigma(z)))=\iota(\mathrm{pr}_{1}(z),\mathrm{pr}_{2}(z))=z$. Hence by (i), $(\mathrm{pr}_{1})_{\#}(\sigma_{\#}\pi)=(\mathrm{pr}_{1}\circ\sigma)_{\#}\pi=(\mathrm{pr}_{2})_{\#}\pi=\nu$ and $(\mathrm{pr}_{2})_{\#}(\sigma_{\#}\pi)=\mu$, so $\sigma_{\#}\pi\in\Pi(\nu,\mu)$; and $\sigma_{\#}(\sigma_{\#}\pi)=(\sigma\circ\sigma)_{\#}\pi=\pi$. By change of variables, $I(\sigma_{\#}\pi)=\int\varphi\circ\sigma\,d\pi$, and $\varphi(\sigma(z))=\lVert\mathrm{pr}_{2}(z)-\mathrm{pr}_{1}(z)\rVert^{2}=\varphi(z)$ because $\mathrm{pr}_{2}(z)-\mathrm{pr}_{1}(z)=(-1)(\mathrm{pr}_{1}(z)-\mathrm{pr}_{2}(z))$ by claims 2 and 3 of Euclidean Space $\mathbb{R}^n$ is a Real Vector Space and claims 2 and 5 of Elementary Identities in a Vector Space, and $\lVert(-1)x\rVert=\lVert x\rVert$ by claim 5 of Properties of the Euclidean Norm on $\mathbb{R}^n$: Square, Distance, Vanishing, Coordinate Bounds, Homogeneity, Triangle Inequalities and Expansion of a Square; so $I(\sigma_{\#}\pi)=I(\pi)$. The map $\pi\mapsto\sigma_{\#}\pi$ thus sends $\Pi(\mu,\nu)$ into $\Pi(\nu,\mu)$, and the same map on $\Pi(\nu,\mu)$ is a two-sided inverse of it, so it is a bijection.

Claim 3. By Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (iv), $\varphi(z)\le2\lVert\mathrm{pr}_{1}(z)\rVert^{2}+2\lVert\mathrm{pr}_{2}(z)\rVert^{2}$ for every $z$. The maps $z\mapsto\lVert\mathrm{pr}_{i}(z)\rVert^{2}$ are Borel, being compositions of $\mathrm{pr}_{i}$ with $x\mapsto\lVert x\rVert^{2}$ (Probability Measures on Euclidean Space and Random Vectors: Standing Notation (ii)), and by change of variables $\int\lVert\mathrm{pr}_{1}(z)\rVert^{2}\,\pi(dz)=\int\lVert x\rVert^{2}\,((\mathrm{pr}_{1})_{\#}\pi)(dx)=M_{2}(\mu)$, and likewise the second integral is $M_{2}(\nu)$, by The Second Moment of a Probability Measure on Euclidean Space and the Probability Measures with Finite Second Moment (i). Claim 1 of the integral theorem gives $I(\pi)\le2M_{2}(\mu)+2M_{2}(\nu)$ in $[0,\infty]$, which is finite when $\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})$. Conversely, Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (iv) with $x=\mathrm{pr}_{2}(z)$ and $y=\mathrm{pr}_{1}(z)$ gives $\lVert\mathrm{pr}_{2}(z)\rVert^{2}\le2\lVert\mathrm{pr}_{1}(z)\rVert^{2}+2\varphi(z)$, using $\lVert\mathrm{pr}_{2}(z)-\mathrm{pr}_{1}(z)\rVert=\lVert\mathrm{pr}_{1}(z)-\mathrm{pr}_{2}(z)\rVert$ as in claim 2; integrating as before, $M_{2}(\nu)\le2M_{2}(\mu)+2I(\pi)$, which is finite under the stated hypotheses, so $\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})$.

Claim 4. The pairing $(S,T)$ is Borel by Probability Measures on Euclidean Space and Random Vectors: Standing Notation (v), so $(S,T)_{\#}\mu\in\mathcal{P}(\mathbb{R}^{d+d})$, and by (i) and (ii) its push-forward by $\mathrm{pr}_{1}$ is $(\mathrm{pr}_{1}\circ(S,T))_{\#}\mu=S_{\#}\mu$, and by $\mathrm{pr}_{2}$ is $T_{\#}\mu$. By change of variables and (ii), $I((S,T)_{\#}\mu)=\int\varphi\circ(S,T)\,d\mu=\int\lVert S(x)-T(x)\rVert^{2}\,\mu(dx)$. For $S=T=\mathrm{id}$, $\mathrm{id}_{\#}\mu=\mu$ and the integrand is $\lVert x-x\rVert^{2}=0$ by claim 2 of Elementary Identities in a Vector Space and claim 3 of Properties of the Euclidean Norm on $\mathbb{R}^n$: Square, Distance, Vanishing, Coordinate Bounds, Homogeneity, Triangle Inequalities and Expansion of a Square, whose integral is $0$ by claim 1 of the integral theorem with the factor $0$.

Claim 5. The maps $T\circ\mathrm{pr}_{2}$ and $S\circ\mathrm{pr}_{1}$ are Borel by Probability Measures on Euclidean Space and Random Vectors: Standing Notation (ii), so the two pairings are Borel and the push-forwards lie in $\mathcal{P}(\mathbb{R}^{d+d})$. By (i) and (ii), the push-forward of $\pi'=(\mathrm{pr}_{1},T\circ\mathrm{pr}_{2})_{\#}\pi$ by $\mathrm{pr}_{1}$ is $(\mathrm{pr}_{1})_{\#}\pi=\mu$ and by $\mathrm{pr}_{2}$ is $(T\circ\mathrm{pr}_{2})_{\#}\pi=T_{\#}((\mathrm{pr}_{2})_{\#}\pi)=T_{\#}\nu$; so $\pi'\in\Pi(\mu,T_{\#}\nu)$, and likewise $(S\circ\mathrm{pr}_{1},\mathrm{pr}_{2})_{\#}\pi\in\Pi(S_{\#}\mu,\nu)$. Now assume $I(\pi)<\infty$ and $\int\lVert T(y)-y\rVert^{2}\,\nu(dy)<\infty$. Regard $(\mathbb{R}^{d+d},\mathcal{B}(\mathbb{R}^{d+d}),\pi)$ as a probability space, so that Borel real functions on $\mathbb{R}^{d+d}$ are random variables on it; let $f=\varphi_{0}$ and $g=\varphi_{0}\circ(\mathrm{pr}_{2},T\circ\mathrm{pr}_{2})$, Borel by Probability Measures on Euclidean Space and Random Vectors: Standing Notation (ii), so that $f(z)=\lVert\mathrm{pr}_{1}(z)-\mathrm{pr}_{2}(z)\rVert$ and $g(z)=\lVert\mathrm{pr}_{2}(z)-T(\mathrm{pr}_{2}(z))\rVert$ by (ii). Then $\int f^{2}\,d\pi=I(\pi)<\infty$, and $\int g^{2}\,d\pi=\int\lVert\mathrm{pr}_{2}(z)-T(\mathrm{pr}_{2}(z))\rVert^{2}\,\pi(dz)=\int\lVert y-T(y)\rVert^{2}\,\nu(dy)<\infty$ by change of variables through $\mathrm{pr}_{2}$, the integrand being $\lVert y-T(y)\rVert^{2}=\lVert T(y)-y\rVert^{2}$ (claim 5 of Properties of the Euclidean Norm on $\mathbb{R}^n$: Square, Distance, Vanishing, Coordinate Bounds, Homogeneity, Triangle Inequalities and Expansion of a Square as in claim 2); so $f$ and $g$ are square-integrable random variables with $\lVert f\rVert_{2}=\sqrt{I(\pi)}$ and $\lVert g\rVert_{2}=\sqrt{\int\lVert T(y)-y\rVert^{2}\,\nu(dy)}$. For every $z$, writing $x=\mathrm{pr}_{1}(z)$ and $y=\mathrm{pr}_{2}(z)$, one has $x-T(y)=(x-y)+(y-T(y))$ by Euclidean Space $\mathbb{R}^n$ is a Real Vector Space, so $\lVert x-T(y)\rVert\le f(z)+g(z)$ by claim 6 of Properties of the Euclidean Norm on $\mathbb{R}^n$: Square, Distance, Vanishing, Coordinate Bounds, Homogeneity, Triangle Inequalities and Expansion of a Square, and $\lVert x-T(y)\rVert^{2}\le(f(z)+g(z))^{2}$ by claim 2 of Monotonicity of Squaring on the Nonnegative Elements of an Ordered Field. By change of variables and (ii), $I(\pi')=\int\varphi\circ(\mathrm{pr}_{1},T\circ\mathrm{pr}_{2})\,d\pi=\int\lVert\mathrm{pr}_{1}(z)-T(\mathrm{pr}_{2}(z))\rVert^{2}\,\pi(dz)\le\int(f+g)^{2}\,d\pi=\lVert f+g\rVert_{2}^{2}$ by monotonicity, and $\lVert f+g\rVert_{2}\le\lVert f\rVert_{2}+\lVert g\rVert_{2}$ by claim 2 of Cauchy-Schwarz and Triangle Inequalities for the Mean-Square Norm; hence $I(\pi')$ is finite and $\sqrt{I(\pi')}\le\lVert f+g\rVert_{2}\le\lVert f\rVert_{2}+\lVert g\rVert_{2}$ by claim 2 of Monotonicity of Squaring on the Nonnegative Elements of an Ordered Field, which is the asserted bound. For the second assertion, one checks from Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (i) that $\sigma\circ(\mathrm{pr}_{1},S\circ\mathrm{pr}_{2})\circ\sigma=(S\circ\mathrm{pr}_{1},\mathrm{pr}_{2})$ pointwise, so by (i) $(S\circ\mathrm{pr}_{1},\mathrm{pr}_{2})_{\#}\pi=\sigma_{\#}\bigl((\mathrm{pr}_{1},S\circ\mathrm{pr}_{2})_{\#}(\sigma_{\#}\pi)\bigr)$; by claim 2, $\sigma_{\#}\pi\in\Pi(\nu,\mu)$ with $I(\sigma_{\#}\pi)=I(\pi)$, so the first assertion applied to $\sigma_{\#}\pi$ and $S$ bounds the cost of $(\mathrm{pr}_{1},S\circ\mathrm{pr}_{2})_{\#}(\sigma_{\#}\pi)$ by the square of $\sqrt{I(\pi)}+\sqrt{\int\lVert S(x)-x\rVert^{2}\,\mu(dx)}$, and claim 2 again shows that the swap preserves this cost.

Claim 6. If $T$ has finite image $F=T(\mathbb{R}^{d})$, then $T^{-1}(\mathbb{R}^{d}\setminus F)=\varnothing$, so $T_{\#}\nu(\mathbb{R}^{d}\setminus F)=\nu(\varnothing)=0$ by Measure, Measure Space, and Probability Measure. Now let $\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})$ and $0<\varepsilon$; put $\eta=\varepsilon^{2}\cdot2^{-1}$ and $\varepsilon'=\varepsilon\cdot2^{-1}$, so that $\eta+\eta=\varepsilon^{2}$, $0<\eta$, $0<\varepsilon'$ by claim 8 of Elementary Order Arithmetic in an Ordered Field, and $\varepsilon'^{2}=\varepsilon^{2}\cdot2^{-1}\cdot2^{-1}\le\eta$ because $2^{-1}<1$ (claim 8 of Elementary Order Arithmetic in an Ordered Field with $\varepsilon=1$, using claim 6 there) and claim 5 of Elementary Arithmetic in an Ordered Field. Let $\kappa:\mathbb{N}\to\mathbb{R}$ be the canonical map. For $n\in\mathbb{N}$ let $K_{n}=\{y:\lVert y\rVert\le\kappa(n)\}$, the complement of $\{y:\lVert y\rVert>\kappa(n)\}$, which is Borel by Measurable Functions and the Lebesgue Integral: Standing Notation (iii) since $y\mapsto\lVert y\rVert$ is Borel; and let $h_{n}=\lVert\cdot\rVert^{2}\mathbf{1}_{K_{n}}$, Borel by claims 1 and 3 of Arithmetic of Measurable Real- and [0,∞]-Valued Maps. As $\kappa$ is increasing (claim 6 of Properties of the Canonical Map from the Natural Numbers to an Ordered Field), $K_{n}\subseteq K_{n+1}$ and $(h_{n}(y))_{n}$ is nondecreasing for each $y$; and by claim 1 of The Archimedean Property of the Real Numbers there is $n$ with $\lVert y\rVert<\kappa(n)$, so $h_{m}(y)=\lVert y\rVert^{2}$ for all $m\ge n$ and $\sup_{n}h_{n}(y)=\lVert y\rVert^{2}$. By Monotone Convergence Theorem, $M_{2}(\nu)=\sup_{n}\int h_{n}\,d\nu$, each $\int h_{n}\,d\nu$ being a real number at most $M_{2}(\nu)<\infty$ by monotonicity. By claim 3 of Approximation Property of the Supremum and the Infimum in $\mathbb{R}$ there is $n$ with $M_{2}(\nu)-\eta<\int h_{n}\,d\nu$. Fix this $n$, put $R=\kappa(n)$ and $K=K_{n}$, and let $g=\lVert\cdot\rVert^{2}\mathbf{1}_{\mathbb{R}^{d}\setminus K}$, Borel likewise. Since $\lVert y\rVert^{2}=h_{n}(y)+g(y)$ for every $y$, claim 1 of the integral theorem gives $M_{2}(\nu)=\int h_{n}\,d\nu+\int g\,d\nu$ with both terms real, whence $\int g\,d\nu=M_{2}(\nu)-\int h_{n}\,d\nu<\eta$.

$K$ is the closed ball of $(\mathbb{R}^{d},d_{E})$ with centre $0_{\mathbb{R}^{d}}$ and radius $R$, since $d_{E}(0_{\mathbb{R}^{d}},y)=\lVert y\rVert$ by claim 2 of Properties of the Euclidean Norm on $\mathbb{R}^n$: Square, Distance, Vanishing, Coordinate Bounds, Homogeneity, Triangle Inequalities and Expansion of a Square and symmetry of the metric; it is compact by claim 2 of A Closed Euclidean Ball is Convex and Compact, the radius $R=\kappa(n)$ being nonnegative by claim 3 of Properties of the Canonical Map from the Natural Numbers to an Ordered Field. The family $(B(c,\varepsilon'))_{c\in K}$ of open balls, which are open by Open Ball in a Metric Space is Open, covers $K$, as $c\in B(c,\varepsilon')$; by Compact Subset Criterion via Open Covers in the Ambient Space there is a finite $J\subseteq K$ with $K\subseteq\bigcup_{c\in J}B(c,\varepsilon')$. Since $0_{\mathbb{R}^{d}}\in K$, $J$ is nonempty, so by Finite Set there are $n'\in\mathbb{N}$ and a bijection $[n']\to J$, $j\mapsto c_{j}$. For $j\in[n']$ let $B_{j}=B(c_{j},\varepsilon')$ and $E_{j}=(K\cap B_{j})\setminus\bigcup_{i\in[n'],\,i<j}B_{i}$, a Borel set (open sets being Borel by Probability Measures on Euclidean Space and Random Vectors: Standing Notation (i), and finite intersections, unions and differences of Borel sets being Borel by σ-Algebra and Measurable Space). For $y\in K$ the set $\{j\in[n']:y\in B_{j}\}$ is nonempty, so it has a least element $j(y)$ by The Natural Numbers Are Well Ordered, and $y\in E_{j(y)}$ while $y\notin E_{j}$ for $j\ne j(y)$; thus the sets $E_{j}$ are pairwise disjoint with union $K$. Define $T:\mathbb{R}^{d}\to\mathbb{R}^{d}$ by $T(y)=c_{j(y)}$ for $y\in K$ and $T(y)=0_{\mathbb{R}^{d}}$ for $y\notin K$. For $B\in\mathcal{B}(\mathbb{R}^{d})$,

$$T^{-1}(B)=\bigcup_{j\in[n'],\,c_{j}\in B}E_{j}\ \cup\ N_{B},\qquad N_{B}=\mathbb{R}^{d}\setminus K\ \text{if }0_{\mathbb{R}^{d}}\in B,\quad N_{B}=\varnothing\ \text{otherwise},$$

a finite union of Borel sets, hence Borel; so $T$ is Borel. Its image is contained in $\{c_{1},\dots,c_{n'}\}\cup\{0_{\mathbb{R}^{d}}\}$, the image of the finite set $[n'+1]$ under the map sending $j\le n'$ to $c_{j}$ and $n'+1$ to $0_{\mathbb{R}^{d}}$, which is finite by claims 1 and 4 of Basic Properties of Finite Sets; so $T(\mathbb{R}^{d})$ is finite by claim 3 there. For $y\in K$, $\lVert T(y)-y\rVert=d_{E}(c_{j(y)},y)<\varepsilon'$ since $y\in B_{j(y)}$, so $\lVert T(y)-y\rVert^{2}\le\varepsilon'^{2}\le\eta$ by claim 2 of Monotonicity of Squaring on the Nonnegative Elements of an Ordered Field; for $y\notin K$, $T(y)-y=(-1)y$ by claims 2 and 3 of Euclidean Space $\mathbb{R}^n$ is a Real Vector Space and Vector Space over a Field, so $\lVert T(y)-y\rVert^{2}=\lVert y\rVert^{2}=g(y)$ by claim 5 of Properties of the Euclidean Norm on $\mathbb{R}^n$: Square, Distance, Vanishing, Coordinate Bounds, Homogeneity, Triangle Inequalities and Expansion of a Square. Hence $\lVert T(y)-y\rVert^{2}\le\eta\mathbf{1}_{K}(y)+g(y)$ for every $y$, and by claim 1 of the integral theorem and The Integral of an Indicator Function is the Measure of the Set,

$$\int\lVert T(y)-y\rVert^{2}\,\nu(dy)\le\eta\,\nu(K)+\int g\,d\nu\le\eta+\eta=\varepsilon^{2},$$

using $\nu(K)\le1$ (claim 2 of Basic Properties of a Measure) and claim 5 of Elementary Arithmetic in an Ordered Field.

Claim 7. Write $F_{+}=\{y\in F:0<\rho(\{y\})\}$, a subset of $F$, finite by claim 3 of Basic Properties of Finite Sets (singletons being Borel by Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (v), so that $\rho(\{y\})$ is defined), and let $w:\mathbb{R}^{d}\to\mathbb{R}$ be $w(y)=\rho(\{y\})^{-1}$ for $y\in F_{+}$ and $w(y)=0$ otherwise, Borel by Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (v) applied to the finite set $F_{+}$. Then $w(y)\rho(\{y\})=\mathbf{1}_{F_{+}}(y)$ for every $y$. Moreover $\rho(\mathbb{R}^{d}\setminus F_{+})=0$: the set $F\setminus F_{+}$ is finite by claim 3 of Basic Properties of Finite Sets, so it is empty or the image of a bijection $[k]\to F\setminus F_{+}$, and then $\rho(F\setminus F_{+})$ is the sum of the $k$ values $\rho(\{y\})=0$ by claim 1 of Basic Properties of a Measure, hence $0$; and $\rho(\mathbb{R}^{d}\setminus F_{+})=\rho(\mathbb{R}^{d}\setminus F)+\rho(F\setminus F_{+})=0$ by the same claim.

Let $q=d+d$ and consider $\mathbb{R}^{q+q}$ with the projections $P_{1}=\mathrm{pr}^{q,q}_{1}$, $P_{2}=\mathrm{pr}^{q,q}_{2}$ and the concatenation $\iota^{q,q}$ of Probability Measures on Euclidean Space and Random Vectors: Standing Notation (v); for $\zeta\in\mathbb{R}^{q+q}$ write $x(\zeta)=\mathrm{pr}_{1}(P_{1}\zeta)$, $y(\zeta)=\mathrm{pr}_{2}(P_{1}\zeta)$, $y'(\zeta)=\mathrm{pr}_{1}(P_{2}\zeta)$ and $v(\zeta)=\mathrm{pr}_{2}(P_{2}\zeta)$, four Borel maps $\mathbb{R}^{q+q}\to\mathbb{R}^{d}$ (Probability Measures on Euclidean Space and Random Vectors: Standing Notation (ii)). Let $\theta_{0}=\pi_{12}\boxtimes\pi_{23}\in\mathcal{P}(\mathbb{R}^{q+q})$ (Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (iii)). The set $D=\{\zeta:y(\zeta)=y'(\zeta)\}$ is Borel: it is the complement of $\{\zeta:\lVert y(\zeta)-y'(\zeta)\rVert^{2}>0\}$, the function $\zeta\mapsto\lVert y(\zeta)-y'(\zeta)\rVert^{2}$ being Borel by Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (iv) and vanishing exactly on $D$ by claim 3 of Properties of the Euclidean Norm on $\mathbb{R}^n$: Square, Distance, Vanishing, Coordinate Bounds, Homogeneity, Triangle Inequalities and Expansion of a Square, claim 3 of Zero Products and Elementary Identities in a Field and claim 2 of Elementary Identities in a Vector Space. Let $h=(w\circ y)\,\mathbf{1}_{D}$, a nonnegative Borel function on $\mathbb{R}^{q+q}$ (claims 1 and 3 of Arithmetic of Measurable Real- and [0,∞]-Valued Maps), and let $\theta$ be the measure with density $h$ with respect to $\theta_{0}$, so that $\int G\,d\theta=\int Gh\,d\theta_{0}$ for every Borel $G:\mathbb{R}^{q+q}\to[0,\infty]$ by claim 3 of that lemma.

Two disintegration identities. Let $g:\mathbb{R}^{q}\to[0,\infty)$ be Borel. We claim the following two identities, referred to below as $(\ast)$:

$$\int g(P_{1}\zeta)\,h(\zeta)\,\theta_{0}(d\zeta)=\int g\,d\pi_{12}\qquad\text{and}\qquad\int g(P_{2}\zeta)\,h(\zeta)\,\theta_{0}(d\zeta)=\int g\,d\pi_{23}.$$

For the first, the integrand $G=(g\circ P_{1})h$ is Borel, so by Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (iii) the function $G\circ\iota^{q,q}$ is measurable for $\mathcal{B}(\mathbb{R}^{q})\otimes\mathcal{B}(\mathbb{R}^{q})$ and $\int G\,d\theta_{0}=\int G\circ\iota^{q,q}\,d(\pi_{12}\otimes\pi_{23})$; here, for $z\in\mathbb{R}^{q}$ and $z'\in\mathbb{R}^{q}$, $G(\iota^{q,q}(z,z'))=g(z)\,w(\mathrm{pr}_{2}z)\,\mathbf{1}[\mathrm{pr}_{2}z=\mathrm{pr}_{1}z']$ by Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (i). By the Tonelli clause of Tonelli's and Fubini's Theorems (both factors being probability measures, hence $\sigma$-finite), this equals $\int_{\mathbb{R}^{q}}\bigl(\int_{\mathbb{R}^{q}}g(z)w(\mathrm{pr}_{2}z)\mathbf{1}[\mathrm{pr}_{2}z=\mathrm{pr}_{1}z']\,\pi_{23}(dz')\bigr)\pi_{12}(dz)$. For fixed $z$, with $y=\mathrm{pr}_{2}z$, the inner integrand is the real constant $g(z)w(y)$ times the indicator of $\mathrm{pr}_{1}^{-1}(\{y\})$, so the inner integral is $g(z)w(y)\pi_{23}(\mathrm{pr}_{1}^{-1}(\{y\}))=g(z)w(y)\rho(\{y\})=g(z)\mathbf{1}_{F_{+}}(y)$ by claim 1 of the integral theorem, The Integral of an Indicator Function is the Measure of the Set and $\pi_{23}\in\Pi(\rho,\nu)$. Thus the first integral in $(\ast)$ equals $\int g(z)\mathbf{1}_{F_{+}}(\mathrm{pr}_{2}z)\,\pi_{12}(dz)$. The integrands $g$ and $g\cdot(\mathbf{1}_{F_{+}}\circ\mathrm{pr}_{2})$ agree off the set $\mathrm{pr}_{2}^{-1}(\mathbb{R}^{d}\setminus F_{+})$, which has $\pi_{12}$-measure $\rho(\mathbb{R}^{d}\setminus F_{+})=0$ as $\pi_{12}\in\Pi(\mu,\rho)$; so they agree $\pi_{12}$-almost everywhere and their integrals coincide by The Lebesgue Integral and Null Sets (iii), proving the first identity. The second is proved in the same way with the order of integration in Tonelli's theorem reversed: for fixed $z'$, with $y'=\mathrm{pr}_{1}z'$, the inner integrand $g(z')w(\mathrm{pr}_{2}z)\mathbf{1}[\mathrm{pr}_{2}z=y']$ equals $g(z')w(y')$ times the indicator of $\mathrm{pr}_{2}^{-1}(\{y'\})$, whose $\pi_{12}$-integral is $g(z')w(y')\rho(\{y'\})=g(z')\mathbf{1}_{F_{+}}(y')$, and the outer integral is $\int g\,d\pi_{23}$ by the same almost-everywhere argument with $\pi_{23}\in\Pi(\rho,\nu)$.

The glued coupling. Taking $g$ the constant $1$ in $(\ast)$ gives $\theta(\mathbb{R}^{q+q})=\int h\,d\theta_{0}=\pi_{12}(\mathbb{R}^{q})=1$, so $\theta\in\mathcal{P}(\mathbb{R}^{q+q})$. Let $\Phi=(x,v):\mathbb{R}^{q+q}\to\mathbb{R}^{d+d}$, the pairing of the Borel maps $x$ and $v$, Borel by Probability Measures on Euclidean Space and Random Vectors: Standing Notation (v), and put $\pi_{13}=\Phi_{\#}\theta\in\mathcal{P}(\mathbb{R}^{d+d})$. For $A\in\mathcal{B}(\mathbb{R}^{d})$, by (i), (ii) and the density formula, $(\mathrm{pr}_{1})_{\#}\pi_{13}(A)=\theta(x^{-1}(A))=\int\mathbf{1}_{x^{-1}(A)}h\,d\theta_{0}$, and $\mathbf{1}_{x^{-1}(A)}=g\circ P_{1}$ with $g=\mathbf{1}_{\mathrm{pr}_{1}^{-1}(A)}$, a Borel map $\mathbb{R}^{q}\to[0,\infty)$; so by $(\ast)$ this equals $\pi_{12}(\mathrm{pr}_{1}^{-1}(A))=\mu(A)$. Likewise $(\mathrm{pr}_{2})_{\#}\pi_{13}(C)=\int(\mathbf{1}_{\mathrm{pr}_{2}^{-1}(C)}\circ P_{2})h\,d\theta_{0}=\pi_{23}(\mathrm{pr}_{2}^{-1}(C))=\nu(C)$. Hence $\pi_{13}\in\Pi(\mu,\nu)$.

The cost bound. Regard $(\mathbb{R}^{q+q},\mathcal{B}(\mathbb{R}^{q+q}),\theta)$ as a probability space, and let $f=\varphi_{0}\circ P_{1}$ and $f'=\varphi_{0}\circ P_{2}$, Borel, so that $f(\zeta)=\lVert x(\zeta)-y(\zeta)\rVert$ and $f'(\zeta)=\lVert y'(\zeta)-v(\zeta)\rVert$ by (ii). By the density formula and $(\ast)$ with $g=\varphi$, $\int f^{2}\,d\theta=\int(\varphi\circ P_{1})h\,d\theta_{0}=\int\varphi\,d\pi_{12}=I(\pi_{12})<\infty$, and similarly $\int f'^{2}\,d\theta=I(\pi_{23})<\infty$; so $f,f'$ are square-integrable random variables on this probability space with $\lVert f\rVert_{2}=\sqrt{I(\pi_{12})}$ and $\lVert f'\rVert_{2}=\sqrt{I(\pi_{23})}$. Let $u(\zeta)=\lVert x(\zeta)-v(\zeta)\rVert^{2}$, Borel by Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets (iv). For $\zeta\in D$ one has $y(\zeta)=y'(\zeta)$, so $x-v=(x-y)+(y'-v)$ at $\zeta$ by Euclidean Space $\mathbb{R}^n$ is a Real Vector Space, whence $\lVert x(\zeta)-v(\zeta)\rVert\le f(\zeta)+f'(\zeta)$ by claim 6 of Properties of the Euclidean Norm on $\mathbb{R}^n$: Square, Distance, Vanishing, Coordinate Bounds, Homogeneity, Triangle Inequalities and Expansion of a Square and $u(\zeta)\le(f(\zeta)+f'(\zeta))^{2}$ by claim 2 of Monotonicity of Squaring on the Nonnegative Elements of an Ordered Field; for $\zeta\notin D$, $h(\zeta)=0$. Hence $uh\le(f+f')^{2}h$ pointwise, and by change of variables, (ii), the density formula, monotonicity and claim 2 of Cauchy-Schwarz and Triangle Inequalities for the Mean-Square Norm,

$$I(\pi_{13})=\int\varphi\circ\Phi\,d\theta=\int u\,d\theta=\int uh\,d\theta_{0}\le\int(f+f')^{2}h\,d\theta_{0}=\int(f+f')^{2}\,d\theta=\lVert f+f'\rVert_{2}^{2}\le\bigl(\lVert f\rVert_{2}+\lVert f'\rVert_{2}\bigr)^{2}.$$

So $I(\pi_{13})$ is finite and $\sqrt{I(\pi_{13})}\le\lVert f\rVert_{2}+\lVert f'\rVert_{2}=\sqrt{I(\pi_{12})}+\sqrt{I(\pi_{23})}$ by claim 2 of Monotonicity of Squaring on the Nonnegative Elements of an Ordered Field.

Claim 8. A Lipschitz $f$ is continuous by A Lipschitz Map is Uniformly Continuous and Continuous, hence Borel by Probability Measures on Euclidean Space and Random Vectors: Standing Notation (ii), and being bounded it is integrable with respect to $\mu$ and to $\nu$ by Probability Measures on Euclidean Space and Random Vectors: Standing Notation (iii). By change of variables in the integrable case, $f\circ\mathrm{pr}_{1}$ and $f\circ\mathrm{pr}_{2}$ are $\pi$-integrable with $\int f\circ\mathrm{pr}_{1}\,d\pi=\int f\,d\mu$ and $\int f\circ\mathrm{pr}_{2}\,d\pi=\int f\,d\nu$. For every $z$, $|f(\mathrm{pr}_{1}(z))-f(\mathrm{pr}_{2}(z))|\le L\,\varphi_{0}(z)$ by the Lipschitz property. By claim 2 of the integral theorem, then monotonicity and homogeneity of the nonnegative integral,

$$\Bigl|\int f\,d\mu-\int f\,d\nu\Bigr|=\Bigl|\int(f\circ\mathrm{pr}_{1}-f\circ\mathrm{pr}_{2})\,d\pi\Bigr|\le\int|f\circ\mathrm{pr}_{1}-f\circ\mathrm{pr}_{2}|\,d\pi\le L\int\varphi_{0}\,d\pi .$$

Finally, on the probability space $(\mathbb{R}^{d+d},\mathcal{B}(\mathbb{R}^{d+d}),\pi)$ the random variable $\varphi_{0}$ is square-integrable with $\lVert\varphi_{0}\rVert_{2}=\sqrt{I(\pi)}$, as $\varphi_{0}^{2}=\varphi$, and the constant $1$ is square-integrable with $\lVert1\rVert_{2}=1$ (The Integral of an Indicator Function is the Measure of the Set with $A=\mathbb{R}^{d+d}$); so $\int\varphi_{0}\,d\pi=\mathbb{E}_{\pi}[\varphi_{0}\cdot1]\le\lVert\varphi_{0}\rVert_{2}\lVert1\rVert_{2}=\sqrt{I(\pi)}$ by claim 1 of Cauchy-Schwarz and Triangle Inequalities for the Mean-Square Norm, the expectation being the nonnegative integral by the identification recorded at the start. Combining, $|\int f\,d\mu-\int f\,d\nu|\le L\sqrt{I(\pi)}$ by claim 5 of Elementary Arithmetic in an Ordered Field.