  For probability measures $\mu,\nu$ and $\rho$ define the cost functionals
\begin{align*} C(\mu,\rho):=\sup_{\pi\in \Pi(\mu,\rho)} \int \langle
x,y\rangle\, \pi(dx,dy),\quad C(\nu,\rho):=\sup_{\pi\in \Pi(\nu,\rho)} \int
\langle x,y\rangle\, \pi(dx,dy), \end{align*} where $\langle\cdot,
\cdot\rangle$ denotes the scalar product and $\Pi(\cdot,\cdot)$ is the set of
couplings. We show that two probability measures $\mu$ and $\nu$ on
$\mathbb{R}^d$ with finite first moments are in convex order (i.e.
$\mu\preceq_c\nu$) iff $C(\mu,\rho)\le C(\nu,\rho)$ holds for all probability
measures $\rho$ on $\mathbb{R}^d$ with bounded support. This generalizes a
result by Carlier. Our proof relies on a quantitative bound for the infimum of
$\int f\,d\nu -\int f\,d\mu$ over all $1$-Lipschitz functions $f$, which is
obtained through optimal transport duality and Brenier's theorem. Building on
this result, we derive new proofs of well-known one-dimensional
characterizations of convex order. We also describe new computational methods
for investigating convex order and applications to model-independent arbitrage
strategies in mathematical finance.
