# Kantorovich vs. Monge: A Numerical Classification of Extremal Multi-Marginal
Mass Transports on Finite State Spaces

Daniela Vögler Faculty of Mathematics, Technische Universität München,
<EMAIL_ADDRESS>+49-89-289-17900 (Tel.), +49-89-289-17932 (Fax)

[2em]2em Abstract. We analyze the validity of Monge’s ansatz regarding the
symmetric multi-marginal Kantorovich optimal transport problem on finite state
spaces, with uniform-marginal constraint.
The class of Monge states reduces the number of unknowns from combinatorial in
both $N$ and $\ell$ to linear in $N$ and $\ell$, where $N$ is the number of
marginals and $\ell$ the number of states. Unfortunately, this low-dimensional
ansatz space is insufficient, i.e., there exist cost functions such that the
corresponding multi-marginal optimal transport problem does not admit a Monge-
type optimizer. We will analyze this insufficiency numerically by utilizing
the convex geometry of the set of admissible trial states for symmetric
respectively pairwise-symmetric cost functions.
We further consider a model problem of optimally coupling $N$ marginals on
three states with respect to a pairwise-symmetric cost function. The
restriction to a state space of three elements allows us to visually compare
Kantorovich’s to Monge’s ansatz space. This visual comparison includes a
consideration of the volumetric ratio.  Keywords: optimal transport, Monge’s
ansatz, N-representability, Birkhoff-von Neumann theorem, density functional
theory
Mathematics Subject Classification: 49K30, 52B12, 52-04, 49S05
Acknowledgments: The author thanks Gero Friesecke and Sören Behr for helpful
discussions.

##  1 Introduction

In general multi-marginal Kantorovich optimal transport problems aim at
coupling $N$ probability measures $\lambda^{(1)},\ldots,\lambda^{(N)}$
optimally with respect to a given cost function $c$ (see (1.1) for a discrete
symmetric optimal transport problem). These problems arise in various fields
of research, ranging from economics [8, 11] through mathematical finance [2,
19] and image processing [1, 32] to electronic structure [13, 6].

One of the central questions in the theory of optimal transportation is: Under
which assumptions exists an optimal coupling that is supported on a graph
(over the first variable)? Such optimizers are then called Monge-solutions
(see (1.5) - (1.7)). In the case of two marginals this question is well
understood; the existence of Monge-solutions is always respectively under very
general conditions guaranteed (see the renowned Birkhoff-von Neumann theorem
[4, 36] regarding finite state spaces respectively, e.g., [35] regarding
continuous state spaces). For multiple marginals the understanding of this
question does not reach the same generality. However, there are isolated
examples for Monge- and non-Monge-solutions. For the former, see [22, 7, 28,
13, 6, 12] as well as the fundamental paper by Gangbo and Święch [20] for an
interesting selection. For the latter, we refer the interested reader to [9,
29, 17, 30, 14, 27, 21, 16] regarding continuous state spaces as well as to
[15, 25, 26, 16] regarding finite state spaces.

The goal of this paper is to investigate the validity of Monge’s ansatz
regarding the symmetric multi-marginal Kantorovich optimal transport problem
on finite state spaces

$\textrm{Minimize
}\int_{X^{N}}c(x_{1},\ldots,x_{N})d\gamma(x_{1},\ldots,x_{N})\\\ \textrm{ over
}\gamma\in\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)\textrm{ subject to
}\gamma\mapsto\overline{\lambda}.$ (1.1)

Here $X$ denotes a finite state space as defined in (1.2),
$c:X^{N}\to\mathbb{R}\cup\\{+\infty\\}$ an arbitrary symmetric cost function,
$\overline{\lambda}$ the uniform marginal as defined in (1.3) and
$\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$ the set of symmetric
probability measures on $X^{N}$, where a probability measure $\gamma$ on
$X^{N}$ is symmetric if

$\gamma\left(A_{1}\times\dots\times
A_{N}\right)=\gamma\left(A_{\sigma(1)}\times\dots\times
A_{\sigma(N)}\right)\\\ \textrm{ for all subsets }A_{1},\ldots,A_{N}\textrm{
of }X\textrm{ and all permutations }\sigma.$

Any $\gamma\in\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$ fulfills
$\gamma\mapsto\overline{\lambda}$ if and only if $\gamma$ has equal one-point
marginals $\overline{\lambda}$, i.e.,

$\gamma\left(X^{k-1}\times A_{k}\times
X^{N-k}\right)=\overline{\lambda}(A_{k})\textrm{ for all subsets
}A_{k}\textrm{ of }X\textrm{ and all }k\in\\{1,\ldots,N\\}.$

Multi-marginal optimal transport problems of form (1.1) were already
considered in [18] and [16]. While [16] discusses the validity of Monge’s
approach in the setting of 3 marginals and 3 sites, [18] introduces a
sufficient ansatz space for problem (1.1) (see Section 2 as well as Remark 3.5
for information about the content of these papers). The present paper
accompanies these previous considerations. In particular, some of the used
nomenclature and notation is already introduced there.

For finite state spaces

$X=\\{a_{1},\ldots,a_{\ell}\\}$ (1.2)

consisting of $\ell$ distinct points $a_{1},\ldots,a_{\ell}$, the uniform
probability measure

$\overline{\lambda}=\sum_{i=1}^{\ell}\frac{1}{\ell}\delta_{a_{i}}$ (1.3)

on $X$ is the prototypical marginal. The corresponding multi-marginal
Kantorovich optimal transport problems, i.e., problems of form (1.1) with
$\mathcal{P}\left(X^{N}\right)$ replacing
$\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$, appear directly as assignment
problems (see [34, 5] for reviews) and arise from continuous problems via
equi-mass discretization [10]. Then, an optimal coupling
$\gamma\in\mathcal{P}\left(X^{N}\right)$ of the $N$ marginals
$\overline{\lambda},\ldots,\overline{\lambda}$ is a Monge-solution if

$\gamma=\sum_{\nu=1}^{\ell}\frac{1}{\ell}\delta_{T_{1}(a_{\nu})}\otimes\dots\otimes\delta_{T_{N}(a_{\nu})}\textrm{
for }N\textrm{ permutations }T_{1},\ldots,T_{N}:X\to X,$ (1.4)

where $T:X\to X$ is a permutation if there exists a permutation of indices
$\tau:\\{1,\ldots,\ell\\}\to\\{1,\ldots,\ell\\}$ such that
$T\left(a_{\nu}\right)=a_{\tau(\nu)}$ for all $\nu\in\\{1,\ldots,\ell\\}$.
Demanding that the $T_{k}$s are permutations ensures that $\gamma$ is indeed a
coupling of $\overline{\lambda},\ldots,\overline{\lambda}$: $T:X\to X$ is a
permutation if and only if it pushes the uniform measure forward to itself,
i.e., $T_{\\#}\overline{\lambda}=\overline{\lambda}$. Here for any probability
measure $\lambda=\sum_{\nu=1}^{\ell}\lambda_{\nu}\delta_{a_{\nu}}$ on $X$ and
any map $T:X\to X$ the push-forward $T_{\\#}\lambda$ of $\lambda$ under $T$ is
defined by
$T_{\\#}\lambda=\sum_{\nu=1}^{\ell}\lambda_{\nu}\delta_{T(a_{\nu})}$. One may
choose $T_{1}=\textrm{id}$, i.e., $T_{1}(a)=a$ for all $a\in X$, by re-
ordering the sum in (1.4).
Regarding (1.1) an admissible trial state $\hat{\gamma}$ is referred to as a
(symmetrized) Monge state if it is the symmetrization of a probability measure
$\gamma$ of form (1.4), i.e.,

$\hat{\gamma}=\sum_{\nu=1}^{\ell}\frac{1}{\ell}S\delta_{T_{1}(a_{\nu})}\otimes\dots\otimes\delta_{T_{N}(a_{\nu})}$
(1.5)

such that

$T_{k}\\#\overline{\lambda}=\overline{\lambda}\textrm{ for all
}k\in\\{1,\ldots,N\\},$ (1.6)

or equivalently,

$T_{1},\ldots,T_{N}:X\to X\textrm{ are permutations.}$ (1.7)

Here $S$ denotes the linear symmetrization operator in $N$ variables as
defined in (2.4). Probability measures on $X^{N}$ of form (1.5)-(1.7) are also
said to be of Monge-type or in Monge-form and restricting the minimization
problem (1.1) to such measures yields the corresponding Monge problem.

A guarantee for the existence of a Monge-solution would allow us to utilize
the Monge approach (1.5)-(1.7) in order to reduce the number of unknowns in
(1.1) from combinatorial in both $N$ and $\ell$ to linear in $N$ and $\ell$.
Unfortunately, for multiple marginals the Monge ansatz is, again, not
sufficient: In [16] an example for a symmetric cost function
$c:X^{3}\to\mathbb{R}\cup\\{+\infty\\}$ such that there exists an unique
minimizer of problem (1.1) that is not of Monge-type is given. In Section 2,
we will quantify this insufficiency of Monge’s approach by utilizing the
convex geometry of the set of admissible trial states.

Note that the restriction to symmetric probability measures in problem (1.1)
is motivated by a physical application. Modeling the electronic structure of a
molecule with $N$ electrons in a discretized setting, is a prototypical
application of multi-marginal optimal transport on finite state spaces. In
this context, $X$ corresponds to a set of $\ell$ discretization points in
$\mathbb{R}^{3}$ and any coupling $\gamma$ of the $N$ marginals
$\overline{\lambda},\ldots,\overline{\lambda}$ describes a joint probability
distribution regarding the electron positions in an $N$-electron molecule.
Then the marginal condition ensures that each discretization point is occupied
equally often and the cost function $c:X^{N}\to\mathbb{R}\cup\\{+\infty\\}$
embodies the electron interaction energy. As electrons are indistinguishable
the considered cost functions are usually symmetric, i.e., invariant under
argument permutation. These symmetric cost functions are ’dual’ to the set of
symmetric probability measures on the product space $X^{N}$ in the following
sense: There always exists an optimal coupling of
$\overline{\lambda},\ldots,\overline{\lambda}$ that is symmetric.

The interaction energy between electrons often displays pairwise structure,
i.e., $c(x_{1},\ldots,x_{N})=\sum_{1\leq i<j\leq N}v(x_{i},x_{j})$, with the
Coulomb cost $\sum_{1\leq i<j\leq N}\frac{1}{|x_{i}-x_{j}|}$ being the
prototypical example. Here $|\cdot|$ is the Euclidean norm in
$\mathbb{R}^{d}$. As discussed in Section 3, this pairwise structure allows us
to reformulate the multi-marginal optimal transport problem (1.1) as

$\textrm{Minimize }\int_{X^{2}}v(x,y)d\mu(x,y)\textrm{ over
}\mu\in\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)\textrm{ subject to
}\mu\mapsto\overline{\lambda}.$ (1.8)

This reformulated problem was initially introduced in [17]. In (1.8),
$\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)\subseteq\mathcal{P}\left(X^{2}\right)$
can be interpreted as a ’reduced version’ of
$\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$.

The above mentioned cost function $c:X^{3}\to\mathbb{R}\cup\\{+\infty\\}$ in
[16], which establishes the insufficiency of Monge’s ansatz regarding problem
(1.1), is of pairwise structure. It immediately follows that there is no
’reduced Monge state’ that solves problem (1.8). In this sense the Monge
ansatz is also insufficient for the reformulated problem (1.8). In Section 3,
we will make use of the convex geometry of the new set of admissible trial
states in order to quantify this insufficiency of Monge’s ansatz.

Finally, in Section 4, we will consider a model problem of optimally coupling
the $N$ marginals $\overline{\lambda},\ldots,\overline{\lambda}$ with respect
to a cost function of pairwise symmetric structure, where the finite state
space $X$ consists only of three states, i.e., $X=\\{a_{1},a_{2},a_{3}\\}$. In
particular $\overline{\lambda}$ is the uniform probability measure on
$\\{a_{1},a_{2},a_{3}\\}$ and the domain of the cost function is given by
$\\{a_{1},a_{2},a_{3}\\}^{N}$. The present setting allows us to draw a visual
comparison between Kantorovich’s and Monge’s ansatz as depicted in Figure 1.
We further compare both optimal transport approaches by volume of (the convex
hull of) the respective set of admissible trial states and establish a
computationally simple upper bound on the optimal value in (1.8).

(a)

(b)

(c)

(d)

Figure 1: The reduced Kantorovich (see (3.4)) respectively Monge (see (3.8))
polytope for $N$ marginals and $3$ states is visualized for $N=3$ (top-left),
$4$ (top-right), $6$ (bottom-left) and $10$ (bottom-right) in green
respectively red. The elements $\left(\mu_{ij}\right)_{i,j=1}^{3}$ of the
polytopes are parametrized by their off-diagonal entries $\mu_{12},\mu_{13}$
and $\mu_{23}$. If a two-dimensional face of the reduced Monge polytope
belongs to the boundary of the reduced Kantorovich polytope the occupied area
is depicted in red.

##  2 Classification of the Extreme Points of a Kantorovich Polytope

Throughout the paper we will consider the finite state space $X$ given by
(1.2). We will denote the set of probability measures on $X$ as
$\mathcal{P}(X)$. Each such probability measure $\lambda\in\mathcal{P}(X)$ can
be canonically identified with a vector in $\mathbb{R}^{\ell}$ in the
following manner

$\lambda_{i}:=\lambda({a_{i}}).$

The vector $(\lambda_{1},\ldots,\lambda_{\ell})$ fulfills
$\sum_{i=1}^{\ell}\lambda_{i}=1$ and $\lambda_{i}\geq 0$ for
$i\in\\{1,\ldots,\ell\\}$ and is therefore an element of the unit simplex. The
probability measure $\lambda$ can then be written in the following manner

$\lambda=\sum_{i=1}^{\ell}\lambda_{i}\delta_{i},$

where here and below we use the shorthand notations

$\delta_{i}:=\delta_{a_{i}},\hskip 56.9055pt\delta_{i_{1}\ldots
i_{N}}:=\delta_{a_{i_{1}}}\otimes\ldots\otimes\delta_{a_{i_{N}}}$ (2.1)

for $a_{i}\in X$ being a single point in the finite state space and
$(a_{i_{1}},\ldots,a_{i_{N}})$ being an element of the product space $X^{N}$.

As announced in the Introduction, we will now take a closer look at the set of
admissible trial states of problem (1.1), i.e., the set

$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right):=\left\\{\gamma\in\mathcal{P}_{\textrm{sym}}\left(X^{N}\right):\gamma\mapsto\overline{\lambda}\right\\}.$
(2.2)

From this point on, we will refer to the elements of this set as symmetric
Kantorovich couplings. The set
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ will be
called (symmetric) Kantorovich polytope for $N$ marginals and $\ell$ states.
As within this paper we focus our attention on the symmetric case, the term
symmetric will be dropped from time to time. It is easy to see that, as a
result of the linearity of the marginal constraint and the finiteness of the
state space $X$,
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ is a compact
and convex set in $\mathbb{R}^{{\ell}^{N}}$ and therefore by Minkowski’s
theorem (see, e.g., [23]) the convex hull of its extreme points.

Recall the following basic definitions and notions of convexity (see, e.g.,
[23, 33]). For $y_{1},\ldots,y_{n}\in\mathbb{R}^{m}$ and
$\lambda_{1},\ldots,\lambda_{n}\geq 0$ such that $\sum_{i=1}^{n}\lambda_{i}=1$

$\lambda_{1}y_{1}+\dots+\lambda_{n}y_{n}=\sum_{i=1}^{n}\lambda_{i}y_{i}$

is called a convex combination of the points $y_{1},\ldots,y_{n}$. A subset
$K\subseteq\mathbb{R}^{m}$ is called convex if for each finite selection of
points in $K$ each possible convex combination of these points is again
contained in $K$. For a subset $V\subseteq\mathbb{R}^{m}$ the convex hull of
$V$, denoted as

$\textrm{conv}(V),$

corresponds to the set of all possible convex combinations of a finite
selection of points in $V$. Obviously a set $K\subseteq\mathbb{R}^{m}$ is
convex if and only if it is equal to its convex hull. Finally an element $k$
of the convex set $K\subseteq\mathbb{R}^{m}$ is called an extreme point if
$k=\lambda_{1}y_{1}+\lambda_{2}y_{2}$ for some $y_{1},y_{2}\in K$ and
$\lambda_{1},\lambda_{2}>0$ such that $\lambda_{1}+\lambda_{2}=1$ implies
$y_{1}=k=y_{2}$. For a considered convex set $K$ the set of extreme points
will from now on be denoted as $\textrm{ext}(K)$.

As $\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ is equal
to the convex hull of its extreme points, we can use the extreme points to
describe the convex structure of the set of symmetric Kantorovich couplings.
Now it follows by a simple contradiction argument that for any given linear
objective function there is always an optimizer that is an extreme point.
Moreover, in our setting of finite states spaces, for any extreme point
$\gamma^{*}$ there is a function $c:X^{N}\to\mathbb{R}$ such that

$\int_{X^{N}}c(x_{1},\ldots,x_{N})d\gamma(x_{1},\ldots,x_{N})>\int_{X^{N}}c(x_{1},\ldots,x_{N})d\gamma^{*}(x_{1},\ldots,x_{N})\\\
\textrm{ for any
}\gamma\in\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right),$
(2.3)

i.e., there is a cost function such that $\gamma^{*}$ is the unique optimizer
of the corresponding optimal transport problem. This is a result of the fact
that $\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ is a
bounded polyhedron, i.e., a polytope, of finite dimension and therefore only
possesses finitely many extreme points each of whom is itself an exposed point
(see, e.g., [33]), i.e., a point in the set
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ that
fulfills (2.3) for some cost function $c:X^{N}\to\mathbb{R}$. As for any cost
function $c:X^{N}\to\mathbb{R}$ there is always an optimizer that is an
extreme point of
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ and vice
versa for any extreme point $\gamma^{*}$ of
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ there is a
cost function $c:X^{N}\to\mathbb{R}$ such that $\gamma^{*}$ is the unique
optimizer, analyzing, how many of the extreme points of
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ are of
Monge-type, seems to be the right approach to investigate the validity of
Monge’s approach. Recall that in the given setting a probability measure
$\gamma\in\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ is
said to be of Monge-type or in Monge-form if there are $N$ permutations
$\tau_{1},\ldots,\tau_{N}:\\{1,\ldots,\ell\\}\to\\{1,\ldots,\ell\\}$ such that

$\gamma=\sum_{i=1}^{\ell}\frac{1}{\ell}S\delta_{\tau_{1}(i)\tau_{2}(i)\ldots\tau_{N}(i)},$

where the symmetrization operator
$S:\mathcal{P}\left(X^{N}\right)\to\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$
is defined by

$\left(S\gamma\right)\left(A_{1}\times\dots\times
A_{N}\right)=\frac{1}{N!}\sum_{\sigma\in
S_{N}}\gamma\left(A_{\sigma(1)}\times\dots\times A_{\sigma(N)}\right)\textrm{
for all }A_{1},\ldots,A_{N}\subseteq X$ (2.4)

with $S_{N}$ being the group of all permutations on the set
$\\{1,\ldots,N\\}$.

As in the given setting it is rather inconvenient and long-winded to check
whether a given symmetric Kantorovich coupling is of Monge-type or not, we
will derive in the following an alternative LP-formulation of problem (1.1),
where Monge-states will correspond exactly to (rescaled) integer points in the
corresponding polytope of admissible trial states. We start by taking a closer
look at the convex geometry of the set of symmetric probability measures on
the product space $X^{N}$.

Note that a probability measure
$\gamma\in\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$ is an extreme point of
$\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$ if and only if it is of the
form

$S\delta_{i_{1}\ldots i_{N}}\textrm{ for some }1\leq i_{1}\leq\dots\leq
i_{N}\leq\ell$ (2.5)

(see [18]). Therefore symmetric Kantorovich couplings, which are of Monge-
type, are an average of $\ell$ not necessarily distinct extremal symmetric
probability measures on the product space $X^{N}$ with respect to the uniform
measure. Below we will elaborate further on this characterization of couplings
in Monge-form, which will be the basis for identifying Monge-states with the
(rescaled) integer points in a certain polytope.

From now on, we will denote the set of extremal symmetric probability
measures, i.e., measures of the form (2.5), as $E_{\textrm{sym}}^{N}$. It was
shown in [18] that $E_{\textrm{sym}}^{N}$ contains $\tbinom{N+\ell-1}{N}$
elements.
As for each pair of these extreme points their support is disjoint, one can
immediately deduce the following result.

######  Proposition 2.1.

$\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$ is a simplex, i.e., the
extremal symmetric probability measures on $X^{N}$ are affinely independent.

Hence, for every $\gamma\in\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$ there
is a unique way to represent $\gamma$ as a convex combination of extremal
symmetric probability measures on $X^{N}$, i.e., there is a unique non-
negative coefficient vector $\alpha\in\mathbb{R}^{|E_{\textrm{sym}}^{N}|}$
fulfilling $\sum\alpha_{i_{1}\ldots i_{N}}=1$ such that

$\gamma=\sum_{1\leq i_{1}\leq\dots\leq i_{N}\leq\ell}\alpha_{i_{1}\ldots
i_{N}}S\delta_{i_{1}\ldots i_{N}}.$ (2.6)

As the extreme points of $\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$ can be
parametrized using their one-point marginal, $\alpha$ can be interpreted as a
probability measure on the set of these one-point marginals.
Given the $k$-point marginal map
$M_{k}:\mathcal{P}\left(X^{N}\right)\to\mathcal{P}\left(X^{k}\right)$ for
$1\leq k\leq N-1$ with

$\left(M_{k}\gamma\right)(A):=\gamma\left(A\times X^{N-k}\right)\textrm{ for
all }A\subseteq X^{k}$ (2.7)

for $\gamma\in\mathcal{P}\left(X^{N}\right)$, with the convention $M_{N}=id$,
note that $M_{1}$ is a bijection from the set of extremal symmetric
probability measures on $X^{N}$, i.e., measures of the form (2.5), to the set
of $\frac{1}{N}$-quantized probability measures

$\mathcal{P}_{\frac{1}{N}}\left(X\right):=\left\\{\lambda\in\mathcal{P}\left(X\right):\lambda(\\{i\\})\in\left\\{0,\frac{1}{N},\ldots,1\right\\}\right\\}$
(2.8)

(see [18]). Hereby the one-point marginal of a measure of form (2.5) is an
empirical measure of the indices $(i_{1},\ldots,i_{N})$, it holds

$M_{1}S\delta_{i_{1}\ldots i_{N}}=\frac{1}{N}\sum_{j=1}^{N}\delta_{i_{j}}.$

In the following $\psi_{N}:\mathcal{P}_{\frac{1}{N}}\left(X\right)\to
E_{\textrm{sym}}^{N}$ will denote the corresponding inverse function.
This parametrization gives rise to the coefficients-to-coupling map
$R:\mathcal{P}\left(\mathcal{P}_{\frac{1}{N}}\left(X\right)\right)\to\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$.
It maps an arbitrary probability measure $\alpha$ on
$\mathcal{P}_{\frac{1}{N}}\left(X\right)$, which, via the underlying
parametrization, corresponds to the coefficients in the representation (2.6),
to the corresponding coupling $\gamma$, i.e., in pedestrian notation

$R\alpha=\sum_{\lambda\in\mathcal{P}_{\frac{1}{N}}\left(X\right)}\alpha_{\lambda}\psi_{N}(\lambda),$

or more elegantly

$R\alpha=\int_{\mathcal{P}_{\frac{1}{N}}\left(X\right)}\psi_{N}(\lambda)d\alpha(\lambda).$

As $\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$ is a simplex, $R$ is
bijective. This enables us to establish the following isomorphic relationship
between two alternative formulations of the set of symmetric Kantorovich
couplings.

######  Lemma 2.2 (isomorphic relationship between couplings and
coefficients).

The coefficients-to-coupling map $R$ maps the polytope

$P_{\textrm{coef}}:=\left\\{\alpha\in\mathbb{R}^{|E_{\textrm{sym}}^{N}|}:A\alpha=\overline{\lambda},\alpha\geq
0\right\\}$ (2.9)

linearly and bijectively to the set of symmetric Kantorovich couplings, i.e.,
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ defined in
(2.2). Here $E_{\textrm{sym}}^{N}$ is the set of extremal symmetric
probability measures on $X^{N}$ and $A$ is the matrix in
$\mathbb{R}^{\ell\times|E_{\textrm{sym}}^{N}|}$, whose columns are given by
the elements of $\mathcal{P}_{\frac{1}{N}}\left(X\right)$, i.e.,

$A:=\begin{pmatrix}\lambda_{1}^{(1)}&\lambda_{1}^{(2)}&&\dots&&\lambda_{1}^{(|E_{\textrm{sym}}^{N}|)}\\\
\vdots&\vdots&&&&\vdots\\\
\lambda_{\ell}^{(1)}&\lambda_{\ell}^{(2)}&&\dots&&\lambda_{\ell}^{(|E_{\textrm{sym}}^{N}|)}\end{pmatrix}.$
(2.10)

The corresponding inverse map is also linear.

###### Proof.

Linearity and injectivity of $R$ as a map from $P_{\textrm{coef}}$ to
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ is an
immediate consequence of the linearity and injective of
$R:\mathcal{P}\left(\mathcal{P}_{\frac{1}{N}}\left(X\right)\right)\to\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$
as introduced above. We further know that any
$\gamma\in\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ is
an element of $\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$. Hence, applying
the parametrization of extremal symmetric probability measures on $X^{N}$ via
their one-point marginals, there exist coefficients
$\alpha\in\mathcal{P}\left(\mathcal{P}_{\frac{1}{N}}\left(X\right)\right)$,
which are non-negative and whose entries sum to $1$ such that

$\gamma=\sum_{\lambda\in\mathcal{P}_{\frac{1}{N}}\left(X\right)}\alpha_{\lambda}\psi_{N}(\lambda)$
(2.11)

and therefore $\gamma=R\alpha$ holds. Applying the linear marginal map $M_{1}$
to (2.11) yields the following.

$\overline{\lambda}=\sum_{\lambda\in\mathcal{P}_{\frac{1}{N}}\left(X\right)}\alpha_{\lambda}\lambda$

Therefore $\alpha$ corresponds to an element of $P_{\textrm{coef}}$. This
implies surjectivity of the considered map $R$. Linearity of the corresponding
inverse map is an immediate consequence of the fact that the extremal
symmetric probability measures on $X^{N}$ of the form (2.5) interpreted as
vectors are linearly independent. ∎

Now it is easy to see that the extreme points of $P_{\textrm{coef}}$
correspond exactly to the extremal symmetric Kantorovich couplings, in the
sense that $R$ maps the corresponding sets of extreme points bijectively to
each other. By standard arguments of polyhedral optimization the extreme
points of $P_{\textrm{coef}}$ have a sparse structure, i.e., any extreme point
of $P_{\textrm{coef}}$ can have at most $\ell$, that is the number of states
in the finite state space $X$, non-zero entries (see, e.g., [3]). In [18] it
was shown that this implies that any extremal Kantorovich coupling is a so
called Quasi-Monge state, i.e., of the form

$\sum_{\nu=1}^{\ell}\alpha_{\nu}S\delta_{T_{1}\left(a_{\nu}\right)}\otimes\dots\otimes\delta_{T_{N}\left(a_{\nu}\right)}$

for $N$ maps $T_{1},\ldots,T_{N}:X\to X$ such that

$\frac{1}{N}\sum_{k=1}^{N}T_{k}\\#\alpha=\overline{\lambda}.$

Here we renounce from using the shorthand notations (2.1) in order to make it
easier to draw a comparison with Monge’s approach (1.5)-(1.7). The ansatz
space of Quasi-Monge states increases the number of unknowns only by
$2\cdot\ell$ compared to the class of symmetrized Monge states and as every
extremal Kantorovich coupling is a Quasi-Monge state, this ansatz space always
contains an optimal coupling, in contrast to Monge’s approach. Note further
that obviously every symmetrized Monge state is a Quasi-Monge state. For
further reading on this sufficient low-dimensional enlargement of the class of
symmetrized Monge states we refer the interested reader to [18]. There also a
characterization of Monge states in the given setting was established. A
probability measure on the product space $X^{N}$ is a symmetrized Monge state
if and only if it is a Quasi-Monge state all of whose site weights
$\alpha_{1},\ldots,\alpha_{\ell}$ are equal to $\frac{1}{\ell}$. In summary,
we get the following corollary.

######  Corollary 2.3.

Extremal symmetric Kantorovich couplings correspond exactly, via the
coefficients-to-coupling map $R$, to the extreme points of
$P_{\textrm{coef}}$. Any of these extreme points of $P_{\textrm{coef}}$ is the
coefficient vector of a coupling in Monge-form if and only if it is an integer
vector scaled by the factor $\frac{1}{\ell}$.

This corollary gives us a numerically-convenient way to compute the set of
extremal Kantorovich couplings and check whether they are of Monge-form or
not. In addition we also want to consider Monge’s approach by itself. For this
purpose we introduce the sets

$\mathcal{P}_{\textrm{sym,Monge}}\left(X^{N}\right):=\left\\{\gamma\in\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right):\gamma\textrm{
is of Monge-form }\eqref{eq:MongeOne}-\eqref{eq:MongeThree}\right\\}$ (2.12)

and

$\mathcal{P}_{\textrm{sym,Monge}}^{\textrm{conv}}\left(X^{N}\right):=\textrm{conv}\left(\mathcal{P}_{\textrm{sym,Monge}}\left(X^{N}\right)\right).$
(2.13)

$\mathcal{P}_{\textrm{sym,Monge}}\left(X^{N}\right)$ is the set of all
symmetrized Monge states. In the following we will refer to
$\mathcal{P}_{\textrm{sym,Monge}}^{\textrm{conv}}\left(X^{N}\right)$ as the
(symmetric) Monge polytope for $N$ marginals and $\ell$ states. For simplicity
we will once again drop the term symmetric from time to time. Note that if
there exists an optimizer of Problem (1.1) which is an element of
$\mathcal{P}_{\textrm{sym,Monge}}^{\textrm{conv}}\left(X^{N}\right)$ then
there exists a Monge-type minimizer.

Having the explanations leading up to Corollary 2.3 in mind, it is easy to see
that $\mathcal{P}_{\textrm{sym,Monge}}\left(X^{N}\right)$ corresponds to the
(scaled by $\frac{1}{\ell}$) integer elements of $P_{\textrm{coef}}$. These
can be for example determined by a simple enumeration of all the ordered
choices of $N-1$ permutations interpreted as coefficient vectors in
$P_{\textrm{coef}}$. Checking which of these scaled integer coefficient
vectors are extremal with respect to the convex hull of them as a whole, gives
us the extremal elements of
$\mathcal{P}_{\textrm{sym,Monge}}^{\textrm{conv}}\left(X^{N}\right)$.

The data in Figure 2 was computed using MATLAB [24] and polymake [31].

Figure 2: The number of extreme points of the symmetric Kantorovich polytope
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ respectively
of the symmetric Monge polytope
$\mathcal{P}_{\textrm{sym,Monge}}^{\textrm{conv}}\left(X^{N}\right)$ is given
in green respectively red. The number of extreme points of
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ that are of
Monge-type (see (1.5)-(1.7)) is depicted in black. Here, as usual, $N$ denotes
the number of marginals and $\ell$ the number of states.

It was already mentioned above, that the extreme points of the polytope
$P_{\textrm{coef}}$ have a sparse structure. In more detail, a coefficient
vector $\alpha\in P_{\textrm{coef}}$ is extremal with respect to the polytope
$P_{\textrm{coef}}$ if and only if its nonzero entries correspond to a
selection of columns of $A$ which are linearly independent (see, e.g., [3]).
That is why, the complexity of computing the extreme points of
$P_{\textrm{coef}}$, and their number, increases faster with the number of
states than with the number of marginals. Suppose you are looking at a setting
where the number of marginals is equal to the number of states. Then, on the
one hand, increasing the number of marginals by $1$ yields
$\tbinom{2N}{N+1}-\tbinom{2N-1}{N}$ more columns in $A$. On the other hand, an
increase in the number of states by $1$ enlarges the number of columns of $A$
by $\tbinom{2N}{N}-\tbinom{2N-1}{N}$. Elementary computations show that in the
second case $A$ has $\frac{1}{N+1}\tbinom{2N}{N}$ more columns than in the
first case. Moreover, in contrast to an increase in the number of marginals,
an increase in the number of states also increases the number of rows of $A$
by 1. Therefore then up to $\ell+1$ columns of $A$ can be linearly
independent. Hence, an increase in the number of states leads to a faster
increasing (compared to an increase in the number of marginals) number of
subsets of linearly independent columns of the constraint matrix $A$ by
yielding a steeper increase in the number of columns of $A$ as well as by
enlarging the dimension of the column space. Each of these subsets corresponds
to an extreme point of $P_{\textrm{coef}}$.

######  Remark 2.4.

This remark lists interpretations and observations regarding Figure 2.

  1. 1)

In the case $N=2$, Figure 2 shows that in the considered cases every extremal
symmetric Kantorovich coupling is of Monge-type. In the given setting this
means that every extreme point of
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ is a
symmetrized permutation matrix. It is easy to see, using the celebrated
Birkhoff-von Neumann theorem [4, 36] as well as the linearity of the
symmetrization operator $S$ (2.4), that this holds true for an arbitrary
number $\ell$ of states. Note, however, that not every symmetrized permutation
matrix is an extreme point of
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$, but only
those symmetrized Monge-states whose corresponding coefficient vectors select
linearly independent columns of $A$.

  2. 2)

In the case $\ell=2$, the number of extremal symmetric Kantorovich couplings
which are of Monge-type increases by 1 each time the marginal number is even.
It is easy to prove that this pattern will continue. Firstly, note that, in
the case $\ell=2$, every symmetric Kantorovich coupling of Monge-type is an
extreme point of
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$. This
follows by a support-argument regarding the corresponding coefficient vectors.
Secondly, we take a look at the symmetrized Monge-states in this setting. We
assume the marginal vectors $\lambda^{(1)},\ldots,\lambda^{(N+1)}$ are sorted
in the columns of $A$ by the first component in decreasing order, i.e.,

$A=\begin{pmatrix}1&\frac{N-1}{N}&\dots&\frac{1}{N}&0\vspace{0.2cm}\\\
0&\frac{1}{N}&\dots&\frac{N-1}{N}&1\end{pmatrix}.$

Then the symmetric Kantorovich couplings of Monge-type are exactly those
couplings with coefficient vectors

$\alpha^{(j)}=\frac{1}{2}e^{(j)}+\frac{1}{2}e^{(N+1-j+1)}$

for $j=1,2,\ldots,\lceil\frac{N+1}{2}\rceil$, where $e^{(i)}$ is the $i$-th
unit vector.

  3. 3)

The setting of $3$ marginals and $3$ sites, i.e., $N=\ell=3$ is the main focus
in [16]. There interested readers can find the $22$ extreme points of the
symmetric Kantorovich polytope explicitly listed including the information
which extremal elements are of Monge-type and which are not. This list also
shows which pairs of permutations (identifiying $T_{1}$ with the identity)
correspond to an extremal symmetric Kantorovich coupling. [16] also visualizes
these $22$ extremal states as molecular packings, where one can identify
irreducible packings with extreme points.

##  3 Classification of the Extreme Points of a Reduced Kantorovich Polytope

In Section 2, we achieved a better understanding of the optimal transport
problem (1.1) by numerically analyzing the convex geometry of the set of
admissible trial states, i.e., the set of symmetric Kantorovich couplings
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$. Motivated
by physical applications we assume from this point on that the given cost
function has pairwise symmetric structure. Then the set of admissible trial
states can be reduced by the linear map $M_{2}$ to a lower-dimensional
polytope thereby decreasing the number of extremal states.

In more detail, we consider the optimal transport problem (1.1) with
$c:X^{N}\to\mathbb{R}$ being a cost function with pairwise symmetric
structure, i.e.,

$c(x_{1},\ldots,x_{N})=\sum_{1\leq i<j\leq N}v(x_{i},x_{j})\textrm{ for all
}(x_{1},\ldots,x_{N})\in X^{N},$ (3.1)

where $v:X^{2}\to\mathbb{R}$ is a symmetric pair-potential, i.e.,
$v(x,y)=v(y,x)$ for all $(x,y)\in X^{2}$. Then the objective function of (1.1)
can be rewritten as

$\int_{X^{N}}c(x_{1},\ldots,x_{N})d\gamma(x_{1},\ldots,x_{N})=\binom{N}{2}\int_{X^{2}}v(x,y)d\left(M_{2}\gamma\right)(x,y),$
(3.2)

where
$\gamma\in\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ is
an arbitrary symmetric Kantorovich coupling. This elementary reformulation was
established in [17]. There also the concept of $N$-representability (see
Definition 3.1) was introduced which we will use in the following to write the
reduced set of admissible trial states in a compact manner.

######  Definition 3.1 ($N$-representability).

A probability measure $\mu\in\mathcal{P}\left(X^{k}\right)$ is called
$N$-representable if there exists a symmetric probability measure $\gamma$ on
the product space $X^{N}$, i.e.,
$\gamma\in\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$, such that $\mu$ is
its $k$-point marginal, i.e.,

$\mu=M_{k}\gamma.$ (3.3)

Any such symmetric probability measure on $X^{N}$ that fulfills (3.3) is then
called a representing measure of $\mu$. In the following the set of
$N$-representable $k$-plans will be denoted by
$\mathcal{P}_{N\textrm{-rep}}\left(X^{k}\right)$.

As we consider pairwise interactions, we will focus our attention on the set
of $N$-representable $2$-point measures, i.e.,
$\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)$. Note, however, that cost
functions $c$ embodying $k$-particle interactions would give rise to a problem
reformulation reducing the set of admissible trial states to a subset of
$\mathcal{P}_{N\textrm{-rep}}\left(X^{k}\right)$. In the case $k=N$, $c$ would
be a symmetric cost which are, as mentioned in the introduction, ’dual’ to the
set of symmetric probability measures on the product space $X^{N}$. In the
same manner, cost functions with symmetric pairwise structure have a dual
relationship with the set of $N$-representable $2$-plans.

By definition of $N$-representability (see Definition 3.1), the set of
$N$-representable $2$-point measures is the image of the set of symmetric
probability measures on $X^{N}$ under the map $M_{2}$, defined in (2.7), i.e.,
$M_{2}\left(\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)\right)=\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)$.
Combining this equality with (3.2) yields that (1.8) is an equivalent
reformulation of the multi-marginal optimal transport problem (1.1) for a cost
function with pairwise symmetric structure (3.1). Here (1.8) can also be
written as

$\min_{\mu\in\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)}\int_{X^{2}}v(x,y)d\mu(x,y),$

where $\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ is
the set of $N$-representable $2$-plans having uniform marginal, i.e.,

$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right):=\left\\{\mu\in\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right):M_{1}(\mu)=\overline{\lambda}\right\\}.$
(3.4)

We will refer to the set
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ as reduced
Kantorovich polytope for $N$ marginals and $\ell$ states. The convex geometry
of this set will be numerically analyzed in the following. Thereby the
validity of Monge’s approach in the given setting, i.e., symmetric multi-
marginal optimal transport on finite state spaces with pairwise symmetric cost
functions, will be tested.

We have seen above that under the assumption of pairwise symmetric cost
functions the optimal transport problem (1.1), where the set of admissible
trial states is given by the high-dimensional set
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$, can be
equivalently formulated as a minimization problem over the lower-dimensional
set $\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ (see
(1.8)). The pairwise symmetric structure implies that any symmetric
Kantorovich coupling influences the value of the objective function of problem
(1.1) only through their respective two-point marginal (see (3.2)). The nature
of this reformulation, applying the two-point marginal map $M_{2}$ on the set
of symmetric Kantorovich couplings, however, entails that the new set of
admissible trial states, i.e., the reduced Kantorovich polytope is only
implicitly known. Only in the two-marginal (N=2) case the reduced Kantorovich
polytope can be understood in a straightforward manner: It corresponds to the
set of symmetric bistochastic matrices scaled by the factor $\frac{1}{\ell}$
(see Remark 3.5 1) below for further consideration of the two-marginal case).
Hence, in the case $N=2$,
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)=\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$
holds. For a better understanding of the multi-marginal ($N>2$) case, we will
in the following, as motivated, view the reduced Kantorovich polytope as the
image of the set of symmetric Kantorovich couplings on $X^{N}$ under the two-
point marginal map, i.e.,

$M_{2}(\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right))=\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right).$
(3.5)

As described in Section 2,
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ corresponds
to the convex hull of its extreme points. Combining this fact with (3.5) and
the linearity of $M_{2}$ yields that the reduced Kantorovich polytope
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ is equal
to the convex hull of the two-point marginals of extremal symmetric
Kantorovich couplings, i.e.,

$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)=\textrm{conv}\left(\left\\{M_{2}\gamma:\gamma\textrm{
is an extreme point of
}\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)\right\\}\right).$
(3.6)

The following proposition is an immediate consequence.

######  Proposition 3.2.

Any extreme point of the reduced Kantorovich polytope for $N$ marginals and
$\ell$ states is the two-point marginal of an extremal symmetric Kantorovich
coupling.

Now the question is whether or not $M_{2}$ represents a bijective relationship
between the sets of extreme points of
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ and
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$. The
following remark sheds light on this issue applying the bijective relationship
between $\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ and
the polytope $P_{\textrm{coef}}$ established in Lemma 2.2 and Corollary 2.3.

######  Remark 3.3.

In Section 2 the extreme points of
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$, i.e., the
set of admissible trial states of problem (1.1), are determined using the
set’s bijective relationship, captured in the coefficients-to-coupling map $R$
introduced in Section 2, to the polytope $P_{\textrm{coef}}$. As explained
above in more detail, the map $R$ identifies any symmetric probability measure
on $X^{N}$ $\gamma$ with a coefficient vector $\alpha$, such that $\gamma$ can
be written as the respective convex combination of the extreme points of
$\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$, i.e., (2.6) holds. These
coefficients are unique due to the disjoint support of the extremal symmetric
probability measures on $X^{N}$. It was proven in [18] that the two-point
marginal map $M_{2}$ is a bijection between the sets of extreme points of
$\mathcal{P}_{\textrm{sym}}\left(X^{N}\right)$ and
$\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)$, respectively. Due to the
linearity of $M_{2}$, given a coefficient vector $\alpha$ and a symmetric
probability measure $\gamma$ on $X^{N}$, such that $\gamma=R\alpha$, i.e.,
(2.6) holds true, then

$M_{2}\gamma=\sum_{1\leq i_{1}\leq\dots\leq
i_{N}\leq\ell}\alpha_{i_{1},\ldots,i_{N}}M_{2}S\delta_{i_{1},\ldots,i_{N}}.$

Only now, these coefficients $\alpha$ representing $M_{2}\gamma$ as a convex
combination of the extreme points of the set of $N$-representable two-point
measures may not be unique, rendering us unable to identify the extreme points
of the reduced Kantorovich polytope with those of the coefficient-polytope
$P_{\textrm{coef}}$.

The remark above illuminates why the extreme points of the set of symmetric
Kantorovich couplings
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ can not be
identified with the extremal elements of the reduced Kantorovich polytope
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ via
$M_{2}$. The two-point marginal map may for example map multiple extreme
points of the set
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$ on a single
point of $\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$;
this point may lie on a face or in the interior of
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ (see [16]
for an well-illustrated example).

Nevertheless, it was established in Proposition 3.2 that every extremal
element of $\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$
has a representing measure that is itself extremal with respect to
$\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$. The extreme
points of this set of symmetric Kantorovich couplings were in Corollary 2.3
identified with the extreme points of $P_{\textrm{coef}}$. In combination with
the in Remark 3.3 established connection between $P_{\textrm{coef}}$ and
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ this leads
us to the following approach to determine the extremal elements of
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$:

  1. 1.

We start by determining the extremal elements of $P_{\textrm{coef}}$. This was
already done within the considerations of Section 2.

  2. 2.

Every such extreme point is multiplied by the matrix
$T\in\mathbb{R}^{\ell^{2}\times|E_{\textrm{sym}}^{N}|}$ which is constructed
as follows. The matrix $A$ as defined in (2.10) lists all the elements of
$\mathcal{P}_{\frac{1}{N}}\left(X\right)$ as columns. It was proven in [18]
that for any element $\lambda$ of $\mathcal{P}_{\frac{1}{N}}\left(X\right)$
the following holds:

$M_{2}\psi_{N}\left(\lambda\right)=\frac{N}{N-1}\lambda\otimes\lambda-\frac{1}{N-1}\left(\textrm{id,id}\right)\\#\lambda,$
(3.7)

where the map $\psi_{N}$ was introduced in Section 2. Note that it was further
established in [18] that measures of form (3.7) for
$\lambda\in\mathcal{P}_{\frac{1}{N}}\left(X\right)$ are exactly the extreme
points of $\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)$. Now we construct
$T$ by replacing any column $\lambda$ of $A$ with
$M_{2}\psi_{N}\left(\lambda\right)$ as given in (3.7) where we canonically
identify matrices with vectors by gluing columns together.

  3. 3.

Finally we check which points of the form

$T\alpha\textrm{ for }\alpha\in\textrm{ext}\left(P_{\textrm{coef}}\right)$

are extremal with respect to
$\textrm{conv}\left(\left\\{T\alpha:\alpha\in\textrm{ext}\left(P_{\textrm{coef}}\right)\right\\}\right)$
and therefore by (3.6) with respect to
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$.

Note that it is computationally more complex to determine the extremal
elements of
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ than those
of $\mathcal{P}_{\textrm{sym},\overline{\lambda}}\left(X^{N}\right)$.

Now, we will incorporate Monge’s approach in the reduced setting.

######  Definition 3.4.

An element of the reduced Kantorovich polytope for $N$ marginals and $\ell$
states is said so be of Monge-type or in Monge-form if it has a representing
measure that is of Monge-form (see (1.5)-(1.7)).

This definition is consistent with our goal to check the validity of Monge’s
approach as any optimizer in Monge-form for problem (1.8) guarantees the
existence of an optimizer in Monge-form for problem (1.1). The set of all
elements of
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ which are
in Monge-form will be denoted as
$\mathcal{P}_{N\textrm{-rep,Monge}}\left(X^{2}\right)$, i.e.,

$\mathcal{P}_{N\textrm{-rep,Monge}}\left(X^{2}\right):=\left\\{M_{2}\gamma:\gamma\textrm{
is of Monge-type }\eqref{eq:MongeOne}-\eqref{eq:MongeThree}\right\\}.$

Analogously to (2.13) we introduce the reduced Monge polytope for $N$
marginals and $\ell$ states
$\mathcal{P}_{N\textrm{-rep,Monge}}^{\textrm{conv}}\left(X^{2}\right)$ as
follows.

$\mathcal{P}_{N\textrm{-rep,Monge}}^{\textrm{conv}}\left(X^{2}\right):=\textrm{conv}\left(\mathcal{P}_{N\textrm{-rep,Monge}}\left(X^{2}\right)\right)$
(3.8)

The extremal elements of the reduced Monge polytope can be determined in the
same manner as those of the reduced Kantorovich polytope (see the description
of the procedure above). Starting point are now the extremal elements of the
Monge polytope $\mathcal{P}_{\textrm{sym,Monge}}\left(X^{N}\right)$
interpreted as coefficient vectors.

Checking which of the extreme points of the reduced Kantorovich polytope
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ correspond
to an extremal element of the reduced Monge polytope
$\mathcal{P}_{N\textrm{-rep,Monge}}^{\textrm{conv}}\left(X^{2}\right)$ tells
us which of the extreme points of
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ are of
Monge-type.

The data in Figure 3 was computed using MATLAB [24] and polymake [31].

Figure 3: The number of extreme points of the reduced Kantorovich polytope
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$
respectively of the reduced Monge polytope
$\mathcal{P}_{N\textrm{-rep,Monge}}^{\textrm{conv}}\left(X^{2}\right)$ is
given in green respectively red. The number of extreme points of
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ that are
of Monge-type (see Definition 3.4) is depicted in black. Here, as usual, $N$
denotes the number of marginals and $\ell$ the number of states.

######  Remark 3.5.

What follows are interpretations and observations regarding Figure 3.

  1. 1)

Combining the convention $M_{N}=\textrm{id}$ with (3.5), it is obvious that
the symmetric Kantorovich polytope for $2$ marginals and $\ell$ states
$\mathcal{P}_{\textrm{sym,}\overline{\lambda}}\left(X^{2}\right)$ coincides
with the reduced Kantorovich polytope for $2$ marginals and $\ell$ states
$\mathcal{P}_{2\textrm{-rep,}\overline{\lambda}}\left(X^{2}\right)$. This fact
was already mentioned above. It was established in Remark 2.4 that every
extreme point of
$\mathcal{P}_{\textrm{sym,}\overline{\lambda}}\left(X^{2}\right)$ is a
symmetrized permutation matrix, i.e., the image of a permutation matrix under
the symmetrization operator (2.4). In the setting of $2$ marginals,
symmetrized permutation matrices exactly correspond to symmetrized Monge
states. See Remark 2.4 for further considerations of the case $N=2$.

  2. 2)

In the case $\ell=2$, Figure 3 depicts that in the considered cases, the
reduced Kantorovich polytope
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ has two
extreme points both of which are in Monge-form. Hence, in any considered case
the line segment
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ coincides
with the respective reduced Monge polytope
$\mathcal{P}_{N\textrm{-rep,Monge}}^{\textrm{conv}}\left(X^{2}\right)$.
One can prove by elementary arguments that this holds true for an arbitrary
number of marginals $N$ in the case of $2$ sites. In a little more detail,
considering the dimension of
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ in the
given case and parametrising the elements of
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ by their
off-diagonal element allows us to deduce that the two extreme points of
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ are given
by

$\mu^{(1)}=M_{2}\left(\frac{1}{2}\psi_{N}\left(\delta_{1}\right)+\frac{1}{2}\psi_{N}\left(\delta_{2}\right)\right)$
(3.9)
$\mu^{(2)}=\begin{cases}M_{2}\left(\psi_{N}\left(\frac{1}{2}\delta_{1}+\frac{1}{2}\delta_{2}\right)\right)&\textrm{
if }N\textrm{ is even}\\\
M_{2}\left(\frac{1}{2}\psi_{N}\left(\frac{N-1}{2N}\delta_{1}+\frac{N+1}{2N}\delta_{2}\right)+\frac{1}{2}\psi_{N}\left(\frac{N+1}{2N}\delta_{1}+\frac{N-1}{2N}\delta_{2}\right)\right)&\textrm{
if }N\textrm{ is odd}\end{cases}$ (3.10)

or in pedestrian notation,

$\mu^{(1)}=\begin{pmatrix}\frac{1}{2}&0\\\ 0&\frac{1}{2}\end{pmatrix}$
$\mu^{(2)}=\begin{cases}\frac{1}{4(N-1)}\begin{pmatrix}N-2&N\\\
N&N-2\end{pmatrix}&\textrm{ if }N\textrm{ is even}\\\
\frac{1}{4N}\begin{pmatrix}N-1&N+1\\\ N+1&N-1\end{pmatrix}&\textrm{ if
}N\textrm{ is odd}.\end{cases}$

As the coefficients in (3.9) and (3.10) are integer multiples of
$\frac{1}{\ell}$ both extreme points are $2$-point marginals of symmetric
Kantorovich couplings in Monge-form and therefore they are themselves elements
of the reduced Kantorovich polytope
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ which are
of Monge-type (see Definition 3.4). Note that $\mu^{(2)}$, which is of Monge-
type and therefore describes a correlated or in other words deterministic
state, converges for $N\to\infty$ to the independent measure
$\lambda\otimes\lambda$ for
$\lambda=\left(\frac{1}{2}\delta_{1}+\frac{1}{2}\delta_{2}\right)$.
These findings coincide with the results in [17], where a model problem of $N$
particles on $2$ sites was considered. There also the set of $N$-representable
$2$-plans $\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)$ for $X$ consisting
of $2$ distinct elements is illustrated. Imposing the here given marginal
condition on these sets leads to the respective line segment
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$.

  3. 3)

The case of $3$ marginals and $3$ sites, i.e., $N=\ell=3$, is a minimal
example of a point in the grid, both with respect to the sum of both
parameters $N+\ell$ and with respect to the minimum of both parameters
$\min\\{N,\ell\\}$, such that not every extremal element of the reduced
Kantorovich polytope
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ is of
Monge-type. In the considered case
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ has $8$
extreme points $5$ of which are in Monge-form. By extension $3$ of them are
not. They are given by

$\frac{1}{2}M_{2}\left(S\delta_{112}\right)+\frac{1}{2}M_{2}\left(S\delta_{233}\right)$
(3.11)

and the two states one generates by imposing the role of the second site on
the first and third site respectively. (3.11) is the unique optimizer of an
optimal transport problem stated in [16]. This problem corresponds to a
molecular packing problem. See [16] for further reading.

##  4 A Model Problem: Optimal Couplings of $N$ marginals on $3$ sites for
pairwise costs

In the following, we focus our attention on symmetric multi-marginal optimal
transport problems (1.1) on $3$ sites, i.e., $X=\\{a_{1},a_{2},a_{3}\\}$. As
in Section 3, we only consider pairwise symmetric costs and therefore are able
to reformulate (1.1) as the lower-dimensional problem (1.8). In particular,
the reduced Kantorovich polytope for $N$ marginals and $3$ states corresponds
to the respective set of admissible trial states. It is easy to see that in
the given setting these polytopes are three-dimensional. As by extension the
reduced Monge polytope is at most three-dimensional, we are able to visually
compare both approaches.

The following visualizations in Figure 4 were generated by extending the above
explained calculations and routines in MATLAB [24] and polymake [31].

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

Figure 4: The reduced Kantorovich respectively Monge polytope for $N$
marginals and $3$ states is visualized for $N=2,3,4,6$ and $10$ in green
respectively red. The plots are arranged in ascending order with respect to
$N$. The elements $\left(\mu_{ij}\right)_{i,j=1}^{3}$ of the polytopes are
parametrized by their off-diagonal entries $\mu_{12},\mu_{13}$ and $\mu_{23}$.
In the case of three marginals the reduced polytopes were initially depicted
in [16]. $V$ indicates the volume of the corresponding polytope. The
volumetric ratio, reduced Monge polytope to reduced Kantorovich polytope, is
depicted in Figure 4.

Note that by definition the reduced Monge polytope is always contained in the
reduced Kantorovich polytope independently of the number of marginals $N$ and
the number of sites $\ell$. Some of the extreme points of the reduced Monge
polytope are also extremal with respect to the reduced Kantorovich polytope
and some lie on faces or in the interior of the latter (see Figure 3 for
specific numbers).

In the setting of $\ell=3$ sites, the numerical analysis of the reduced
setting discussed in Section 3 yields that in the case of $N=2,3,\ldots,9$ and
$10$ marginals there are always $5$ prominent extreme points of the reduced
Kantorovich polytope that are of Monge-type. This is also indicated by Figure
1 as well as Figure 4. In the illustrations they can be identified with the
four extreme points on the co-ordinate axes, including the origin, as well as
the ’peak’ in the front of the polytopes. In formulas these extreme points can
be written as depicted in Table 1.

Nomenclature | Abstract Notation | Matrix Notation
---|---|---
|  | $\frac{1}{3}M_{2}\psi_{N}(\delta_{1})$ | $\begin{pmatrix}\frac{1}{3}&0&0\\\ 0&\frac{1}{3}&0\\\ 0&0&\frac{1}{3}\\\ \end{pmatrix}$
$EA^{(N)}$ | $EA^{(N)}$ | $+\frac{1}{3}M_{2}\psi_{N}(\delta_{2})$
|  | $+\frac{1}{3}M_{2}\psi_{N}(\delta_{3})$
| $ER^{(3m)}$ | $M_{2}\psi_{N}\left(\frac{1}{3}\delta_{1}+\frac{1}{3}\delta_{2}+\frac{1}{3}\delta_{3}\right)$ | $\begin{pmatrix}\frac{N-3}{9(N-1)}&\frac{N}{9(N-1)}&\frac{N}{9(N-1)}\\\ \frac{N}{9(N-1)}&\frac{N-3}{9(N-1)}&\frac{N}{9(N-1)}\\\ \frac{N}{9(N-1)}&\frac{N}{9(N-1)}&\frac{N-3}{9(N-1)}\\\ \end{pmatrix}$
|  | $\frac{1}{3}M_{2}\psi_{N}\left(\frac{m}{N}\delta_{1}+\frac{m}{N}\delta_{2}+\frac{m+1}{N}\delta_{3}\right)$ | $\begin{pmatrix}\frac{N-2}{9N}&\frac{N+1}{9N}&\frac{N+1}{9N}\\\ \frac{N+1}{9N}&\frac{N-2}{9N}&\frac{N+1}{9N}\\\ \frac{N+1}{9N}&\frac{N+1}{9N}&\frac{N-2}{9N}\\\ \end{pmatrix}$
$ER^{(N)}$ | $ER^{(3m+1)}$ | $+\frac{1}{3}M_{2}\psi_{N}\left(\frac{m}{N}\delta_{1}+\frac{m+1}{N}\delta_{2}+\frac{m}{N}\delta_{3}\right)$
|  | $+\frac{1}{3}M_{2}\psi_{N}\left(\frac{m+1}{N}\delta_{1}+\frac{m}{N}\delta_{2}+\frac{m}{N}\delta_{3}\right)$ | $\begin{pmatrix}\frac{N-2}{9N}&\frac{N+1}{9N}&\frac{N+1}{9N}\\\ \frac{N+1}{9N}&\frac{N-2}{9N}&\frac{N+1}{9N}\\\ \frac{N+1}{9N}&\frac{N+1}{9N}&\frac{N-2}{9N}\\\ \end{pmatrix}$
|  | $\frac{1}{3}M_{2}\psi_{N}\left(\frac{m}{N}\delta_{1}+\frac{m+1}{N}\delta_{2}+\frac{m+1}{N}\delta_{3}\right)$
| $ER^{(3m+2)}$ | $+\frac{1}{3}M_{2}\psi_{N}\left(\frac{m+1}{N}\delta_{1}+\frac{m}{N}\delta_{2}+\frac{m+1}{N}\delta_{3}\right)$
|  | $+\frac{1}{3}M_{2}\psi_{N}\left(\frac{m+1}{N}\delta_{1}+\frac{m+1}{N}\delta_{2}+\frac{m}{N}\delta_{3}\right)$
| $E12^{(2m)}$ | $\frac{2}{3}M_{2}\psi_{N}\left(\frac{1}{2}\delta_{1}+\frac{1}{2}\delta_{2}\right)+\frac{1}{3}M_{2}\psi_{N}\delta_{3}$ | $\begin{pmatrix}\frac{N-2}{6(N-1)}&\frac{N}{6(N-1)}&0\\\ \frac{N}{6(N-1)}&\frac{N-2}{6(N-1)}&0\\\ 0&0&\frac{1}{3}\\\ \end{pmatrix}$
$E12^{(N)}$
| $\frac{1}{3}M_{2}\psi_{N}\left(\frac{N-1}{2N}\delta_{1}+\frac{N+1}{2N}\delta_{2}\right)$ | $\begin{pmatrix}\frac{N-1}{6N}&\frac{N+1}{6N}&0\\\ \frac{N+1}{6N}&\frac{N-1}{6N}&0\\\ 0&0&\frac{1}{3}\\\ \end{pmatrix}$
| $E12^{(2m+1)}$ | $+\frac{1}{3}M_{2}\psi_{N}\left(\frac{N+1}{2N}\delta_{1}+\frac{N-1}{2N}\delta_{2}\right)$
|  | $+\frac{1}{3}M_{2}\psi_{N}(\delta_{3})$
Table 1: The extreme points $EA^{(N)}$, $ER^{(N)}$ and $E12^{(N)}$ of the
reduced Kantorovich polytope are depicted in abstract and matrix notation.
Hereby $N$ corresponds to the number of marginals and $m\in\mathbb{N}_{0}$ is
a non-negative integer, allowing us to distinguish between the various cases
regarding $N$. As all the coefficients in the ’Abstract Notation’-column are
integer multiples of $\frac{1}{3}$, these extreme points are of Monge-type.

In Figure 4, $EA^{(N)}$ corresponds to the origin, $ER^{(N)}$ to the ’peak’
and $E12^{(N)}$ to the non-origin extreme point on the $\mu_{12}$-axis.
$E12^{(N)}$ assumes an exemplary role in Table 1. The corresponding extreme
points on the $\mu_{13}$\- respectively $\mu_{23}$-axis can be expressed
analogously in abstract as well as matrix notation and will be denoted by
$E13^{(N)}$ respectively $E23^{(N)}$.

So far we know by numerical analysis that $EA^{(N)}$, $ER^{(N)}$, $E12^{(N)}$,
$E13^{(N)}$ and $E23^{(N)}$ are extreme points of the reduced Kantorovich
polytope $\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$
in the cases of $N=2,3,\ldots,9$ and $10$ marginals. One can prove that this
holds true for a general number $N\geq 2$ of marginals. For $E12^{(N)}$,
$E13^{(N)}$ and $E23^{(N)}$ one can show this by following the same approach
taken in Remark 3.5 2). In case of $EA^{(N)}$ and $ER^{(N)}$ it is an
immediate consequence of Theorem 4.1. In the following $d:X\times
X\to\mathbb{R}$ will denote the discrete metric defined by

$d(x,y):=\begin{cases}1&\textrm{ if }x\neq y\\\ 0&\textrm{ if
}x=y.\end{cases}$

######  Theorem 4.1.

We consider the reduced multi-marginal optimal transport problem (1.8) for
$N\geq 2$ marginals and $\ell=3$ sites.

  1. a)

For the attractive cost function $d:X\times X\to\mathbb{R}$ the unique
minimizer is given by $EA^{(N)}$.

  2. b)

For the repulsive cost function $c_{R}:X\times X\to\mathbb{R}$ given by

$c_{R}(x,y):=\begin{cases}\frac{1}{d(x,y)}&\textrm{ if }x\neq y\\\ B&\textrm{
if }x=y\end{cases}$ (4.1)

for some constant $B>1$, the unique minimizer is given by $ER^{(N)}$.

###### Proof.

In the following the elements of the reduced Kantorovich polytope will always
be interpreted as matrices. Along those lines $D$ respectively $CR$
corresponds to the matrix notation of $d$ respectively $c_{R}$, i.e.,
$D_{ij}:=d(a_{i},a_{j})$ respectively $CR_{ij}:=c_{R}(a_{i},a_{j})$ for
$i,j\in\\{1,2,\ldots,\ell\\}$, and $\langle\cdot,\cdot\rangle$ denotes the
standard matrix scalar product.

  1. a)

Note that by non-negativity of the cost function $d$ the objective value of an
arbitrary admissible state
$\mu\in\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ is
non-negative, i.e., $\langle D,\mu\rangle\geq 0$. As $EA^{(N)}$ is admissible
and yields an objective value of $0$, i.e., $\langle D,EA^{(N)}\rangle=0$, it
is an optimizer of the corresponding problem (1.8). Positivity of $D$ in its
off-diagonal entries and the marginal constraint ensure that $EA^{(N)}$ is the
unique optimizer.

  2. b)

To prove the second assertion, we drop the marginal constraint in a
reformulated version of the considered problem (1.8) and calculate the
extremal elements of $\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)$, which
solve the new optimization problem. There will be a unique convex combination
of the optimal extreme points of
$\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)$, namely $ER^{(N)}$, that lies
in $\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$. This
state then corresponds to the unique minimizer of problem (1.8).

We consider the problem

$\min_{\mu\in\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)}\langle
CR,\mu\rangle.$ (4.2)

Subsequently changing the objective function to
$\langle{CR}-{\mathbb{1}},\cdot\rangle$, where all the entries of
${\mathbb{1}}\in\mathbb{R}^{3\times 3}$ are given by $1$, and plugging in the
marginal constraint allows us to reformulate (4.2) as follows

$\max_{\mu\in\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)}\mu_{12}+\mu_{13}+\mu_{23}.$
(4.3)

By dropping the marginal constraint, further restricting the admissible set to
the extreme points (3.7) of $\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)$
and rescaling the new admissible set by $N^{2}$ leads to the new optimization
problem

$\max_{\lambda\in
N\mathcal{P}_{\frac{1}{N}}\left(X\right)}\lambda_{1}\lambda_{2}+\lambda_{1}\lambda_{3}+\lambda_{2}\lambda_{3}=\lambda_{1}\lambda_{2}+(N-\lambda_{3})\lambda_{3}.$
(4.4)

Assuming $\left(\lambda_{1}^{*},\lambda_{2}^{*},\lambda_{3}^{*}\right)$ is an
optimizer of problem (4.4), then elementary calculations show that
$\left(\lambda_{1}^{*},\lambda_{2}^{*}\right)$ fulfills

$\left(\lambda_{1}^{*},\lambda_{2}^{*}\right)\in\begin{cases}\left\\{\left(\frac{r}{2},\frac{r}{2}\right)\right\\}&\textrm{
if }r\textrm{ is even}\\\
\left\\{\left(\frac{r-1}{2},\frac{r+1}{2}\right),\left(\frac{r+1}{2},\frac{r-1}{2}\right)\right\\}&\textrm{
if }r\textrm{ is odd,}\end{cases}$ (4.5)

where $r:=N-\lambda_{3}^{*}$. Otherwise
$\left(\lambda_{1}^{*},\lambda_{2}^{*},\lambda_{3}^{*}\right)$ would not be
optimal. Note that (4.4) always admits a maximizer as
$\mathcal{P}_{\frac{1}{N}}\left(X\right)$ is finite.
This allows us to identify problem (4.4) with the one-parameter optimization
problem

$\max\left\\{\max_{r\in 2\mathbb{N}_{0},r\leq
N}-\frac{3}{4}r^{2}+Nr,\max_{r\in 2\mathbb{N}_{0}+1,r\leq
N}-\frac{3}{4}r^{2}+Nr-\frac{1}{4}\right\\}.$ (4.6)

Elementary calculations reveal that $r\in\\{0,1,\ldots,N\\}$ is optimal
regarding (4.6) if and only if

$r\in\begin{cases}\\{2m\\}&\textrm{ if }N=3m\textrm{ for
}m\in\mathbb{N}_{0}\\\ \\{2m,2m+1\\}&\textrm{ if }N=3m+1\textrm{ for
}m\in\mathbb{N}_{0}\\\ \\{2m+1,2m+2\\}&\textrm{ if }N=3m+2\textrm{ for
}m\in\mathbb{N}_{0}.\end{cases}$

It immediately follows that $\lambda\in
N\mathcal{P}_{\frac{1}{N}}\left(X\right)$ is optimal with respect to problem
(4.4) if and only if

$\lambda\in\begin{cases}N\left\\{\left(\frac{m}{N},\frac{m}{N},\frac{m}{N}\right)\right\\}&\textrm{if
}N=3m\textrm{ for }m\in\mathbb{N}_{0}\\\
N\left\\{\left(\frac{m}{N},\frac{m}{N},\frac{m+1}{N}\right),\left(\frac{m}{N},\frac{m+1}{N},\frac{m}{N}\right),\left(\frac{m+1}{N},\frac{m}{N},\frac{m}{N}\right)\right\\}&\textrm{if
}N=3m+1\textrm{ for }m\in\mathbb{N}_{0}\\\
N\left\\{\left(\frac{m}{N},\frac{m+1}{N},\frac{m+1}{N}\right),\left(\frac{m+1}{N},\frac{m}{N},\frac{m+1}{N}\right),\left(\frac{m+1}{N},\frac{m+1}{N},\frac{m}{N}\right)\right\\}&\textrm{if
}N=3m+2\textrm{ for }m\in\mathbb{N}_{0}.\end{cases}$ (4.7)

Recall that (3.7) allows us (after dropping the factor $N$ in (4.7)) to
identify the given maximizers with exactly those extremal elements of
$\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)$ that maximize the sum of
their off-diagonal entries. As by Minkowski’s theorem any element of
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ can be
written as convex combination of the extreme points of
$\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)$ and in any of the considered
cases in (4.7) there are unique coefficients given by

$\displaystyle\alpha_{1}=1$ $\displaystyle\hskip 28.45274pt\textrm{if
}N=3m\textrm{ for }m\in\mathbb{N}_{0}$
$\displaystyle\alpha_{1}=\alpha_{2}=\alpha_{3}=\frac{1}{3}$ else

that allow us to write $\overline{\lambda}$ as a convex combination of the
respective optimizers in (4.7), it is easy to see that $ER^{(N)}$ as defined
in Table 1 is the unique optimizer of (4.3) and thereby (4.2) for any $N\geq
2$.

∎

######  Remark 4.2.

It is an immediate consequence of Theorem 4.1 that

$\gamma_{GS}:=\frac{1}{3}S\delta_{11\ldots 1}+\frac{1}{3}S\delta_{22\ldots
2}+\frac{1}{3}S\delta_{33\ldots 3}$

respectively

$\gamma_{C}:=\frac{1}{3}S\delta_{1\tau(1)\ldots\tau^{(N-1)}(1)}+\frac{1}{3}S\delta_{2\tau(2)\ldots\tau^{(N-1)}(2)}+\frac{1}{3}S\delta_{3\tau(3)\ldots\tau^{(N-1)}(3)},$

where $\tau:\\{1,2,3\\}\to\\{1,2,3\\}$ is the cyclic permutation defined by
$\tau(1)=2$, $\tau(2)=3$, $\tau(3)=1$ and $\tau^{(i)}$ denotes the $i$-th
composition of $\tau$ with itself, is a solution to the optimal transport
problem (1.1) for the Gangbo-Święch cost function $c_{GS}:X^{N}\to\mathbb{R}$
defined by

$c_{GS}\left(x_{1},\ldots,x_{N}\right):=\sum_{1\leq i<j\leq N}d(x_{i},x_{j})$
(4.8)

respectively the Coulomb cost function $c_{C}:X^{N}\to\mathbb{R}$ defined by

$c_{C}\left(x_{1},\ldots,x_{N}\right):=\sum_{1\leq i<j\leq
N}c_{R}(x_{i},x_{j}).$

Here (4.8) is a discretization of the pair-cost considered in [20]. Note that
one could replace $d(\cdot,\cdot)$ in (4.8) with $d(\cdot,\cdot)^{p}$ for any
$p>1$ without changing $c_{GS}$.

Next, we examine the behavior of the sequences $\left(EA^{(N)}\right)_{N\geq
2}$, $\left(ER^{(N)}\right)_{N\geq 2}$, $\left(E12^{(N)}\right)_{N\geq 2}$,
$\left(E13^{(N)}\right)_{N\geq 2}$ and $\left(E23^{(N)}\right)_{N\geq 2}$ for
$N$ tending to $\infty$. Taking a look at the right column in Table 1, it is
easy to see that the following holds true.

$\displaystyle
EA^{(N)}\xrightarrow[]{N\to\infty}\begin{pmatrix}\frac{1}{3}&0&0\\\
0&\frac{1}{3}&0\\\ 0&0&\frac{1}{3}\end{pmatrix}$
$\displaystyle=:EA^{(\infty)},ER^{(N)}\xrightarrow[]{N\to\infty}\begin{pmatrix}\frac{1}{9}&\frac{1}{9}&\frac{1}{9}\\\
\frac{1}{9}&\frac{1}{9}&\frac{1}{9}\\\
\frac{1}{9}&\frac{1}{9}&\frac{1}{9}\end{pmatrix}=:ER^{(\infty)},$
$\displaystyle E12^{(N)}$
$\displaystyle\xrightarrow[]{N\to\infty}\begin{pmatrix}\frac{1}{6}&\frac{1}{6}&0\\\
\frac{1}{6}&\frac{1}{6}&0\\\ 0&0&\frac{1}{3}\end{pmatrix}=:E12^{(\infty)}$

Here $E12^{(\infty)}$ assumes again an exemplary role and $E13^{(\infty)}$ as
well as $E23^{(\infty)}$ are defined in an analogous manner. One can express
these ’limit extreme points’ in a more probabilistic manner
$EA^{(\infty)}=\frac{1}{3}\delta_{11}+\frac{1}{3}\delta_{22}+\frac{1}{3}\delta_{33}$,
$ER^{(\infty)}=\left(\frac{1}{3}\delta_{1}+\frac{1}{3}\delta_{2}+\frac{1}{3}\delta_{3}\right)\otimes\left(\frac{1}{3}\delta_{1}+\frac{1}{3}\delta_{2}+\frac{1}{3}\delta_{3}\right)$
as well as
$E12^{(\infty)}=\frac{2}{3}\left(\frac{1}{2}\delta_{1}+\frac{1}{2}\delta_{2}\right)\otimes\left(\frac{1}{2}\delta_{1}+\frac{1}{2}\delta_{2}\right)+\frac{1}{3}\delta_{3}$
corresponding to the ’Abstract Notation’-column in Table 1 via (3.7).

In the following, $\mathcal{D}_{\infty\textrm{-rep},\overline{\lambda}}$ will
denote the convex hull of these ’limit extreme points’, i.e.,

$\mathcal{D}_{\infty\textrm{-rep},\overline{\lambda}}=\textrm{conv}\left(\left\\{EA^{(\infty)},ER^{(\infty)},E12^{(\infty)},E13^{(\infty)},E23^{(\infty)}\right\\}\right).$
(4.9)

For an illustration of $\mathcal{D}_{\infty\textrm{-rep},\overline{\lambda}}$
see Figure 4.

figureThe diamond-shaped polytope
$\mathcal{D}_{\infty\textrm{-rep},\overline{\lambda}}$, as defined in (4.9),
is depicted in blue. The elements $(\mu_{ij})^{3}_{i,j=1}$ of the polytope are
parametrized by their off-diagonal entries $\mu_{12},\mu_{13}$ and $\mu_{23}$.
The volume of $\mathcal{D}_{\infty\textrm{-rep},\overline{\lambda}}$ is
indicated in the upper-right corner.

figureThe volumetric ratio, reduced Monge polytope to reduced Kantorovich
polytope for $N$ marginals and $3$ states, is depicted in dependency of the
number of marginals $N$.

It was proven in [17] that $N$-representability becomes an increasingly
stringent condition as $N$ grows, in more detail,
$\mathcal{P}_{N\textrm{-rep}}\left(X^{2}\right)\subseteq\mathcal{P}_{\hat{N}\textrm{-rep}}\left(X^{2}\right)$
for any $N\geq\hat{N}\geq 2$. It follows immediately that the reduced
Kantorovich polytope for $N$ marginals and $3$ states
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ is
contained in the reduced Kantorovich polytope
$\mathcal{P}_{\hat{N}\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ for
$\hat{N}$ marginals and $3$ states. As
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ is closed
and convex, $\mathcal{D}_{\infty\textrm{-rep},\overline{\lambda}}$ is a subset
of the reduced Kantorovich polytope
$\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)$ for any
number $N\geq 2$ of marginals and $3$ sites. In summary we get the following
chain of inequalities

$\min_{\mu\in\mathcal{P}_{N\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)}V[\mu]\leq\min_{\mu\in\mathcal{P}_{(N+1)\textrm{-rep},\overline{\lambda}}\left(X^{2}\right)}V[\mu]\leq\dots\leq\min_{\mu\in\mathcal{D}_{\infty\textrm{-rep},\overline{\lambda}}}V[\mu]\leq
V[ER^{(\infty)}],$ (4.10)

where $V[\mu]:=\int_{X^{2}}v(x,y)d\mu(x,y)$. The inequalities (4.10) show that
for any number of marginals $N\geq 2$ we can find an upper bound of the
optimal value in (1.8) by computing the objective value of the ’attractive
limit extreme point’ $EA^{(\infty)}$, the ’repulsive limit extreme point’
$ER^{(\infty)}$ and the ’axis limit extreme points’ $E12^{(\infty)}$,
$E13^{(\infty)}$ as well as $E23^{(\infty)}$ and choosing the smallest one.
Note that this improves the in physics common mean field approximation
$V[ER^{(\infty)}]$, where one usually considers repulsive pair-costs
$v:X\times X\to\mathbb{R}$.

Finally, we note that the volume portion of the reduced Kantorovich polytope
that is occupied by the reduced Monge polytope exhibits oscillatory behavior
with decreasing amplitude when interpreted as a function of the number of
marginals $N$, see Figure 4. The considered volumetric ratio oscillates around
a value above 0.9 where even marginals when directly compared to the odd
marginals produce a higher ratio. In the sense that an optimizer in the
occupied volume yields the existence of a Monge-solution, the Monge ansatz
seems to be ’better’ for an even number of marginals.

## References

  * [1] Martial Agueh and Guillaume Carlier.  Barycenters in the Wasserstein Space.  SIAM J. Math. Anal., 43(2):904–924, 2011.
  * [2] Mathias Beiglböck, Pierre Henry-Labordère, and Friedrich Penkner.  Model-independent bounds for option prices - a mass transport approach.  Finance Stoch., 17(3):477–501, 7 2013.
  * [3] Dimitri P. Bertsekas.  Convex Optimization Theory.  Athena Scientific, Belmont, MA, 1 edition, 2009.
  * [4] Garrett Birkhoff.  Tres observaciones sobre el algebra lineal.  Universidad Nacional de Tucuman Revista, Serie A, 5(3):147–151, 1946.
  * [5] Rainer Burkard, Mauro Dell’Amico, and Silvano Martello.  Assignment Problems. Revised reprint. Society for Industrial and Applied Mathematics, 2012.
  * [6] Giuseppe Buttazzo, Luigi De Pascale, and Paola Gori-Giorgi.  Optimal-transport formulation of electronic density-functional theory.  Phys. Rev. A, 85:062502, 6 2012.
  * [7] Guillaume Carlier.  On a Class of Multidimensional Optimal Transportation Problems.  J. Convex Anal., 10(2):517–529, 02 2003.
  * [8] Guillaume Carlier and Ivar Ekeland.  Matching for Teams.  Econom. Theory, 42(2):397–418, 2 2010.
  * [9] Guillaume Carlier and Bruno Nazaret.  Optimal transportation for the determinant.  ESAIM Control Optim. Calc. Var., 14(4):678–698, 2008.
  * [10] Huajie Chen, Gero Friesecke, and Christian B. Mendl.  Numerical Methods for a Kohn-Sham Density Functional Model Based on Optimal Transport.  Journal of Chemical Theory and Computation, 10(10):4360–4368, 2014\.  PMID: 26588133.
  * [11] Pierre-André Chiappori, Robert J. McCann, and Lars P. Nesheim.  Hedonic price equilibria, stable matching, and optimal transport: equivalence, topology, and uniqueness.  Econom. Theory, 42(2):317–354, Feb 2010.
  * [12] Maria Colombo, Luigi De Pascale, and Simone Di Marino.  Multimarginal Optimal Transport Maps for One-dimensional Repulsive Costs.  Canad. J. Math., 67:350–368, 2013.
  * [13] Codina Cotar, Gero Friesecke, and Claudia Klöppelberg.  Density Functional Theory and Optimal Transportation with Coulomb Cost.  Comm. Pure Appl. Math., 66(4):548–599, 2013.
  * [14] Codina Cotar, Gero Friesecke, and Brendan Pass.  Infinite-body optimal transport with Coulomb cost.  Calc. Var. Partial Differential Equations, 54(1):717–742, 9 2015\.
  * [15] J. Csima.  Multidimensional Stochastic Matrices and Patterns.  J. Algebra, 14(2):194 – 202, 1970.
  * [16] Gero Friesecke.  A simple counterexample to the Monge ansatz in multi-marginal optimal transport, convex geometry of the set of Kantorovich plans, and the Frenkel-Kontorova model.  ArXiv e-prints, 8 2018.
  * [17] Gero Friesecke, Christian B. Mendl, Brendan Pass, Codina Cotar, and Claudia Klüppelberg.  N-density representability and the optimal transport limit of the Hohenberg-Kohn functional.  The Journal of Chemical Physics, 139(16):164109, 2013.
  * [18] Gero Friesecke and Daniela Vögler.  Breaking the Curse of Dimension in Multi-Marginal Kantorovich Optimal Transport on Finite State Spaces.  SIAM J. Math. Anal., 50(4):3996–4019, 2018.
  * [19] Alfred Galichon, Pierre Henry-Labordère, and Nizar Touzi.  A stochastic control approach to no-arbitrage bounds given marginals, with an application to lookback options.  Ann. Appl. Probab., 24(1):312–336, 2014.
  * [20] Wilfrid Gangbo and Andrzej Święch.  Optimal maps for the multidimensional Monge-Kantorovich problem.  Comm. Pure Appl. Math., 51(1):23–45, 1998.
  * [21] Augusto Gerolin, Anna Kausamo, and Tapio Rajala.  Non-existence of optimal transport maps for the multi-marginal repulsive harmonic cost.  ArXiv e-prints, 5 2018.
  * [22] Henri Heinich.  Probléme de Monge pour probabilités.  C. R. Math. Acad. Sci. Paris, 334(9):793–795, 12 2002.
  * [23] Lars Hoermander.  Notions of Convexity.  Birkhäuser, Boston, Massachusetts, 1994.
  * [24] The MathWorks Inc.  MATLAB and Statistics Toolbox R2018b (MATLAB 9.5).  Natick, Massachusetts, United States, 2018.
  * [25] V. M. Kravtsov.  Combinatorial properties of noninteger vertices of a polytope in a three-index axial assignment problem.  Cybernet. Systems Anal., 43(1):25–33, Jan 2007.
  * [26] Nathan Linial and Zur Luria.  On the Vertices of the d-Dimensional Birkhoff Polytope.  Discrete Comput. Geom., 51(1):161–170, Jan 2014.
  * [27] Abbas Moameni and Brendan Pass.  Solutions to multi-marginal optimal transport problems concentrated on several graphs.  ESAIM Control Optim. Calc. Var., 23(2):551–567, 2017.
  * [28] Brendan Pass.  Uniqueness and Monge Solutions in the Multimarginal Optimal Transportation Problem.  SIAM J. Math. Anal., 43(6):2758–2775, 2011.
  * [29] Brendan Pass.  On the local structure of optimal measures in the multi-marginal optimal transportation problem.  Calc. Var. Partial Differential Equations, 43(3):529–536, Mar 2012\.
  * [30] Brendan Pass.  Remarks on the semi-classical Hohenberg-Kohn functional.  Nonlinearity, 26(9):2731, 2013.
  * [31] The polymake team.  polymake 3.2.  Discrete Mathematics/Geometry Institut für Mathematik der Technischen Universität Berlin, Straße des 17. Juni 136, Berlin, Deutschland, 2018.
  * [32] Julien Rabin, Gabriel Peyré, Julie Delon, and Marc Bernot.  Wasserstein Barycenter and Its Application to Texture Mixing.  Scale Space and Variational Methods in Computer Vision, pages 435–446, 2012.
  * [33] R. Tyrrell Rockafellar.  Convex Analysis.  Princeton University Press, New Jersey, 1997.
  * [34] Frits C. R. Spieksma.  Multi Index Assignment Problems: Complexity, Approximation, Applications, pages 1–12.  Springer US, Boston, MA, 2000.
  * [35] Cédric Villani.  Optimal Transport: Old and New.  Springer Verlag, Berlin Heidelberg, 2009.
  * [36] John von Neumann.  A certain zero-sum two person game equivalent to the optimal assignment problem.  Contributions to the Theory of Games, 11:5–12, 1953.