###### Abstract

The Sinkhorn algorithm is the most popular method for solving the entropy
minimization problem called the Schrödinger problem: in the non-degenerate
cases, the latter admits a unique solution towards which the algorithm
converges linearly. Here, motivated by recent applications of the Schrödinger
problem with respect to structured stochastic processes (such as increasing
ones), we study the Sinkhorn algorithm in degenerate cases where it might
happen that no solution exist at all. We show that in this case, the algorithm
ultimately alternates between two limit points. Moreover, these limit points
can be used to compute the solution of a relaxed version of the Schrödinger
problem, which appears as the $\Gamma$-limit of a problem where the marginal
constraints are replaced by asymptotically large marginal penalizations,
exactly in the spirit of the so-called unbalanced optimal transport. Finally,
our work focuses on the support of the solution of the relaxed problem, giving
its typical shape and designing a procedure to compute it quickly. We showcase
promising numerical applications related to a model used in cell biology.

CONVERGENCE OF THE SINKHORN ALGORITHM WHEN THE SCHRÖDINGER PROBLEM HAS NO
SOLUTION

Authors: Aymeric Baradat1, Elias Ventre2.
1 - Université Claude Bernard Lyon 1, CNRS UMR 5208, Institut Camille Jordan,
Villeurbanne, France.
2 - University of British Columbia, Vancouver, BC, Canada.
Corresponding author<EMAIL_ADDRESS>
Keywords: Schrödinger problem, the Sinkhorn algorithm, matrix scaling

MSC code: 65F35, 65B99, 92-08

##  1 Introduction

The Schrödinger problem has been introduced by Schrödinger himself in the 30’s
[27, 28] in the context of statistical mechanics. It is one of these problems
in mathematics for which there is periodically a resurgence of interest, as
witnessed by the numerous works which it was the object of for almost 100
years, among which [10, 7, 34, 9, 24]. The last of these resurgences, over the
past twenty years, occurred because of its close links with optimal transport.
On the one hand, the Schrödinger problem, that comes with a temperature
parameter in its classical formulation, converges towards optimal transport
when the temperature goes to zero [21, 17, 18]. On the other hand, it is often
much easier to compute the solutions of the Schrödinger problem than the ones
of the optimal transport problem [8, 23], thanks to the so-called Sinkhorn
algorithm [29]. This algorithm converges exponentially fast (_i.e._ at a
linear rate, following the usual terminology in the field), at least when it
is applied to a reference matrix whose entries are all below bounded by a
positive number.

It is known in the theory of matrix scaling that when the reference matrix has
nonnegative but possibly cancelling entries, the data in the Schrödinger
problem may be chosen in such a way that the latter admits no solution. This
is the so-called _non-scalable_ case. Also, when the data are located at the
boundary of those for which there is a solution, the so-called _approximately
scalable_ case, the Schrödinger problem has a solution but the convergence of
the Sinkhorn algorithm is not linear anymore.

In this paper, we want to study the Sinkhorn algorithm in the degenerate case
where the Schrödinger problem has no solution. Our main finding is that for
such problem, the Sinkhorn algorithm leads to exactly two limit points, each
of them being the solution of a Schrödinger problem with modified data, that
we characterize themselves as solutions of auxiliary optimization problems.
Also, we show that these limit points are related to a problem where the
marginal constraints of the original problem are replaced by marginal
penalizations. Moreover, the Schrödinger problem related to the modified data
is seen to belong to the approximately scalable case in general. We therefore
provide a new outlook on the question of the support of the solution in this
case, allowing to design an approximate method for improving the Sinkhorn
algorithm’s convergence both in the approximately scalable and non-scalable
cases.

For simplicity and because it fits with the context of our numerical
explorations and needs, we decided to work in finite spaces, even though some
of the results might be generalizable.

####  The Schrödinger problem in finite spaces

Let $\mathcal{D}=\\{x_{1},\dots,x_{N}\\}$ and
$\mathcal{F}=\\{y_{1},\dots,y_{M}\\}$ be two nonempty finite spaces and
$R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$ be a nonnegative measure
on $\mathcal{D}\times\mathcal{F}$. Of course, we can identify $R$ with a
matrix $R=(R_{ij})\in\mathbb{R}_{+}^{N\times M}$ by setting
$R_{ij}:=R(\\{(x_{i},y_{j})\\})$. Assuming that $R$ models the coupling
between the initial and final positions of the particles of a large system, we
interpret $R_{ij}$ as the sum of the masses of all the particles being in
$x_{i}$ at the initial time, and in $y_{j}$ at the final time.

Let us choose $\mu\in\mathcal{M}_{+}(\mathcal{D})$ and
$\nu\in\mathcal{M}_{+}(\mathcal{F})$. Once again, we see $\mu=(\mu_{i})$ and
$\nu=(\nu_{j})$ as vectors of $\mathbb{R}_{+}^{N}$ and $\mathbb{R}_{+}^{M}$
respectively.

We call $\Pi(\mu,\nu)$ the subset of
$\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$ consisting of all those
matrices $\bar{R}$ whose row and column sums give $\mu$ and $\nu$
respectively, that is, such that

$\forall
i=1,\dots,N,\quad\sum_{j}\bar{R}_{ij}=\mu_{i},\qquad\mbox{and}\qquad\forall
j=1,\dots M,\quad\sum_{i}\bar{R}_{ij}=\nu_{j}.$

In our interpretation, it means that for the system described by $\bar{R}$,
the sum of the masses of all the particles being in $x_{i}\in\mathcal{D}$ at
the initial time is $\mu_{i}$, and the sum of the masses of all the particles
being in $y_{j}\in\mathcal{F}$ at the final time is $\nu_{j}$. In particular,
for $\Pi(\mu,\nu)$ to be nonempty, $\mu$ and $\nu$ need to share their total
mass.

######  Remark 1.

Let us point out to the readers acquainted with the notations used in the
Optimal Transport literature that calling
$X:\mathcal{D}\times\mathcal{F}\to\mathcal{D}$ and
$Y:\mathcal{D}\times\mathcal{F}\to\mathcal{F}$ the canonical projections and
denoting by # the push forward operation on measures, the measure $\bar{R}$
belongs to $\Pi(\mu,\nu)$ provided $X{}_{\\#}\bar{R}=\mu$ and
$Y{}_{\\#}\bar{R}=\nu$. Actually, we will not use these notations, and prefer
to define $\mu^{\bar{R}}:=X{}_{\\#}\bar{R}$ and
$\nu^{\bar{R}}:=Y{}_{\\#}\bar{R}$, see formula 3.

We call the Schrödinger problem w.r.t. $R$ between $\mu$ and $\nu$ the convex
optimization problem consisting in minimizing among $\Pi(\mu,\nu)$ the
relative entropy w.r.t $R$:

$\operatorname{Sch}(R;\mu,\nu):=\min\Big{\\{}H(\bar{R}|R)\,\Big{|}\,\bar{R}\in\Pi(\mu,\nu)\Big{\\}},$

where for all $\bar{R}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$, the
relative entropy of $P$ w.r.t. $R$ is defined by

$H(\bar{R}|R):=\sum_{ij}\Big{\\{}\bar{R}_{ij}\log\frac{\bar{R}_{ij}}{R_{ij}}+R_{ij}-\bar{R}_{ij}\Big{\\}},$

taking the conventions $a\log\frac{a}{0}=+\infty$ if $a>0$, and $0\log
0=0\log\frac{0}{0}=0$. Notice that if $H(\bar{R}|R)<+\infty$, then for all
$i,j$ such that $R_{ij}=0$, we also have $\bar{R}_{ij}=0$, _i.e._ ,
$\bar{R}\ll R$ in the sense of measures.

######  Remark 2.

Of course, as before, for a solution $R^{*}$ to exist, $\mu$ and $\nu$ need to
have the same total mass, and $R^{*}$ will then have the same total mass as
$\mu$ and $\nu$.

By strict convexity of the relative entropy as a function of $\bar{R}$, when
there is a solution, the latter is unique. Also, the relative entropy being
lower semicontinuous w.r.t. $\bar{R}$ and $\Pi(\mu,\nu)$ being compact, the
existence of a solution $R^{*}$ for $\operatorname{Sch}(R;\mu,\nu)$ is
equivalent to the existence of a $\bar{R}\in\Pi(\mu,\nu)$ satisfying
$H(\bar{R}|R)<+\infty$. In what follows, such an $\bar{R}$ is called a
competitor for $\operatorname{Sch}(R;\mu,\nu)$.

Heuristically, we seek for the measure $R^{*}$ that is the closest possible to
$R$ in the entropic sense while imposing its first and second marginals.

In virtue of the Sanov theorem [25], this problem has an interpretation in
terms of large deviations. It is also known to be connected to optimal
transport problems, see [17, 18, 11, 5]: if for all $i,j$, $c_{ij}$ models the
cost to transport a unit of mass from $x_{i}$ to $y_{j}$, and
$R_{ij}\propto\exp(-c_{ij}/\varepsilon)$ for some small $\varepsilon>0$, then
the solution of $\operatorname{Sch}(R;\mu,\nu)$ is a good approximation of a
solution of the optimal transport problem between $\mu$ and $\nu$, of cost
$(c_{ij})$.

####  The Sinkhorn algorithm

When the solution of $\operatorname{Sch}(R;\mu,\nu)$ exists, it is well known
for a very long time that this solution turns out to be the limit of the
sequences $(P^{n})_{n\in\mathbb{N}^{*}}$ and $(Q^{n})_{n\in\mathbb{N}^{*}}$
appearing in the following so-called the Sinkhorn algorithm, also called IPFP
for _iterative proportional fitting procedure_ [29, 30, 14, 22]:

$\left\\{\begin{gathered}Q^{0}:=R,\\\ \begin{aligned} \forall n\geq 0,\quad
P^{n+1}&:=\operatorname*{arg\,min}\Big{\\{}H(P|Q^{n}),&&X{}_{\\#}P=\mu\Big{\\}},\\\
\forall n\geq 0,\quad
Q^{n+1}&:=\operatorname*{arg\,min}\Big{\\{}H(Q|P^{n+1}),&&Y{}_{\\#}Q=\nu\Big{\\}}.\end{aligned}\end{gathered}\right.$
(1)

This formulation is implicit as it involves minimization problems. In fact,
easy results concerning these problems, detailed in Corollary 6 below, give
access to an explicit and easily computable version, which takes the following
form, when expressed in terms of the so called _dual variables_ or
_potentials_
$(a^{n})_{n\in\mathbb{N}^{*}}\in(\mathbb{R}_{+}^{N})^{\mathbb{N}}$ and
$(b^{n})_{n\in\mathbb{N}}\in(\mathbb{R}_{+}^{M})^{\mathbb{N}}$:

$\left\\{\begin{gathered}\forall j,\quad b^{0}_{j}:=1,\\\\[10.0pt]
\begin{aligned} &\forall n\geq 0,\quad\forall
i,j,&a^{n+1}_{i}&:=\frac{\mu_{i}}{\displaystyle\sum_{j^{\prime}}b^{n}_{j^{\prime}}R_{ij^{\prime}}},&P^{n+1}_{ij}&:=a^{n+1}_{i}b^{n}_{j}R_{ij},\\\
&\forall n\geq 0,\quad\forall
i,j,&b^{n+1}_{j}&:=\frac{\nu_{j}}{\displaystyle\sum_{i^{\prime}}a^{n+1}_{i^{\prime}}R_{i^{\prime}j}},&Q_{ij}^{n+1}&:=a^{n+1}_{i}b^{n+1}_{j}R_{ij}.\end{aligned}\end{gathered}\right.$
(2)

A reason for the popularity of this algorithm is that in a lot of contexts,
the sequences of potentials $(a^{n})$ and $(b^{n})$, and hence the sequence of
couplings $(P^{n})$ and $(Q^{n})$ converge at a linear rate, and the limit of
$(P^{n})$ and $(Q^{n})$ coincide with the unique solution of
$\operatorname{Sch}(R;\mu,\nu)$. For this reason, the Sinkhorn algorithm is
nowadays the most efficient way to compute approximate solutions of optimal
transport problems [8, 2, 23].

Observe that _a priori_ , the existence of a solution for the Schrödinger
problem is not necessary to give a meaning to the Sinkhorn algorithm.
Actually, we will see that there are lots of situations where the Schrödinger
problem has no solution, and yet the Sinkhorn algorithm is perfectly well
defined. These are the cases that we want to study in this text.

####  A degenerate case

As we just said, our aim is to study the Sinkhorn algorithm in the cases where
the existence of a solution of the Schrödinger problem is either false, or at
least nontrivial. This includes the case where $\mu$ and $\nu$ do not have the
same total mass, see Remark 2. However, this is not the main new situation
that we want to encompass, since the Sinkhorn algorithm behaves trivially
under normalization. More interestingly, we will give a detailed study of the
case where some entries of $R$ cancel, or in optimal transport terms, when the
cost function takes the value $+\infty$.

In that situation, it can be hard to exhibit a competitor, since the natural
candidate that is usually chosen, namely, the product measure of $\mu$ and
$\nu$, is not absolutely continuous w.r.t. $R$ in general. In fact, there are
cases where it is easy to see that no competitor exists. We give in Appendix A
an explicit and simple example of such a case. To illustrate our findings, we
also describe the behaviour of the Sinkhorn algorithm applied to this example.

Note that beyond the theoretical interest, there are practical motivations for
studying cases where the problem has no solution. Indeed, the Schrödinger
problem can be used as follows. Suppose that $\mu$ and $\nu$ are some observed
densities of a random phenomenon at two different timepoints, obtained for
instance by building the empirical distributions associated to some collected
data. Suppose also that we have at our disposal a good model for this
phenomenon, that is, a reference stochastic process chosen based on our
knowledge of the system prior to the observations of $\mu$ and $\nu$. Let us
call $R$ the coupling of this process between the two studied timepoints. If
we believe enough in our model and in our data, but still the marginals of $R$
are not $\mu$ and $\nu$, then it is reasonable to try to improve our model by
looking for the coupling that is the closest to $R$ (for instance in the
entropic sense), but which is compatible with the data: this means solving
$\operatorname{Sch}(R;\mu,\nu)$.

Now imagine that there is no lower-bound for the coupling $R$, which can be
perfectly justified (think for instance of a nondecreasing process, like the
size of some randomly growing phenomenon). Then, small measurement errors due
to imprecision of the devices or even to too restricted samplings may result
in the non-existence of any coupling with marginals $\mu$ and $\nu$ being
absolutely continuous with respect to $R$: the Schrödinger problem would thus
have no solution. In that case, we would like to be able to find a coupling
which explains the best the data while being entropically close to $R$. Some
methods are available for doing so, like for instance algorithms solving the
so-called unbalanced problem [6], but at the cost of introducing a new
parameter quantifying the balance between the proximity to the data and to the
reference coupling, whose value will often be arbitrarily chosen. We show in
this article that interestingly, the Sinkhorn algorithm allows to overcome
this choice in the specific situation where the data are more trustworthy than
the model.

In particular, we were motivated by an application of the Sinkhorn algorithm
related to systems biology, and more specifically to the treatment of single-
cell data. The quick progresses of acquisition methods for such data raises
the hope of a better understanding of the cell-differentiation process, which
would in turn pave the way for major medical breakthroughs. In the seminal
papers [26, 16], Schiebinger and his coauthors suggest to analyse the
collected data through an approach based on optimal transport and more
specifically on the Schrödinger problem.

In this field, the unknown is the law of the evolution of the quantity of mRNA
molecules in the cells through time: this evolution cannot be followed, as our
techniques of measurement destroy the cells. Hence, to study it between two
timepoints, the approach consists in:

  1. (i)

choosing a reference theoretical model $R$, where for all $i,j$, $R_{ij}$ is
the expected quantity of cells whose mRNA levels are given by the vector
$x_{i}$ at the initial time, and by $y_{j}$ at the final time;

  2. (ii)

measuring the mRNA levels of samples of cells at the initial and final times
to get approximate distributions $\mu$ and $\nu$ of these levels among the
population of cells under study;

  3. (iii)

solving the Schrödinger problem $\operatorname{Sch}(R;\mu,\nu)$ to get a law
$R^{*}$ that is close to our model $R$, but which explains the data.

In the case of Schiebinger, $R$ is the coupling produced by a Brownian motion
between two time points, and therefore admits a below bound. In a separated
work [31], the second author argues that a more realistic model would be
obtained by replacing the Brownian motion by a _piecewise deterministic Markov
process_ as described in [13]. For such models, dynamical constraints
involving mRNAs half-life times lead to a degenerate $R$ and the corresponding
Schrödinger problem could thus have no solution, not because of a lack in the
model, but because of inaccuracies in the measurements. Our results show that
the Sinkhorn algorithm can still be used in this situation, without any pre-
treatment of the data. We refer once again to Appendix A for a further
discussion on this topic.

####  Contributions

In this article, we work with a potentially degenerate $R$, and our main
contributions are the following.

  * •

We show that the two sequences $(P^{n})_{n\in\mathbb{N}^{*}}$ and
$(Q^{n})_{n\in\mathbb{N}^{*}}$ defined in (1) converge towards two possibly
different matrices $P^{*}$ and $Q^{*}$, each of them being the solution of a
Schrödinger problem with modified marginals. More precisely, the matrice
$P^{*}$ is the solution of the problem $\operatorname{Sch}(R;\mu,\nu^{*})$,
where $\nu^{*}$ minimizes the relative entropy w.r.t. $\nu$ within the set of
marginals $\bar{\nu}$ for which the Schrödinger problem
$\operatorname{Sch}(R;\mu,\bar{\nu})$ admits a solution, and a similar
statement holds for $Q^{*}$. This result, stated at Theorem 11, is the main
result of Section 3.

  * •

We show in Section 4 that the Sinkhorn algorithm enables to compute the
solution of a modified Schrödinger problem where the marginal _constraints_
are replaced by marginal _penalizations_ : as shown at Theorem 17, the limit
of the solution of the problem

$\min\Big{\\{}H(\bar{R}|R)+\lambda\big{(}H(\mu^{\bar{R}}|\mu)+H(\nu^{\bar{R}}|\nu)\big{)}\,\Big{|}\,\bar{R}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})\Big{\\}},$

(where once again, $\mu^{\bar{R}}$ and $\nu^{\bar{R}}$ are the first and
second marginal of $\bar{R}$, see Remark 1) converges towards the
componentwise geometric mean of the two limits $P^{*}$ and $Q^{*}$ of the
Sinkhorn algorithm as $\lambda\to+\infty$.

  * •

In Section 5, we recall a well known necessary and sufficient condition on
$R$, $\mu$ and $\nu$ for $\operatorname{Sch}(R;\mu,\nu)$ to admit a solution.
Using this condition, we develop at Proposition 28 a procedure to find the
(common) support $\mathcal{S}$ of $P^{*}$ and $Q^{*}$ without computing them.
As explained in Subsection 5.2, our motivation is that the convergence rate of
the Sinkhorn algorithm is linear if and only if $\mathcal{S}$ coincides with
the support of $R$. When it is not the case, as often when
$\operatorname{Sch}(R;\mu,\nu)$ has no solution, we can therefore improve the
speed of convergence by first computing $\mathcal{S}$, and then by applying
the Sinkhorn algorithm to
$\operatorname{Sch}(\mathds{1}_{\mathcal{S}}R;\mu,\nu)$ instead of
$\operatorname{Sch}(R;\mu,\nu)$, which does not change the limits $P^{*}$ and
$Q^{*}$.

  * •

Section 6 is an application of the developments made at Section 5. We
implement an approximate but fast algorithm, usable in practice, allowing to
recover an estimate of the support $\mathcal{S}$. We then compare the Sinkhorn
algorithm and the technique coming from [6] with our method consisting in
first computing $\mathcal{S}$ with our approximate algorithm and then applying
the Sinkhorn algorithm to
$\operatorname{Sch}(\mathds{1}_{\mathcal{S}}R;\mu,\nu)$. We also detail the
regimes in which our method is a significant improvement of the other
techniques.

Some of the results of this paper can be generalized by replacing
$\mathcal{D}$ and $\mathcal{F}$ by general Polish spaces without much effort.
This is the reason why we will often write $H(P|R)<+\infty$ instead of $P\ll
R$: these are equivalent in the finite case, but not in the continuous one. In
the latter case, we often need the stronger entropic assumption. Even if we
decided to stick to the finite case in order to stress the key arguments that
make everything work in practice, we believe that the continuous case is also
interesting, and we wish to study it in a further work.

Before coming up with our contributions, we recall a few facts about the
relative entropy functional at Section 2.

##  2 Notations, properties of the entropy and terminology

In this preliminary section, we introduce some notations, provide well known
elementary results concerning the entropy, and recall the terminology usually
used in the theory of matrix scaling.

###  2.1 Notations

Let us first give a few notations that will be used systematically in this
work. Most of them were already given in the introduction.

  * •

Whenever $I$ is a finite set of labels and $\mathcal{E}=\\{u_{k},\ k\in I\\}$
is a finite set indexed by $I$, we denote by $\mathcal{M}_{+}(\mathcal{E})$
the set of nonnegative measures on $\mathcal{E}$. This set is identified with
$\mathbb{R}_{+}^{I}$ through the the correspondence

$\mathsf{r}\in\mathcal{M}_{+}(\mathcal{E})\leftrightsquigarrow(\mathsf{r}_{k}:=\mathsf{r}(\\{u_{k}\\}))_{k\in
I}\in\mathbb{R}^{I}.$

For all $\mathsf{r}\in\mathcal{M}_{+}(\mathcal{E})$, we denote by
$\mathsf{M}(\mathsf{r}):=\sum_{k}\mathsf{r}_{k}$ its total mass. If
$\mathsf{M}(\mathsf{r})=1$, we say that $\mathsf{r}$ is a probability measure
on $\mathcal{E}$, and we write $\mathsf{r}\in\mathcal{P}(\mathcal{E})$. The
topology considered on $\mathcal{M}_{+}(\mathcal{E})$ is nothing but the one
of $\mathbb{R}^{I}$.

  * •

In the same way, we identify the set $\mathcal{F}(\mathcal{E};\mathbb{R})$ of
real functions $Z$ on $\mathcal{E}$ with $\mathbb{R}^{I}$ through the
correspondence

$Z\in\mathcal{F}(\mathcal{E};\mathbb{R})\leftrightsquigarrow(Z_{k}:=Z(u_{k}))_{k\in
I}\in\mathbb{R}_{+}^{I}.$

Depending on the context, we will either call such functions $Z$ _test
functions_ , or _random variables_ , thinking of $\mathcal{E}$ as a measurable
set. The random variables that we will consider will actually often be
slightly more general, and be allowed to take the value $-\infty$, in which
case we will tell it explicitly.

  * •

Through our identifications, the duality between
$\mathcal{M}_{+}(\mathcal{E})$ and $\mathcal{F}(\mathcal{E};\mathbb{R})$ is
nothing but the usual scalar product on $\mathbb{R}^{I}$, and denoted for all
$Z\in\mathcal{F}(\mathcal{E};\mathbb{R})$ and
$\mathsf{r}\in\mathcal{M}_{+}(\mathcal{E})$ by

$\langle Z,\mathsf{r}\rangle:=\sum_{k}Z_{k}\mathsf{r}_{k}.$

When $Z$ possibly takes the value $-\infty$, we always choose by convention
$-\infty\times 0=0$.

  * •

In the context of the introduction, when $\mathcal{D}=\\{x_{1},\dots,x_{N}\\}$
and $\mathcal{F}=\\{y_{1},\dots,y_{M}\\}$ are two nonempty finite spaces and
$\mathcal{E}=\mathcal{D}\times\mathcal{F}$, then the corresponding $I$ is the
product space $\\{1,\dots,N\\}\times\\{1,\dots,M\\}$, and
$\bar{R}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$, is seen as a
matrix. We define its marginals $\mu^{\bar{R}}\in\mathcal{M}_{+}(\mathcal{D})$
and $\nu^{\bar{R}}\in\mathcal{M}_{+}(\mathcal{F})$ by the formulas

$\forall
i=1,\dots,N,\quad\mu^{\bar{R}}_{i}:=\sum_{j}\bar{R}_{ij},\quad\mbox{and}\quad\forall
j=1,\dots M,\quad\nu^{\bar{R}}_{j}:=\sum_{i}\bar{R}_{ij}.$ (3)

Of course, $\mu^{\bar{R}}$ and $\nu^{\bar{R}}$ have the same total mass as
$\bar{R}$, that is:

$\mathsf{M}(\bar{R})=\mathsf{M}(\mu^{\bar{R}})=\mathsf{M}(\nu^{\bar{R}}).$ (4)

In particular, if $R$ is a probability measure, its marginals are probability
measures as well.

  * •

As before, if $\mu\in\mathcal{M}_{+}(\mathcal{D})$ and
$\nu\in\mathcal{M}_{+}(\mathcal{F})$, we call $\Pi(\mu,\nu)$ the set of those
$R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$ such that $\mu^{R}=\mu$
and $\nu^{R}=\nu$.

  * •

For the sake of simplicity, we do not use different notations for the same
functions applied in different context. For instance, notations for the total
mass $\mathsf{M}$ or the relative entropy $H$ (see Definition 3 below) might
be applied to different sets $\mathcal{E}$ namely $\mathcal{D}$, $\mathcal{F}$
and $\mathcal{D}\times\mathcal{F}$.

###  2.2 First properties of the relative entropy

This subsection only contains easy and very well known results concerning the
relative entropy that will be useful in the sequel. We stick to the finite
case as this is the one studied in this paper, and we provide some proofs for
the readers who are not acquainted with this notion of entropy, but all the
properties given here are known in a much wider context, see for instance
[19].

As already said in the introduction, the relative entropy is defined as
follows.

######  Definition 3.

Let $\mathcal{E}=\\{u_{k},\ k\in I\\}$ be a finite set and
$\mathsf{r}=(\mathsf{r}_{k})\in\mathcal{M}_{+}(\mathcal{E})$. For all
$\bar{\mathsf{r}}=(\bar{\mathsf{r}}_{k})\in\mathcal{M}_{+}(\mathcal{E})$, the
relative entropy of $\bar{\mathsf{r}}$ w.r.t $\mathsf{r}$ is the value in
$[0,+\infty]$ given by

$H(\bar{\mathsf{r}}|\mathsf{r}):=\sum_{k}\Big{\\{}\bar{\mathsf{r}}_{k}\log\frac{\bar{\mathsf{r}}_{k}}{\mathsf{r}_{k}}+\mathsf{r}_{k}-\bar{\mathsf{r}}_{k}\Big{\\}}=\sum_{k}\bar{\mathsf{r}}_{k}\log\frac{\bar{\mathsf{r}}_{k}}{\mathsf{r}_{k}}+\mathsf{M}(\mathsf{r})-\mathsf{M}(\bar{\mathsf{r}}),$

with convention $a\log\frac{a}{0}=+\infty$ for all $a>0$, and $0\log
0=0\log\frac{0}{0}=0$.

First, this definition provides a convex function with good continuity
properties. We state them in the following proposition, for which we omit the
straightforward proof.

######  Proposition 4.

Let $\mathcal{E}$ be a finite set and
$\mathsf{r}\in\mathcal{M}_{+}(\mathcal{E})$. The functional

$\bar{\mathsf{r}}\in\mathcal{M}_{+}(\mathcal{E})\mapsto
H(\bar{\mathsf{r}}|\mathsf{r})\in[0,+\infty]$

is strictly convex, lower semicontinuous, and continuous on its domain, which
is the closed set
$\\{\bar{\mathsf{r}}\ll\mathsf{r}\\}\subset\mathcal{M}_{+}(\mathcal{E})$.

For a given $\bar{\mathsf{r}}\in\mathcal{M}_{+}(\mathcal{E})$, the functional

$\mathsf{r}\in\mathcal{M}_{+}(\mathcal{E})\mapsto
H(\bar{\mathsf{r}}|\mathsf{r})\in[0,+\infty]$

is convex and continuous for the canonical topology of $[0,+\infty]$. Its
domain is the open set
$\\{\mathsf{r}\gg\bar{\mathsf{r}}\\}\subset\mathcal{M}_{+}(\mathcal{E})$.

The most useful property of the relative entropy is the computation of its
Legendre transform. This property can be stated as follows.

######  Theorem 5.

Let $\mathcal{E}=\\{u_{k},\ k\in I\\}$ be a finite set, and
$\mathsf{r}\in\mathcal{M}_{+}(\mathcal{E})$. For all test function $Z$
possibly taking the value $-\infty$ on $\mathcal{E}$ and all nonnegative
measure $\bar{\mathsf{r}}$ on $\mathcal{E}$, we have

$\langle Z,\bar{\mathsf{r}}\rangle\leq H(\bar{\mathsf{r}}|\mathsf{r})+\langle
e^{Z}-1,\mathsf{r}\rangle,$ (5)

with conventions $e^{-\infty}=0$, $-\infty\times 0=0$ and
$+\infty-\infty=+\infty$.

Moreover, equality in $\mathbb{R}$ holds if and only if
$\bar{\mathsf{r}}\ll\mathsf{r}$ and for all $k\in I$,

$Z_{k}=\log\frac{\bar{\mathsf{r}}_{k}}{\mathsf{r}_{k}}\in[-\infty,+\infty)$
(6)

with convention $\log\frac{0}{a}=-\infty$ for all $a\geq 0$.

###### Proof.

Let $\mathsf{r},\bar{\mathsf{r}}$ and $Z$ be as in the statement of the
theorem. If $H(\bar{\mathsf{r}}|\mathsf{r})=+\infty$, there is nothing to
prove, and we assume $\bar{\mathsf{r}}\ll\mathsf{r}$.

By direct real computations, with the same conventions as in the statement of
the theorem, we find that for all $k\in I$:

$Z_{k}\bar{\mathsf{r}}_{k}\leq\bar{\mathsf{r}}_{k}\log\frac{\bar{\mathsf{r}}_{k}}{\mathsf{r}_{k}}+\mathsf{r}_{k}-\bar{\mathsf{r}}_{k}+\big{(}e^{Z_{k}}-1\big{)}\mathsf{r}_{k},$

with equality if and only if $\mathsf{r}_{k}=\bar{\mathsf{r}}_{k}=0$ or
$\mathsf{r}_{k}>0$ and

$Z_{k}=\log\frac{\bar{\mathsf{r}}_{k}}{\mathsf{r}_{k}}\in[-\infty,+\infty).$

We find (5) and (6) by summing this inequality over $k$. ∎

This theorem will be useful as such, but also implies the following corollary
which gives a full understanding of one step in the Sinkhorn algorithm (1).

######  Corollary 6.

Let $\mathcal{D}$ and $\mathcal{F}$ be two finite sets, and
$\bar{R},R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$. With the
notations of (3), we have

$H(\mu^{\bar{R}}|\mu^{R})\leq H(\bar{R}|R)\qquad\mbox{and}\qquad
H(\nu^{\bar{R}}|\nu^{R})\leq H(\bar{R}|R).$ (7)

In the case where $H(\bar{R}|R)$ is finite, equality holds if and only if for
all $i,j$, respectively:

$\bar{R}_{ij}=\frac{\mu^{\bar{R}}_{i}}{\mu^{R}_{i}}R_{ij}\qquad\mbox{and}\qquad\bar{R}_{ij}=\frac{\nu^{\bar{R}}_{j}}{\nu^{R}_{j}}R_{ij},$

with convention $\frac{0}{0}=0$.

In particular, given $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$ and
$\mu\in\mathcal{M}_{+}(\mathcal{D})$, the problem

$\min\Big{\\{}H(P|R)\,\Big{|}\,\mu^{P}=\mu\Big{\\}}$ (8)

admits a solution if and only if $H(\mu|\mu^{R})<+\infty$, and in this case,
this solution $P$ is unique and satisfies for all $i,j$

$P_{ij}=\frac{\mu_{i}}{\mu^{R}_{i}}R_{ij}$ (9)

with convention $\frac{0}{0}=0$. Moreover, $H(P|R)=H(\mu|\mu^{R})$.

Similarly, given $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$ and
$\nu\in\mathcal{M}_{+}(\mathcal{F})$, the problem

$\min\Big{\\{}H(Q|R)\,\Big{|}\,\nu^{Q}=\nu\Big{\\}}$

admits a solution if and only if $H(\nu|\nu^{R})<+\infty$, and in this case,
this solution $Q$ is unique and satisfies for all $i,j$

$Q_{ij}=\frac{\nu_{j}}{\nu^{R}_{j}}R_{ij}$

with convention $\frac{0}{0}=0$. Moreover, $H(Q|R)=H(\nu|\nu^{R})$.

###### Proof.

The first inequality in (7) is a direct application of (5) with
$\mathsf{r}=R$, $\bar{\mathsf{r}}=\bar{R}$ and for all $i,j$,

$Z_{ij}=\log\frac{\mu^{\bar{R}}_{i}}{\mu^{R}_{i}}.$

The second inequality is proved in the same way, and the equality case is a
consequence of (6).

For the second part of the statement, let us observe that for all $P$
satisfying the constraint in (8), because of (7), $H(P|R)\geq H(\mu|\mu^{R})$,
which – by the equality case – is attained if and only if (9) holds. The
problem involving the second marginal is treated in the same way. ∎

###  2.3 The Schrödinger problem: assumptions and terminology

Let $\mathcal{D}=\\{x_{1},\dots,x_{N}\\}$ and
$\mathcal{F}=\\{y_{1},\dots,y_{M}\\}$ be two nonempty finite sets, and let us
choose a reference measure
$R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$. Given
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$,
the Schrödinger problem, already defined in the introduction, rewrites with
the notations of Subsection 2.1:

$\operatorname{Sch}(R;\mu,\nu):=\min\Big{\\{}H(\bar{R}|R)\,\Big{|}\,\bar{R}\in\mathcal{M}_{+}(\mu,\nu)\mbox{
such that }\mu^{\bar{R}}=\mu\mbox{ and }\nu^{\bar{R}}=\nu\Big{\\}}.$ (10)

######  Remark 7.

Here, we define $\operatorname{Sch}(R;\mu,\nu)$ as the optimal value of our
problem. However, with an abusive terminology, we will refer to the minimizer
of the r.h.s. of (10) as ”the solution of $\operatorname{Sch}(R;\mu,\nu)$”.
More generally, we will call ”the problem $\operatorname{Sch}(R;\mu,\nu)$” the
optimization problem consisting in computing the value
$\operatorname{Sch}(R;\mu,\nu)$.

As we will see in Theorem 11, the Sinkhorn algorithm (1) associated with the
problem $\operatorname{Sch}(R;\mu,\nu)$ is well defined if and only if the
following assumption holds.

######  Assumption 8.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$,
and let us call

$\mathcal{E}:=\Big{\\{}(x_{i},y_{j})\in\mathcal{D}\times\mathcal{F}\mbox{ such
that }R_{ij}>0,\,\mu_{i}>0\mbox{ and }\nu_{j}>0\Big{\\}}.$ (11)

We say that the triple $(R;\mu,\nu)$ satisfies Assumption 8 provided
$R^{0}:=\mathds{1}_{\mathcal{E}}\cdot R$ is such that:

$H(\mu|\mu^{R^{0}})<+\infty\qquad\mbox{and}\qquad H(\nu|\nu^{R^{0}})<+\infty.$
(12)

This assumption is easily seen to be necessary for
$\operatorname{Sch}(R;\mu,\nu)$ to admit a solution. Under Assumption 8 either
$\mathsf{M}(\mu)=\mathsf{M}(\nu)=0$, or none of them is $0$. In the second
case, up to replacing $\mathcal{D}$ by $\mathcal{D}^{\prime}$, the support of
$\mu$, $\mathcal{F}$ by $\mathcal{F}^{\prime}$, the support of $\nu$, and $R$
by its restriction (or equivalently of the one of $R^{0}$) on
$\mathcal{D}^{\prime}\times\mathcal{F}^{\prime}$, we end up with the following
assumption, that will often be used in this paper.

######  Assumption 9.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$.
We say that the triple $(R;\mu,\nu)$ satisfies Assumption 9 provided the
support of $\mu$ and $\mu^{R}$ is $\mathcal{D}$ and the support of $\nu$ and
$\nu^{R}$ is $\mathcal{F}$.

The Schrödinger problem (10) consists in minimizing a convex function under
linear constraints. Therefore, the functional
$(\mu,\nu)\in\mathcal{M}_{+}(\mathcal{D})\times\mathcal{M}_{+}(\mathcal{F})\mapsto\operatorname{Sch}(R;\mu,\nu)\in[0,+\infty]$
is convex.

In the case where Assumption 9 holds, following the usual terminology of the
matrix scaling theory (except for the last item which is more exotic), see
[14], we say that:

  * •

The problem is _scalable_ if $(\mu,\nu)$ is in the relative interior of the
domain of $\operatorname{Sch}(R;\cdot)$. In this case,
$\mathsf{M}(\mu)=\mathsf{M}(\nu)$, the Schrödinger problem admits a unique
solution $R^{*}$, $R^{*}\sim R$ in the sense of measures, and the Sinkhorn
algorithm converges towards $R^{*}$, at a linear rate. In Lemma 24, we recall
an explicit necessary and sufficient condition on $R$, $\mu$, $\nu$ for
$\operatorname{Sch}(R;\mu,\nu)$ to be scalable.

  * •

The problem is _approximately scalable_ if $(\mu,\nu)$ is at the relative
boundary of the domain of $\operatorname{Sch}(R;\cdot)$. In this case,
$\mathsf{M}(\mu)=\mathsf{M}(\nu)$, the Schrödinger problem admits a unique
solution $R^{*}$, and the Sinkhorn algorithm converges towards $R^{*}$.
However, in this case, the support of $R^{*}$ is strictly included in the
support of $R$ (else, we easily see that we are in the scalable case), and the
rate cannot be linear anymore: as proved in [1], a linear rate of convergence
for the Sinkhorn algorithm is not compatible with the appearance of new zero
entries at the limit. We recall at Theorem 23 a necessary and sufficient
condition on $R$, $\mu$ and $\nu$ for $\operatorname{Sch}(R;\mu,\nu)$ to be at
least approximately scalable, that is, either approximately scalable or
scalable.

  * •

The problem is _non-scalable_ if $\mathsf{M}(\mu)=\mathsf{M}(\nu)$, but the
Schrödinger problem $\operatorname{Sch}(R;\mu,\nu)$ does not admit a solution.
This is the case when the condition of Theorem 23 does not hold. This case is
the main case of interest in this work.

  * •

The problem is _unbalanced_ if $\mathsf{M}(\mu)\neq\mathsf{M}(\nu)$. Calling
$\mu^{\prime}:=\mu/\mu(\mathcal{D})$ and $\nu^{\prime}:=\nu/\nu(\mathcal{F})$
their normalized versions, we will say that $\operatorname{Sch}(R;\mu,\nu)$ is
respectively unbalanced scalable, unbalanced approximately scalable and
unbalanced non-scalable whenever
$\operatorname{Sch}(R;\mu^{\prime},\nu^{\prime})$ is scalable, approximately
scalable or non-scalable.

Yet, with an abuse of terminology, we will often refer to the non-scalable
case for results that are true in _any_ situation, including the balanced and
unbalanced non-scalable ones, which are often the most difficult.

##  3 The Sinkhorn algorithm in the non-scalable case

In this section, we consider
$R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$
that we identify respectively with a matrix and two vectors, as before.

The goal of this section is to show that under obvious necessary assumptions,
then the algorithm given in (1) is well defined, and that the sequences
$(P^{n})_{n\in\mathbb{N}^{*}}$ and $(Q^{n})_{n\in\mathbb{N}^{*}}$ that it
provides converge separately towards matrices $P^{*}$ and $Q^{*}$ that we
define now. It will be obvious from their definition that these matrices
coincide if and only if the problem $\operatorname{Sch}(R;\mu,\nu)$ defined in
(10) admits a solution, that is, if it is at least approximately scalable.
Hence our proof recovers the classical fact that the Sinkhorn algorithm
converges towards the solution of the Schrödinger problem as soon as the
latter exists.

The first step to define $P^{*}$ and $Q^{*}$ is to define a pair of new
marginals $\mu^{*}\in\mathcal{M}_{+}(\mathcal{D})$ and
$\nu^{*}\in\mathcal{M}_{+}(\mathcal{F})$ as solutions of the following
optimization problem:

$\begin{gathered}\mu^{*}:=\operatorname*{arg\,min}\Big{\\{}H(\bar{\mu}|\mu)\,\Big{|}\,\bar{\mu}=\mu^{Q}\mbox{
for some }Q\mbox{ with }H(Q|R)<+\infty\mbox{ and }\nu^{Q}=\nu\Big{\\}},\\\
\nu^{*}:=\operatorname*{arg\,min}\Big{\\{}H(\bar{\nu}|\nu)\,\Big{|}\,\bar{\nu}=\nu^{P}\mbox{
for some }P\mbox{ with }H(P|R)<+\infty\mbox{ and
}\mu^{P}=\mu\Big{\\}}.\end{gathered}$ (13)

The question of existence of $\mu^{*}$ and $\nu^{*}$ is treated in Theorem 11
below. Of course, if the problem $\operatorname{Sch}(R;\mu,\nu)$ admits a
competitor, then $\mu^{*}=\mu$ and $\nu^{*}=\nu$.

######  Remark 10.

In the unbalanced case, notice that the total mass of $\nu^{*}$ is the one of
$\mu$, and the total mass of $\mu^{*}$ is the one of $\nu$, that is,
$\mathsf{M}(\nu^{*})=\mathsf{M}(\mu)$ and
$\mathsf{M}(\mu^{*})=\mathsf{M}(\nu)$.

Then $P^{*}$ and $Q^{*}$ are simply defined as the solutions of the
Schrödinger problems $\operatorname{Sch}(R;\mu,\nu^{*})$ and
$\operatorname{Sch}(R;\mu^{*},\nu)$ respectively, that is:

$P^{*}:=\operatorname*{arg\,min}\Big{\\{}H(P|R)\,\Big{|}\,P\in\Pi(\mu,\nu^{*})\Big{\\}}\quad\mbox{and}\quad
Q^{*}:=\operatorname*{arg\,min}\Big{\\{}H(Q|R)\,\Big{|}\,Q\in\Pi(\mu^{*},\nu)\Big{\\}}.$
(14)

Of course, if the problem $\operatorname{Sch}(R;\mu,\nu)$ admits a competitor,
and hence a solution, then both $P^{*}$ and $Q^{*}$ coincide with this
solution.

Our convergence theorem can be stated as follows.

######  Theorem 11.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$
satisfy Assumption 8. The sequences $(P^{n})_{n\in\mathbb{N}^{*}}$ and
$(Q^{n})_{n\in\mathbb{N}^{*}}$ from (1), the marginals $\mu^{*}$ and $\nu^{*}$
from (13) and the matrices $P^{*}$ and $Q^{*}$ from (14) are well defined, and

$P^{n}\underset{n\to+\infty}{\longrightarrow}P^{*}\qquad\mbox{and}\qquad
Q^{n}\underset{n\to+\infty}{\longrightarrow}Q^{*}.$

######  Remark 12.

  * •

Assumption 8 is necessary: it is straightforward to check that if $Q^{1}$ from
(1) is well defined, then $H(Q^{1}|R^{0})<+\infty$. In particular, projecting
on the second marginal, we conclude that $H(\nu|\nu^{R^{0}})<+\infty$. Arguing
in the same way with $P^{2}$ in place of $Q^{1}$ and the second marginal in
place of the first one, we see that if $P^{2}$ is well defined, then
$H(\mu|\mu^{R^{0}})<+\infty$. In particular, there is nothing to check before
starting the algorithm: if the algorithm is able to compute $P^{2}$, then it
means that our assumption is satisfied and that the convergence holds.

  * •

Note that the topology for the convergence stated in the theorem does not
matter since we are working in finite dimensional spaces. However, we believe
that the result is still true replacing $\mathcal{D}$ and $\mathcal{F}$ by
general Polish spaces. In this case, the convergence needs to be understood in
the sense of the narrow topology, a topology for which the sequences $(P^{n})$
and $(Q^{n})$ can be proved to be compact due to the properties of their
marginals.

  * •

Remarkably, we will be able to prove this theorem without deriving the
optimality conditions for $\mu^{*}$ and $\nu^{*}$. However, these optimality
conditions will be needed in the next section, and hence written at
Proposition 19.

  * •

As developed in [7], there exists a strong analogy between the relative
entropy and square Euclidean distances, and this in spite of the lack of
symmetry of the first. In particular, following this analogy, the Sinkhorn
algorithm (1) consists in iteratively othogonaly projecting on the convex sets
of measures absolutely continuous w.r.t. $R$ satisfying the first and second
marginal constraint respectively.

With this picture in mind, we can give in Figure 1 a visual representation of
the scalable and non-scalable case. In the scalable case, the two convex sets
intersect, and the sequences $(P^{n})_{n\in\mathbb{N}^{*}}$ and
$(Q^{n})_{n\in\mathbb{N}^{*}}$ converge towards the point of the intersection
that is the closest to $R$. In the non-scalable case, the two convex sets do
not intersect. However, the sequences $(P^{n})_{n\in\mathbb{N}^{*}}$ and
$(Q^{n})_{n\in\mathbb{N}^{*}}$ still converge respectively to $P^{*}$ and
$Q^{*}$, the two extreme points of the shortest line segment connecting both
sets. Theorem 11 indeed justifies this type of behaviour for the Sinkhorn
algorithm.

One should still keep in mind that this analogy and our drawings are only
sketchy. In reality, the projections are not orthogonal, and the convex sets
have polygonal borders.

Figure 1:  Sketchy representation of the Sinkhorn algorithm in the scalable
case (to the left) and nonscalable case (to the right).

###### Proof.

Step 1: All the objects are well defined.

Let us first show that under the assumption of the theorem, the sequences
$(P^{n})_{n\in\mathbb{N}^{*}}$ and $(Q^{n})_{n\in\mathbb{N}^{*}}$ are well
defined. We start with $P^{1}$. As by assumption
$\mu\ll\mu^{R^{0}}\ll\mu^{R}$, Corollary 6 shows that $P^{1}$ is well defined,
and that for all $i,j$,

$P^{1}_{ij}=\frac{\mu_{i}}{\mu^{R}_{i}}R_{ij},$

with convention $\frac{0}{0}=0$. Clearly, $R^{0}\ll P^{1}$, as the support of
the latter is

$\Big{\\{}(x_{i},y_{j})\in\mathcal{D}\times\mathcal{F}\mbox{ s.t.\
}R_{ij}>0\mbox{ and }\mu_{i}>0\Big{\\}}\supset\mathcal{E}.$

Therefore, $\nu\ll\nu^{R^{0}}\ll\nu^{P^{1}}$. So once again, Corollary 6 shows
that $Q^{1}$ is well defined, and that for all $i,j$,

$Q^{1}_{ij}=\frac{\nu_{j}}{\nu^{P^{1}}_{j}}P^{1}_{ij},$

with convention $\frac{0}{0}=0$. The support of $Q^{1}$ is

$\displaystyle\Big{\\{}(x_{i},y_{j})\in\mathcal{D}\times\mathcal{F}$
$\displaystyle\mbox{ s.t.\ }P^{1}_{ij}>0\mbox{ and }\nu_{j}>0\Big{\\}}$
$\displaystyle=\Big{\\{}(x_{i},y_{j})\in\mathcal{D}\times\mathcal{F}\mbox{
s.t. }R_{ij}>0\mbox{ and }\mu_{i}>0\mbox{ and
}\nu_{j}>0\Big{\\}}=\mathcal{E}.$

Then, a direct induction argument relying on the following formulas holding
for all $n\in\mathbb{N}$ and all $i,j$:

$P^{n+1}_{ij}=\frac{\mu_{i}}{\mu^{Q^{n}}_{i}}Q^{n}_{ij}\qquad\mbox{and}\qquad
Q^{n+1}_{ij}=\frac{\nu_{i}}{\nu^{P^{n+1}}_{j}}P^{n+1}_{ij},$ (15)

with convention $\frac{0}{0}=0$ show that for all $n\geq 2$, $P^{n}$ and
$Q^{n}$ are well defined and admit $\mathcal{E}$ as their common support.

Let us now show that $\mu^{*}$ and $\nu^{*}$ are well defined. Their role are
symmetric, so we just need to show that $\mu^{*}$ is well defined. First
$Q^{1}$ satisfies $\nu^{Q^{1}}=\nu$ and $H(Q^{1}|R)<+\infty$. Therefore, the
problem

$\inf\Big{\\{}H(\bar{\mu}|\mu),\,\bar{\mu}=\mu^{Q}\mbox{ for some }Q\mbox{
with }H(Q|R)<+\infty\mbox{ and }\nu^{Q}=\nu\Big{\\}}$

consists in minimizing the continuous (on its domain) and strictly convex
function $\bar{\mu}\mapsto H(\bar{\mu}|\mu)$ over the nonempty compact convex
set

$\Big{\\{}\bar{\mu}=\mu^{Q}\mbox{ for some }Q\mbox{ with }H(Q|R)<+\infty\mbox{
and }\nu^{Q}=\nu\Big{\\}}.$

Hence, it admits a unique solution $\mu^{*}$.

Finally, let us show that $P^{*}$ and $Q^{*}$ are well defined. Once again,
their role are symmetric, so we only show the existence of $Q^{*}$. We already
saw that $\mu^{*}$ is well defined. By definition of the latter, there exists
$\bar{Q}$ with $\nu^{\bar{Q}}=\nu$, $\mu^{\bar{Q}}=\mu^{*}$ and
$H(\bar{Q}|R)<+\infty$. So $\operatorname{Sch}(R;\mu^{*},\nu)$ consists in
minimizing the continuous (on its domain) and strictly convex function
$Q\mapsto H(Q|R)$ on the nonempty compact convex set

$\Pi(\mu^{*},\nu)\cap\Big{\\{}Q\in\mathcal{P}(\mathcal{D}\times\mathcal{F})\mbox{
such that }H(Q|R)<+\infty\Big{\\}}.$

So it admits a unique solution $Q^{*}$.

Step 2: A formula for $H(Q|R)$, for all $Q\in\Pi(\mu^{*},\nu)$ with
$H(Q|R)<+\infty$.

Recalling $Q^{0}=R$, we infer from (15) that for all
$(x_{i},y_{j})\in\mathcal{E}$ and $n\in\mathbb{N}^{*}$,

$Q^{n}_{ij}=\frac{\nu_{j}}{\nu^{P^{n}}_{j}}\times\frac{\mu_{i}}{\mu^{Q^{n-1}}_{i}}\times\dots\times\frac{\nu_{j}}{\nu^{P^{1}_{j}}}\times\frac{\mu_{i}}{\mu^{Q^{0}}_{i}}\times
R_{ij}.$ (16)

Observe that in the product in the r.h.s., because we assumed that
$(x_{i},y_{j})\in\mathcal{E}$, the common support of all the iterates of the
Sinkhorn algorithm, all the factors are positive.

In addition, for all $Q\in\Pi(\mu^{*},\nu)$ with finite entropy w.r.t. $R$,
the support of $Q$ is included in $\mathcal{E}$. This is because
$H(\mu^{*}|\mu)<+\infty$, and thereby $\mu^{*}\ll\mu$. Therefore, we deduce
that $Q\ll Q^{n}$.

So as a consequence of (16), for all $i,j$ in the support of $Q$,

$\displaystyle\log\frac{Q_{ij}}{R_{ij}}$
$\displaystyle=\log\frac{Q_{ij}}{Q^{n}_{ij}}+\sum_{k=1}^{n}\Big{\\{}\log\frac{\nu_{j}}{\nu^{P^{k}}_{j}}+\log\frac{\mu_{i}}{\mu^{Q^{k-1}}_{i}}\Big{\\}}$
$\displaystyle=\log\frac{Q_{ij}}{Q^{n}_{ij}}+\sum_{k=1}^{n}\Big{\\{}\log\frac{\nu_{j}}{\nu^{P^{k}}_{j}}-\log\frac{\mu^{*}_{i}}{\mu_{i}}+\log\frac{\mu^{*}_{i}}{\mu^{Q^{k-1}}_{i}}\Big{\\}}.$

Let us multiply this equality by $Q_{ij}$, and sum over $i,j$. We get

$\sum_{i,j}Q_{ij}\log\frac{Q_{ij}}{R_{ij}}=\sum_{i,j}Q_{ij}\log\frac{Q_{ij}}{Q^{n}_{ij}}+\sum_{k=1}^{n}\left\\{\sum_{i,j}Q_{ij}\log\frac{\nu_{j}}{\nu^{P^{k}}_{j}}-\sum_{i,j}Q_{ij}\log\frac{\mu^{*}_{i}}{\mu_{i}}+\sum_{i,j}Q_{ij}\log\frac{\mu^{*}_{i}}{\mu^{Q^{k-1}}_{i}}\right\\}.$

Fix $k$ in $\\{1,\dots,n\\}$, and consider the term:

$\sum_{i,j}Q_{ij}\log\frac{\nu_{j}}{\nu^{P^{k}}_{j}}=\sum_{j}\left(\sum_{i}Q_{ij}\right)\log\frac{\nu_{j}}{\nu^{P^{k}}_{j}}.$

As the second marginal of $Q$ is $\nu$, we find

$\sum_{i,j}Q_{ij}\log\frac{\nu_{j}}{\nu^{P^{k}}_{j}}=\sum_{j}\nu_{j}\log\frac{\nu_{j}}{\nu^{P^{k}}_{j}}.$

As the first marginal of $Q$ is $\mu^{*}$, we can reason in the same way for
the other terms of the same type, and find

$\sum_{i,j}Q_{ij}\log\frac{Q_{ij}}{R_{ij}}=\sum_{i,j}Q_{ij}\log\frac{Q_{ij}}{Q^{n}_{ij}}+\sum_{k=1}^{n}\left\\{\sum_{j}\nu_{j}\log\frac{\nu_{j}}{\nu^{P^{k}}_{j}}-\sum_{i}\mu^{*}_{i}\log\frac{\mu^{*}_{i}}{\mu_{i}}+\sum_{i}\mu^{*}_{i}\log\frac{\mu^{*}_{i}}{\mu^{Q^{k-1}}_{i}}\right\\}.$

We now use the Definition 3 of the relative entropy to find that this identity
means:

$\displaystyle H(Q|R)-\mathsf{M}(R)+\mathsf{M}(Q)$
$\displaystyle=H(Q|Q^{n})-\mathsf{M}(Q^{n})+\mathsf{M}(Q)$
$\displaystyle\qquad+\sum_{k=1}^{n}\Big{\\{}H(\nu|\nu^{P^{k}})-\mathsf{M}(\nu)+\mathsf{M}(\nu^{P^{k}})-H(\mu^{*}|\mu)+\mathsf{M}(\mu^{*})-\mathsf{M}(\mu)$
$\displaystyle\hskip
147.0pt+H(\mu^{*}|\mu^{Q^{k-1}})-\mathsf{M}(\mu^{*})+\mathsf{M}(\mu^{Q^{k-1}})\Big{\\}},$

or, simplifying the masses $\mathsf{M}(Q)$ and $\mathsf{M}(\mu^{*})$ appearing
several times,

$\displaystyle
H(Q|R)-\mathsf{M}(R)=H(Q|Q^{n})-\mathsf{M}(Q^{n})\\!+\\!\sum_{k=1}^{n}\\!\Big{\\{}H(\nu|\nu^{P^{k}})-\mathsf{M}(\nu)+$
$\displaystyle\mathsf{M}(\nu^{P^{k}})-H(\mu^{*}|\mu)-\mathsf{M}(\mu)$
$\displaystyle+H(\mu^{*}|\mu^{Q^{k-1}})+\mathsf{M}(\mu^{Q^{k-1}})\Big{\\}}.$

Let us check how the masses simplify. By (4) as $Q$ and $Q^{k}$,
$k\in\mathbb{N}^{*}$ admit $\nu$ as their second marginals, they have the same
total masses, and it coincides with the one of their first marginals. Namely,

$\forall k\geq
1,\qquad\mathsf{M}(\nu)=\mathsf{M}(Q^{k})=\mathsf{M}(Q)=\mathsf{M}(\mu^{Q^{k}})=\mathsf{M}(\mu^{*}).$

In the same way,

$\forall k\geq
1,\qquad\mathsf{M}(\mu)=\mathsf{M}(P^{k})=\mathsf{M}(\nu^{P^{k}}).$

And finally, as $Q^{0}=R$, we also have

$\mathsf{M}(R)=\mathsf{M}(Q^{0})=\mathsf{M}(\mu^{Q^{0}}).$

Coming back to our entropy identity, we can simplify more and get

$H(Q|R)=H(Q|Q^{n})+\sum_{k=1}^{n}\Big{\\{}H(\nu|\nu^{P^{k}})-H(\mu^{*}|\mu)\Big{\\}}+\sum_{k=1}^{n}H(\mu^{*}|\mu^{Q^{k-1}}).$
(17)

(The only subtlety is that for $k=1$ and for $k=1$ only,
$\mathsf{M}(\mu^{Q^{k-1}})$ does not simplify with $\mathsf{M}(\nu)$. But then
$\mathsf{M}(\mu^{Q^{k-1}})=\mathsf{M}(\mu^{Q^{0}})$ simplifies with
$\mathsf{M}(R)$, and $\mathsf{M}(\nu)$ simplifies with $\mathsf{M}(Q^{n})$.)

Now, we claim that every term in the first sum in the r.h.s. is nonnegative,
_i.e._ that $H(\nu|\nu^{P^{k}})\geq H(\mu^{*}|\mu)$. First, we know that
$\nu^{Q^{k}}=\nu$, so that:

$\displaystyle H(\nu|\nu^{P^{k}})$
$\displaystyle=\sum_{j}\nu_{j}\log\frac{\nu_{j}}{\nu^{P^{k}}_{j}}+\mathsf{M}(\nu^{P^{k}})-\mathsf{M}(\nu)$
$\displaystyle=\sum_{ij}Q^{k}_{ij}\log\frac{\nu_{j}}{\nu^{P^{k}}_{j}}+\mathsf{M}(\nu^{P^{k}})-\mathsf{M}(\nu)$
$\displaystyle=\sum_{ij}Q^{k}_{ij}\log\frac{Q^{k}_{ij}}{P^{k}_{ij}}+\mathsf{M}(P^{k})-\mathsf{M}(Q^{k})=H(Q^{k}|P^{k}),$

where we used at the third line that from (15), we know that for all $j$ and
$1\leq k\leq n$, $\nu_{j}/\nu^{P^{k}}_{j}=Q^{k}_{ij}/P^{k}_{ij}$. Second, by
Corollary 6,

$H(\nu|\nu^{P^{k}})=H(Q^{k}|P^{k})\geq
H(\mu^{Q^{k}}|\mu^{P^{k}})=H(\mu^{Q^{k}}|\mu).$

Finally, $Q^{k}$ has finite entropy w.r.t. $R$ (use for instance (16) with
$n=k$) and admits $\nu$ as a first marginal. So by optimality of $\mu^{*}$,
$H(\mu^{Q^{k}}|\mu)\geq H(\mu^{*}|\mu)$. Our claim follows.

Step 3: Consequence of (17), convergence of the marginals.

As a consequence of Step 2, both sums in the r.h.s. of (17) are bounded sums
of nonnegative terms. Therefore, they converge as $n\to+\infty$, and their
terms tend to $0$ as $k\to+\infty$. We deduce in particular that

$H(\mu^{*}|\mu^{Q^{n}})\underset{n\to+\infty}{\longrightarrow}0.$

In particular, by continuity of $H$ w.r.t. its second variable as stated in
Proposition 4, and by compactness of
$\\{\bar{\mu}\in\mathcal{M}_{+}(\mathcal{D})\mbox{ s.t.
}\mathsf{M}(\bar{\mu})=\mathsf{M}(\nu)\\}$, $\mu^{Q^{n}}\to\mu^{*}$. So now
let us pick $\bar{Q}$ any limit point of $(Q^{n})$. Such a limit point exist
by compactness of $\\{Q\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})\mbox{
s.t. }\mathsf{M}(Q)=\mathsf{M}(\nu)\\}$. It follows from
$\mu^{Q^{n}}\to\mu^{*}$ that $\mu^{\bar{Q}}=\mu^{*}$.

Step 4: $\bar{Q}=Q^{*}$.

Let us show that $\bar{Q}=Q^{*}$, so that actually the whole sequence
$(Q^{n})$ converges towards $Q^{*}$. On the one hand, passing to the limit
$n\to+\infty$ along the subsequences generating $\bar{Q}$ in (17) and using
the continuity of $H$ w.r.t. the second variable as stated in Proposition 4,
we find

$H(Q|R)=H(Q|\bar{Q})+\sum_{k=1}^{+\infty}\Big{\\{}H(\nu|\nu^{P^{k}})-H(\mu^{*}|\mu)\Big{\\}}+\sum_{k=1}^{+\infty}H(\mu^{*}|\mu^{Q^{k-1}}).$
(18)

On the other hand, as for all $n\in\mathbb{N}^{*}$, $Q^{n}\ll R$, this is also
true for $\bar{Q}$. In particular, $H(\bar{Q}|R)<+\infty$, and as
$\bar{Q}\in\Pi(\mu^{*},\nu)$, we can apply (18) with $\bar{Q}$ in place of
$Q$, and find

$H(\bar{Q}|R)=\sum_{k=1}^{+\infty}\Big{\\{}H(\nu|\nu^{P^{k}})-H(\mu^{*}|\mu)\Big{\\}}+\sum_{k=1}^{+\infty}H(\mu^{*}|\mu^{Q^{k-1}}).$
(19)

Now it remains to apply (18) with $Q=Q^{*}$ and to plug the previous equality
to find

$H(Q^{*}|R)=H(Q^{*}|\bar{Q})+H(\bar{Q}|R).$

As by optimality of $R^{*}$, $H(\bar{Q}|R)\geq H(Q^{*}|R)$, we can conclude
that $H(Q^{*}|Q)=0$. Therefore, $\bar{Q}=Q^{*}$, as announced.

The proof of $P^{n}\to P^{*}$ follows the same lines. ∎

As a free output of the proof of Theorem 1, we can show that we could have
swapped $\mu$ and $\bar{\mu}$, and $\nu$ and $\bar{\nu}$ in the definitions
(13) of $\mu^{*}$ and $\nu^{*}$ respectively. This is justified in the
following remark.

######  Remark 13.

Observe the following optimization problem, where $R$, $\mu$ and $\nu$ are
given, and where the competitor is $\bar{\nu}$:

$\min\Big{\\{}H(\nu|\bar{\nu})\,\Big{|}\,\bar{\nu}=\nu^{P},\mbox{ for some
}P\mbox{ with }H(P|R)<+\infty\mbox{ and }\mu^{P}=\mu\Big{\\}}.$ (20)

This problem is almost the same as the one defining $\nu^{*}$ in (13), except
from the fact that $\nu$ and $\bar{\nu}$ are swapped in the relative entropy.
In this remark, we justify that the solution of this problem is $\nu^{*}$ as
well, and that the corresponding optimal value is $H(\mu^{*}|\mu)$.

Provided there exists a competitor $\bar{\nu}$ for this problem with
$H(\nu|\bar{\nu})<+\infty$, we can find $P$ such that $H(P|R)<+\infty$ and
$P\in\Pi(\mu,\bar{\nu})$, and $Q\in\mathcal{P}(\mathcal{D}\times\mathcal{F})$
defined for all $i,j$ by

$Q_{ij}:=\frac{\nu_{j}}{\bar{\nu}_{j}}P_{ij},$

which is legitimate since $H(\nu|\bar{\nu})<+\infty$. We have then
$H(Q|R)<+\infty$ and $\nu^{Q}=\nu$. Hence, using the definition (13) of
$\mu^{*}$, we have

$H(\nu|\bar{\nu})=H(Q|P)\geq H(\mu^{Q}|\mu)\geq H(\mu^{*}|\mu),$

where the first equality is a direct computation, and where the first
inequality is obtained using Corollary 6.

On the other hand, as soon as the assumption of Theorem 11 holds, $\nu^{*}$ is
a competitor for the problem in (20), and so in particular $H(\nu|\nu^{*})\geq
H(\mu^{*}|\mu)$. But because the terms of the first series in (18) tend to $0$
and $\nu^{P^{k}}\to\nu^{*}$, we conclude that actually,
$H(\nu|\nu^{*})=H(\mu^{*}|\mu)$ and $\nu^{*}$ is a solution of (20). Finally,
it is easy to see that a solution $\bar{\nu}$ of (20) must satisfy
$\bar{\nu}\ll\nu$ (because conditioning on the support of $\nu$ reduces the
entropy), and by strict convexity of $\bar{\nu}\mapsto H(\nu|\bar{\nu})$ on
the set $\\{\bar{\nu}\ll\nu\\}$, under the assumption of Theorem 11, the
problem (20) admits $\nu^{*}$ as its unique solution, so that (20) can be used
as an alternative definition of $\nu^{*}$.

Of course, we could argue in the same way to provide an alternative definition
of $\mu^{*}$, and we have the following equalities:

$H(\nu|\nu^{*})=H(\mu^{*}|\mu)\qquad\mbox{and}\qquad
H(\mu|\mu^{*})=H(\nu^{*}|\nu).$

In particular, $\mu^{*}\sim\mu$ and $\nu^{*}\sim\nu$ in the sense of measures.

We also give another remark concerning the generalization of Theorem 11 to
Polish spaces.

######  Remark 14.

We crucially use the fact that $\mathcal{D}$ and $\mathcal{F}$ are finite in
order to obtain (18) and (19). In the continuous case, as $H$ is not more than
lower semicontinuous w.r.t the second variable, identity (18) becomes an
inequality, where $=$ is replaced by $\geq$, which is the good direction for
the proof. The difficulty is then to find an equality sign in (19).

##  4 $\mathbf{\Gamma}$-convergence in the marginal penalization problem

In this section, we want to show that when $R$, $\mu$ and $\nu$ are such that
the Schrödinger problem $\operatorname{Sch}(R;\mu,\nu)$ has no solution, then
the limit points $P^{*}$ and $Q^{*}$ given by Theorem 11 are relevant in view
of the possible applications of the Sinkhorn algorithm.

To do so, let us think of $R$ as an imperfect theoretical model describing the
coupling between the initial and final positions of the particles of a large
system. Also, let us imagine that $\mu$ and $\nu$ are data obtained by
measuring the positions of the particles of the actual system that $R$ is
supposed to describe, at the initial and final time. In this situation, if
$\operatorname{Sch}(R;\mu,\nu)$ has a solution $R^{*}$, this solution is
interpreted as the model that is the closest to $R$ that can explain the data.

However, even when $R$ is a rather good model, and when $\mu$ and $\nu$ are
rather precise measurements, it is possible that
$\operatorname{Sch}(R;\mu,\nu)$ has no solution for several reasons:

  * •

The first reason could be that our modeling does not take into account some
physical phenomena. For instance, in Subsection 4.1, we will consider the case
where the true system allows creation or annihilation of mass with very small
probability, whereas the modeling does not.

  * •

Another reason could be that $\mu$ and $\nu$ are only approximations of the
real marginals. This can result from imprecise or biased measurements, or from
a restricted amount of collected data. This will be considered in Subsection
4.2.

In both cases, it is very natural to relax the marginal constraints in (10) by
introducing a fitting term in the value functional, that cancels when the
constraints are satisfied, but which remains finite otherwise.

The main result of this section asserts that in these two situations, that are
actually very close, the limit points $P^{*}$ and $Q^{*}$ of the Sinkhorn
algorithm allow to compute the solution of the relaxed problem when the new
fitting term takes the form of an entropy, in the limit where the level of
marginal penalization tends to $+\infty$. The second case is a direct
consequence of the first one, but that we wanted to keep separated because it
does not have the same interpretation.

###  4.1 Unbalanced problems

In this subsection, we give ourselves
$R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$
as before, and we study the following optimization problem, which is a
reasonable modification of $\operatorname{Sch}(R;\mu,\nu)$ where the marginal
constraints are replaced with marginal penalizations:

$\min\Big{\\{}H(\bar{R}|R)+\lambda\big{(}H(\mu^{\bar{R}}|\mu)+H(\nu^{\bar{R}}|\nu)\big{)}\,\Big{|}\,\bar{R}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})\Big{\\}},$
(21)

where $\lambda>0$ parametrizes the level of penalization.

This approach is extremely reminiscent of the idea introduced by Liero, Mielke
and Savaré in [20] to deal with unbalanced data, that is, when
$\mathsf{M}(\mu)\neq\mathsf{M}(\nu)$, in optimal transport problems. This was
the starting point of the theory of _unbalanced optimal transport_ , also
discovered independently by other teams [15, 6].

More precisely, we will study the limit of the problem in (21) as
$\lambda\to+\infty$. In this limit, it is actually more convenient to call
$\varepsilon=1/\lambda$ and to multiply the value functional by $\varepsilon$,
to find the problem that we call
$\operatorname{Sch}^{\varepsilon}(R;\mu,\nu)$:

$\operatorname{Sch}^{\varepsilon}(R;\mu,\nu):=\min\Big{\\{}\varepsilon
H(\bar{R}|R)+H(\mu^{\bar{R}}|\mu)+H(\nu^{\bar{R}}|\nu)\,\Big{|}\,\bar{R}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})\Big{\\}}.$

As we want to study the behavior of this problem in the limit $\varepsilon\to
0$, we define the following functionals:

$\displaystyle\Lambda^{\varepsilon}:\bar{R}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$
$\displaystyle\mapsto\varepsilon
H(\bar{R}|R)+H(\mu^{\bar{R}}|\mu)+H(\nu^{\bar{R}}|\nu),$
$\displaystyle\Lambda:\bar{R}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$
$\displaystyle\mapsto\chi_{H(\bar{R}|R)<+\infty}+H(\mu^{\bar{R}}|\mu)+H(\nu^{\bar{R}}|\nu),$

where $\chi_{H(\bar{R}|R)<+\infty}$ is the convex indicatrix taking value $0$
on the set

$\Big{\\{}\bar{R}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})\mbox{ such
that }H(\bar{R}|R)<+\infty\Big{\\}},$

and $+\infty$ elsewhere.

The following proposition follows from standard arguments in the theory of
$\Gamma$-convergence, see for instance [3, Theorem 1.47], and from the strict
convexity of the relative entropy w.r.t. its first variable. We omit the
proof.

######  Proposition 15.

We have:

$\Gamma-\lim_{\varepsilon\to 0}\Lambda^{\varepsilon}=\Lambda.$

In particular, assuming that $\Lambda$ is not uniformly infinite, let us call
$R_{\mathrm{opt}}$ one of its minimizers, $\mu^{g}:=\mu^{R_{\mathrm{opt}}}$
and $\nu^{g}:=\nu^{R_{\mathrm{opt}}}$. The marginals $\mu^{g}$ and $\nu^{g}$
do not depend on the choice of $R_{\mathrm{opt}}$, and as $\varepsilon\to 0$,
the unique solution $R^{\varepsilon}$ of
$\operatorname{Sch}^{\varepsilon}(R;\mu,\nu)$ exists and converges towards the
solution of $\operatorname{Sch}(R;\mu^{g},\nu^{g})$.

######  Remark 16.

In the notations $\mu^{g}$ and $\nu^{g}$, the $g$ stands for _geometric_. This
is because as shown in Theorem 17, $\mu^{g}$ and $\nu^{g}$ are respectively
the componentwise geometric means of $\mu$ and $\mu^{*}$, and of $\nu$ and
$\nu^{*}$.

Therefore, studying the behavior of
$\operatorname{Sch}^{\varepsilon}(R;\mu,\nu)$ in the limit $\varepsilon\to 0$
reduces to the study of the Schrödinger problem with modified marginals
$\mu^{g}$ and $\nu^{g}$. The following theorem shows the link between $R^{*}$
– the solution of $\operatorname{Sch}(R;\mu^{g},\nu^{g})$ – on the one hand,
and $P^{*}$ and $Q^{*}$ from Theorem 11 on the other hand.

######  Theorem 17.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$
satisfy Assumption 8. Then the functional $\Lambda$ is not uniformly infinite.
Moreover, considering $P^{*}$ and $Q^{*}$ as given by Theorem 11, and
$\mu^{g}$ and $\nu^{g}$ as given by Proposition 15, the solution of
$\operatorname{Sch}(R;\mu^{g},\nu^{g})$ is the componentwise geometric mean of
$P^{*}$ and $Q^{*}$, that is, the matrix $R^{*}$ defined for all $i,j$ by

$R^{*}_{ij}:=\sqrt{P^{*}_{ij}Q^{*}_{ij}}.$ (22)

Also, if $\mu^{*}$ and $\nu^{*}$ are defined by (13), $\mu^{g}$ and $\nu^{g}$
are the componentwise geometric means of $\mu^{*}$ and $\mu$ for the first
one, and of $\nu^{*}$ and $\nu$ for the second one. In other terms, we have
for all $i,j$,

$\mu^{g}_{i}=\sqrt{\mu^{*}_{i}\mu_{i}}\qquad\mbox{and}\qquad\nu^{g}_{j}=\sqrt{\nu^{*}_{j}\nu_{j}}.$
(23)

######  Remark 18.

  * •

Having in mind the approach of [20], we can give the following interpretation
of the matrix $R^{*}$: In the degenerate case where the Schrödinger problem
has no solution, it is necessary to allow creation and annihilation of mass to
find solutions. Following [20], we can do this by replacing the balanced
problem $\operatorname{Sch}(R;\mu,\nu)$ by the unbalanced problem
$\operatorname{Sch}^{\varepsilon}(R;\mu,\nu)$. Following this analogy,
$\lambda=\frac{1}{\varepsilon}$ parametrizes the cost of creating particles.
The matrix $R^{*}$ from Theorem 17 is therefore the limit of these solutions
when the cost of creating or destroying matter tends to $+\infty$.

  * •

A small adaptation of the proof shows that given $\alpha\in[0,1]$, if we
replace the problem in (21) by

$\min\Big{\\{}H(\bar{R}|R)+\lambda\Big{(}(1-\alpha)H(\mu^{\bar{R}}|\mu)+\alpha
H(\nu^{\bar{R}}|\nu)\Big{)}\,\Big{|}\,\bar{R}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})\Big{\\}},$

and if we call $R^{\alpha,\lambda}$ its solution, then as $\lambda\to+\infty$,
we have for all $i,j$:

$R^{\alpha,\lambda}_{ij}\underset{\lambda\to+\infty}{\longrightarrow}\big{(}P^{*}_{ij}\big{)}^{1-\alpha}\big{(}Q^{*}_{ij}\big{)}^{\alpha}.$

To prove this theorem, we will need to study carefully the optimality
conditions for $\mu^{*}$ and $\nu^{*}$. This could be done writing the Karush-
Kuhn-Tucker conditions for the corresponding optimalization problems. We will
rather adopt a more hand by hand approach, that is more likely to be
generalizable in the continuous case. This is done in the following
proposition.

######  Proposition 19.

Assume that the conditions of Theorem 11 are fulfilled. For all $i,j$, we have

$P^{*}_{ij}=\frac{\mu_{i}}{\mu^{*}_{i}}Q^{*}_{ij}\qquad\mbox{and}\qquad
Q^{*}_{ij}=\frac{\nu_{j}}{\nu^{*}_{j}}P^{*}_{ij},$ (24)

with convention $\frac{0}{0}=0$. In particular, $P^{*}$ and $Q^{*}$ are
equivalent, and we call $\mathcal{S}$ their common support. Also, recall the
definition of $\mathcal{E}$ in (11). Of course
$\mathcal{S}\subset\mathcal{E}$. Finally, we call for all $i,j$

$\varphi_{i}:=\log\frac{\mu^{*}_{i}}{\mu_{i}}\qquad\mbox{and}\qquad\psi_{j}:=\log\frac{\nu^{*}_{j}}{\nu_{j}}.$
(25)

For all $(i,j)\in\mathcal{E}$, $\varphi_{i}$ and $\psi_{j}$ are well defined
in $\mathbb{R}$, and:

$\left\\{\begin{aligned} \varphi_{i}+\psi_{j}&=0,&&\mbox{if
}(i,j)\in\mathcal{S},\\\ \varphi_{i}+\psi_{j}&\geq 0,&&\mbox{if
}(i,j)\in\mathcal{E}.\end{aligned}\right.$ (26)

###### Proof of Proposition 19.

To get (24), it suffices to let $n$ tend to $+\infty$ in (15). The fact that
$\mathcal{S}\subset\mathcal{E}$ relies on the closed property of $P^{n}$ and
$Q^{n}$ defined in (1) to have its support included in $\mathcal{E}$ for
$n\geq 2$. If $(i,j)\in\mathcal{E}$, let us check that $\varphi_{i}$ and
$\psi_{j}$ are well defined. On the one hand, by definition of $\mathcal{E}$,
$i$ is in the support of $\mu$ and $j$ is in the support of $\nu$. On the
other hand, as observed in Remark 13, $\mu^{*}\sim\mu$ and $\nu^{*}\sim\nu$.
Our claim follows.

Now, let $(i,j)\in\mathcal{S}$. A consequence of (24) is

$P^{*}_{ij}=\frac{\mu^{*}_{i}}{\mu_{i}}\frac{\nu^{*}_{j}}{\nu_{j}}P^{*}_{ij}=\exp(\varphi_{i}+\psi_{j})P^{*}_{ij}.$

As $(i,j)$ is in the support of $P^{*}$ by definition of $\mathcal{S}$, we
conclude that $\varphi_{i}+\psi_{j}=0$.

Finally, it remains to prove that for all $(i,j)\in\mathcal{E}$,
$\varphi_{i}+\psi_{j}\geq 0$. For this we use the optimality of
$H(\mu^{*}|\mu)=H(\mu^{Q^{*}}|\mu)$ over all $Q$ satisfying $\nu^{Q}=\nu$. So
let us take $(i,j)\in\mathcal{E}$. As $\nu_{j}>0$, there exists $i^{\prime}$
such that $(i^{\prime},j)\in\mathcal{S}$, that is, such that
$Q^{*}_{i^{\prime}j}>0$. Let us define for $\varepsilon>0$

$Q^{\varepsilon}=Q^{*}+\varepsilon\delta_{ij}-\varepsilon\delta_{i^{\prime}j},$

where $\delta_{ij}$ is the matrix whose only nonzero coefficient is a one at
position $(i,j)$, and similarly for $\delta_{i^{\prime}j}$. If $\varepsilon$
is sufficiently small,
$Q^{\varepsilon}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\nu^{Q^{\varepsilon}}=\nu$ and with obvious notations,
$\mu^{Q^{\varepsilon}}=\mu^{*}+\varepsilon\delta_{i}-\varepsilon\delta_{i^{\prime}}$.
Therefore, for such $\varepsilon$,

$H(\mu^{Q^{\varepsilon}}|\mu)\geq H(\mu^{*}|\mu).$

derivating to the right this inequality at $\varepsilon=0$, we find

$\log\frac{\mu_{i}^{*}}{\mu_{i}}-\log\frac{\mu_{i^{\prime}}^{*}}{\mu_{i^{\prime}}}\geq
0,$

which rewrites $\varphi_{i}-\varphi_{i^{\prime}}\geq 0$. But
$(i^{\prime},j)\in\mathcal{S}$ so $\varphi_{i^{\prime}}=-\psi_{j}$, and so
$\varphi_{i}+\psi_{j}\geq 0$. ∎

With this proposition at hand, we can prove Theorem 17.

###### Proof of Theorem 17.

The fact that under Assumption 8, $\Lambda$ is not uniformly infinite follows
from observing that $\Lambda(R^{0})<+\infty$, where $R_{0}$ was defined
Assumption 8. Now we reason in two steps. First we will prove using
Proposition 19 that $R^{*}$ defined by (22) is an optimizer of $\Lambda$, and
then that it is the solution of the Schrödinger problem between its marginals.

Step 1: $R^{*}$ is an optimizer of $\Lambda$.

To see that $R^{*}$ is an optimizer of $\Lambda$, we first give a formula
relating the vectors $\varphi$ and $\psi$ as defined by formula (25) and the
marginals $\mu^{R^{*}}$ and $\nu^{R^{*}}$ of $R^{*}$. Using (24) and the
definition (22) of $R^{*}$, we see that for all $i,j$,

$R^{*}_{ij}=\sqrt{\frac{\nu_{j}}{\nu^{*}_{j}}}P^{*}_{ij}=\sqrt{\frac{\mu_{i}}{\mu^{*}_{i}}}Q^{*}_{ij}.$
(27)

Summing respectively these identities w.r.t. $i$ and $j$, we deduce that for
all $i,j$,

$\mu^{R^{*}}_{i}=\sqrt{\mu_{i}^{*}\mu_{i}}=\sqrt{\frac{\mu^{*}_{i}}{\mu_{i}}}\mu_{i}\qquad\mbox{and}\qquad\nu^{R^{*}}_{j}=\sqrt{\nu_{j}^{*}\nu_{j}}=\sqrt{\frac{\nu^{*}_{j}}{\nu_{j}}}\nu_{j}.$

Let us define for all $i,j$:

$Z^{\mu}_{i}:=\log\frac{\mu^{R^{*}}_{i}}{\mu_{i}}=\frac{1}{2}\varphi_{i}\qquad\mbox{and}\qquad
Z^{\nu}_{j}:=\log\frac{\nu^{R^{*}}_{j}}{\nu_{j}}=\frac{1}{2}\psi_{j}.$

Note that for all $(x_{i},y_{j})\in\mathcal{E}$, $Z^{\mu}_{i}$ and
$Z^{\nu}_{j}$ are well defined in $\mathbb{R}$.

Now let $\bar{R}$ be such that $\Lambda(\bar{R})<+\infty$. Using inequality
(5) to bound from below each relative entropy, we have

$\displaystyle\Lambda(\bar{R})$
$\displaystyle=H(\mu^{\bar{R}}|\mu)+H(\nu^{\bar{R}}|\nu)$
$\displaystyle\geq\langle Z^{\mu},\mu^{\bar{R}}\rangle-\langle
e^{Z^{\mu}}-1,\mu\rangle+\langle Z^{\nu},\nu^{\bar{R}}\rangle-\langle
e^{Z^{\nu}}-1,\nu\rangle$
$\displaystyle=\frac{1}{2}\langle\varphi,\mu^{\bar{R}}\rangle+\frac{1}{2}\langle\psi,\nu^{\bar{R}}\rangle-\sum_{i}\big{\\{}\mu_{i}^{R^{*}}-\mu_{i}\big{\\}}-\sum_{j}\big{\\{}\nu_{j}^{R^{*}}-\nu_{j}\big{\\}}$
$\displaystyle=\frac{1}{2}\langle\varphi\oplus\psi,\bar{R}\rangle+\mathsf{M}(\mu)+\mathsf{M}(\nu)-2\mathsf{M}(R^{*}),$

where $\varphi\oplus\psi$ is the matrix defined for all $i,j$ by
$(\varphi\oplus\psi)_{ij}:=\varphi_{i}+\psi_{j}$. Now, because of the second
line of (26), as the support of $\bar{R}$ is easily seen to be a subset of
$\mathcal{E}$, we get

$\Lambda(\bar{R})\geq\mathsf{M}(\mu)+\mathsf{M}(\nu)-2\mathsf{M}(R^{*}).$

On the other hand, by definition of $Z^{\mu}$ and $Z^{\nu}$,

$\displaystyle\Lambda(R^{*})$
$\displaystyle=H(\mu^{R^{*}}|\mu)+H(\nu^{R^{*}}|\nu)$ $\displaystyle=\langle
Z^{\mu},\mu^{R^{*}}\rangle+\mathsf{M}(\mu)-\mathsf{M}(R^{*})+\langle
Z^{\nu},\nu^{R^{*}}\rangle+\mathsf{M}(\nu)-\mathsf{M}(R^{*})$
$\displaystyle=\frac{1}{2}\langle\varphi\oplus\psi,R^{*}\rangle+\mathsf{M}(\mu)+\mathsf{M}(\nu)-2\mathsf{M}(R^{*}).$

But now, as the support of $R^{*}$ is precisely $\mathcal{S}$, by the first
line of (26), we get

$\Lambda(R^{*})=\mathsf{M}(\mu)+\mathsf{M}(\nu)-2\mathsf{M}(R^{*}).$

We deduce that $\Lambda(\bar{R})\geq\Lambda(R^{*})$ and $R^{*}$ is indeed an
optimizer of $\Lambda$. In particular, $\mu^{g}=\mu^{R^{*}}$ and
$\nu^{g}=\nu^{R^{*}}$, which proves (23).

Step 2: $R^{*}$ is the solution of $\operatorname{Sch}(R;\mu^{g},\nu^{g})$.

To show that $R^{*}$ solves the Schrödinger problem between its marginals, we
consider another $\bar{R}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$
such that $\bar{R}\in\Pi(\mu^{g},\nu^{g})$ and $H(\bar{R}|R)<+\infty$. Then,
for $\varepsilon>0$, we define

$P^{\varepsilon}:=P^{*}+\varepsilon(\bar{R}-R^{*}).$

As $R^{*}\ll P^{*}$ (see (22)), whenever $\varepsilon$ is sufficiently small,
$P^{\varepsilon}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$, and in
addition, we easily check that $P^{\varepsilon}\in\Pi(\mu,\nu^{*})$. So by
definition (14) of $P^{*}$,

$H(P^{*}|R)\leq H(P^{\varepsilon}|R).$

Derivating this inequality to the right at $\varepsilon=0$, we find

$\sum_{ij}R^{*}_{ij}\log\frac{P^{*}_{ij}}{R_{ij}}\leq\sum_{ij}\bar{R}_{ij}\log\frac{P^{*}_{ij}}{R_{ij}},$

with convention $\frac{0}{0}=0$, $0\log 0=0$ and $a\log 0=-\infty$ for all
$a>0$. In particular, we deduce that $\bar{R}\ll P^{*}\sim R^{*}$, and our
inequality rewrites

$H(R^{*}|R)+H(\bar{R}|P^{*})-H(R^{*}|P^{*})\leq H(\bar{R}|R).$

The last thing to observe is that because of (27), $R^{*}$ is the solution of
the Schrödinger problem $\operatorname{Sch}(P^{*};\mu^{g},\nu^{g})$: a direct
application of (5) with
$Z_{ij}=\log\frac{R^{*}_{ij}}{P^{*}_{ij}}=-\bar{\psi}_{j}$ (which is well
defined on the support of $P^{*}$, and so on the support of $\bar{R}$)
provides

$\displaystyle H(\bar{R}|P^{*})$ $\displaystyle\geq\langle
Z,\bar{R}\rangle-\langle e^{Z}-1,P^{*}\rangle$
$\displaystyle=-\langle\bar{\psi},\nu^{g}\rangle+\sum_{ij}P^{*}_{ij}-R^{*}_{ij}$
$\displaystyle=\langle
Z,R^{*}\rangle+\sum_{ij}P^{*}_{ij}-R^{*}_{ij}=H(R^{*}|P^{*}).$

The result follows. ∎

######  Remark 20.

In Step 2, we used a particular case of the following more general result that
is proved in the same way:

######  Lemma 21.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu,\mu^{\prime}\in\mathcal{M}_{+}(\mathcal{D})$ and
$\nu,\nu^{\prime}\in\mathcal{M}_{+}(\mathcal{F})$. Assume that
$\operatorname{Sch}(R;\mu,\nu)$ admits a solution $P$ and that
$\operatorname{Sch}(P;\mu^{\prime},\nu^{\prime})$ admits a solution $Q$. Then
the unique solution of $\operatorname{Sch}(R;\mu^{\prime},\nu^{\prime})$
exists: it is $Q$.

###  4.2 Balanced version

In the last subsection, we interpreted the fact that
$\operatorname{Sch}(R;\mu,\nu)$ has no solution by the fact that our model
does not incorporate the ability of the real system to create or destroy mass.
In that case, the total mass of $R^{*}$ is not the same as the one of $\mu$
and $\nu$ in general, even when the latter two coincide. Therefore, $R^{*}$
cannot be interpreted directly as a joint law for the initial and final
positions of the particles. Following the lines of [20], we see that its
interpretation is actually rather complicated.

In this subsection, we want to consider the case where the real system under
study is truly balanced, that is, no creation of annihilation of mass is
possible at all. In this situation, whatever the way we are obtaining the
data, $\mu$ and $\nu$ must have the same mass, and up to renormalizing, we can
assume that they are probability measures. We want to interpret the fact that
$\operatorname{Sch}(R;\mu,\nu)$ has no solution by the fact that $\mu$ and
$\nu$ are imperfect measurements of the true marginals, and we want to find a
_probability_ measure $\bar{R}^{*}$ that is entropically close to $R$ while
having its marginals entropically close to $\mu$ and $\nu$, that can be
interpreted as a joint law.

Therefore, we introduce the following problem that is a slight modification of
$\operatorname{Sch}^{\varepsilon}$ where the competitor $\bar{R}$ needs to be
a probability measure: for all
$R\in\mathcal{P}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{P}(\mathcal{D})$ and $\nu\in\mathcal{P}(\mathcal{F})$,

$\overline{\operatorname{Sch}}{}^{\varepsilon}(R;\mu,\nu):=\min\Big{\\{}\varepsilon
H(\bar{R}|R)+H(\mu^{\bar{R}}|\mu)+H(\nu^{\bar{R}}|\nu)\,\Big{|}\,\bar{R}\in\mathcal{P}(\mathcal{D}\times\mathcal{F})\Big{\\}}.$

The following theorem states the behaviour of this optimization problem as
$\varepsilon\to 0$, and is a direct adaptation of Theorem 17 to the balanced
case.

######  Theorem 22.

Let $R\in\mathcal{P}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{P}(\mathcal{D})$ and $\nu\in\mathcal{P}(\mathcal{F})$ satisfy
the conditions of Assumption 8, and call

$\mathcal{Z}:=\sum_{ij}\sqrt{P^{*}_{ij}Q^{*}_{ij}},$

where $P^{*}$ and $Q^{*}$ are given by Theorem 11. Then for all
$\varepsilon>0$, the solution $\bar{R}^{\varepsilon}$ of
$\overline{\operatorname{Sch}}{}^{\varepsilon}(R;\mu,\nu)$ exists, is unique,
and satisfies for all $i,j$:

$\bar{R}^{\varepsilon}_{ij}\underset{\varepsilon\to
0}{\longrightarrow}\bar{R}^{*}_{ij}:=\frac{\sqrt{P^{*}_{ij}Q^{*}_{ij}}}{\mathcal{Z}}.$

Its marginals are given for all $i,j$ by

$\mu^{\bar{R}^{*}}_{i}=\frac{\sqrt{\mu^{*}_{i}\mu_{i}}}{\mathcal{Z}}\qquad\mbox{and}\qquad\nu^{\bar{R}^{*}}_{j}=\frac{\sqrt{\nu^{*}_{j}\nu_{j}}}{\mathcal{Z}}.$

###### Proof.

Theorem 22 is a direct consequence of Theorem 17 once noticed the following
fact: If $R,\mu,\nu$ are as in the statement of the theorem, if
$\varepsilon>0$ and if $R^{\varepsilon}$ is the solution of
$\operatorname{Sch}^{\varepsilon}(R;\mu,\nu)$, then
$R^{\varepsilon}/\mathsf{M}(R^{\varepsilon})$ is the solution of
$\overline{\operatorname{Sch}}{}^{\varepsilon}(R;\mu,\nu)$. To see this,
consider $R^{\prime}\in\mathcal{P}(\mathcal{D}\times\mathcal{F})$. Direct
computations imply

$\displaystyle
H(R^{\prime}|R)=\frac{H\Big{(}\mathsf{M}(R^{\varepsilon})R^{\prime}\Big{|}R\Big{)}}{\mathsf{M}(R^{\varepsilon})}+\log\frac{1}{\mathsf{M}(R^{\varepsilon})}+1-\frac{1}{\mathsf{M}(R^{\varepsilon})},$
$\displaystyle
H\bigg{(}\frac{R^{\varepsilon}}{\mathsf{M}(R^{\varepsilon})}\bigg{|}R\bigg{)}=\frac{H(R^{\varepsilon}|R)}{\mathsf{M}(R^{\varepsilon})}+\log\frac{1}{\mathsf{M}(R^{\varepsilon})}+1-\frac{1}{\mathsf{M}(R^{\varepsilon})}.$

By optimality of $R^{\varepsilon}$,
$H(\mathsf{M}(R^{\varepsilon})R^{\prime}|R)\geq H(R^{\varepsilon}|R)$, and
therefore $H(R^{\prime}|R)\geq
H(R^{\varepsilon}/\mathsf{M}(R^{\varepsilon})|R)$. Our claims follows, and
hence the theorem as $\mathsf{M}$ is a continuous functional and
$\mathcal{Z}=\mathsf{M}(R^{*})$, where $R^{*}$ is given by Theorem 17. ∎

##  5 Existence and support of the solutions to Schrödinger problems

In this section, our goal is to give a detailed study of the support of the
solution of $\operatorname{Sch}(R;\mu,\nu)$ when the latter exists, or of the
common one of $P^{*}$, $Q^{*}$ and $R^{*}$ from Theorems 11 and 17 in the non-
scalable case. This study will rely on a new interpretation of the well known
existence conditions for the Schrödinger problem in finite spaces, for which
we refer to [4, 14].

We start with our new formulation of these conditions of existence, which is
very close to the ones introduced by Brualdi [4], but has the advantage of
helping understanding the shape of the support of the optimizers seen as a
bipartite graph.

In the second part of the section, we provide a theoretical procedure allowing
to get the support of the optimizers, both in the approximately scalable and
non-scalable cases, without using the Sinkhorn algorithm. This procedure will
be used in the next section as a preliminary step, before launching the
Sinkhorn algorithm, in order to recover a linear rate for the latter.

###  5.1 A necessary and sufficient condition of existence for the Schrödinger
problem in finite spaces

Let us state a necessary and sufficient condition on $R$, $\mu$ and $\nu$ for
the existence of a solution $R^{*}$ of $\operatorname{Sch}(R;\mu,\nu)$, that
is, for $\operatorname{Sch}(R;\mu,\nu)$ to be scalable or approximately
scalable. In order to do so, we need to give a few definitions. First, we
endow the set $\mathcal{D}\cup\mathcal{F}$ with a bipartite graph structure
related to $R$: we set

$\forall i=1,\dots,N\mbox{ and }j=1,\dots,M,\qquad x_{i}\triangle
y_{j}\Leftrightarrow R_{ij}>0.$

We have $x_{i}\triangle y_{j}$ whenever it is possible to travel from $x_{i}$
to $y_{j}$ under $R$. We write indifferently $x_{i}\triangle y_{j}$ or
$y_{j}\triangle x_{i}$.

With this structure in hand, we are able to push forward or pull backward
subsets of $\mathcal{D}$ and $\mathcal{F}$, that is, we define:

$\displaystyle\forall A\subset\mathcal{D},$ $\displaystyle F_{R}(A)$
$\displaystyle:=\Big{\\{}y\in\mathcal{F}\,|\,\exists x\in A\mbox{ s.t.\
}x\triangle y\Big{\\}},$ (28) $\displaystyle\forall B\subset\mathcal{F},$
$\displaystyle D_{R}(B)$ $\displaystyle:=\Big{\\{}x\in\mathcal{D}\,|\,\exists
y\in B\mbox{ s.t.\ }x\triangle y\Big{\\}}.$

Heuristically, for all $A\subset\mathcal{D}$, $F_{R}(A)$ is the set of all
possible final positions of particles starting from $A$, under $R$.
Correspondingly, for all $B\subset\mathcal{F}$, $D_{R}(B)$ is the set of all
possible initial positions of particles arriving in $B$ under $R$. Notice the
explicit mention of $R$ in the notations: in the following, we will allow
ourselves to replace $R$ by any other measure
$\bar{R}\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$.

The main result of this section is the following.

######  Theorem 23.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$.
The three following assertions are equivalent:

  1. (a)

$\mathsf{M}(\mu)=\mathsf{M}(\nu)$ and for all $A\subset\mathcal{D}$,
$\mu(A)\leq\nu(F_{R}(A))$.

  2. (b)

$\mathsf{M}(\mu)=\mathsf{M}(\nu)$ and for all $B\subset\mathcal{F}$,
$\nu(B)\leq\mu(D_{R}(B))$.

  3. (c)

$\operatorname{Sch}(R;\mu,\nu)$ is scalable or approximately scalable.

Note that the implications (c) $\Rightarrow$ (a) and (c) $\Rightarrow$ (b) are
straightforward, and that only the reverse implications are challenging. Also,
we already noticed in Subsection 2.3 that (c) implies Assumption 8. Hence, it
is also the case for (a) and (b).

The proof relies on the following Lemma 24, which gives a necessary and
sufficient condition on $R$, $\mu$ and $\nu$ ensuring $R^{*}$ to have the same
support as $R$, that is, to be in the scalable case. In this statement, we use
the notations $\mu^{R}$ and $\nu^{R}$ as defined in (3), and we work under
Assumption 9, which is always possible under Assumption 8 up to considering
subspaces of $\mathcal{D}$ and $\mathcal{F}$, see Subsection 2.3.

######  Lemma 24.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$,
satisfying Assumption 9. The three following assertions are equivalent:

  1. (a’)

$\mathsf{M}(\mu)=\mathsf{M}(\nu)$ and for all $A\subset\mathcal{D}$,
$\mu(A)\leq\nu(F_{R}(A))$, with a strict inequality whenever
$\mu^{R}(A)<\nu^{R}(F_{R}(A))$.

  2. (b’)

$\mathsf{M}(\mu)=\mathsf{M}(\nu)$ and for all $B\subset\mathcal{F}$,
$\nu(B)\leq\mu(D_{R}(B))$, with a strict inequality whenever
$\nu^{R}(B)<\mu^{R}(D_{R}(B))$.

  3. (c’)

$\operatorname{Sch}(R;\mu,\nu)$ is scalable.

In plain words, it highlights the difference between the approximately
scalable and scalable cases, by showing that the scalable case consists in
assuming as much strict inequalities in (a) or in (b) as possible. Although
both Theorem 23 and Lemma 24 can be directly deduced from the work of Brualdi
[4], we provide in Appendix B a short and independent proof based on
topological arguments.

###  5.2 Theoretical construction of the support

In the scalable case, the Sinkhorn algorithm is known to have a linear rate of
convergence. On the other hand, in the approximately scalable case, the
algorithm still converges, but the (unknown) convergence rate cannot be linear
[1].

In this subsection, we study the support of the solution of the Schrödinger
problem in the approximately scalable and non-scalable cases for the following
reason. Take $R$, $\mu$ and $\nu$ such that $\operatorname{Sch}(R;\mu,\nu)$ is
approximately scalable, $R^{*}$ the solution of this problem, and
$\mathcal{S}$ the support of $R^{*}$. Then
$\operatorname{Sch}(\mathds{1}_{\mathcal{S}}R;\mu,\nu)$ is scalable and its
solution is $R^{*}$. In particular, the Sinkhorn algorithm applied to this
problem has a linear rate of convergence. Interestingly, a similar reasoning
is valid in the non-scalable case, as we show in Proposition 25 below.

Without loss of generality, and for the sake of simplicity, in the whole
subsection, we work under Assumption 9. By Remark 13, if $\mu^{*}$ and
$\nu^{*}$ are defined by (13), we have $\nu\sim\nu^{*}$ and $\mu\sim\mu^{*}$.
So under Assumption 9, they have a full support as well.

######  Proposition 25.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$
satisfying Assumption 9. Let us call $\mathcal{S}$ the common support of
$P^{*}$ and $Q^{*}$ from Theorem 11, and $R^{*}$ from Theorem 17. Let
$(P^{n})_{n\in\mathbb{N}^{*}}$ and $(Q^{n})_{n\in\mathbb{N}^{*}}$ be given by
(2) applied to $\operatorname{Sch}(\mathds{1}_{\mathcal{S}}R;\mu,\nu)$. They
converge respectively towards $P^{*}$ and $Q^{*}$, both of them at a linear
rate.

###### Proof.

Let $(P^{n})_{n\in\mathbb{N}^{*}}$ and $(Q^{n})_{n\in\mathbb{N}^{*}}$ be given
by the equivalent formulations (1) and (2) applied to
$\operatorname{Sch}(\mathds{1}_{\mathcal{S}}R;\mu,\nu)$. Let us show that
$(P^{n})$ converges towards $P^{*}$ at a linear rate. The case of $(Q^{n})$
follows the same arguments. The idea is that if
$(\tilde{P}^{n})_{n\in\mathbb{N}^{*}}$ and
$(\tilde{Q}^{n})_{n\in\mathbb{N}^{*}}$ are given by (1) and (2) applied to
$\operatorname{Sch}(\mathds{1}_{\mathcal{S}}R;\mu,\nu^{*})$, then for all
$n\in\mathbb{N}^{*}$, $P^{n}=\tilde{P}^{n}$. As the problem
$\operatorname{Sch}(\mathds{1}_{\mathcal{S}}R;\mu,\nu^{*})$ is scalable (its
solution, $P^{*}$, has the same support as $\mathds{1}_{\mathcal{S}}R$), the
rate of convergence of $(\tilde{P}^{n})$ towards $P^{*}$ is linear, and the
result follows.

So let us prove by induction that for all $n\in\mathbb{N}^{*}$,
$P_{n}=\tilde{P}_{n}$. According to (1), $P^{1}$ and $\tilde{P}^{1}$ are
solutions to the same problem, and therefore coincide. Let us now consider
$n\in\mathbb{N}^{*}$ such that $P^{n}=\tilde{P}^{n}$ and show that
$P^{n+1}=\tilde{P}^{n+1}$. By construction, the support of $P^{n+1}$ and
$\tilde{P}^{n+1}$ is $\mathcal{S}$, so we just need to check that for all
$(x_{j},y_{j})\in\mathcal{S}$, $P^{n+1}_{ij}=\tilde{P}^{n+1}_{ij}$. By the
first line of (26), for all $(x_{i},y_{j})\in\mathcal{S}$, we have:

$\nu_{j}=\frac{\mu^{*}_{i}\nu^{*}_{j}}{\mu_{i}}.$

Hence, for all $(x_{i},y_{j})\in\mathcal{S}$:

$\displaystyle
P^{n+1}_{ij}=\frac{\mu_{i}}{\mu^{Q^{n}}_{i}}Q^{n}_{ij}=\frac{\mu_{i}}{\displaystyle\sum_{j^{\prime}}\frac{\nu_{j^{\prime}}}{\nu^{P^{n}}_{j^{\prime}}}P^{n}_{ij^{\prime}}}\times\frac{\nu_{j}}{\nu^{P^{n}}_{j}}P^{n}_{ij}=\frac{\mu_{i}}{\displaystyle\frac{\mu_{i}^{*}}{\mu_{i}}\sum_{j^{\prime}}\frac{\nu^{*}_{j^{\prime}}}{\nu^{\tilde{P}^{n}}_{j^{\prime}}}\tilde{P}^{n}_{ij^{\prime}}}\times\frac{\mu_{i}^{*}}{\mu_{i}}\frac{\nu^{*}_{j}}{\nu^{\tilde{P}^{n}}_{j}}\tilde{P}^{n}_{ij}=\frac{\mu_{i}}{\mu^{\tilde{Q}^{n}}_{i}}\tilde{Q}^{n}_{ij}=\tilde{P}^{n+1}_{ij},$

where the change from $P^{n}$ to $\tilde{P}^{n}$ in the middle coming from the
induction assumption $P^{n}=\tilde{P}^{n}$. The result follows. ∎

Therefore, even in the non-scalable case, a way to improve the Sinkhorn
algorithm consists in first finding $\mathcal{S}$, and then computing the
solution of a scalable problem. We propose in this subsection a theoretical
procedure allowing to get this support without using the Sinkhorn algorithm in
both the approximately and non-scalable cases, and we will propose an
approximate method for achieving this task numerically at Section 6.

To detail our procedure, we introduce a class of subsets of $\mathcal{D}$
associated with a triple $(R;\mu,\nu)$.

######  Definition 26.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$
satisfying Assumption 9. Let us consider $R^{*}$ from Theorem 17. We say that
a subset $A\subset\mathcal{D}$ is the _source of an isolated scalable problem_
(or for short that $A$ is a SISP set) for $(R;\mu,\nu)$ if $A\neq\emptyset$
and:

  * •

The set $(\mathcal{D}\backslash A)\times F_{R}(A)$ is $R^{*}$-negligible,
_i.e._

$R^{*}\Big{(}(\mathcal{D}\backslash A)\times F_{R}(A)\Big{)}=0.$ (29)

  * •

For all $x_{i}\in A$ and $y_{j}\in\mathcal{F}$,

$R^{*}_{ij}>0\quad\Leftrightarrow\quad R_{ij}>0.$ (30)

We show at Figure 2 an illustration of what a SISP is.

Of course, as $P^{*}$ and $Q^{*}$ from Theorem 11 are equivalent to $R^{*}$ in
the sense of measures, we could have replaced $R^{*}$ in the previous
definition by one of them.

Figure 2:  If $A\subset\mathcal{D}$, up to reordering the lines, we can assume
that it corresponds to the last lines. Up to reordering the columns, we can
assume that $F_{R}(A)$ corresponds to the last columns. Then, $R$ has the form
given in the picture. In this situation, $A$ is a SISP set for $(R;\mu,\nu)$
if $R^{*}$ cancels on the block $\mathfrak{B}$ and if the supports of $R$ and
$R^{*}$ coincide on the block $\mathfrak{C}$.

On the one hand, SISP sets always exist, at least under Assumption 9, as
announced in the following lemma. Its proof is our main task in this part of
our work, and is given at the end of the subsection.

######  Lemma 27.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$
satisfying Assumption 9. Then there exists a SISP set for $(R;\mu,\nu)$.

On the other hand, once we know how to find SISP sets, an iterative procedure
consisting in finding SISP sets for a sequence of more and more restricted
problems makes is possible to reconstruct the whole subset $\mathcal{S}$.

######  Proposition 28.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$
satisfying Assumption 9. Let us call $\mathcal{S}$ the common support of
$P^{*}$ and $Q^{*}$ from Theorem 11, and $R^{*}$ from Theorem 17.

We define by inference $(R^{n})_{n\in\mathbb{N}}$ a sequence in
$\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$(\mathcal{D}^{n})_{n\in\mathbb{N}}$ a nonincreasing sequence of subsets of
$\mathcal{D}$ and $(\mathcal{F}^{n})_{n\in\mathbb{N}}$ a nonincreasing
sequence of subsets of $\mathcal{F}$ in the following way:

  * •

For $n=0$, we set $R^{0}:=R$, $\mathcal{D}^{0}:=\mathcal{D}$ and
$\mathcal{F}^{0}:=\mathcal{F}$;

  * •

For all $n\in\mathbb{N}$, if $\mathcal{D}^{n}$ and $\mathcal{F}^{n}$ are
nonempty and
$(R^{n}\llcorner_{\mathcal{D}^{n}\times\mathcal{F}^{n}};\mu\llcorner_{\mathcal{D}^{n}},\nu\llcorner_{\mathcal{F}^{n}})$
satisfies Assumption 9, we pick $M_{n}$ a SISP set as given by Lemma 27, and
we set:

$\displaystyle\mathcal{D}^{n+1}:=\mathcal{D}^{n}\backslash
M_{n},\qquad\mathcal{F}^{n+1}:=\mathcal{F}^{n}\backslash
F_{R^{n}\llcorner_{\mathcal{D}^{n}\times\mathcal{F}^{n}}}(M_{n}),$
$\displaystyle\forall i,j,\quad R^{n+1}_{ij}:=\left\\{\begin{aligned}
&0,&&\mbox{if }y_{j}\in
F_{R^{n}\llcorner_{\mathcal{D}^{n}\times\mathcal{F}^{n}}}(M_{n})\mbox{ and
}x_{i}\in\mathcal{D}^{n+1},\\\\[5.0pt]
&R^{n}_{ij},&&\mbox{otherwise}.\end{aligned}\right.$

Otherwise, we set $R^{n+1}:=R^{n}$, $\mathcal{D}^{n+1}:=\mathcal{D}^{n}$ and
$\mathcal{F}^{n+1}:=\mathcal{F}^{n}$.

With this construction, the sequence
$(R^{n},\mathcal{D}^{n},\mathcal{F}^{n})_{n\in\mathbb{N}}$ is stationary. More
precisely, there exists $N\in\mathbb{N}^{*}$ such that for all $n\geq N$,

$\mathcal{D}^{n}=\emptyset,\qquad\mathcal{F}^{n}=\emptyset,\qquad
R^{n}=\mathds{1}_{\mathcal{S}}R.$

Figure 3:  In the situation of Figure 2 where we have reordered the lines and
columns, our procedure consists in recursively add zeros to $R$ at positions
where we know thanks to (29) that $R^{*}$ admits a zero. We know that we did
not forget any zero in $M\times F_{R}(M)$ thanks to (35).

An illustration of the procedure at each iteration, is provided in Figure 3.
An illustration of the full procedure in a specific non-scalable case is
provided in Figure 4.

Figure 4: Illustration of the procedure of Proposition 28 when the matrix $R$
is upper diagonal and $R^{*}$ a staircase matrix (see Appendix A for more
details). Only $A_{0}$ is a SISP set for $(R;\mu,\nu)$, but $A_{1}$ and
$A_{2}$ are SISP sets for the restricted problems at the iterations 2 and 3.
In this example, the procedure is stationary after 3 steps. In this example,
the SISP set at each iteration is the unique maximal $\theta$-set existing for
the reduced problem, as presented in Definition 30. We remark that we can also
build $\nu^{*}$, the second marginal of $P^{*}$ defined in Theorem 11, on the
successive SISP sets obtained along the procedure, thanks to the second step
of the proof of Proposition 27 which ensures that the ratio
$\frac{\nu}{\nu^{*}}$ is constant inside the maximal $\theta$-sets.

###### Proof.

In this proof, in order to lighten the notations, we call
$R^{n}_{r}:=R^{n}\llcorner_{\mathcal{D}^{n}\times\mathcal{F}^{n}}$. We will
prove by inference the following facts. For all $n\in\mathbb{N}$:

  1. 1.

Calling $\mathcal{S}^{n}$ the support of $R^{n}$, and therefore
$\mathcal{S}^{0}$ the support of $R$,

$\mathcal{S}^{n}\cap\Big{(}(\mathcal{D}\times\mathcal{F})\backslash(\mathcal{D}^{n}\times\mathcal{F}^{n})\Big{)}=\mathcal{S}\cap\Big{(}(\mathcal{D}\times\mathcal{F})\backslash(\mathcal{D}^{n}\times\mathcal{F}^{n})\Big{)}=\mathcal{S}^{0}\cap\left(\bigcup_{k=0}^{n-1}M_{k}\times
F_{R^{k}_{r}}(M_{k})\right),$ (31)

$\mathcal{S}\subset\mathcal{S}^{n}$, and
$\mathcal{S}^{n}\cap(\mathcal{D}^{n}\times\mathcal{F}^{n})=\mathcal{S}^{0}\cap(\mathcal{D}^{n}\times\mathcal{F}^{n})$.

  2. 2.

$\mathcal{D}^{n}$ is empty if and only if $\mathcal{F}^{n}$ is empty.

  3. 3.

If $\mathcal{D}^{n}$ and $\mathcal{F}^{n}$ are not empty,
$(R^{n}_{r};\mu\llcorner_{\mathcal{D}^{n}},\nu\llcorner_{\mathcal{F}^{n}})$
satisfies Assumption 9 and the matrices $P^{n,*}$, $Q^{n,*}$ and $R^{n,*}$
associated with
$(R^{n}_{r};\mu\llcorner_{\mathcal{D}^{n}},\nu\llcorner_{\mathcal{F}^{n}})$
through Theorems 11 and 17 are the restrictions of $P^{*}$, $Q^{*}$ and
$R^{*}$ to $\mathcal{D}^{n}\times\mathcal{F}^{n}$.

This is enough to prove the proposition: if the conclusion of the inference is
true, then by the third point and Lemma 27, as long as $\mathcal{D}^{n}$ and
$\mathcal{F}^{n}$ are nonempty,
$(R^{n}_{r};\mu\llcorner_{\mathcal{D}^{n}},\nu\llcorner_{\mathcal{F}^{n}})$
admits a SISP set $M_{n}$, which is not empty by definition. Therefore,
$(\mathcal{D}^{n})$ is strictly decreasing in the sense of inclusion as long
as it is not empty, so it has to reach $\emptyset$ at a certain rank $N$. At
this rank, because of the first point, we also have
$\mathcal{F}^{N}=\emptyset$, and because of (31),
$\mathcal{S}^{N}=\mathcal{S}$, so that the conclusion follows. So let us prove
the inference.

At rank $0$, everything is clear, so let us assume that the conclusions of
points one, two and three hold at rank $n$, and prove them at rank $n+1$.
First, if $\mathcal{D}^{n}$ is empty, by assumption $\mathcal{F}^{n}$ is empty
as well, so we have reached a stationary point, and everything is still true
at rank $n+1$. So we can assume without loss of generality that
$\mathcal{D}^{n}$, and hence $\mathcal{F}^{n}$, are nonempty. By assumption,
$(R^{n}_{r};\mu\llcorner_{\mathcal{D}^{n}},\nu\llcorner_{\mathcal{F}^{n}})$
satisfies Assumption 9, and by Lemma 27, we can find a SISP set $M_{n}$. In
this context, let us check the points one by one at rank $n+1$.

First point. Observing that $\mathcal{D}^{n}$ is the disjoint union of $M_{n}$
and $\mathcal{D}^{n+1}$, and that $\mathcal{F}^{n}$ is the disjoint union of
$F_{R^{n}_{r}}(M_{n})$ and $\mathcal{F}^{n+1}$, we have

$\displaystyle(\mathcal{D}\times\mathcal{F})$
$\displaystyle\backslash(\mathcal{D}^{n+1}\times\mathcal{F}^{n+1})$
$\displaystyle=\Big{(}(\mathcal{D}\times\mathcal{F})\backslash(\mathcal{D}^{n}\times\mathcal{F}^{n})\Big{)}\cup\big{(}\mathcal{D}^{n+1}\times
F_{R^{n}_{r}}(M_{n})\big{)}\cup\big{(}M_{n}\times\mathcal{F}^{n+1}\big{)}\cup\big{(}M_{n}\times
F_{R^{n}_{r}}(M_{n})\big{)}.$

So in order to prove (31) at rank $n+1$, we need to show that

$\displaystyle\mathcal{S}^{n+1}\cap\Big{(}(\mathcal{D}\times\mathcal{F})\backslash(\mathcal{D}^{n}\times\mathcal{F}^{n})\Big{)}=\mathcal{S}^{n}\cap\Big{(}(\mathcal{D}\times\mathcal{F})\backslash(\mathcal{D}^{n}\times\mathcal{F}^{n})\Big{)},$
(32)
$\displaystyle\mathcal{S}^{n+1}\cap\big{(}M_{n}\times\mathcal{F}^{n+1}\big{)}=\mathcal{S}\cap\big{(}M_{n}\times\mathcal{F}^{n+1}\big{)}=\emptyset,$
(33) $\displaystyle\mathcal{S}^{n+1}\cap\big{(}\mathcal{D}^{n+1}\times
F_{R^{n}_{r}}(M_{n})\big{)}=\mathcal{S}\cap\big{(}\mathcal{D}^{n+1}\times
F_{R^{n}_{r}}(M_{n})\big{)}=\emptyset,$ (34)
$\displaystyle\mathcal{S}^{n+1}\cap\big{(}M_{n}\times
F_{R^{n}_{r}}(M_{n})\big{)}=\mathcal{S}\cap\big{(}M_{n}\times
F_{R^{n}_{r}}(M_{n})\big{)}=\mathcal{S}^{0}\cap\big{(}M_{n}\times
F_{R^{n}_{r}}(M_{n})\big{)}.$ (35)

To prove these equalities, the main tool is the following formula which is a
direct consequence of the construction:

$\mathcal{S}^{n+1}=\mathcal{S}^{n}\backslash\big{(}\mathcal{D}^{n+1}\times
F_{R^{n}_{r}}(M_{n})\big{)}.$ (36)

With this formula at hand, we see that (32) follows from
$\mathcal{D}^{n+1}\times
F_{R^{n}_{r}}(M_{n})\subset\mathcal{D}^{n}\times\mathcal{F}^{n}$. We also
deduce very easily that
$\mathcal{S}^{n+1}\cap(\mathcal{D}^{n+1}\times\mathcal{F}^{n+1})=\mathcal{S}^{n}\cap(\mathcal{D}^{n+1}\times\mathcal{F}^{n+1})=\mathcal{S}^{0}\cap(\mathcal{D}^{n+1}\times\mathcal{F}^{n+1})$,
where the last equality follows from the first point at rank $n$.

Then, to prove (33), as both $\mathcal{S}$ (by assumption) and
$\mathcal{S}^{n+1}$ (by (36)) are included in $\mathcal{S}^{n}$, it suffices
to show that $\mathcal{S}^{n}\cap(M_{n}\times\mathcal{F}^{n+1})=\emptyset$.
But that last assertion follows from the definition of
$\mathcal{F}^{n+1}=\mathcal{F}^{n}\backslash F_{R^{n}_{r}}(M_{n})$: these are
precisely the columns where $R^{n}_{r}$ has only zero entries on the
intersection with the lines $M_{n}$.

To prove (34), let us observe that the equality
$\mathcal{S}^{n+1}\cap(\mathcal{D}^{n+1}\times
F_{R^{n}_{r}}(M_{n}))=\emptyset$ is a direct consequence of (36). The other
equality, namely, $\mathcal{S}\cap(\mathcal{D}^{n+1}\times
F_{R^{n}_{r}}(M_{n}))=\emptyset$ follows from the fact that $M_{n}$ is a SISP
set for
$(R^{n}_{r};\mu\llcorner_{\mathcal{D}^{n}},\nu\llcorner_{\mathcal{F}^{n}})$,
so that (29) applies with $M_{n}$ instead of $A$, $R^{n}_{r}$ instead of $R$,
$\mathcal{D}^{n}$ instead of $\mathcal{D}$ and
$R^{n,*}=R^{*}\llcorner_{\mathcal{D}^{n}\times\mathcal{F}^{n}}$ instead of
$R^{*}$ (here, we use the point three at rank $n$). Notice that as
$\mathcal{S}\cap(\mathcal{D}^{n+1}\times F_{R^{n}_{r}}(M_{n}))=\emptyset$ and
$\mathcal{S}\subset\mathcal{S}^{n}$, by (36), we have also proved that
$\mathcal{S}\subset\mathcal{S}^{n+1}$.

Finally, to prove (35), as
$\mathcal{S}\subset\mathcal{S}^{n+1}\subset\mathcal{S}^{0}$, we just need to
prove that $\mathcal{S}\cap(M_{n}\times
F_{R^{n}_{r}}(M_{n}))=\mathcal{S}^{0}\cap(M_{n}\times F_{R^{n}_{r}}(M_{n}))$.
But this is a direct consequence of the fact that $M_{n}$ is a SISP set for
$(R^{n}_{r};\mu\llcorner_{\mathcal{D}^{n}},\nu\llcorner_{\mathcal{F}^{n}})$,
so that (30) applies with $R^{n}_{r}$ instead of $R$, $R^{n,*}$ instead of
$R^{*}$, $M_{n}$ instead of $A$ and $\mathcal{F}^{n}$ instead of
$\mathcal{F}$.

Second point. By definition, $\mathcal{D}^{n+1}$ is empty if and only if
$M_{n}=\mathcal{D}^{n}$. So if $\mathcal{D}^{n+1}=\emptyset$, then by
Assumption 9 applied to
$(R^{n}_{r};\mu\llcorner_{\mathcal{D}^{n}},\nu\llcorner_{\mathcal{F}^{n}})$,
we clearly have $F^{R^{n}_{r}}(M_{n})=\mathcal{F}^{n}$ and so
$\mathcal{F}^{n+1}=\emptyset$. On the other hand, if
$\mathcal{D}^{n+1}\neq\emptyset$ we have

$\displaystyle 0<\mu(\mathcal{D}^{n+1})$
$\displaystyle=P^{*}\Big{(}\mathcal{D}^{n+1}\times\mathcal{F}\Big{)}=P^{*}\Big{(}\mathcal{D}^{n+1}\times\mathcal{F}^{n}\Big{)}$
$\displaystyle=P^{*}\Big{(}\mathcal{D}^{n+1}\times
F_{R^{n}_{r}}(M_{n})\Big{)}+P^{*}\Big{(}\mathcal{D}^{n+1}\times\mathcal{F}^{n+1}\Big{)}$
$\displaystyle=P^{*}\Big{(}\mathcal{D}^{n+1}\times\mathcal{F}^{n+1}\Big{)},$

where the inequality comes from Assumption 9, the second equality is an easy
consequence of (31) at rank $n$, and the last one comes from the definition of
SISP sets (see (29)) and from the fact that $P^{*}$ and $R^{*}$ has the same
support. So $\mathcal{F}^{n+1}$ cannot be empty, as announced.

Third point. To check that
$(R^{n+1}_{r};\mu\llcorner_{\mathcal{D}^{n+1}},\nu\llcorner_{\mathcal{F}^{n+1}})$
satisfies Assumption 9, we need only need to show that the support of
$\mu^{R^{n+1}}\llcorner_{\mathcal{D}^{n+1}}$ is $\mathcal{D}^{n+1}$ and the
support of $\nu^{R^{n+1}}\llcorner_{\mathcal{F}^{n+1}}$ is
$\mathcal{F}^{n+1}$. But as we already proved that
$\mathcal{S}\subset\mathcal{S}^{n+1}$, we know that $P^{*}\ll R^{n+1}$ and
$Q^{*}\ll R^{n+1}$, so the conclusion follows from the fact that $\mu$ and
$\nu$ have full support by Assumption 9 applied to $(R;\mu,\nu)$. The last
thing to check, namely that the matrices $P^{n+1,*}$, $Q^{n+1,*}$ and
$R^{n+1,*}$ associated with
$(R^{n+1}_{r};\mu\llcorner_{\mathcal{D}^{n+1}},\nu\llcorner_{\mathcal{F}^{n+1}})$
through Theorems 11 and 17 are the restrictions of $P^{*}$, $Q^{*}$ and
$R^{*}$ to $\mathcal{D}^{n+1}\times\mathcal{F}^{n+1}$ is a direct consequence
of Proposition 29 below (that we wanted to separate to the rest of the proof
because we will use it again later), and of (34). ∎

######  Proposition 29.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$
satisfy Assumption 8. Let $P^{*}$, $Q^{*}$ and $R^{*}$ be the matrices
associated with the problem $\operatorname{Sch}(R;\mu,\nu)$ by Theorems 11 and
17. Finally, let $A\subset\mathcal{D}$ be such that

$R^{*}\Big{(}(\mathcal{D}\backslash A)\times F_{R}(A)\Big{)}=0.$ (37)

Then (with slightly sloppy notations), $P^{*}\llcorner_{A\times F_{R}(A)}$,
$Q^{*}\llcorner_{A\times F_{R}(A)}$ and $R^{*}\llcorner_{A\times F_{R}(A)}$
are the matrices associated with the restricted problem
$\operatorname{Sch}(R\llcorner_{A\times
F_{R}(A)},\mu\llcorner_{A},\nu\llcorner_{F_{R}(A)})$ by Theorem 11 and 17.

Similarly, calling $A^{\prime}:=\mathcal{D}\backslash A$ and
$F^{\prime}:=\mathcal{F}\backslash F_{R}(A)$,
$P^{*}\llcorner_{A^{\prime}\times F^{\prime}}$,
$Q^{*}\llcorner_{A^{\prime}\times F^{\prime}}$ and
$R^{*}\llcorner_{A^{\prime}\times F^{\prime}}$ are the matrices associated
with the restricted problem $\operatorname{Sch}(R\llcorner_{A^{\prime}\times
F^{\prime}},\mu\llcorner_{A^{\prime}},\nu\llcorner_{F^{\prime}})$ by Theorem
11 and 17.

###### Proof.

We show the result in the case of $P^{*}$, related to the ”block” $A\times
F_{R}(A)$. The case of $Q^{*}$ is similar, the case of $R^{*}$ easily follows
from the two previous ones, and the similar results on $A^{\prime}\times
F^{\prime}$ follow the same lines. Let $\nu^{*}$ be defined by (13). The first
thing to prove is

$\displaystyle\nu^{*}\llcorner_{F_{R}(A)}$ $\displaystyle\hskip
5.0pt=\operatorname*{arg\,min}\Big{\\{}H(\bar{\nu}|\nu\llcorner_{F_{R}(A)})\,\Big{|}\,\bar{\nu}=\nu^{P}\mbox{
for some }P\mbox{ with }H(P|R_{A\times F_{R}(A)})<+\infty\mbox{ and
}\mu^{P}=\mu\llcorner_{A}\Big{\\}}.$

The measure $\nu^{*}\llcorner_{F_{R}(A)}$ is a competitor for the problem in
the r.h.s. because it corresponds to $P:=P^{*}\llcorner_{A\times F_{R}(A)}$.
Let us show that it is the optimizer. To do this, we call $\bar{\nu}$ the
optimizer, and we show that $\nu^{*}_{F_{R}(A)}=\bar{\nu}^{*}$. Let us
consider $\bar{P}$ a $P$ corresponding to $\bar{\nu}^{*}$ in the problem
above, $\bar{P}^{*}$ the matrix obtained by replacing the entries of $P^{*}$
on $A\times F_{R}(A)$ by the entries of $\bar{P}$, and
$\bar{\nu}^{*}:=\nu^{\bar{P}^{*}}$. We have

$H(\bar{\nu}^{*}|\nu)=H(\bar{\nu}|\nu\llcorner_{F_{R}(A)})+H(\nu^{*}\llcorner_{F^{\prime}}|\nu\llcorner_{F^{\prime}})\leq
H(\nu^{*}\llcorner_{F_{R}(A)}|\nu\llcorner_{F_{R}(A)})+H(\nu^{*}\llcorner_{F^{\prime}}|\nu\llcorner_{F^{\prime}})=H(\nu^{*}|\nu),$

where the inequality, being a consequence of the optimality of
$\bar{\nu}^{*}$, is an equality if and only if
$\nu^{*}_{F_{R}(A)}=\bar{\nu}^{*}$. But by optimality of $\nu^{*}$ in (13)
this inequality is indeed an equality, and therefore
$\nu^{*}_{F_{R}(A)}=\bar{\nu}^{*}$.

It remains to show that $P^{*}\llcorner_{A\times F_{R}(A)}$ is the solution of
$\operatorname{Sch}(R\llcorner_{A\times
F_{R}(A)},\mu\llcorner_{A},\nu^{*}\llcorner_{F_{R}(A)})$. For this, let us
consider $\bar{P}$ the solution of $\operatorname{Sch}(R\llcorner_{A\times
F_{R}(A)},\mu\llcorner_{A},\nu^{*}\llcorner_{F_{R}(A)})$, and $\bar{P}^{*}$
the matrix obtained by replacing the entries of $P^{*}$ on $A\times F_{R}(A)$
by the entries of $\bar{P}$. Because of (37), we have

$\displaystyle H(\bar{P}^{*}|R)$ $\displaystyle=H(\bar{P}|R\llcorner_{A\times
F_{R}(A)})+H(P^{*}\llcorner_{A^{\prime}\times
F^{\prime}}|R\llcorner_{A^{\prime}\times F^{\prime}})$ $\displaystyle\leq
H(P^{*}\llcorner_{A\times F_{R}(A)}|R\llcorner_{A\times
F_{R}(A)})+H(P^{*}\llcorner_{A^{\prime}\times
F^{\prime}}|R\llcorner_{A^{\prime}\times F^{\prime}})=H(P^{*}|R^{*}),$

where the inequality is a consequence of the optimality of $\bar{P}$, and is
an equality if and only if $\bar{P}=P^{*}\llcorner_{A\times F_{R}(A)}$. But by
optimality of $P^{*}$, this inequality is indeed an equality, and we conclude
that $\bar{P}=P^{*}\llcorner_{A\times F_{R}(A)}$. The proposition is proved. ∎

Now, we want to prove Lemma 27. To do this, we introduce a new class of
subsets of $\mathcal{D}$, associated with a triple $(R;\mu,\nu)$.

######  Definition 30.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$
satisfying Assumption 9. The maximal $\theta$ associated to $(R;\mu,\nu)$ is
defined by:

$\theta_{m}:=\max_{\begin{subarray}{c}A\subset\mathcal{D}\\\
A\neq\emptyset\end{subarray}}\frac{\mu(A)}{\nu(F_{R}(A))}.$ (38)

We say that $A\subset\mathcal{D}$ is a maximal $\theta$-set for $(R;\mu,\nu)$
if $A$ is a maximizer of (38). We say that it is a _smallest_ maximal
$\theta$-set if in addition, it is a minimal element in the sense of inclusion
among all maximal $\theta$-sets associated with $(R;\mu,\nu)$.

As maximal $\theta$-sets are optimizers of a finite function (thanks to
Assumption 9) on a finite set (the set of all nonempty subsets of
$\mathcal{D}$), any triple $(R;\mu,\nu)$ satisfying Assumption 9 admits at
least one maximal $\theta$-set. The set of all maximal $\theta$-sets being
itself finite, we know that there exists at least one minimal element in this
set, so that smallest maximal $\theta$-sets always exist under Assumption 9.
Hence, Lemma 27 is an obvious consequence of the following proposition, whose
proof heavily relies on the optimality conditions stated in Proposition 19.

######  Proposition 31.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$
satisfying Assumption 9. A smallest maximal $\theta$-set for $(R;\mu,\nu)$ is
a SISP set for $(R;\mu,\nu)$.

###### Proof.

Let $\mu^{*}$ and $\nu^{*}$ be defined by (13). We first define the two
following quantities

$\theta^{\mathcal{D}}:=\max_{i\in\mathcal{D}}\frac{\mu_{i}}{\mu^{*}_{i}},\qquad\theta^{\mathcal{F}}:=\max_{j\in\mathcal{F}}\frac{\nu_{j}^{*}}{\nu_{j}}.$

(Recall that under Assumption 9, by Remark 13, $\mu^{*}$ and $\nu^{*}$ defined
by (13) have full support.) Then, we define the two following sets, that are
nonempty subsets of $\mathcal{D}$ and $\mathcal{F}$ respectively:

$\overline{M}:=\Big{\\{}x_{i}\in\mathcal{D}\mbox{ s.t.\
}\frac{\mu_{i}}{\mu^{*}_{i}}=\theta^{\mathcal{D}}\Big{\\}},\qquad\overline{F}:=\Big{\\{}y_{j}\in\mathcal{F}\mbox{
s.t.\ }\frac{\nu^{*}_{j}}{\nu_{j}}=\theta^{\mathcal{F}}\Big{\\}}.$

The main argument of the proof consists in showing that $\theta^{\mathcal{D}}$
and $\theta^{\mathcal{F}}$ coincide with $\theta_{m}$, the maximal $\theta$
for $(R;\mu,\nu)$. Even if $\overline{M}$ is not a smallest maximal-$\theta$
set in general (more precisely, it is a maximal $\theta$-set that is not
minimal in general), we show at Step 3 below how this information allows to
conclude. As before, $P^{*}$ is the matrix defined by (14).

Step 1: $\theta^{\mathcal{D}}=\theta^{\mathcal{F}}$.

Let $x_{i}\in\overline{M}$ and $y_{j}\in\mathcal{F}$ be such that
$P^{*}_{ij}>0$ (such a $j$ exists thanks to Assumption 9). By the first line
of (26), we have

$\frac{\mu^{*}_{i}}{\mu_{i}}\frac{\nu^{*}_{j}}{\nu_{j}}=1,$ (39)

which implies that $\nu_{j}^{*}/\nu_{j}=\theta^{\mathcal{D}}$, and hence that
$\theta^{\mathcal{F}}\geq\theta^{\mathcal{D}}$. The other inequality is proved
in the same way, and the result follows. From now on, we call

$\overline{\theta}:=\theta^{\mathcal{D}}=\theta^{\mathcal{F}}.$

Step 2: $\overline{\theta}=\theta_{m}$.

First, $\overline{\theta}\geq\theta_{m}$. Indeed, for any
$A\subset\mathcal{D}$, by Theorem (23), as $\operatorname{Sch}(R;\mu,\nu^{*})$
is at least approximately scalable, we have $\mu(A)\leq\nu^{*}(F_{R}(A))$. But
on the other hand, by definition of $\overline{\theta}$, we have
$\nu^{*}\leq\overline{\theta}\nu$, so that actually,
$\mu(A)\leq\overline{\theta}\nu(F_{R}(A))$, and hence
$\overline{\theta}\geq\theta_{m}$.

Also, $\overline{\theta}\leq\theta_{m}$. To see this, let us first observe
that
$P^{*}((\mathcal{D}\backslash\overline{M})\times\overline{F})=P^{*}(\overline{M}\times(\mathcal{F}\backslash\overline{F}))=0$.
This is because if $x_{i},y_{j}$ are such that $P^{*}_{ij}>0$, still by (39),
$\mu_{i}/\mu^{*}_{i}=\overline{\theta}$ if and only if
$\nu^{*}_{j}/\nu_{j}=\overline{\theta}$, so that $x_{i}\in\overline{M}$ if and
only if $y_{j}\in\overline{F}$. Therefore, on the one hand,
$F_{R}(\overline{M})\subset\overline{F}$, and on the other hand, projecting on
both marginals:

$\displaystyle\mu(\overline{M})$
$\displaystyle=P^{*}(\overline{M}\times\mathcal{F})=P^{*}(\overline{M}\times\overline{F})+\underbrace{P^{*}(\overline{M}\times(\mathcal{F}\backslash\overline{F}))}_{=0}$
$\displaystyle=P^{*}(\overline{M}\times\overline{F})+\underbrace{P^{*}((\mathcal{D}\backslash\overline{M})\times\overline{F})}_{=0}=P^{*}(\mathcal{D}\times\overline{F})=\nu^{*}(\overline{F}).$

As by definition of $\overline{F}$ and $\overline{\theta}$,
$\nu^{*}(\overline{F})=\overline{\theta}\nu(\overline{F})$, we conclude that

$\overline{\theta}\nu(F_{R}(\overline{M}))\leq\overline{\theta}\nu(\overline{F})=\mu(\overline{M}),$

so that $\overline{\theta}\leq\theta_{m}$, as announced.

Step 3: Conclusion.

We are now in position to conclude. Let $A$ be a smallest maximal
$\theta$-set. As $\operatorname{Sch}(R;\mu,\nu^{*})$ is at least approximately
scalable, we know that $\mu(A)\geq\nu^{*}(F_{R}(A))$. On the other hand, as
$A$ is a maximal $\theta$-set, we know that
$\mu(A)=\theta_{m}\nu(F_{R}(A))=\overline{\theta}\nu(F_{R}(A))$. But by
definition of $\overline{\theta}$, we know that
$\overline{\theta}\nu\geq\nu^{*}$, so that $\nu^{*}(F_{R}(A))\leq\mu(A)$. We
conclude that $\nu^{*}(F_{R}(A))=\mu(A)$, that
$\nu^{*}\llcorner_{A}=\overline{\theta}\nu\llcorner_{A}=\theta_{m}\nu\llcorner_{A}$,
and hence that

$P^{*}\Big{(}(\mathcal{D}\backslash A)\times
F_{R}(A)\Big{)}=P^{*}(\mathcal{D}\times F_{R}(A))-P^{*}(A\times
F_{R}(A))=\nu^{*}(F_{R}(A))-\mu(A)=0,$

so that (29) holds.

In addition, by Proposition 29 the measure $P^{*}\llcorner_{A\times F_{R}(A)}$
is the solution of the problem $\operatorname{Sch}(R\llcorner_{A\times
F_{R}(A)};\mu\llcorner_{A},\nu^{*}\llcorner_{F_{R}(A)})$. So in order to prove
(35), it suffices to prove that this problem is scalable. For this purpose, we
will use Lemma 24. We call $R_{r}:=R\llcorner_{A\times F_{R}(A)}$. Let $B$ be
a nonempty strict subset of $A$. As $A$ is a minimal element in the set of
maximal $\theta$-sets for $(R;\mu,\nu)$, we know that
$\mu(B)<\theta_{m}\nu(F_{R}(B))$. Then, as $F_{R}(B)\subset F_{R}(A)$, we have
$F_{R}(B)=F_{R_{r}}(B)$, so $\mu(B)<\theta_{m}\nu(F_{R_{r}}(B))$. Finally, as
$\nu^{*}\llcorner_{A}=\theta_{m}\nu\llcorner_{A}$, we have
$\mu(B)<\nu^{*}(F_{R_{r}}(B))$. So Lemma 24 applies, and
$\operatorname{Sch}(R\llcorner_{A\times
F_{R}(A)};\mu\llcorner_{A},\nu^{*}\llcorner_{F_{R}(A)})$ is scalable, which
concludes the proof. ∎

We close this section with a remark concerning the stability with respect to
union of SISP sets.

######  Remark 32.

It is easy to check that SISP sets associated with a triple $(R;\mu,\nu)$ are
stable by union. Therefore, there exists an upper bound in the set of all SISP
sets for $(R;\mu,\nu)$, that we call the _largest_ SISP set. If we want the
procedure described in Proposition 28 to be as fast as possible, it is logical
to look for SISP sets that are as large as possible, in order to minimize the
rank $N$ at which the procedure reaches its stationary point. This is what we
are going to do in the next section.

##  6 Numerical applications

A simple consequence of the theoretical procedure described in the previous
section is that in a lot of cases, if the problem
$\operatorname{Sch}(R;\mu,\nu)$ is non-scalable, the matrices $P^{*}$ and
$Q^{*}$ from Theorem 11, and $R^{*}$ from Theorem 17 have more zero entries
than $R$. For instance, if the problem is balanced (_i.e._
$\mathsf{M}(\mu)=\mathsf{M}(\nu)$), and if the bipartite graph of $R$ is
connected (that is, $\mu(A)=\nu(F_{R}(A))$ only holds for $A=\emptyset$ or
$A=\mathcal{D}$, which is a reasonable assumption in a lot of contexts), we
can check that the matrix $R^{1}$ from Proposition 28 cannot coincide with
$R$.

Therefore, typically, the Sinkhorn algorithm in the non-scalable case does not
converge linearly. We are going to detail an approximate algorithm allowing to
find the common support of $P^{*}$, $Q^{*}$ and $R^{*}$, and therefore to
recover a linear rate of convergence for the Sinkhorn algorithm by Proposition
25.

###  6.1 Stopping criterion

Before any numerical application, we need to define a stopping criterion for
the Sinkhorn algorithm when the Schrödinger problem is non-scalable. When the
problem is scalable, the classical criterion that is used is the duality gap
estimated at each step $n\in\mathbb{N}^{*}$ of the Sinkhorn algorithm:

$\displaystyle{SC}^{n}=H(P^{n}|R)-\langle\log(a^{n}),\mu\rangle-\langle\log(b^{n-1}),\nu\rangle,$
(40)

where $P^{n}$, $a^{n}$ and $b^{n-1}$ are defined by the relations (2). Indeed,
it is known that this quantity is always positive when the relation
$P^{n}_{ij}=a^{n}_{i}b^{n-1}_{j}R_{ij}$ holds for all $i,j$, and the Fenchel-
Rockafellar duality ensures that ${SC}_{n}\to 0$ as $n\to\infty$ when the
problem is scalable, i.e. when $(a^{n})$ and $(b^{n})$ converge.

In the approximately scalable case, numerical instabilities may appear when
$n\to\infty$ because $(a^{n})$ and $(b^{n})$ do not converge, but this
criterion may remain useful if the error that is tolerated is not too small.
However, in the non-scalable case, this criterion does not hold as the problem
$\operatorname{Sch}(R;\mu,\nu)$ has no solution. The results presented in the
previous sections allow nevertheless to define an approximate criterion.
Indeed, it has been shown in [6] that for a given $\lambda>0$, the problem
defined with the notations of Subsection 2.1 by

$\textrm{Schu}_{\lambda}(R;\mu,\nu)=\min\Big{\\{}H(P|R)+\lambda
H(\nu^{P}|\nu)\,\Big{|}\,P\mbox{ s.t.\ }\mu^{P}=\mu\Big{\\}},$ (41)

can be solved numerically with a generalization of the Sinkhorn algorithm.
More precisely, the duality gap defined for all $n\in\mathbb{N}^{*}$ by

${SCu}_{\lambda}^{n}=H(P^{n}|R)+\lambda\left(H(\nu^{P^{n}}|\nu)-\Big{\langle}1-(1/b^{n-1})^{\frac{1}{\lambda}},\nu\Big{\rangle}\right)-\langle\log(a^{n}),\mu\rangle,$
(42)

converges to $0$ whenever $P^{n}_{ij}:=a^{n}_{i}b^{n-1}_{j}R_{ij}$ for all
$i,j$, and with $a^{n},b^{n}$ defined for all $n\in\mathbb{N}^{*}$ by the
relations:

$\left\\{\begin{gathered}\forall j,\quad b^{0}_{j}:=1,\\\\[10.0pt]
\begin{aligned} &\forall n\geq 0,\quad\forall
i,&a^{n+1}_{i}&:=\frac{\mu_{i}}{\displaystyle\sum_{j}b^{n}_{j}R_{ij}},\\\
&\forall n\geq 0,\quad\forall
j,&b^{n+1}_{j}&:=\Bigg{(}\frac{\nu_{j}}{\displaystyle\sum_{i}a^{n+1}_{i}R_{ij}}\Bigg{)}^{\frac{\lambda}{1+\lambda}}.\end{aligned}\end{gathered}\right.$
(43)

On the other hand, a slight modification of our $\Gamma$-convergence result of
Proposition 15 asserts that $P^{*}$, from Theorem 11, is the limit of the
solution of (41) as $\lambda\to+\infty$. So if we now define
$SCu_{\lambda}^{n}$ by the formula (42) where $P^{n}$, $a^{n}$ and $b^{n-1}$
are computed with the standard Sinkhorn algorithm (2), instead of the modified
one (43), we conclude that for all $\varepsilon>0$, there exists a threshold
$\lambda^{\varepsilon}$ such that for all $\lambda\geq\lambda^{\varepsilon}$,

$\limsup_{n\to+\infty}SCu_{\lambda}^{n}\leq\varepsilon.$

Therefore, the stopping criterion (42) can still be used for the sequence
$(P^{n})_{n\in\mathbb{N}}$ generated by the classical Sinkhorn algorithm (2),
as long as $\lambda^{\varepsilon}$ is chosen to be sufficiently large w.r.t.
$\varepsilon$. In practice, we observe that taking
$\lambda^{\varepsilon}=\frac{1}{\varepsilon}$ works well. This is what we are
going to do in the following section, considering a level of error
$\varepsilon:=10^{-3}$.

###  6.2 An approximate numerical method for constructing the support of
$R^{*}$

An interesting application of the theoretical procedures described in
Subsection 5.2 is the construction of an approximate algorithm allowing the
identification of the support $\mathcal{S}$ of $R^{*}$ w.r.t. $R$, when the
problem $\operatorname{Sch}(R;\mu,\nu)$ is approximately scalable or non-
scalable. Our motivation is double. First, the fact that knowing the support
allows to recover a linear rate for the Sinkhorn algorithm, thanks to
Proposition 25, suggests that if this algorithm is simple enough, the full
procedure to obtain $R^{*}$ (preprocessing to know the $\mathcal{S}$ and then
the Sinkhorn algorithm applied to
$\operatorname{Sch}(\mathds{1}_{\mathcal{S}}R;\mu,\nu)$ is excepted to be
faster than the Sinkhorn algorithm applied directly to
$\operatorname{Sch}(R;\mu,\nu)$. Second, even if we have shown that the
Sinkhorn algorithm converges in the non-scalable case, the Schrödinger
potentials appearing in the Sinkhorn procedure are likely to be too high to be
computed numerically before the algorithm has converged: thus, finding the
support before to run the Sinkhorn algorithm can be an advantage, as long as
this preprocessing prevents this potential’s explosion and even if it does not
accelerate the procedure.

We recall that we identified, at each iteration $n$ of the theoretical
procedure described in Proposition 28, a SISP set, denoted $M_{n}$, for a
reduced problem
$\operatorname{Sch}(R\llcorner_{\mathcal{D}^{n}\times\mathcal{F}^{n}};\mu\llcorner_{\mathcal{D}^{n}},\nu\llcorner_{\mathcal{F}^{n}})$.
Regarding the proof of Lemma 27, a natural choice for $M_{n}$ is an union of
smallest maximal $\theta$-sets for
$\operatorname{Sch}(R\llcorner_{\mathcal{D}^{n}\times\mathcal{F}^{n}};\mu\llcorner_{\mathcal{D}^{n}},\nu\llcorner_{\mathcal{F}^{n}})$,
that we introduced in Definition 30. The approximate procedure that we design
for finding $M_{n}$ at each iteration consists in two steps:

  1. 1.

Find the largest set $M_{n}^{\prime}$ (in the sens of inclusion) such that the
supports of $R$ and $R^{*}$ coincide in $M_{n}^{\prime}\times\mathcal{F}$;

  2. 2.

Find $M_{n}\subset M_{n}^{\prime}$, an union of smallest maximal $\theta$-set
for $(R\llcorner_{M_{n}^{\prime}\times
F_{R}(M_{n}^{\prime})};\mu\llcorner_{M_{n}^{\prime}},\nu\llcorner_{F_{R}(M_{n}^{\prime})})$.

Indeed, finding $M_{n}^{\prime}$ in the first step can be achieved with a
slightly modified Sinkhorn procedure applied to
$\operatorname{Sch}(R\llcorner_{\mathcal{D}^{n}\times\mathcal{F}^{n}},\mu\llcorner_{\mathcal{D}^{n}},\nu\llcorner_{\mathcal{F}^{n}})$.
More precisely, we initialize $U=\mathcal{D}^{n}$ and $V=\mathcal{F}^{n}$, and
at each step of this Sinkhorn procedure, for all $x_{i}\in\mathcal{D}^{n}$, if
there exists $y_{j}\in V$ such that the coupling obtained at this step is
smaller than a given threshold, we set $U=U\backslash x_{i}$. Since Lemma 27
ensures that there exists at least one SISP set, $M_{n}^{\prime}$ is not empty
and thus we do not set all the lines to $0$ (given that the threshold is
sufficiently small). For all $y_{j}\in V$, if $y_{j}\notin
F_{R\llcorner_{\mathcal{D}^{n}\times\mathcal{F}^{n}}}(U)$, we set
$V=V\backslash y_{j}$.

After several steps of this modified Sinkhorn procedure, the elements
remaining in $U$ and $V$ coincide respectively with $M_{n}^{\prime}$ and
$F_{R\llcorner_{\mathcal{D}^{n}\times\mathcal{F}^{n}}}(M_{n}^{\prime})$. Since
the Sinkhorn algorithm restricted to $U\times V$ converges linearly (thanks to
Proposition 25), we then rapidly observe the convergence of the procedure.

For detailing how to obtain $M_{n}$ at the second step, we need to introduce
the notion of connected components:

######  Definition 33.

Let $U\subset\mathcal{D}$ and $V\subset\mathcal{F}$. We say that $A\times B$
is a connected component of the graph $(U\cup V,\triangle)$ (using the
notation of 28) if:

  * •

$R\llcorner_{U\times V}(A\times B^{c})=R\llcorner_{U\times V}(A^{c}\times
B)=0$;

  * •

$(A\cup B,\triangle)$ is connected.

Finding connected components for such undirected graph is a classical task in
Graph theory, for which there exists ready-to-use algorithms [12]. The
following proposition justifies that we can define at the second step of the
procedure

$\displaystyle
M_{n}=\bigcup\left\\{U_{i}\,\bigg{|}\,\frac{\mu(U_{i})}{\nu(V_{i})}=\max\limits_{j=1,\cdots,C}\frac{\mu(U_{j})}{\nu(V_{j})}\right\\},$
(44)

where $U_{1}\times V_{1},\cdots,U_{C}\times V_{C}$ are the connected
components of the graph $(M_{n}^{\prime}\cup
F_{R\llcorner_{\mathcal{D}^{n}\times\mathcal{F}^{n}}}(M_{n}^{\prime}),\triangle)$.

######  Proposition 34.

Let $R\in\mathcal{M}_{+}(\mathcal{D}\times\mathcal{F})$,
$\mu\in\mathcal{M}_{+}(\mathcal{D})$ and $\nu\in\mathcal{M}_{+}(\mathcal{F})$
satisfying Assumption 9, and let us denote $U$ the largest subset of
$\mathcal{D}$ (in the sens of inclusion) such that the supports of $R$ and
$R^{*}$ coincide in $U\times\mathcal{F}$, and $V=F_{R}(U)$. Let us denote
$(U_{1}\times V_{1},\cdots,U_{C}\times V_{C})$ the $C$ connected components of
the graph $(U\cup V,\triangle)$. Then, for any $U_{i}\times V_{i}$ such that

$\displaystyle\frac{\mu(U_{i})}{\nu(V_{i})}=\max\limits_{j=1,\cdots,C}\frac{\mu(U_{j})}{\nu(V_{j})},$
(45)

$U_{i}$ is a smallest maximal $\theta$-sets for $(R;\mu,\nu)$.

As we have already shown that for all $n\leq 0$, Assumption 9 holds for
$(R\llcorner_{\mathcal{D}^{n}\times\mathcal{F}^{n}},\mu\llcorner_{\mathcal{D}^{n}},\nu\llcorner_{\mathcal{F}^{n}})$,
we can apply Proposition 34 to this triple.

###### Proof.

The main point consists in showing that for any smallest maximal $\theta$-set
for $(R;\mu,\nu)$, that we denote $A$, $A\times F_{R}(A)$ is a connected
component of $(U\cup V,\triangle)$. Indeed, assuming this claim, as there
always exists such smallest maximal $\theta$-set, any connected component
$U_{i}\times V_{i}$ maximizing (45), $U_{i}$ is a maximal $\theta$-set:
$U_{i}$ is then necessarily a smallest one because it cannot contain any other
connected component, and then any smallest maximal $\theta$-set neither. We
now prove the claim.

As $A$ is a SISP set (thanks to the proof of Proposition 27), it is in $U$ and
we have necessarily $R^{*}(A^{c}\times F_{R}(A))=R^{*}(A\times
F_{R}(A)^{c})=0$. The supports of $R$ and $R^{*}$ being the same in $U\times
V$ by construction of $U$ and $V$, we have thus $R\llcorner_{U\times
V}(A^{c}\times F_{R}(A))=R\llcorner_{U\times V}(A\times F_{R}(A)^{c})=0$.
Moreover, $A\times F_{R}(A)$ cannot contain any connected component of $(U\cup
V,\triangle)$: if it was the case, this connected component would characterize
a maximal $\theta$-set for $(R;\mu,\nu)$, which would contradict the
minimality of $A$. Thus, $A\times F_{R}(A)$ is necessarily a connected
component of the graph $(U\cup V,\triangle)$. ∎

We provide in Algorithm 1 the pseudo-code of this iterative method.

Algorithm 1 Find the support Supp of $R^{*}$: return Supp

0: $\bullet$ A set of minimal factors: $\\{m_{i},\,i=1,\cdots,N\\}$,$\bullet$
A stopping criterion: $stop(a,b,R,\mu,\nu)$.We set $A=\mathcal{D}$,
$B=\mathcal{F}$, $\textrm{Supp}=\textrm{Support}(R)$.

while $A\neq\emptyset$ do

$\bar{R}=R\llcorner_{A\times B}$, $\bar{\mu}=\mu\llcorner_{A}$,
$\bar{\nu}=\nu\llcorner_{B}$

$b=\mathds{1}_{B}$, $a=\mathds{1}_{A}$

$U=A$, $V=B$

while $stop(a,b,\bar{R},\bar{\mu},\bar{\nu})\neq 1$ do

for $x_{i}\in U$ do

if $\sum\limits_{y_{j}\in V}\bar{R}_{ij}b_{j}=0$ then

$U=U\backslash\\{x_{i}\\}$

else

$a_{i}=\displaystyle{\frac{\bar{\mu}_{i}}{\sum\limits_{y_{j}\in
V}\bar{R}_{ij}b_{j}}}$

if $\min\limits_{y_{j}\in V}a_{i}\times b_{j}<m_{i}$ then

$U=U\backslash\\{x_{i}\\}$

end if

end if

end for

$\bar{R}=\bar{R}\llcorner_{U\times V}$, $\bar{\mu}=\bar{\mu}\llcorner_{U}$,
$a=a\llcorner_{U}$

for $y_{j}\in V$ do

if $\sum\limits_{x_{i}\in A}\bar{R}_{ij}a_{i}=0$ then

$V=V\backslash\\{y_{j}\\}$

else

$b_{j}=\displaystyle{\frac{\bar{\nu}_{j}}{\sum\limits_{x_{i}\in
A}\bar{R}_{ij}a_{i}}}$

end if

end for

$\bar{R}=\bar{R}\llcorner_{U\times V}$, $\bar{\nu}=\bar{\nu}\llcorner_{V}$,
$b=b\llcorner_{V}$

end while

$(U_{1}\times V_{1},\cdots,U_{C}\times V_{C})=\textrm{connected components of
the graph}\,(U\cup V,\triangle)$

$U\times V=\bigcup\big{\\{}U_{i}\times
V_{i},\,|\,\frac{\mu(U_{i})}{\nu(V_{i})}=\max\limits_{j=1,\cdots,C}\frac{\mu(U_{j})}{\nu(V_{j})}\big{\\}}$

$A=A\backslash U$, $B=B\backslash V$

$\textrm{Supp}=\textrm{Supp}\backslash\big{(}(A\backslash U)\times V\big{)}$

end while

Let us make a few comment on this algorithm. The support of $R^{*}$ only
depends on the support of $R$ and not of its values, so when identifying
$\mathcal{S}$, we can equivalently consider the problem
$\operatorname{Sch}(R^{\prime};\mu,\nu)$, where $R^{\prime}=\mathds{1}_{R\neq
0}$. This explains why we consider the minimum $\min\limits_{y_{j}\in
V}a_{i}\times b_{j}$ rather than $\min\limits_{y_{j}\in V}a_{i}\times
b_{j}\times R_{ij}$ in Algorithm 1.

The stopping criterion corresponds to the criterion (42) detailed in the
previous section, which has to be smaller than a certain threshold
$\varepsilon$ to be satisfied.

Choosing in an appropriate way the set of minimal factors
$\\{m_{i},\,i=1,\cdots,N\\}$, is crucial: it determines the level of
approximation that is considered as acceptable, i.e. the minimal value at
which we can consider that the algorithm should create a new zero entries. In
practice, we observe that

$\displaystyle m_{i}:=\frac{1}{n}\frac{\mu_{i}}{\mu_{i}^{R}},$ (46)

seems to be a good tradeoff between efficiency and security in most of the
cases that we explored.

######  Remark 35.

  * •

This method can be seen as an improvement over the naive approximate method
which consists, at each iteration of the Sinkhorn algorithm applied to
$\operatorname{Sch}(R;\mu,\nu)$, to set all the entries of $R$ that are
smaller to a certain threshold to zero. With our method, we do not have to
identify all the zero entries one by one but line by line, which can avoid
numerical errors which may appear for some entries converging slowly to zero.
However, this is done at the cost of identifying the connected components of
some subgraphs of $(\mathcal{D}\cup\mathcal{F},\triangle)$ at every iteration
of the procedure, which slows down the algorithm in cases where these
subgraphs are large. Note that in typical cases where the matrix $R$ is
structured, we do not expect to find more than one connected component at each
iteration as in Figure 4.

  * •

We emphasize the fact that this algorithm is only approximate, and this for
two reasons. The first one occurs when the set of thresholds
$\\{m_{i},\,i=1,\dots,N\\}$ are too large. Then, we can set to zero lines
which should not be, just because some of the entries of $R^{*}$ should be
small on this line. This must be avoided as then the algorithm cannot converge
towards $R^{*}$.

The other case where our algorithm does not identify $\mathcal{S}$ exactly is
either when the threshold $\varepsilon$ of the stopping criterion is large, or
when the thresholds $\\{m_{i},\,i=1,\dots,N\\}$ are small. Then, the Sinkhorn
algorithm can satisfy the stopping criterion before all the zeroes have been
identified. This is not a big problem, since it means that the algorithm
converges well without having to identify the additional zeroes of $R^{*}$.

With these observations, we conclude that the thresholds
$\\{m_{i},\,i=1,\dots,N\\}$ need to be taken rather small w.r.t. the threshold
$\varepsilon$ of the stopping criterion, even though of course, if they are
taken too small, efficiency is lost since then the algorithm just behaves as
Sinkhorn without any improvement.

We illustrate in Figure 5 the efficiency of this procedure. As emphasized at
the beginning of this section, the number of iterations of a Sinkhorn-like
algorithm is a better indicator than the computation time of the full
procedure because our main motivation to find the support of $R^{*}$ before
running the Sinkhorn algorithm is that we want the Schrödinger potentials not
to explode before reaching convergence, which could be the case for non-
scalable problems. We thus represent the number of iterations needed for the
Sinkhorn algorithm to converge as a function of the number of additional zero
entries in $R^{*}$ w.r.t. to $R$, and compare when we apply or not the
preprocessing described in Algorithm 1. For varying the number of zero
entries, we take $R$ upper-diagonal, build $\mu$ and $\nu$ similarly to what
we described in Figure 4, and then vary the number of blocks from $1$