[proofplan]
[Brenier's theorem](/theorems/7477) gives existence of at least one optimal quadratic transport plan induced by the gradient of a convex potential. To prove uniqueness, take two optimal plans and average them; the average is still optimal because the cost functional is affine in the plan. Quadratic [Kantorovich duality](/theorems/6799) places the averaged optimal plan on the subdifferential graph of a convex potential, and disintegration then shows that almost every conditional measure lives inside one subdifferential. Since convex functions are differentiable Lebesgue-a.e. and $\mu \ll \mathcal{L}^n$, those subdifferentials are singletons $\mu$-a.e.; hence the averaged plan is graph-induced, forcing both original optimal plans to be the same graph plan.
[/proofplan]
[step:Use Brenier's theorem to obtain an optimal graph plan]
Let $\operatorname{id}_{\mathbb{R}^n}:\mathbb{R}^n \to \mathbb{R}^n$ denote the identity map. By Brenier's theorem for quadratic cost, applied to $\mu,\nu \in \mathcal{P}_2(\mathbb{R}^n)$ with $\mu \ll \mathcal{L}^n$, there exists a proper lower semicontinuous convex function $\varphi_0:\mathbb{R}^n \to (-\infty,\infty]$, finite and differentiable $\mu$-a.e., such that the Borel map $T_0:\mathbb{R}^n \to \mathbb{R}^n$, defined by $T_0(x)=\nabla\varphi_0(x)$ on the differentiability set of $\varphi_0$ and arbitrarily on its $\mu$-null complement, satisfies
\begin{align*}
\gamma_0 := (\operatorname{id}_{\mathbb{R}^n},T_0)_\#\mu \in \Pi(\mu,\nu).
\end{align*}
Moreover, $\gamma_0$ is optimal for the quadratic Kantorovich problem. Hence the infimum is attained.
[/step]
[step:Average two optimal plans and keep optimality]
Let $\pi_1,\pi_2 \in \Pi(\mu,\nu)$ be optimal plans. Define the averaged plan
\begin{align*}
\pi := \frac{1}{2}\pi_1+\frac{1}{2}\pi_2.
\end{align*}
Since marginals are linear under convex combinations, $\pi \in \Pi(\mu,\nu)$. Let $m$ be the optimal value:
\begin{align*}
m := \int_{\mathbb{R}^n \times \mathbb{R}^n} |x-y|^2 \, d\pi_1(x,y) = \int_{\mathbb{R}^n \times \mathbb{R}^n} |x-y|^2 \, d\pi_2(x,y).
\end{align*}
By linearity of integration with respect to finite measures,
\begin{align*}
\int_{\mathbb{R}^n \times \mathbb{R}^n} |x-y|^2 \, d\pi(x,y) = \frac{1}{2}m+\frac{1}{2}m = m.
\end{align*}
Thus $\pi$ is also optimal.
[/step]
[step:Place the averaged optimal plan on a convex subdifferential graph]
Since $\mu,\nu\in\mathcal{P}_2(\mathbb{R}^n)$, every admissible coupling has finite quadratic cost bounded by the second moments of $\mu$ and $\nu$. Therefore quadratic Kantorovich duality and [complementary slackness](/theorems/2559) apply to the optimal plan $\pi$. In Brenier form, this yields a proper lower semicontinuous convex function
\begin{align*}
\varphi:\mathbb{R}^n \to (-\infty,\infty]
\end{align*}
such that the averaged optimal plan $\pi$ is concentrated on the subdifferential graph
\begin{align*}
\Gamma_\varphi := \{(x,y)\in \mathbb{R}^n \times \mathbb{R}^n : y \in \partial\varphi(x)\}.
\end{align*}
That is,
\begin{align*}
\pi(\Gamma_\varphi)=1.
\end{align*}
[guided]
The reason for passing to the averaged plan is that optimality is preserved by convex combinations, while duality gives a rigid geometric description of every optimal plan. We first check that the duality theorem is applicable. The hypotheses $\mu,\nu\in\mathcal{P}_2(\mathbb{R}^n)$ imply that the quadratic cost has finite integral on admissible couplings; for example, $|x-y|^2\leq 2|x|^2+2|y|^2$, and the two right-hand moments are finite under any $\rho\in\Pi(\mu,\nu)$. Thus the quadratic Kantorovich problem is finite, and optimality of $\pi$ permits applying quadratic Kantorovich duality with complementary slackness.
In the quadratic-cost Brenier formulation of duality, complementary slackness provides a proper lower semicontinuous convex function
\begin{align*}
\varphi:\mathbb{R}^n \to (-\infty,\infty]
\end{align*}
such that $\pi$ is concentrated on the dual equality set. For the cost $|x-y|^2$, the usual conversion from the quadratic dual potentials to a convex Brenier potential rewrites that equality set exactly as the subdifferential relation $y\in\partial\varphi(x)$. Therefore, if we define
\begin{align*}
\Gamma_\varphi := \{(x,y)\in \mathbb{R}^n \times \mathbb{R}^n : y \in \partial\varphi(x)\},
\end{align*}
then complementary slackness gives
\begin{align*}
\pi(\Gamma_\varphi)=1.
\end{align*}
This is the point at which optimality becomes a pointwise geometric constraint: the second coordinate $y$ is not arbitrary once the first coordinate $x$ is fixed; it must lie in the convex set $\partial\varphi(x)$.
[/guided]
[/step]
[step:Disintegrate the averaged plan over its first marginal]
Since $\mathbb{R}^n$ is a standard Borel space and the first marginal of $\pi$ is $\mu$, the disintegration theorem gives a $\mu$-a.e. uniquely determined family of Borel probability measures $(\pi_x)_{x\in\mathbb{R}^n}$ on $\mathbb{R}^n$ such that, for every nonnegative Borel function
\begin{align*}
F:\mathbb{R}^n \times \mathbb{R}^n \to [0,\infty],
\end{align*}
one has
\begin{align*}
\int_{\mathbb{R}^n \times \mathbb{R}^n} F(x,y)\,d\pi(x,y) = \int_{\mathbb{R}^n}\left(\int_{\mathbb{R}^n} F(x,y)\,d\pi_x(y)\right)d\mu(x).
\end{align*}
Apply this identity to the indicator function $\mathbb{1}_{\Gamma_\varphi^c}:\mathbb{R}^n \times \mathbb{R}^n \to \{0,1\}$. Since $\pi(\Gamma_\varphi^c)=0$,
\begin{align*}
0 = \int_{\mathbb{R}^n} \pi_x(\mathbb{R}^n \setminus \partial\varphi(x))\,d\mu(x).
\end{align*}
The integrand is nonnegative, so
\begin{align*}
\pi_x(\partial\varphi(x))=1
\end{align*}
for $\mu$-a.e. $x\in\mathbb{R}^n$.
[/step]
[step:Use differentiability of convex functions to make the conditional measures Dirac masses]
From the previous step, $\pi_x(\partial\varphi(x))=1$ for $\mu$-a.e. $x$, so $\partial\varphi(x)\neq\varnothing$ for $\mu$-a.e. $x$. In particular, $\varphi(x)<\infty$ for $\mu$-a.e. $x$. Convex functions are differentiable $\mathcal{L}^n$-a.e. on their effective domain, and $\mu \ll \mathcal{L}^n$ transfers this null exceptional set to a $\mu$-null set. Let $D_\varphi \subset \mathbb{R}^n$ denote the Borel set of points at which $\varphi$ is finite and differentiable. Then $\mu(D_\varphi)=1$, and for every $x\in D_\varphi$,
\begin{align*}
\partial\varphi(x)=\{\nabla\varphi(x)\}.
\end{align*}
Combining this singleton identity with $\pi_x(\partial\varphi(x))=1$, we obtain
\begin{align*}
\pi_x = \delta_{\nabla\varphi(x)}
\end{align*}
for $\mu$-a.e. $x\in\mathbb{R}^n$.
[guided]
The disintegration step reduced uniqueness to the size of the set $\partial\varphi(x)$. If $\partial\varphi(x)$ contains many points, then the conditional measure $\pi_x$ could split mass among them. The absolute continuity hypothesis prevents this for almost every $x$.
The previous disintegration step tells us that $\pi_x(\partial\varphi(x))=1$ for $\mu$-a.e. $x$. Since $\pi_x$ is a probability measure, this forces $\partial\varphi(x)\neq\varnothing$ for $\mu$-a.e. $x$. A point with nonempty convex subdifferential lies in the effective domain of $\varphi$, so $\varphi(x)<\infty$ for $\mu$-a.e. $x$.
Now we apply the differentiability theorem for convex functions: a convex function is differentiable $\mathcal{L}^n$-a.e. on its effective domain. Because $\mu \ll \mathcal{L}^n$, every $\mathcal{L}^n$-null set is also $\mu$-null. Therefore the set
\begin{align*}
D_\varphi := \{x\in\mathbb{R}^n : \varphi \text{ is finite and differentiable at } x\}
\end{align*}
satisfies $\mu(D_\varphi)=1$.
At each $x\in D_\varphi$, convex analysis gives the exact identification of the subdifferential:
\begin{align*}
\partial\varphi(x)=\{\nabla\varphi(x)\}.
\end{align*}
From the previous step we already know that $\pi_x(\partial\varphi(x))=1$ for $\mu$-a.e. $x$. Hence, for $\mu$-a.e. $x\in D_\varphi$,
\begin{align*}
\pi_x(\{\nabla\varphi(x)\})=1.
\end{align*}
A probability measure assigning mass $1$ to a singleton is the Dirac mass at that point, so
\begin{align*}
\pi_x = \delta_{\nabla\varphi(x)}
\end{align*}
for $\mu$-a.e. $x\in\mathbb{R}^n$.
[/guided]
[/step]
[step:Identify the averaged plan as a graph pushforward]
Since the gradient of a convex function is Borel on its differentiability set, choose a Borel extension
\begin{align*}
T:\mathbb{R}^n \to \mathbb{R}^n
\end{align*}
such that $T(x)=\nabla\varphi(x)$ on $D_\varphi$ and define it arbitrarily on $\mathbb{R}^n\setminus D_\varphi$. The disintegration formula and the identity $\pi_x=\delta_{T(x)}$ for $\mu$-a.e. $x$ give, for every nonnegative Borel function
\begin{align*}
F:\mathbb{R}^n \times \mathbb{R}^n \to [0,\infty],
\end{align*}
the equality
\begin{align*}
\int_{\mathbb{R}^n \times \mathbb{R}^n} F(x,y)\,d\pi(x,y) = \int_{\mathbb{R}^n} F(x,T(x))\,d\mu(x).
\end{align*}
This is precisely
\begin{align*}
\pi = (\operatorname{id}_{\mathbb{R}^n},T)_\#\mu.
\end{align*}
[/step]
[step:Force both original optimal plans to equal the same graph plan]
Let
\begin{align*}
G_T := \{(x,y)\in \mathbb{R}^n \times \mathbb{R}^n : y=T(x)\}
\end{align*}
be the graph of $T$. Since $\pi=(\operatorname{id}_{\mathbb{R}^n},T)_\#\mu$, we have $\pi(G_T)=1$. From
\begin{align*}
\pi=\frac{1}{2}\pi_1+\frac{1}{2}\pi_2,
\end{align*}
we get, for every Borel set $A\subset \mathbb{R}^n \times \mathbb{R}^n$,
\begin{align*}
\pi_i(A)\leq 2\pi(A)
\end{align*}
for $i\in\{1,2\}$. Applying this to $A=G_T^c$ gives
\begin{align*}
\pi_i(G_T^c)=0
\end{align*}
for $i\in\{1,2\}$.
Now let $\rho\in\Pi(\mu,\nu)$ be any plan concentrated on $G_T$. For Borel sets $B,C\subset\mathbb{R}^n$,
\begin{align*}
\rho(B\times C)=\rho((B\times C)\cap G_T)=\mu(B\cap T^{-1}(C)).
\end{align*}
The right-hand side is exactly the value of $(\operatorname{id}_{\mathbb{R}^n},T)_\#\mu$ on $B\times C$. Since Borel rectangles generate the product Borel $\sigma$-algebra, $\rho=(\operatorname{id}_{\mathbb{R}^n},T)_\#\mu$. Applying this to $\rho=\pi_1$ and $\rho=\pi_2$ yields
\begin{align*}
\pi_1=\pi_2=(\operatorname{id}_{\mathbb{R}^n},T)_\#\mu.
\end{align*}
Thus the optimal plan is unique.
[/step]
[step:Conclude uniqueness of the Brenier gradient]
The first step produced at least one optimal Brenier graph plan, and the previous step shows that every optimal plan equals the unique plan
\begin{align*}
\pi_*=(\operatorname{id}_{\mathbb{R}^n},T)_\#\mu.
\end{align*}
Therefore $\pi_*=(\operatorname{id}_{\mathbb{R}^n},T)_\#\mu$, where $T$ is a Borel representative of $\nabla\varphi$ and $\varphi$ is a convex Brenier potential. Equivalently, $\pi_*=(\operatorname{id}_{\mathbb{R}^n},\nabla\varphi)_\#\mu$ with this representative understood.
Finally, let $\psi:\mathbb{R}^n \to (-\infty,\infty]$ be another convex Brenier potential, and let $S:\mathbb{R}^n \to \mathbb{R}^n$ be a Borel representative of $\nabla\psi$ defined $\mu$-a.e. Define
\begin{align*}
G_S := \{(x,y)\in \mathbb{R}^n \times \mathbb{R}^n : y=S(x)\}
\end{align*}
to be the graph of $S$. If $(\operatorname{id}_{\mathbb{R}^n},S)_\#\mu$ is optimal, then uniqueness of the optimal plan gives
\begin{align*}
(\operatorname{id}_{\mathbb{R}^n},S)_\#\mu=(\operatorname{id}_{\mathbb{R}^n},T)_\#\mu.
\end{align*}
The common measure is concentrated on both graphs $G_S$ and $G_T$. Hence
\begin{align*}
\mu(\{x\in\mathbb{R}^n:S(x)=T(x)\})=1.
\end{align*}
Since $S=\nabla\psi$ $\mu$-a.e. and $T=\nabla\varphi$ $\mu$-a.e., this proves
\begin{align*}
\nabla\psi=\nabla\varphi
\end{align*}
$\mu$-a.e. This completes the proof.
[/step]
