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 (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 §pushforward. 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 §functions, 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 §cost, each such $\pi$ being a coupling of its own push-forwards by Couplings of Two Probability Measures on Euclidean Space and Their Quadratic Cost §coupling. 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 §projections, 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 Function and the 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 §pairs and Pairs of Euclidean Points: Coordinate Projections, Pairings, the Product Measure on a Euclidean Space, Borel Norm Functions and Finite Sets §product (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 §coupling.

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 §pairing, 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 §projections, $\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 Elementary Properties of the Euclidean Norm on $\mathbb{R}^n$; 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 §functions, $\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 §borel-maps), 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 §moment. 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 §functions 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 §pairs, 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 Elementary Properties of the Euclidean Norm on $\mathbb{R}^n$, 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 §borel-maps, 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 §borel-maps, 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 Elementary Properties of the Euclidean Norm on $\mathbb{R}^n$ 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 Elementary Properties of the Euclidean Norm on $\mathbb{R}^n$, 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 §projections 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 Measure Spaces and the Lebesgue Integral: Standing Notation §measurable 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, Absolute Values, and Pointwise Limits of Measurable Real-Valued Functions. 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 Elementary Properties of the Euclidean Norm on $\mathbb{R}^n$ 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 §spaces, and finite intersections, unions and differences of Borel sets being Borel by Sigma-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})$,

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 Elementary Properties of the Euclidean Norm on $\mathbb{R}^n$. 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,

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 §finite-sets, 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 §finite-sets 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 §pairs; 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 §borel-maps). 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 §product). 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 §functions and vanishing exactly on $D$ by claim 3 of Elementary Properties of the Euclidean Norm on $\mathbb{R}^n$, 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, Absolute Values, and Pointwise Limits of Measurable Real-Valued Functions), 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)$:

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 §product 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 §projections. By the Tonelli clause of Tonelli and Fubini 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: Almost-Everywhere Comparison, Markov's Inequality, and Dominated Convergence Almost Everywhere §comparison, 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 §pairs, 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 §functions. 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 Elementary Properties of the Euclidean Norm on $\mathbb{R}^n$ 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,

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, hence Borel by Probability Measures on Euclidean Space and Random Vectors: Standing Notation §borel-maps, and being bounded it is integrable with respect to $\mu$ and to $\nu$ by Probability Measures on Euclidean Space and Random Vectors: Standing Notation §measures. 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,

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.