<figure><img src="image_1.png" /></figure>

HAL
open science

# The football model, stochastic ordering and martingale transport

Gaoyue Guo, Nicolas Juillet, Wenpin Tang

<figure><img src="image_2.png" /></figure>

## To cite this version:

Gaoyue Guo, Nicolas Juillet, Wenpin Tang. The football model, stochastic ordering and martingale transport. 2025. hal-04983005v2

HAL Id: hal-04983005

https://hal.science/hal-04983005v2

Preprint submitted on 15 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
THE FOOTBALL MODEL, STOCHASTIC ORDERING AND
MARTINGALE TRANSPORT

GAOYUE GUO, NICOLAS JUILLET, AND WENPIN TANG

ABSTRACT. Tournaments are competitions between a number of teams, the outcome of which determines the relative strength or rank of each team. In many cases, the strength of a team in the tournament is given by a score. Perhaps the most fundamental mathematical result in the theory of random tournaments is Moon's theorem, which provides a necessary and sufficient condition for a feasible score sequence via majorization. To give a probabilistic interpretation of Moon's result, Aldous and Kolesnik introduced the football model, the existence of which gives a short proof of Moon's theorem. However, the proof of Aldous and Kolesnik is "noncanonical", leading to the question of a canonical construction of the football model. The purpose of this paper is to provide explicit constructions of the football model with an additional stochastic ordering constraint, which can be formulated by martingale transport. Two solutions are given: one is by solving an entropy optimization problem via Sinkhorn's algorithm, and the other relies on the idea of shadow couplings. It turns out that both constructions yield the property of strong stochastic transitivity. Nontransitive version of the football model is also considered.

*Key words*: Entropy optimization, martingale transport, pairwise comparison, score sequence, Sinkhorn's algorithm, shadow coupling, football model, stochastic ordering, strong stochastic transitivity, tournament.

AMS 2010 Mathematics Subject Classification: 60E15, 62F07, 60G42, 49Q22.

# 1. INTRODUCTION

A tournament refers to pairwise competitions between a number of teams (or players), all participating in a sport or game, in order to determine the winner or to produce a ranking of the teams. There has been a long history of research in tournaments and their rankings, including psychology [68, 69], game and economic theory [6, 42, 43], combinatorial and graph theory [29, 30, 51], applied probability [2, 3, 67], statistics [25, 38, 47], and more recently, large language models via direct preference optimization [22, 59]. In everyday language, a tournament often means a *single-elimination* or *knockout tournament*. This paper focuses on the *n*-team *round-robin tournament*, where each team competes against every other team. This often refers to *pairwise comparisons* in the statistics literature. A common metric for evaluating the strengths of the teams in the tournament is the *score sequence*, or simply, the *score* [41, 46]. We will consider the football model [4], which is an alternative to the famous Bradley-Terry model [17]. The main contribution of this work is to provide two explicit constructions of the football model that satisfies a stochastic ordering constraint, exploiting various probabilistic techniques such as stochastic ordering, martingale transport and entropy optimization. This answers a question of Aldous and Kolesnik (see Subsection 1.5 for details).

Date: October 15, 2025.
1.1. **Generalized tournament matrices and Moon's theorem.** To provide context, let us first introduce the notion of tournament matrices and scores. Consider an n-team tournament where each team competes N times against every other opponent. For $i \neq j$, let $p_{ij} \in [0, 1]$ be a number representing the relative strength between the teams i and j. It can be interpreted as:

* the proportion of the number of wins of team i over team j in N games;

* the probability that team i wins over its opponent j in a random game.

**Definition 1.1** (Generalized tournament matrices and scores). Let $n \ge 1$.

(1) Denote by $G_n$ the set of matrices $P = (p_{ij})_{1 \le i \ne j \le n}$, where $p_{ij} \in [0, 1]$ and $p_{ij} + p_{ji} = 1$ for each $i \neq j$. Such matrices are called the generalized tournament matrices of dimension n. (Note that the diagonals are undetermined.) Moreover, if $p_{ij} \in \{0, 1\}$ for each $i \neq j$, then P is a tournament matrix.

(2) Denote by $G'_n$ the set of matrices $P = (p_{ij})_{1 \le i,j \le n}$, where $p_{ii} = 1/2$ for each $i = 1, \dots, n$ and $(p_{ij})_{1 \le i \ne j \le n} \in G_n$.

(3) For $\mathbf{x} = (x_1, \dots, x_n) \in [0, +\infty)^n$, the subset $G_n(\mathbf{x}) \subset G_n$ (resp. $G'_n(\mathbf{x}) \subset G'_n$) denotes the collection of matrices $P$ such that

$$ \sum_{j \neq i} p_{ij} = x_i \quad \left( \text{resp. } \sum_{j=1}^{n} p_{ij} = x_i \right), $$

for each $i = 1, \dots, n$. The vector $\mathbf{x}$ is called the (generalized tournament) score.

In the probabilistic setting, for every $P \in G_n(\mathbf{x})$, $\mathbf{x}$ can be interpreted as the vector of expected number of wins because $x_i = \sum_{j \neq i} p_{ij}$ for each $i$. By convention, for $P \in G'_n(\mathbf{x})$, a half win is artificially added by setting $p_{ii} = 1/2$ (this would correspond to a match of a team against a copy of itself). In the deterministic case, $Nx_i$ is the number of games that team i wins. A classical theorem of [50] provides a necessary and sufficient condition for the score $\mathbf{x}$ so that $G_n(\mathbf{x})$ (or $G'_n(\mathbf{x})$) is not empty. To state Moon's result, we need the notion of majorization.

**Definition 1.2** (Majorization). Let $n \ge 1$. On the set of vectors $\mathbf{x} = (x_1, \dots, x_n)$ with increasing coordinates (i.e., $x_1 \le \dots \le x_n$), the partial order of majorization is defined by

$$ \mathbf{x} \preceq \mathbf{y} \quad \text{if and only if} \quad \sum_{i=1}^{n} x_i = \sum_{i=1}^{n} y_i \quad \text{and} \quad \sum_{i=1}^{k} x_i \ge \sum_{i=1}^{k} y_i, \quad \text{for each } k < n. \quad (1.1) $$

This order<fnref n="1" /> extends to $\mathbb{R}^n$ by $\mathbf{x} \preceq \mathbf{y}$ if and only if $\tilde{\mathbf{x}} \preceq \tilde{\mathbf{y}}$, where the vectors $\tilde{\mathbf{x}}, \tilde{\mathbf{y}}$ are the rearranged versions of $\mathbf{x}, \mathbf{y}$ (i.e., $\tilde{\mathbf{x}} = (x_{\sigma(1)}, \dots, x_{\sigma(n)})$, $\tilde{\mathbf{y}} = (y_{\sigma'(1)}, \dots, y_{\sigma'(n)})$ for some permutations $\sigma, \sigma'$) such that the coordinates of $\tilde{\mathbf{x}}, \tilde{\mathbf{y}}$ are nondecreasing.

**Theorem 1.3** ([50]). Let $\mathbf{x} \in \mathbb{R}^n$.

(1) $G_n(\mathbf{x})$ is nonempty if and only if $\mathbf{x} \preceq (0, 1, \dots, n-1)$.

(2) $G'_n(\mathbf{x})$ is nonempty if and only if $\mathbf{x} \preceq (\frac{1}{2}, \frac{3}{2}, \dots, n-\frac{1}{2})$.

<footnote n="1">¹On $\mathbb{R}^n$, it is only a preorder: the reflexivity fails. See Appendix A for reminders.</footnote>
Note that a result similar to Theorem 1.3 for tournament matrices (see Definition 1.1 (1)) and $\mathbf{x} \in \mathbb{Z}^n$ was proved by H. G. [41]. See also [39], [8, 23, 27, 64] for recent developments on the enumeration of score sequences.

1.2. **Zermelo(-Bradley-Terry) model.** Definition 1.1 is quite general. Several parametric models have been developed to provide further structures to $p_{ij}$'s in pairwise comparisons. The most popular and well-studied example is the Zermelo model [72], which is widely known as the Bradley-Terry model [17].

[72] is arguably the first to consider the inference problem in tournaments. In his model, each team $i$'s strength is specified by a positive number $u_i$, which is called the “force” (Spieltärke in German), with $\sum_{i=1}^n u_i = 1$. For some $k \ge 1$, the outcome of the $k^{th}$ game between the teams $i$ and $j$ can be represented by a Bernoulli variable $B_{ij}^k$ with parameter $\frac{u_i}{u_i+u_j}$, where $\{B_{ij}^k = 1\}$ means that team $i$ beats team $j$ in their $k^{th}$ game, and $\{B_{ij}^k = 0\}$ indicates the other way around. Hence, the generalized tournament matrix is specified by $p_{ij} := \frac{u_i}{u_i+u_j}$. Assuming that the random variables $(B_{ij}^k)_{i,j,k}$ are all independent, the maximum likelihood estimate (MLE) is used to infer the vector parameter $\mathbf{u} = (u_1, \dots, u_n)$. Under an irreducibility condition (see Definition 3.1 below), the MLE is uniquely determined by a system of $n$ equations<fnref n="2" />

$$x_i = \sum_{j \neq i} \frac{u_i}{u_i + u_j}, \quad \text{for } 1 \le i \le n. \qquad (1.2)$$

The Zermelo model was rediscovered in [17, 28]. Since then, there have been various extensions such as the Plackett-Luce model [45, 58] and the Mallows model [47]. See [16, 63] for historical notes and further references.

The Zermelo model can be reparametrized by $v_i = \log u_i$, which yields the (generalized) linear model [49, Section 7.5]:

$$p_{ij} = \frac{u_i}{u_i + u_j} = \frac{1}{1 + e^{v_j - v_i}}, \quad \text{for } 1 \le i \neq j \le n. \qquad (1.3)$$

**Remark 1.4.** Here we provide some explanations of the reparametrization $v_i = \log u_i$. The idea, following Thurston and Mosteller [68, 69, 53], is to associate each player $i$ with an independent random variable $X_i$, and the generalized tournament matrix is specified by $p_{ij} = P(X_i \ge X_j)$ for $1 \le i, j \le n$. One choice is that $X_i = v_i + V$, where $v_i$ is the indicator of strength, and $V$ is the noise. Hence, $p_{ij} = P(W \le v_i - v_j)$, where $W$ is distributed as $V - V'$, with $V'$ an independent copy of $V$. In practice, $W$ is often only required to be a (symmetric) random variable. Specializing to the case where $W$ has the cumulative distribution function $F(t) = \frac{1}{1+e^{-t}}$ recovers the model (1.3). Alternatively, we can take $X_i$ to be an independent exponential random variable with parameter $u_i^{-1}$ to recover the model (1.3).

It is easy to check that the model (1.3) enjoys the property of strong stochastic transitivity (SST):

$$p_{ij} \ge 1/2 \text{ and } p_{jk} \ge 1/2 \quad \implies \quad p_{ik} \ge \max(p_{ij}, p_{jk}), \qquad (1.4)$$

because the left side of (1.4) is equivalent to $v_i \ge v_j \ge v_k$. This property will be the center of our study in Sections 2 and 5.

<footnote n="2">²Though the likelihood is homogeneous of degree zero, the uniqueness is guaranteed by the constraint $\sum_{i=1}^n u_i = 1$.</footnote>
1.3. **The football model.** Recently, [4] introduced a new parametric model for tournaments – the football model. The idea is in the same spirit to Remark 1.4 by associating each team with a random variable, or equivalently a probability distribution, for paired comparisons. It is described as follows.

The model is parametrized by

$$ \Theta_n := \left\{ (\mu_1, \dots, \mu_n) \in \mathcal{P}(\{0, \dots, n-1\})^n : \sum_{i=1}^n \mu_i = \sum_{k=0}^{n-1} \delta_k \right\}, \quad (1.5) $$

where $\mathcal{P}(\{0, \dots, n-1\})$ denotes the set of probability measures on $\{0, \dots, n-1\}$, and $\delta_k$ is the Dirac measure on $k$. Each team $i$ is assigned a probability measure $\mu_i$. For every $k \ge 1$, the outcome of the $k^{th}$ game between the teams $i$ and $j$ is determined as follows: Let $X_i^k \sim \mu_i$ and $X_j^k \sim \mu_j$ be independent.

* If $X_i^k > X_j^k$, then team $i$ beats team $j$.

* If $X_i^k < X_j^k$, then team $j$ beats team $i$.

* If $X_i^k = X_j^k$, then each team is granted a probability $\frac{1}{2}$ (by external randomization) to win the game.

Here the random variables $X_i^k$ and $X_j^k$ can be interpreted as the number of goals scored by the team $i$ and $j$, which explains the name “football”. When the scores are equal, a tie-breaking method such that penalty shoot-out or coin tossing will take place and decide which team wins the game with equal probability. Specifically, let $Z_{ij}^k$ be a Bernoulli variable with parameter $\frac{1}{2}$, independent of $(X_i^k, X_j^k)$. Set

$$ B_{ij}^k := \mathbb{1}_{\{X_i^k > X_j^k\}} + \mathbb{1}_{\{X_i^k = X_j^k, Z_{ij}^k = 1\}}, $$

so $\{B_{ij}^k = 1\}$ means that team $i$ beats team $j$ in their $k^{th}$ game, and $\{B_{ij}^k = 0\}$ indicates the other way around. The corresponding tournament matrix is specified by

$$ p_{ij} := \mathbb{P}(X_i > X_j) + \frac{1}{2}\mathbb{P}(X_i = X_j), \quad (1.6) $$

where $X_i \sim \mu_i$ for $1 \le i \le n$, and $(X_i)_{1 \le i \le n}$ are pairwise independent. The constraint $\sum_{i=1}^n \mu_i = \sum_{k=0}^n \delta_k$ in the definition (1.5) implies that the random number of goals of a random team in a match is uniform on $\{0, 1, \dots, n-1\}$. This, of course, is not realistic when $n$ is large. However, it is a central element in the short probabilistic proof of Aldous and Kolesnik that will be recalled in Section 1.4.

1.4. **Proof of Moon's theorem via the football model.** As pointed out in [4], a remarkable property of this model is that for $(\mu_1, \dots, \mu_n) \in \Theta_n$, the score $x_i = \sum_{j \ne i} p_{ij}$ of team $i$ after exactly one game against every other opponent is the expected number of goals (or points) $\mathbb{E}X_i = \int x d\mu_i(x)$. To see this, let $\chi(\mu_i, \mu_j)$ be the probability for team $i$ to win over team $j$. The map $\chi(\cdot, \cdot)$ can be extended into a linear map for signed measures with compact support. In fact, $\chi(\alpha, \beta) := \int f(y-x)d(\alpha \otimes \beta)(x, y)$, where $f(z) := \mathbb{1}_{\{z>0\}} + \frac{1}{2}\mathbb{1}_{\{z=0\}}$. Denoting by $\lambda := \sum_{k=0}^{n-1} \delta_k$, we have

$$ \chi(\mu_i, \lambda) = \sum_{k=0}^{n} \chi(\mu_i, \delta_k) = \frac{1}{2} + \mathbb{E}X_i, $$
on one hand, and

$$ \chi(\mu_i, \lambda) = \chi(\mu_i, \mu_i) + \sum_{j \neq i} \chi(\mu_i, \mu_j) = \frac{1}{2} + \sum_{j \neq i} p_{ij}, $$

on the other hand.

As a consequence, the more difficult implication in Moon's theorem (Theorem 1.3) can easily be proved by using a famous result commonly attributed to [65] for the set of probability measures with a finite first moment. Earlier versions for discrete measures with finite support such as in Muirhead's inequality (see e.g., [31]) are sufficient for this purpose.

**Theorem 1.5** ([31, 20, 65]). Let $\rho$ and $\mu$ be two probability measures on $\mathbb{R}^d$ having a finite first moment. Then the following two conditions are equivalent:

(1) The measures are in convex order $\rho \preceq_C \mu$, i.e., for any convex function $\varphi : \mathbb{R}^d \to \mathbb{R}$, we have $\int \varphi d\rho \le \int \varphi d\mu$.

(2) There exists a pair of random variables $(X, Y)$ such that $X \sim \rho$, $Y \sim \mu$, and $\mathbb{E}(Y|X = x) = x$ for $\rho$-almost every $x$.

Moreover, for $d=1$, if $\rho$ and $\mu$ are uniform measures on $\{x_1 \le \dots \le x_n\}$ and $\{y_1 \le \dots \le y_n\}$ respectively, the conditions (1) and (2) are satisfied if and only if $\mathbf{x} \preceq \mathbf{y}$. In this case, the condition (2) is translated as follows: For $1 \le i \le n$, let $\mu_i$ be the conditional distribution of $Y$ given $X = x_i$. We have $\sum_{i=1}^n \mu_i = \sum_{j=1}^n \delta_{y_j}$ and $\int y d\mu_i(y) = x_i$ for each $i$.

The following proof of Theorem 1.3 by Aldous and Kolesnik is hard to beat.

*Proof of Theorem 1.3.* The implication $\mathcal{G}_n(\mathbf{x})$ is nonempty $\Rightarrow \mathbf{x} \preceq (0, \dots, n-1)$ follows from the fact that the expected number of wins of the $k$-weakest teams must be no less than $1 + \dots + (k-1) = \binom{k}{2}$.

For the more difficult implication, assume that $\mathbf{x} \preceq (0, \dots, n-1)$. By Theorem 1.5, there exists $(\mu_1, \dots, \mu_n) \in \Theta_n$ such that $\int y d\mu_i(y) = x_i$ for each $i$. In the football model, $\int y d\mu_i(y)$ is the score of team $i$. So the matrix $P = (p_{ij})_{1 \le i \ne j \le n}$ defined by (1.6) is a generalized tournament matrix in $\mathcal{G}_n(\mathbf{x})$. $\square$

**1.5. Motivation and guideline.** The above short proof relies on the observation that for $\mathbf{x} \preceq (0, \dots, n-1)$, the set

$$ \Theta_n(\mathbf{x}) := \left\{ (\mu_1, \dots, \mu_n) \in \Theta_n : \int y d\mu_i(y) = x_i \text{ for } 1 \le i \le n \right\}, \quad (1.7) $$

is nonempty. This fact follows from a soft argument by applying Theorem 1.5 that is non-constructive. It is natural to call for an explicit construction of $(\mu_1, \dots, \mu_n) \in \Theta_n(\mathbf{x})$, which was also asked in [4]. We quote them twice. First from the abstract:

In particular, our proof of Moon's theorem on mean score sequences seems more constructive than previous proofs. This provides a comparatively concrete introduction to a longstanding mystery, the lack of a canonical construction for a joint distribution in the representation theorem for convex order.

Quotation from the final “Discussion” section in [4]:
To us, the most interesting part of the bigger picture surrounding convex order is that there is apparently no “canonical” choice of joint distribution in (3), (8): proofs may be constructive but they involve rather arbitrary choices and the resulting joint distributions are not easily described. Recent literature on peacocks [7]<fnref n="3" /> studies continuous-parameter processes increasing in convex order, via many different constructions, and ideas from that literature might be relevant in our context.

Also note that the generic dimension of $\Theta_n(\mathbf{x})$ is $(n-1)(n-2)$, which is larger than the dimension $\frac{1}{2}(n-1)(n-2)$ of $G_n(\mathbf{x})$, and $n-1$ of the Zermelo model. The fact that the football model has more degrees of freedom is an advantage because it allows one more flexibility in modeling, e.g., to fit nontransitive situations. It is also a weakness because the system of equations that identifies $\Theta_n(\mathbf{x})$ is underdetermined, which is not the case for the Zermelo model. It is worth noting that [40] recently introduced *Coxeter tournaments*, where players may collaborate and play solitaire games. They developed alternative versions of the Bradley-Terry model and the football model, with proofs that are more algebraic and geometric in nature.

The main objective of this paper is to provide explicit constructions of $(\mu_1, \dots, \mu_n) \in \Theta_n(\mathbf{x})$ with the SST as an additional constraint. This approach provides “canonical” representations of the football model, thereby responding to the invitation mentioned above. Our approach is based on martingale optimal transport, which is a topic closely related to the peacocks (*Processus Croissant pour l'Ordre Convexe*). Two different algorithmic solutions are given: one is obtained by solving an entropic martingale transport problem via Sinkhorn’s algorithm [56, 61], and the other is related to the concept of a *shadow* that was introduced to define a class of interesting martingale transport plans [12, 13, 10].

**Organization of the paper:** The remainder of the paper is organized as follows. In Section 2, we prove (in an abstract way) that $\Theta_n(\mathbf{x})$ contains an element that yields the SST. Two explicit constructions of $(\mu_1, \dots, \mu_n) \in \Theta_n(\mathbf{x})$ with the stochastic ordering constraint are presented in Sections 3 and 4. In an opposite direction, we illustrate in Section 5 that $\Theta_n(\mathbf{x})$ also has nontransitive solutions. Finally, general theoretic clarifications on the partial orders induced by a generalized tournament matrix are given in Appendix A.

## 2. EXISTENCE OF AN SST SOLUTION IN THE FOOTBALL MODEL

We have seen in Subsection 1.2 that the Zermelo model (1.3) enjoys the property of SST (1.4). Here we show that the SST also appears in the football model. In fact, for $\mathbf{x} \preceq (0, 1, \dots, n-1)$, $\Theta_n(\mathbf{x})$ contains an element that satisfies an additional constraint $\mu_1 \preceq_{\text{sto}} \dots \preceq_{\text{sto}} \mu_n$ (see Proposition 2.2, with the definition of $\preceq_{\text{sto}}$ being recalled at the end of this introduction). With this intermediate constraint, the corresponding generalized tournament matrix enjoys the SST (see Proposition 2.4). The proof of Proposition 2.2 relies on a result of [54] (see also [55, Section 2.6]) concerning the existence of conditional martingale kernels that are increasing in stochastic order. This result is equivalent to the existence of 1-Lipschitz martingale transport plans (see e.g., [11] before Lemma 3.3), which is well known in the field of peacocks and martingale transport, and has been discovered several times independently. For instance, [37] proved the existence in a nonconstructive way based on Choquet theory;

<footnote n="3">³It is the book by Hirsch, Profeta, Roynette and Yor referred to as [32] in the present paper.</footnote>
[44] is based on Hobson's approach to the Skorokhod embedding problem (SEP); [11] relies on Root's solution that is also not constructive; and the sunset coupling [13] is also related to the SEP, and to Kellerer's solution for which it gives a more explicit construction.

Here two novel approaches are developed in the context of the football model. The first method is on entropy optimization in the space of martingale transport plans, which will be detailed in Section 3. This method was proposed by [34], but was applied to $P = (p_{ij})_{1\le i\ne j\le n}$ in $G_n(\mathbf{x})$ instead of $\Theta_n(\mathbf{x})$.<fnref n="4" /> The second method is a direct construction using shadows, and will be explained in Section 4.

Recall that two probability measures $\mu$ and $\mu'$ on $\mathbb{R}$ are in stochastic order $\mu \preceq_{\text{sto}} \mu'$ if and only if their cumulative distribution functions satisfy $F_{\mu} \ge F_{\mu'}$. There are many equivalent criterions, e.g., there exists a coupling $(X, X')$ such that $X \sim \mu, X' \sim \mu'$ and $\mathbb{P}(X \le X') = 1$ (see [60, Chapter 1]).

**Lemma 2.1.** Let $\mu$ and $\nu$ be two probability measures on $\mathbb{R}$ such that $\mu \preceq_{\text{sto}} \nu$. Then for $X \sim \mu$ and $Y \sim \nu$ independent,

$$ \mathbb{P}(X > Y) + \frac{1}{2}\mathbb{P}(X = Y) \ge 1/2. $$

Moreover, if $\mu \ne \nu$, then the equality is strict.

*Proof.* Consider a coupling $(X, X')$ such that $X \sim \mu, X' \sim \nu$ and $\mathbb{P}(X' \ge X) = 1$. Let $Y \sim \nu$ be independent of $(X, X')$ so that $(X, Y) \sim \mu \times \nu$. Since $(Y, X') \sim \nu \times \nu$, we get

$$ \frac{1}{2} = \mathbb{P}(Y > X') + \frac{1}{2}\mathbb{P}(Y = X') = \mathbb{E}[f(Y - X')], $$

where

$$ f(z) := \mathbb{1}_{\{z>0\}} + \frac{1}{2}\mathbb{1}_{\{z=0\}}. \qquad (2.1) $$

Observing that $Y - X' \le Y - X$ almost surely, we have:

$$ \frac{1}{2} = \mathbb{E}[f(Y - X')] \le \mathbb{E}[f(Y - X)] = \mathbb{P}(Y > X) + \frac{1}{2}\mathbb{P}(Y = X). \qquad (2.2) $$

Furthermore, if $\mu \ne \nu$, the event $\{X < X'\}$ has nonzero probability. Note that $Y$ is independent from $(X, X')$, and has the same law as $X'$. Therefore, $\mathbb{P}(X < Y \le X') > 0$. It follows $\mathbb{P}(f(Y - X') \le 1/2, f(Y - X) = 1) > 0$, and the inequality in (2.2) is strict. $\square$

The following proposition is an improved version of Strassen's theorem, following [54].

**Proposition 2.2.** Let $\mathbf{x}, \mathbf{y} \in \mathbb{R}^n$ be such that $x_1 \le \dots \le x_n$ and $y_1 \le \dots \le y_n$ and assume $\mathbf{x} \preceq \mathbf{y}$. There exists $(\mu_1, \dots, \mu_n)$ such that $\sum_{i=1}^n \mu_i = \sum_{k=1}^n \delta_{y_k}$, $\int x d\mu_i = x_i$ for each $1 \le i \le n$ and $(\mu_i)_{1\le i\le n}$ is increasing in $\preceq_{\text{sto}}$.

If $\mathbf{y} = (0, 1, \dots, n-1)$, it can be formulated that there exists $(\mu_1, \dots, \mu_n) \in \Theta_n(\mathbf{x})$ such that $(\mu_i)_{1\le i\le n}$ is increasing in $\preceq_{\text{sto}}$.

*Proof.* As seen in the introduction, Theorem 1.5 has an improved version when $\rho, \mu$ are supported on $\mathbb{R}$. In addition to $\mathbb{E}(Y|X = x) = x$ for $\rho$-almost every $x$, it was proved in [54] that there exists a family $(\mu_x)_{x \in \mathbb{R}}$ of (regular) conditional laws (i.e., $\mathbb{E}(f(Y)|X = x) = \int f(y) d\mu_x(y)$ for every positive $f$) that is increasing in stochastic order: $x \le x'$ implies

<footnote n="4">⁴[4] also promoted Joe's approach. One main motivation of [4] is to make a connection between Joe's result and Moon's theorem.</footnote>
$\mu_x \preceq_{\text{sto}} \mu_{x'}$. Exactly as in the second part of Theorem 1.5, this result translates in the discrete setting into the statement of the proposition. $\square$

**Remark 2.3.** In Proposition 2.2, if $\mu_i \preceq_{\text{sto}} \mu_{i+1}$, then $\mu_i = \mu_{i+1}$ is equivalent to $x_i = x_{i+1}$. This can be seen by the coupling $(X_i, X_{i+1})$ such that $X_i \sim \mu_i$, $X_{i+1} \sim \mu_{i+1}$ and $X_i \le X_{i+1}$ almost surely. Then $\mathbb{E}X_i = \mathbb{E}X_{i+1}$ is equivalent to $X_i = X_{i+1}$ almost surely.

Combining Lemma 2.1, Proposition 2.2 and Remark 2.3 yields that for each $\mathbf{x} \preceq (0, \dots, n-1)$, there exists $(\mu_1, \dots, \mu_n) \in \Theta_n(\mathbf{x})$ such that the corresponding generalized tournament matrix $P$ satisfies:

* $p_{ij} \ge 1/2$ if and only if $x_i \ge x_j$.

* $p_{ij} = 1/2$ if and only if $x_i = x_j$.

The next proposition proves the SST (1.4).

**Proposition 2.4.** For $\mathbf{x} \preceq (0, \dots, n-1)$, let $(\mu_1, \dots, \mu_n) \in \Theta_n(\mathbf{x})$ be specified as in Proposition 2.2. Then $(\mu_i)_{1 \le i \le n}$ is increasing in $\preceq_{\text{sto}}$, and hence, the generalized tournament matrix defined by (1.6) satisfies the SST.

*Proof.* Take $(\mu_1, \dots, \mu_n) \in \Theta_n(\mathbf{x})$ as in Proposition 2.2, and assume that $p_{kj} \ge 1/2$ and $p_{ji} \ge 1/2$. As explained just above, we have $x_i \le x_j \le x_k$. Since $\mu_i \preceq_{\text{sto}} \mu_j$, there exist $X_i \sim \mu_i$, $X_j \sim \mu_j$ and $X_k \sim \mu_k$ such that $X_i \le X_j$ almost surely and $X_k$ is independent of $(X_i, X_j)$. With the function $f$ defined in (2.1), we have $f(X_k - X_i) \ge f(X_k - X_j)$ almost surely. Taking the expectation on both sides we get $p_{ki} \ge p_{kj}$. Similarly, $\mu_j \preceq_{\text{sto}} \mu_k$ yields $p_{ki} \ge p_{ji}$. Thus, $p_{ki} \ge \max(p_{ji}, p_{kj})$. $\square$

**Remark 2.5.** In view of Propositions 2.2 and 2.4, one may construct tournament matrices satisfying the SST via the (unique) strong Markov martingale coupling, denoted by $\pi$, between the discrete measures:

$$ \sum_{i=1}^{n} \delta_{x_i} \quad \text{and} \quad \sum_{i=1}^{n} \delta_{i-1}. $$

A common approach to constructing $\pi$ relies on the fact that every discrete martingale can be embedded into a continuous-time martingale. More precisely, the so-called Bass martingale $(X_t)_{0 \le t \le 1}$ satisfies Law$(X_0, X_1) = \pi$, see, e.g., [7]. This martingale can be identified by solving a stochastic control problem, see [24, 1]. However, since

$$ \sum_{i=1}^{n} \delta_{x_i} \quad \text{and} \quad \sum_{i=1}^{n} \delta_{i-1}. $$

are both discrete measures, convergence of the numerical schemes suggested in the aforementioned references is not a priori guaranteed, and approximation of this stochastic control problem becomes more technical. We propose two alternative constructions that are specifically adapted to the discrete setting.

### 3. THE ENTROPIC CONSTRUCTION

In this section, we provide a construction of $(\mu_1, \dots, \mu_n) \in \Theta_n(\mathbf{x})$ satisfying the stochastic ordering constraint $\mu_1 \preceq_{\text{sto}} \dots \preceq_{\text{sto}} \mu_n$ by solving an entropy optimization problem. An iterative algorithm is given. For ease of presentation, we identify $\Theta_n$ and $\Theta_n(\mathbf{x})$, whose
elements are $(\mu_1, \dots, \mu_n) \in \mathcal{P}(\{0, \dots, n-1\})^n$, with the set of matrices $M = (m_{ij})_{1 \le i,j \le n}$ defined by $m_{ij} := \mu_i(\{j-1\})$. More precisely,

$$ \Theta_n = \{M \in \mathcal{M}_n(\mathbb{R}) : M \text{ is a doubly stochastic matrix}\}, $$

$$ \Theta_n(\mathbf{x}) = \left\{ M \in \Theta_n : \sum_{j=1}^{n} m_{ij}(j-1) = x_i \text{ for } 1 \le i \le n \right\}, $$

where the measure $\mu_i$ is encoded as a vector of $n$ nonnegative coefficients corresponding to the $i$-th row of $M$ by $\mu_i = \sum_{j=1}^n m_{ij}\delta_{j-1}$ for $1 \le i \le n$.

To proceed further, we need the following irreducibility definition, which is implicit in the results established in [28, 52, 72].

**Definition 3.1** (Irreducibility condition). Let $P \in \mathcal{G}_n(\mathbf{x})$, and $\tilde{\mathbf{x}}$ be the rearranged version of $\mathbf{x}$. The following conditions are equivalent:

(1) For each (nontrivial) partition $I \cup I^c$ of $\{1, \dots, n\}$, there exist $i \in I$ and $j \in I^c$ such that $p_{ij} > 0$.

(2) For each $i, j \in \{1, \dots, n\}$, there exist $r \le n$ and a chain of coefficients $(i_\ell)_{\ell=1, \dots, r}$ such that $i_1 = i$, $i_r = j$, and $p_{i_\ell i_{\ell+1}} > 0$ for each $\ell < r$. This condition can be interpreted as the strong connectivity in the oriented graph with vertices $\{1, \dots, n\}$ and edges $s \to t$ open if and only if $p_{st} > 0$.

(3) There exists $r \ge 1$ such that all the entries of $P^r = \underbrace{P \cdots P}_{r \text{ times}}$ are strictly positive.

$$ (4) \sum_{i=1}^{n} \tilde{x}_i = \frac{n(n-1)}{2} \text{ and } \sum_{i=1}^{k} \tilde{x}_i < \frac{k(k-1)}{2} \text{ for } k < n. $$

If one of the conditions (1)-(3) is satisfied, we say that $P$ is an irreducible generalized tournament matrix. If the condition (4) is satisfied, $\mathbf{x}$ is called an irreducible score.

The following two conditions for reducible (i.e., non irreducible) generalized tournament matrices and scores are equivalent:

(1) There exists a partition $I \cup I^c$ such that $p_{ij} = 0$ for each $(i, j) \in I \times I^c$.

(2) There exists $1 \le k < n$ such that $(\tilde{x}_1, \dots, \tilde{x}_k) \preceq (0, \dots, k-1)$ and $(\tilde{x}_{k+1}, \dots, \tilde{x}_n) \preceq (k, \dots, n-1)$.

In particular, if $P \in \mathcal{G}_n(\mathbf{x})$ is irreducible (resp. reducible), then so is every $Q \in \mathcal{G}_n(\mathbf{x})$.

For $\mathbf{x} \preceq (0, \dots, n-1)$, define the entropy by

$$ H : M \in \Theta_n(\mathbf{x}) \mapsto \sum_{i,j=1}^{n} m_{ij} \log(m_{ij}). \qquad (3.1) $$

First, we show in Subsection 3.1 that if $\mathbf{x} \preceq (0, \dots, n-1)$ is an irreducible score, then $H$ has a unique minimizer whose entries are strictly positive. In Subsection 3.2, we prove that $P \in \mathcal{G}_n(\mathbf{x})$ corresponding to this unique minimizer satisfies the SST. Next in Subsection 3.3, we provide an algorithm and prove its convergence to the minimizer of $H$. Finally, we extend the results to the reducible case in Subsection 3.4.

<footnote n="5">⁵Notice $\frac{k(k-1)}{2} = 1 + 2 + \dots + k$ and look at Remark 3.4 for a similar condition where $x \preceq y \neq (1, \dots, n)$.</footnote>
**3.1. Minimizer of H.** Recall the definition of irreducible scores from Definition 3.1. The following proposition is useful in proving that $H$ has a unique minimizer, as well as the convergence of Sinkhorn's algorithm in Subsection 3.3. As pointed out by a referee, this result appears to be well-known as it corresponds to the main case in [18, Theorem 2], which was further developed by [19, 21]. However, the following proof is different, and seems to be original.

**Proposition 3.2.** Let $\mathbf{x} \preceq (0, \dots, n-1)$ be an irreducible score (so that $G_n(\mathbf{x})$ and $\Theta_n(\mathbf{x})$ are nonempty). Then there exists $M \in \Theta_n(\mathbf{x})$ such that $m_{ij} > 0$ for all $(i, j)$.

*Proof.* We consider the partial order on $\Theta_n(\mathbf{x})$, for which $M = (m_{ij})_{ij}$ dominates $N = (n_{ij})_{ij}$ if and only if $m_{ij} > 0 \Rightarrow n_{ij} > 0$ for all $(i, j)$. Let $M$ be maximal for this order. Here we gather a few observations on $M$:

* For each row $i$, the set of indices $j$ with $m_{ij} > 0$ is a nonempty (discrete) interval. Assume by contradiction that $m_{ij} = 0$, $m_{ij_0} > 0$ and $m_{ij_1} > 0$ with $j_0 < j < j_1$. Then there exists $i' \neq i$ such that $m_{i'j} > 0$. For $\lambda \in (0, 1)$ satisfying $\lambda y_{j_1} + (1-\lambda)y_{j_0} = y_j$, consider $M + hM'$, where $M'$ is a matrix with all entries zero except the six following $m'_{ij_0} = -m'_{i'j_0} = \lambda$, $m'_{ij_1} = m'_{i'j_1} = 1-\lambda$ and $m'_{i'j} = -m'_{ij} = 1$. For $h < \min(m_{ij}, m_{i'j_0}/\lambda, m_{i'j_1}/(1-\lambda))$, we have $M+hM' \in \Theta_n(\mathbf{x})$, which yields a contradiction to the maximality of $M$.

* If the four entries $m_{ij}$, $m_{ij'}$, $m_{i'j}$ and $m_{i'j'}$ are positive, then for the two rows $i$ and $i'$ we have $m_{ik} > 0 \Leftrightarrow m_{i'k} > 0$ (the positive entries are the same for the two rows). The proof follows the same argument as in the first point.

For row $i$, let $J_i$ be the set of indices $j$ such that $m_{ij} > 0$. With the two observations above and the fact that $\sum_j m_{ij}y_j = x_i$, we see that $M$ is a block matrix with the blocks going from left to right. Assume by contradiction that there are at least two blocks. Then $J_1 = \dots = J_k = \{1, \dots, \ell\}$ for some $k < n$, and $J_{k+1} \neq J_k$. In fact, $J_{k+1} \cap J_k = \emptyset$ or $\{\ell\}$. Consider two cases $k \ge \ell$ and $k \le \ell-1$: For $k \ge \ell$, note that $(i \le k$ and $j > \ell) \Rightarrow m_{ij} = 0$. Thus, the sum of all the entries is larger than that of all entries $m_{ij}$ with $i \le k$ or $j \ge \ell+1$, which is $k + (n-\ell)$. (Here $k$ corresponds to the sum over the $k$ first rows, and $n-\ell$ is the sum over the $n-\ell$ last columns.) Note that the sum of all the entries of $\Theta_n$ is $n$. By analyzing the equality case, we get $k = \ell$ and $J_{k+1} \cap J_k = \emptyset$, which contradicts the fact that $\mathbf{x}$ is irreducible because the probability that a team of the rows $0, \dots, k$ defeats a team indexed by $k+1, \dots, n$ would be zero (see the second part of Definition 3.1). For $k \le \ell-1$, we consider the lower left of $M$: $(i > k$ and $j < \ell) \Rightarrow m_{ij} = 0$. The sum of all the entries is larger or equal to $(n-k) + (\ell-1) \ge n$, which corresponds to the $n-k$ last rows plus the $\ell-1$ first columns. Again this sum is $n$, so $m_{1\ell} = \dots = m_{k\ell} = 0$ which contradicts the fact that $J_1 = \dots = J_k = \{1, \dots, \ell\}$. $\square$

The following result shows that if $\mathbf{x}$ is irreducible, the entropy function $H$ on $\Theta_n(\mathbf{x})$ has a unique minimum point.

**Proposition 3.3.** Let $\mathbf{x} \preceq (0, \dots, n-1)$ be an irreducible score. Then the function $H$ defined by (3.1) has a unique minimizer, whose coefficients $m_{ij}$'s are all strictly positive.

*Proof.* The function $H$ is continuous on the compact set $\Theta_n(\mathbf{x})$ (with the convention $0 \times \log 0 = 0$). Moreover, it is strictly convex so $H$ has a unique minimizer denoted by $M^{(0)}$.
Suppose by contradiction that one of the entries of $M^{(0)}$ is zero. By Proposition 3.2, let $M^{(1)} \in \Theta_n(\mathbf{x})$ with all strictly positive entries. For each $0 \le \lambda \le 1$, $M^{(\lambda)} := \lambda M^{(1)} + (1 - \lambda)M^{(0)} \in \Theta_n(\mathbf{x})$. Because all the coefficients of $M^{(1)}$ are strictly positive, the derivative of $\lambda \in [0, 1] \mapsto H(M^{(\lambda)})$ at $\lambda = 0^+$ is $-\infty$. This contradicts the fact that $M^{(0)}$ minimizes $H$. $\square$

**Remark 3.4.**

(1) Minimizing $P \in \mathcal{G}_n(\mathbf{x}) \mapsto \sum_{ij} p_{ij} \log(p_{ij})$ (instead of $H: M \mapsto \sum_{ij} m_{ij} \log(m_{ij})$) is used by [34] to construct a generalized tournament matrix of the Zermelo-Bradley-Terry model. Similar to Proposition 3.2, we can show that if $\mathbf{x}$ is irreducible, then $\mathcal{G}_n(\mathbf{x})$ contains some $P$ with all strictly positive entries. The argument is as follows: let $P \in \mathcal{G}_n(\mathbf{x})$ with the maximal number of nonzero entries. Assume by contradiction that $p_{i_1i_0} = 0$, so $p_{i_0i_1} = 1$. By Definition 3.1, there exists a chain $p_{i_1i_2}, \dots, p_{i_{N-1}i_N}, p_{i_N,i_0} > 0$. Now we operate as follows to get $p_{i_1i_0} > 0$ and guarantee $P$ stays in $\mathcal{G}_n(\mathbf{x})$: for sufficiently small $h$, replace $p_{ij_{i_1i_2}}$ with $p_{ij_{i_1i_2}} - h$, with the convention $i_N+1 = i_0$. Hence, $p_{ij_{i_1i_2}}$ is replaced with $p_{ij_{i_1i_2}} + h$. It suffices to use the infinite derivative of $p \mapsto p \log p$ at $p = 0^+$ as in Proposition 3.3 to conclude. Note that this argument also works for other functions $f$, provided that the derivative of $f$ at $0^+$ is $-\infty$.

(2) As was pointed out by a referee, both the entropic argument and Joe's theorem can be generalized to non round-robin tournaments. This has been done in [33, Theorem 7], where the complete graph case is equivalent to Joe's theorem.

(3) Proposition 3.3 can be easily generalized to irreducible pairs $\mathbf{x} \preceq \mathbf{y}$, where $\mathbf{y} = (k, k+1, \dots, \ell-1, \ell)$ is different from $(0, \dots, n-1)$. Definition 3.1 (4) can be adapted as follows: $\sum_{i=1}^{n'} \tilde{x}_{k-1+i} = \sum_{i=1}^{n'} \tilde{y}_{k-1+i}$ and $\sum_{i=1}^{k} \tilde{x}_i < \sum_{i=1}^{k} \tilde{y}_i$ for $k < n'$. (with $n' = \ell - k + 1$). In this case, $\Theta_n(\mathbf{x})$ is replaced by the set of double stochastic matrices with $\sum_{j=1}^{n'} m_{ij} y_{k-1+j} = x_{k-1+i}$ for each $1 \le i \le n'$.

**3.2. Property of SST.** The following theorem shows that for an irreducible $\mathbf{x}$, the minimizer of $H$ on $\Theta_n(\mathbf{x})$ satisfies the stochastic ordering constraint, and hence, the corresponding generalized tournament matrix enjoys the property of SST.

**Theorem 3.5.** Let $\mathbf{x} \preceq (0, \dots, n-1)$ be an irreducible score, and $M$ be the (unique) minimizer of $H$ on $\Theta_n(\mathbf{x})$. Then $(\mu_i)_{1 \le i \le n}$ is increasing in stochastic order, and the generalized tournament matrix corresponding to $M$ (defined by (1.6)) satisfies the SST.

*Proof.* Since $\mathbf{x}$ is irreducible, $H$ has a unique minimizer on $\Theta_n(\mathbf{x})$ by Proposition 3.3. Moreover, the minimizer is in the interior (of the affine linear space spanned by $\Theta_n(\mathbf{x})$). By the Karush-Kuhn-Tucker theorem [15, Section 5.5.3], we can neglect the constraints $m_{ij} \ge 0$, and only consider those on the rows with multipliers $a_i$'s, on the columns with multipliers $b_j$'s, and on the barycenters with multipliers $c_i$'s. Differentiating the Lagrangian, we get:

$$ (1 + \log(m_{ij})) + a_i + b_j + c_i y_j = 0, \quad \text{where } y_j = j - 1. $$

So we have $m_{ij} = e^{-a_i - 1} e^{-b_j} e^{-c_i y_j}$, and the probability vector $\mu_i = (m_{i1}, \dots, m_{in})$ is specified by the vector $(e^{-b_1}, \dots, e^{-b_n})$ modulated by entrywise product with the shape vector $(e^{-c_i y_1}, \dots, e^{-c_i y_n})$ and normalized by the constant $e^{-a_i}$. For ease of presentation, we fix
$(b_1, \dots, b_n)$, let $\bar{\mu} = (e^{-b_1}, \dots, e^{-b_n})$, and define $\bar{\mu}(s)$ by

$$ \bar{\mu}(s) = (\bar{\mu}_1(s), \dots, \bar{\mu}_n(s)) = (\bar{\mu}_1 e^{sy_1}, \dots, \bar{\mu}_n e^{sy_n}). $$

(Note that $\bar{\mu}(0) = \bar{\mu}$, and recall $y_j = j - 1$.) Let $C(s) = \sum_{j=1}^n \bar{\mu}_j e^{sy_j}$, so that $\mu_i = C(-c_i)^{-1} \bar{\mu}(-c_i)$ for $1 \le i \le n$. Now for $s < t$ and $1 \le k < n$, we have:

$$ \frac{\sum_{j=1}^{k} \bar{\mu}_j(t)}{\sum_{j=1}^{k} \bar{\mu}_j(s)} \le e^{(t-s)y_k} < e^{(t-s)y_{k+1}} \le \frac{\sum_{j=k+1}^{n} \bar{\mu}_j(t)}{\sum_{j=k+1}^{n} \bar{\mu}_j(s)}. $$

Therefore,

$$ \frac{\sum_{j=1}^{k} \bar{\mu}_j(t)}{\sum_{j=k+1}^{n} \bar{\mu}_j(t)} < \frac{\sum_{j=1}^{k} \bar{\mu}_j(s)}{\sum_{j=k+1}^{n} \bar{\mu}_j(s)} $$

so that

$$ g_k : s \mapsto C(s)^{-1} \sum_{j=1}^{k} \bar{\mu}_j(s) = \frac{\sum_{j=1}^{k} \bar{\mu}_j(s)}{\sum_{j=1}^{k} \bar{\mu}_j(s) + \sum_{j=k+1}^{n} \bar{\mu}_j(s)} = \left( 1 + \frac{\sum_{j=k+1}^{n} \bar{\mu}_j(s)}{\sum_{j=1}^{k} \bar{\mu}_j(s)} \right)^{-1} $$

is strictly increasing. Comparing the cumulative distribution functions of

$$ \mu(s) = \sum_{j=1}^{n} C(s)^{-1} \bar{\mu}_j(s) \delta_{y_j} \quad \text{and} \quad \mu(t) = \sum_{j=1}^{n} C(t)^{-1} \bar{\mu}_j(t) \delta_{y_j}, \qquad (3.2) $$

we see that $\mu(s) \preceq_{\text{sto}} \mu(t)$. Recall that $\mu_i = C(-c_i)^{-1} \mu(-c_i)$ for some $c_i \in \mathbb{R}$. As a result, the probability measures $\mu_1, \dots, \mu_n$ are totally ordered (in stochastic order), and are in the same order as their barycenters $x_i = \int y d\mu_i(y)$. Hence, $(\mu_i)_{1 \le i \le n}$ is increasing in stochastic order. The property of SST follows readily from Proposition 2.4. $\square$

**3.3. Sinkhorn's algorithm and convergence.** In this subsection, we provide a numerical scheme to approximate the minimizer of the entropy in the football model with the stochastic ordering constraint (and the SST) as in Theorem 3.5. This scheme is closely related to the one for the entropic optimal transport with discrete marginals, which was initiated in the pioneering work of Benamou, Carlier, Cuturi, Nenna and Peyré [14]. Our situation is different due to the martingale constraint. Thus, we need to adapt the existing theory carefully, depending on whether $x$ is irreducible or not (see Subsection 3.4).$^6$ Our presentation closely follows the one by [9].

Set $C := C_1 \cap C_2 \cap C_3 \in \mathbb{R}^{n \times n}$, where

$$ C_1 := \left\{ M \in \mathbb{R}_+^{n \times n} : \sum_{j=1}^{n} m_{ij} = 1, \text{ for all } i = 1, \dots, n \right\}, \qquad (3.3) $$

$$ C_2 := \left\{ M \in \mathbb{R}_+^{n \times n} : \sum_{i=1}^{n} m_{ij} = 1, \text{ for all } j = 1, \dots, n \right\}, \qquad (3.4) $$

$$ C_3 := \left\{ M \in \mathbb{R}_+^{n \times n} : \sum_{j=1}^{n} (j-1)m_{ij} = x_i, \text{ for all } i = 1, \dots, n \right\}. \qquad (3.5) $$

<footnote>$^6$We point out that some similar but different algorithms have been considered for martingale optimization problems satisfying either marginal or expectation constraints [5, 26].</footnote>
By Proposition 3.2, $C \cap (0, \infty)^{n \times n} \neq \emptyset$.

Let $h: \mathbb{R} \to (-\infty, \infty]$ be defined by $h(x) = x \log x$ for $x \ge 0$, and $\infty$ otherwise. The goal is to solve the convex optimization problem:

$$ \min_{M \in C} \left\{ H(M) := \sum_{i,j=1}^{n} h(m_{ij}) \right\}. \qquad (3.6) $$

To this end, we use Bregman's (pseudo-)distance $D: \mathbb{R}^{n \times n} \times \text{idom}(H) \to [0, \infty]$ defined by

$$ \begin{aligned} D(L, M) &:= H(L) - H(M) - \nabla H(M) \cdot (L - M) \\ &= \sum_{i,j=1}^{n} \left[ l_{ij} \log l_{ij} - m_{ij} \log m_{ij} - (1 + \log m_{ij})(l_{ij} - m_{ij}) \right], \end{aligned} \qquad (3.7) $$

where $\text{dom}(H) = \mathbb{R}_+^{n \times n}$ denotes the domain of $H$, and its interior and boundary are denoted by $\text{idom}(H) = (0, \infty)^{n \times n}$ and $\text{bdom}(H) = \{M \in \mathbb{R}_+^{n \times n} : m_{ij} = 0 \text{ for some } (i, j)\}$, respectively.

Set $M^0 := (m_{ij}^0 \equiv 1) \in \text{idom}(H)^7$, so $D(L, M^0) = H(L) - \sum_{i,j=1}^n l_{ij} + K$, where $K \in \mathbb{R}$ depends only on $M^0$. Because $\sum_{i,j=1}^n l_{ij} = n$ for $L \in C$, we have:

$$ \arg \min_{L \in C} H(L) = \arg \min_{L \in C} D(L, M^0) \quad \text{and} \quad \min_{L \in C} H(L) = \min_{L \in C} D(L, M^0) + n - K. $$

For $k \ge 4$, let $C_k := C_k \mod 3$. For $k \ge 1$, define $M_k$ as the Bregman projection of $M_{k-1}$ on $C_k$:

$$ M^k := \arg \min_{L \in C_k} D(L, M^{k-1}), \qquad (3.8) $$

where $M^k$ is uniquely determined. The following theorem establishes the convergence of the sequence $(M^k)_{k \ge 0}$, where $M^0$ is given previously and $M^k$ is uniquely defined. Finally, the convergence of our scheme is summarized as below.

**Theorem 3.6.** Let $\mathbf{x} \preceq (0, \dots, n-1)$ be an irreducible score, so that $C \cap (0, \infty)^{n \times n} \neq \emptyset$. Then

$$ \lim_{k \to \infty} M^k = \arg \min_{L \in C} D(L, M^0) = \arg \min_{L \in C} H(L). \qquad (3.9) $$

*Proof.* Note that $\text{dom}(h) = \mathbb{R}_+$, $\text{idom}(h) = (0, \infty)$ and $\text{bdom}(h) = \{0\}$. It is easy to check:

(1) $h$ is a proper convex function. Moreover, $h$ is closed because $\{x \in \text{dom}(h) : h(x) \le \alpha\}$ is closed for all $\alpha \in \mathbb{R}$.

(2) $h$ is Legendre because

* $h$ is differentiable on $\text{idom}(h)$;

* $\lim_{t \to 0+} h'(x + t(y - x))(y - x) = -\infty$ for all $x \in \text{bdom}(f)$ and $y \in \text{idom}(h)$;

* $h$ is strictly convex on $\text{idom}(h)$.

(3) $h$ is co-finite because $\lim_{t \to \infty} \frac{h(tx)}{t} = \infty$ for all $x \ne 0$.

(4) $h$ is very strictly convex because $h''(x) > 0$ for all $x \in \text{idom}(h)$.

Since $C_1, C_2, C_3$ are all affine subsets, it suffices to apply [9, Theorem 4.3] to conclude. $\square$

<footnote>$^7$The initialization $M^0 := (m_{ij}^0 \equiv 1) \in \text{idom}(H)$ is chosen so that the entropy optimization problem $\min_{L \in C} H(L)$ is equivalent to the problem $\min_{L \in C} D(L, M^0)$, which can be provably solved by Sinkhorn's algorithm.</footnote>
See also [57] for recent developments on the convergence rate of Bregman’s iteration under further technical assumptions, which we do not pursue here.

Next we propose a computational scheme inspired by Theorem 3.6. The key is to compute numerically, for each $M \in (0, \infty)^{n \times n}$, its Bregman projection on $C_k$. We distinguish three cases:

* $k=1$: We introduce the Lagrange multiplier $\lambda = (\lambda_1, \dots, \lambda_n) \in \mathbb{R}^n$, and set $\Phi(L, \lambda) := D(L, M) + \sum_{i=1}^n \lambda_i (\sum_{j=1}^n l_{ij} - 1)$. Differentiating $\Phi$ with respect to $m_{ij}$ and $\lambda_i$ yields $\partial_{l_{ij}} \Phi = \log(l_{ij}) - \log(m_{ij}) + \lambda_i$ and $\partial_{\lambda_i} \Phi = \sum_{j=1}^n l_{ij} - 1$. By setting these to zero, we get

$$ \arg \min_{L \in C_1} D(L, M) = \left( \frac{m_{ij}}{\sum_{k=1}^{n} m_{ik}} \right)_{1 \le i, j \le n} . \qquad (3.10) $$

* $k=2$: The same reasoning as in the previous case yields:

$$ \arg \min_{L \in C_2} D(L, M) = \left( \frac{m_{ij}}{\sum_{k=1}^{n} m_{ik}} \right)_{1 \le i, j \le n} . \qquad (3.11) $$

* $k=3$: Define $\Phi(L, \lambda) := D(L, M) + \sum_{i=1}^n \lambda_i (\sum_{j=1}^n (j-1)l_{ij} - x_i)$, and differentiate $\Phi$ with respect to $m_{ij}$ and $\lambda_i$ yields $\partial_{l_{ij}} \Phi = \log(l_{ij}) - \log(m_{ij}) + \lambda_i(j-1)$ and $\partial_{\lambda_i} \Phi = \sum_{j=1}^n (j-1)l_{ij} - x_i$. Setting these to zero yields

$$ \arg \min_{L \in C_3} D(L, M) = \left( m_{ij} r_{ij}^{j-1} \right)_{1 \le i, j \le n}, \qquad (3.12) $$

where $r_{ij}$ is the unique positive root to the polynomial equation $\sum_{j=1}^n (j-1)m_{ij}r_{ij}^{j-1} - x_i = 0$. Let $f(r) := \sum_{j=1}^n (j-1)m_{ij}r_{ij}^{j-1} - x_i$. It is easy to see that $f$ is strictly increasing on $[0, \infty)$ with $f(0) \le 0$ and $\lim_{r \to \infty} f(r) = +\infty$. Thus, it is easy (and quick) to find a numerical root of $f$ on $[0, \infty)$ by Newton's method.

To summarize, we normalize the rows ($k=1$), the columns ($k=2$) and the barycenters of the rows ($k=3$) sequentially. The two first operations correspond to the standard steps in Sinkhorn's algorithm. The third one is not exactly a barycenter normalization of the rows $\mu_i$, because the total mass of every row changes from 1 to another value in (3.12). From a theoretical viewpoint, we can merge the steps 1 and 3 to a minimization problem in two variables. It can be checked that the minimum is specified by the proper value of $s$ for $\mu(s)$ in (3.2). Nevertheless, the current presentation is more tractable from a numerical viewpoint.

**3.4. The reducible case.** Without loss of generality, we assume $x = \tilde{x}$, i.e., $x_1 \le \dots \le x_n$. The reducible case can be treated similarly on each of the irreducible components, i.e., for score subsequences $(x_{k_{r+1}}, \dots, x_{k_{r+1}}) \preceq (k_r, k_r+1, \dots, k_{r+1}-1)$, where all inequalities in (1.1) are strict except the equality of the two total sums. More precisely, put $k_0 := 0$, and define recursively $k_{r+1}$ whenever $k_r < n$ by

$$ k_{r+1} := \min \left\{ k \in \mathbb{N} : k_r < k \le n \text{ and } \sum_{i=1}^{k} x_i = \frac{k(k-1)}{2} \right\}. $$

Let $R \in \mathbb{N}$ be such that $k_R = n$. Obviously, $R \ge 2$. We see from Definition 3.1 that for the reducible case, there are different leagues with the probability for a team in a lower league
to defeat a team of an upper league being zero. As a result, an element $M \in \Theta_n(\mathbf{x})$ must take the form:

$$M = \begin{pmatrix} M_1 & 0 & \cdots & 0 \\ 0 & M_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & M_R \end{pmatrix},$$

where $M_r$ has $k_r - k_{r-1}$ rows, the number of teams in the $r$th league. Multiplying $M$ on the left with the row vectors $(1, \dots, 1, 0, \dots, 0)$ with $k_r$ entries 1, and on the right with the column vector $(y_1, \dots, y_n) = (0, 1, \dots, n-1)$, permits us to prove that the $M_r$'s are square matrices. Conversely, every block diagonal matrix made of $R$ doubly stochastic matrices that satisfy

$$\sum_{j=1}^{k_r - k_{r-1}} m_{ij} (k_{r-1} + j - 1) = x_{k_{r-1}+1} \quad \text{for each } i \le k_r - k_{r-1},$$

is clearly an element of $\Theta_n(\mathbf{x})$. Also note that the constraints in the definition of $\Theta_n(\mathbf{x})$ are separable.<fnref n="8" />

Now we see that in the minimization problem, the function $h : m \mapsto m \log m$ is zero on the off-diagonal blocks, and $H(M) = \sum_{r=1}^R H_r(M_r)$, where each $H_r$ is the entropy function defined on $M_{k_r-k_{r-1}}(\mathbb{R})$. Since the constraints on each block are separable, we have $R$ separate problems, each of which has a unique minimizer.<fnref n="9" /> Numerically, we can just solve the $R$ problems separately, which requires us to adapt the step 3 properly by replacing $j-1$ by the proper $y_j$. However, the $R$ problems can be merged into one with the same steps as in the irreducible case. It suffices to set the initial matrix with ones on the diagonal blocks, and zeros on the off-diagonal blocks.

### 4. THE SHADOW CONSTRUCTION

In this section, we provide another construction of $M \in \Theta_n(\mathbf{x})$ for $\mathbf{x} \preceq (0, \dots, n-1)$ with the stochastic ordering constraint. We present it in the form of a martingale coupling, i.e., a probability measure $\pi^*$ on $\mathbb{R}^2$ with marginals the uniform measures $\{x_1, \dots, x_n\}$ and $\{0, \dots, n-1\}$, and a martingale constraint corresponding to (3.12). Then $M = (\mu_1, \dots, \mu_n) \in \Theta_n(\mathbf{x})$ is defined by $\mu_i(A) = n\pi^*(\{x_i\} \times A)$ for each $1 \le i \le n$ and Borel $A \subset \mathbb{R}$, or $m_{ij} = \pi^*(\{(x_i) \times \{j-1\}\})$ for each $1 \le i, j \le n$.

We start by recalling the terminology in [12]. For any closed set $E \subset \mathbb{R}$, denote by $\mathfrak{M}(E)$ the set of measures on $E$. Next we generalize the convex order (see Theorem 1.5) to the extended convex order, denoted by $\preceq_E$. For $\mu, \nu \in \mathfrak{M}(\mathbb{R})$, we say $\mu \preceq_E \nu$ if and only if for any nonnegative convex function $f : \mathbb{R} \to \mathbb{R}$,

$$\int_{\mathbb{R}} f(x)\mu(dx) \le \int_{\mathbb{R}} f(x)\nu(dx).$$

<footnote n="8">⁸The set $\Theta_n(\mathbf{x})$ is the Cartesian product of factors $\Theta_{k_r-k_{r-1}}(\mathbf{x}_r, \mathbf{y}_r)$ with $\mathbf{x}_k = (x_{k_{r-1}+1}, \dots, x_{k_r})$ and $\mathbf{y}_k = (k_{r-1}+1, \dots, k_r)$, and the appropriate definition concerning the constraints.</footnote>

<footnote n="9">⁹This is a slight extension of Proposition 3.3, which is mentioned in Remark 3.4. The difference is that the vector $\mathbf{y} = (0, 1, \dots, n-1)$ is replaced by $(k_{r-1}+1, \dots, k_r)$.</footnote>
Furthermore, for $\mu, \nu \in \mathfrak{M}(\mathbb{R})$, we say $\mu \preceq_{+} \nu$ if and only if for any nonnegative function $f: \mathbb{R} \to \mathbb{R}$,

$$ \int_{\mathbb{R}} f(x)\mu(dx) \leq \int_{\mathbb{R}} f(x)\nu(dx). $$

Define the notion of shadow: Let $\mu, \nu \in \mathfrak{M}(\mathbb{R})$ such that $\mu \preceq_E \nu$. Then there exists a measure $S^\nu(\mu)$, called the shadow of $\mu$ in $\nu$, such that

(1) $S^\nu(\mu) \preceq_{+} \nu$.

(2) $\mu \preceq_C S^\nu(\mu)$.

(3) If $\eta$ is another measure satisfying (1) and (2), then we have $S^\nu(\mu) \preceq_C \eta$.

Now we use shadows to construct a martingale coupling $\pi^*$. Denoting $\frac{1}{k} \sum_{i=1}^k \delta_{z_i}$ by $U_{(z_1, \dots, z_k)}$ we first fix

$$ \mu := U_{(x_1, \dots, x_n)}, \quad \nu := U_{(0, \dots, n-1)}, \qquad (4.1) $$

and for each permutation $\sigma$ we denote $(x_{\sigma(k)})_{1 \leq k \leq n}$ by $(x_k^\sigma)_{1 \leq k \leq n}$. By [12, Lemma 4.13 and Example 4.20], we can carry out the following algorithm because we have $\frac{1}{n}\delta_{x_k^\sigma} \preceq_E \nu_{k-1}$ at each step $k$.

(1) Let $\eta_0^\sigma := 0$ and $\nu_0^\sigma := \nu$ (we have $\eta_0^\sigma + \nu_0^\sigma = \nu$, and $\eta_k^\sigma + \nu_k^\sigma = \nu$ for each $1 \leq k \leq n$ with $\eta_0^\sigma \leq \eta_1^\sigma \leq \dots$ and $\nu_0^\sigma \geq \nu_1^\sigma \geq \dots$).

(2) For $k = 1, \dots, n$, define recursively:

$$ \eta_k^\sigma := \eta_{k-1}^\sigma + S^{\nu_{k-1}^\sigma} \left( \frac{1}{n} \delta_{x_k^\sigma} \right), \quad \nu_k^\sigma := \nu_{k-1}^\sigma - S^{\nu_{k-1}^\sigma} \left( \frac{1}{n} \delta_{x_k^\sigma} \right). \qquad (4.2) $$

(In particular, $\eta_k^\sigma + \nu_k^\sigma = \nu$ for each $1 \leq k \leq n$.)

(3) Define $\pi^\sigma \in \mathcal{P}(\mathbb{R}^2)$ (the set of probability measures on $\mathbb{R}^2$) by

$$ \pi^\sigma(dx, dy) := \sum_{k=1}^{n} \delta_{x_k^\sigma}(dx) \otimes S^{\nu_{k-1}^\sigma} \left( \frac{1}{n} \delta_{x_k^\sigma} \right) (dy). \qquad (4.3) $$

(4) Define $\pi^* \in \mathcal{P}(\mathbb{R}^2)$ by

$$ \pi^* := \frac{1}{n!} \sum_{\sigma \in \mathfrak{S}(n)} \pi^\sigma. \qquad (4.4) $$

Recall that the output measure $\pi^*$ can be used to construct an element of $\Theta_n(\mathbf{x})$ by defining $\mu_i(A) = n\pi^*(\{x_i\} \times A)$, or equivalently, $m_{ij} = n\pi^*(\{x_i\} \times \{j-1\})$.

Figure 1 illustrates the algorithm and the associated shadows for a small example ($n=3$).

**Remark 4.1.** Here we give a probabilistic interpretation of the above algorithm. If we want to know how many goals are scored by team $i$, we pick uniformly at random a sequence $x_{\sigma(1)}, x_{\sigma(2)}, \dots$ until $\sigma(k) = i$. During this process, the measure $\nu$ is filled with

$$ \eta_{k-1}^\sigma = S^\nu \left( \frac{1}{n} \delta_{x_{\sigma(1)}} \right) + S^{\nu_1^\sigma} \left( \frac{1}{n} \delta_{x_{\sigma(2)}} \right) + \dots + S^{\nu_{k-2}^\sigma} \left( \frac{1}{n} \delta_{x_{\sigma(k-1)}} \right) \preceq_{+} \nu, \qquad (4.5) $$

and it remains $\nu_{k-1}^\sigma = \nu - \eta_{k-1}^\sigma$. Now we embed $\frac{1}{n}\delta_{x_{\sigma(k)}}$ in $\nu_{k-1}^\sigma$, and obtain (up to a factor of $1/n$) a probability distribution that underlies the random number of goals. Finally, note that this conditional distribution is obtained by conditioning on $\sigma(1), \dots, \sigma(k)$ rather than the entire $\sigma$.