lemmatheorem lemma corollarytheorem corollary propositiontheorem proposition
definitiontheorem definition definitionPropositiontheorem
definitionProposition remarktheorem remark

#  Quantitative contraction rates for Sinkhorn’s algorithm: beyond bounded
costs and compact marginals

Giovanni Conforti  École Polytechnique Département de Mathématiques Appliqués,
Palaiseau, France<EMAIL_ADDRESS>,  Alain Oliviero
Durmus  École Polytechnique Département de Mathématiques Appliqués, Palaiseau,
France<EMAIL_ADDRESS>and  Giacomo Greco  Eindhoven
University of Technology Department of Mathematics and Computer Science,
5600MB Eindhoven, Netherlands<EMAIL_ADDRESS>
###### Abstract.

We show non-asymptotic exponential convergence of Sinkhorn iterates to the
Schrödinger potentials, solutions of the quadratic Entropic Optimal Transport
problem on $\mathbb{R}^{d}$. Our results hold under mild assumptions on the
marginal inputs: in particular, we only assume that they admit an
asymptotically positive log-concavity profile, covering as special cases log-
concave distributions and bounded smooth perturbations of quadratic
potentials. More precisely, we provide exponential $\mathrm{L}^{1}$ and
pointwise convergence of the iterates (resp. their gradient and Hessian) to
Schrödinger potentials (resp. their gradient and Hessian) for large enough
values of the regularization parameter. As a corollary, we establish
exponential convergence of Sinkhorn plans with respect to a symmetric relative
entropy and in $\mathbf{W}_{1}$. Up to the authors’ knowledge, these are the
first results which establish exponential convergence of Sinkhorn’s algorithm
in a general setting without assuming bounded cost functions or compactly
supported marginals.

GC acknowledges funding from the grant SPOT (ANR-20-CE40- 0014)

AD would like to thank the Isaac Newton Institute for Mathematical Sciences
for support and hospitality during the programme _The mathematical and
statistical foundation of future data-driven engineering_ when work on this
paper was undertaken. This work was supported by: EPSRC grant number
EP/R014604/1.

GG thanks École Polytechnique for its hospitality, where this research has
been carried out, and NDNS+ for funding his visit there (NDNS/2023.004 PhD
Travel Grant). GG is supported from the NWO Research Project 613.009.111
"Analysis meets Stochastics: Scaling limits in complex systems"

###### Contents

  1. 1 Introduction
  2. 2 Sinkhorn’s algorithm and its geometric convergence
    1. 2.1 A first (informal) statement of our results and assumptions
    2. 2.2 Main results and proof strategy
    3. 2.3 Literature review
  3. 3 Integrated semiconvexity propagation along Sinkhorn
  4. 4 Wasserstein distance w.r.t a measure with log-concave profile
  5. 5 Proof of the main results
    1. 5.1 Exponential convergence of the gradients and $\mathbf{W}_{1}$ along Sinkhorn’s algorithm
    2. 5.2 Exponential pointwise and entropic convergence of Sinkhorn’s algorithm
    3. 5.3 Exponential convergence of the Hessians along Sinkhorn’s algorithm
  6. A Explicit rates for strongly log-concave marginals
  7. B Explicit construction and contractive properties of $\mathbf{W}_{f^{\nu}_{x}}$
  8. C Technical results and moment bounds
    1. C.1 Proof of

##  1\. Introduction

Given two distributions $\mu,\nu\in\mathcal{P}(\mathbb{R}^{d})$, we consider
in this paper the problem of finding a solution to the entropic optimal
transport problem: for $T\geq 0$,

(EOT) $\text{minimize
}\int\frac{|x-y|^{2}}{2}\,\mathrm{d}\pi+T\,\mathscr{H}(\pi|\mu\otimes\nu)\text{
under the constraint }\pi\in\Pi(\mu,\nu)\;,$

where $\mathscr{H}$ denotes relative entropy (or Kullback-Leibler divergence)
and $\Pi(\mu,\nu)$ is the set of couplings of $\mu$ and $\nu$. In EOT, the
case $T=0$ corresponds to the Monge-Kantorovich or Optimal Transport (OT)
problem and $T$ has to be understood as a regularization parameter.
Remarkably, the variational problem EOT, whose study in statistical machine
learning is motivated by its computational advantages over the Monge-
Kantorovich problem, is equivalent to the Schrödinger problem (SP), that is an
old statistical mechanics problem whose origin can be traced back to a
question posed by Schrödinger [58] about the most likely behaviour of a cloud
of diffusive particles conditionally on the observation of its configuration
at an initial and final time. We refer to [42, 47] for surveys on SP and EOT.
Both SP and EOT stand nowadays at the crossroads of different research areas
in both pure and applied mathematics including stochastic optimal control [13,
14, 45], geometry of optimal transport and functional inequalities [3, 19, 17,
29, 31, 32, 50], abstract metric spaces [34, 35, 46], numerics for PDEs and
gradient flows [5, 6, 4, 51].

Recently, optimal transport and its entropic version have drawn interest in
statistics and machine learning [52] and find a wealth of applications such as
inverse problems [1], classification [40], generative modeling [21, 2],
computer vision [25, 61] and signal processing [39]. However, in most cases,
the numerical resolution of OT is challenging while EOT turns to be simpler to
tackle computationally for moderately small parameter $T$ and its solutions
are used as proxies for solutions of OT [20, 52]. This is in particular
justified by [48, 53] which show convergence of solutions of EOT to OT as
$T\to 0$. As just mentioned, EOT is simpler to solve compared to OT mainly due
to Sinkhorn’s algorithm also called Iterative Proportional Fitting Procedure,
which can be used at the basis of numerical schemes aiming to get approximate
solutions for EOT .

To introduce Sinkhorn’s algorithm, that is the object of study of this work,
we begin by recalling that under mild assumptions on the marginals $\mu,\nu$
(see for instance [16, Proposition 2.2]), (EOT) admits a unique minimizer
$\pi^{\star}\in\Pi(\mu,\nu)$, referred to as the Schrödinger plan, and there
exist two functions $\varphi^{\star}\in\mathrm{L}^{1}(\mu)$ and
$\psi^{\star}\in\mathrm{L}^{1}(\nu)$, called Schrödinger potentials, such that

(1) $\pi^{\star}(\mathrm{d}x\mathrm{d}y)=(2\pi
T)^{-d/2}\,\exp\biggl{(}-\frac{|x-y|^{2}}{2T}-\varphi^{\star}(x)-\psi^{\star}(y)\biggr{)}\mathrm{d}x\,\mathrm{d}y\;.$

Moreover, the couple $(\varphi^{\star},\psi^{\star})$ is unique up to constant
translations $a\mapsto(\varphi^{\star}+a,\psi^{\star}-a)$. If we suppose that
the marginals admit densities of the form

(2)
$\mu(\mathrm{d}x)=\exp(-U_{\mu}(x))\mathrm{d}x\,,\qquad\nu(\mathrm{d}x)=\exp(-U_{\nu}(x))\mathrm{d}x\;,$

then, imposing that a noisy transport plan of the form (1) belongs to
$\Pi(\mu,\nu)$ one finds that $\varphi^{\star},\psi^{\star}$ solve the
Schrödinger system

(3) $\begin{cases}\varphi^{\star}=U_{\mu}+\log P_{T}\exp(-\psi^{\star})\\\
\psi^{\star}=U_{\nu}+\log P_{T}\exp(-\varphi^{\star})\;,\end{cases}$

where $(P_{t})_{t\geq 0}$ is the Markov semigroup generated by the standard
$d$-dimensional Brownian motion $(B_{t})_{t\geq 0}$:
$P_{t}f(x)=\mathbb{E}[f(x+B_{t})]$ for any non-negative measurable function
$f:\mathbb{R}^{d}\to\mathbb{R}$.

Then, starting from a given initialization
$\psi^{0},\varphi^{0}\colon\mathbb{R}^{d}\to\mathbb{R}$, Sinkhorn’s algorithm
solves (3) as a fixed point problem by generating two sequences of potentials
$\\{\varphi^{n},\psi^{n}\\}_{n\in\mathbb{N}}$, called Sinkhorn potentials
according to the following recursive scheme:

(4) $\begin{cases}\varphi^{n+1}=U_{\mu}+\log P_{T}\exp(-\psi^{n})\\\
\psi^{n+1}=U_{\nu}+\log P_{T}\exp(-\varphi^{n+1})\;.\end{cases}$

The main goal of the present paper is to provide a new convergence analysis
for Sinkhorn iterates to the Schrödinger potentials $\varphi^{\star}$ and
$\psi^{\star}$, and hence the optimal plan $\pi^{\star}$. When working on
discrete spaces, convergence of the algorithm is well-known and we refer to
the book [52] for an extensive overview. However, let us mention that in the
discrete setting, convergence of Sinkhorn’s algorithm dates back at least to
[59] and [60]. More recently, [30, 9] show that Sinkhorn iterates are
equivalent to a sequence of iterates associated to an appropriate contraction
in the Hilbert projective metrics, therefore proving the convergence of the
algorithm boils down to studying a fixed-point problem, whose (exponential)
convergence can be deduced from Birckoff’s Theorem. In the case of continuous
state spaces, the use of Birckoff’s Theorem for the Hilbert metrics has also
been employed by [12], in order to establish the exponential convergence of
Sinkhorn’s algorithm under the condition that the state space is compact or
that the cost function is bounded.

When considering non-compact continuous state spaces and unbounded costs,
results are scarcer; see Section 2.3 for a detailed literature overview. In
addition, to the best of our knowledge and understanding, exponential
convergence of Sinkhorn iterates for unbounded costs and marginals $\mu,\nu$
with unbounded support is still an open problem. In this paper, we address
this important question by showing convergence of Sinkhorn iterates and its
derivatives in the landmark example of the quadratic cost for a large class of
unbounded marginal distributions and large enough values of the regularization
parameter $T$.

More precisely, our first main contribution is to establish exponential
$L^{1}$-convergence bounds for the gradients of Sinkhorn iterates
$\\{\nabla\varphi^{n},\nabla\psi^{n}\\}_{n\in\mathbb{N}^{*}}$ to the gradients
of the Schrödinger potentials $\nabla\varphi^{\star}$ and
$\nabla\psi^{\star}$. Our second main contribution is to prove quantitative
exponential convergence of Sinkhorn’s plans towards the Schrödinger plan in
Wasserstein distance and, under more stringent assumptions, also in relative
entropy. In contrast to our result, previous works only established linear
convergence rates in the unbounded cost setting, see [33]. Finally, we also
establish exponential convergence of the Hessians of the Sinkhorn iterates
$\\{\nabla^{2}\varphi^{n},\nabla^{2}\psi^{n}\\}_{n\in\mathbb{N}^{*}}$ to
$\nabla^{2}\varphi^{\star}$ and $\nabla^{2}\psi^{\star}$. These are the first
convergence result we are aware of concerning second order derivatives, though
in the bounded cost setting one could bootstrap convergence of first and
second derivatives from convergence of iterates using the results of [22].

Organisation. The rest of the paper is organized as follows. We start Section
2 introducing Sinkhorn’s algorithm in more detail, then we state informally
our main results in Section 2.1. We give in Section 2.2 the complete
statements of our results in their general form and provide a sketch of our
proof strategy for the readers’ convenience. Section 2.3 is devoted to a
review of the existing literature on the convergence of Sinkhorn’s algorithm.
In Section 3 and Section 4, we develop important tools that we use to obtain
our main results. In particular, in Section 3 we prove integrated
semiconvexity estimates for Sinkhorn iterates and Schrödinger potentials. In
Section 4 we show how (synchronous and reflection) coupling techniques can be
employed in order to bound the Wasserstein distance between two probability
measures with asymptotically positive log-concavity profile. This result is
stated in a general setting as it is of independent interest. In Section 5 we
carry out the proof of the main results assembling together the results of the
previous sections. In Appendix A we derive explicit expression for the rates
appearing in Section 2.2 in the strongly log-concave setting. Lastly, we
gather at Appendices B and C some technical results and moment bounds employed
in Section 5.

Notation. The set $\mathbb{R}^{d}$ is endowed with the standard Euclidean
metric and we denote by $\left|\cdot\right|$ and
$\left\langle\cdot,\cdot\right\rangle$ the corresponding norm and scalar
product. For $a\in\mathbb{R}$, we denote by $a_{+}=0\vee a$ and
$a_{-}=0\vee(-a)$. We denote by $\mathcal{B}(\mathbb{R}^{d})$ the Borel
$\sigma$-field of $\mathbb{R}^{d}$, by $\mathcal{P}(\mathbb{R}^{d})$ the space
of probability measures defined on
$(\mathbb{R}^{d},\mathcal{B}(\mathbb{R}^{d}))$ and by
$\mathcal{P}_{2}(\mathbb{R}^{d})$ the subset of probability measures with
finite second moment. For any probability measure
$\mu\in\mathcal{P}(\mathbb{R}^{d})$ we denote with
$M_{k}(\mu)=\int\left|x\right|^{k}\mathrm{d}\mu(x)$ its $k^{th}$ moment and
define its covariance matrix as $\mathrm{Cov}(\mu)\coloneqq\int xx^{\sf
T}\mathrm{d}\mu-(\int x\mathrm{d}\mu)(\int x\mathrm{d}\mu)^{\sf T}$. For two
distributions $\mu_{1},\mu_{2}$ on
$(\mathbb{R}^{d},\mathcal{B}(\mathbb{R}^{d}))$, $\Pi(\mu_{1},\mu_{2})$ denotes
the set of couplings between $\mu_{1}$ and $\mu_{2}$, i.e.,
$\xi\in\Pi(\mu_{1},\mu_{2})$ if and only if $\xi$ is a probability measure on
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ and
$\xi(\mathsf{A}\times\mathbb{R}^{d})=\mu_{1}(\mathsf{A})$ and
$\xi(\mathbb{R}^{d}\times\mathsf{A})=\mu_{2}(\mathsf{A})$. On
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ we consider the projection operators
$\mathrm{proj}_{x}(a,b)=a$ and $\mathrm{proj}_{y}(a,b)=b$. For any couple of
functions $f,g\colon\mathbb{R}^{d}\to\mathbb{R}$ let us denote by $f\oplus
g\colon\mathbb{R}^{2d}\to\mathbb{R}$ the function defined as $f\oplus
g(x,y)\coloneqq f(x)+g(y)$. For $f:\mathbb{R}^{d}\to\mathbb{R}^{d}$ measurable
and $\mu\in\mathcal{P}(\mathbb{R}^{d})$, we denote by $f_{\sharp}\mu$ the
pushforward measure of $\mu$ by $f$:
$f_{\sharp}\mu(\mathsf{A})=\mu(f^{-1}(\mathsf{A}))$ for any
$\mathsf{A}\in\mathcal{B}(\mathbb{R}^{d})$. For any function
$f\colon\mathbb{R}^{d}\to\mathbb{R}$, denote by
$\left\|f\right\|_{\infty}\coloneqq\sup_{x\in\mathbb{R}^{d}}|f(x)|$ its
supremum norm. Denote by $\rm\mathrm{L}^{1}(\mu)$ the set of functions
integrable with respect to $\mu$. For two (non-negative) measures
$\lambda_{1},\lambda_{2}$, we define the relative entropy (or Kullback Leibler
divergence) between $\lambda_{1}$ and $\lambda_{2}$ as
$\mathscr{H}(\lambda_{1}\mid\lambda_{2})=\int_{\mathbb{R}^{d}}\log(\mathrm{d}\lambda_{1}/\mathrm{d}\lambda_{2})\mathrm{d}\lambda_{1}$
if $\lambda_{1}\ll\lambda_{2}$ and
$\int_{\mathbb{R}^{d}}[\log(\mathrm{d}\lambda_{1}/\mathrm{d}\lambda_{2})]_{-}\mathrm{d}\lambda_{1}<+\infty$,
and $\mathscr{H}(\lambda_{1}\mid\lambda_{2})=+\infty$ otherwise. We denote by
$\mathrm{Leb}$ the Lebesgue measure and define
$\mathrm{Ent}(\lambda_{1})=\int\log(\mathrm{d}\lambda_{1}/\mathrm{d}\mathrm{Leb})\mathrm{d}\mathrm{Leb}$
if $\lambda_{1}\ll\mathrm{Leb}$ and $+\infty$ otherwise.

##  2\. Sinkhorn’s algorithm and its geometric convergence

To gain a deeper understanding and build intuition on how Sinkhorn’s
algorithms works, it is worth recalling some basic results on the Schrödinger
problem: this interlude allows us to naturally introduce objects and concepts
that will be used extensively throughout our convergence analysis. The
Schrödinger problem (SP) associated with EOT , is an equivalent formulation of
EOT inspired by Schrödinger’s seminal work [58] about the behaviour of
Brownian particles conditionally to some observations; see [42, Sec. 6] for a
more detailed discussion on the statistical mechanics motivation behind
Schrödinger’s question.

Let $\mu,\,\nu$ be two probability measures on $\mathbb{R}^{d}$ satisfying
$\mathrm{Ent}(\mu)+\mathrm{Ent}(\nu)<+\infty$, where $\mathrm{Ent}$ denotes
the relative entropy with respect to the Lebesgue measure. Given a horizon
$T>0$, SP consists in solving the optimization problem

(SP) $\text{minimize }\mathscr{H}(\pi|\mathrm{R}_{0,T})\text{ under the
constraint }\pi\in\Pi(\mu,\nu)\;,$

where $\mathscr{H}(\pi|\mathrm{R}_{0,T})$ is the relative entropy between the
coupling $\pi$ and $\mathrm{R}_{0,T}$ defined by

(5) $\mathrm{R}_{0,T}(\mathrm{d}x\mathrm{d}y)=(2\pi
T)^{-d/2}\,\exp\biggl{(}-\frac{|x-y|^{2}}{2T}\biggr{)}\mathrm{d}x\,\mathrm{d}y\;.$

Indeed, the identity

(6)
$T\mathscr{H}(\pi|\mathrm{R}_{0,T})-T\mathrm{Ent}(\mu)-T\mathrm{Ent}(\nu)-\frac{dT}{2}\log(2\pi
T)=\int\frac{|x-y|^{2}}{2}\mathrm{d}\pi+T\mathscr{H}(\pi|\mu\otimes\nu)\,\text{
for any}\;\pi\in\Pi(\mu,\nu)$

immediately implies the equivalence between SP and EOT .

Based on the formulation (SP) it has been pointed out [5] that Sinkhorn’s
algorithm is a special case of the Bregman’s iterated projection algorithm for
relative entropy. Indeed, in the current setup Bregman’s iterated projection
algorithm produces two sequence of plans
$(\pi^{n,n},\pi^{n+1,n})_{n\in\mathbb{N}^{*}}$ starting from a positive
measure $\pi^{0,0}$ according to the following recursion:

(7)
$\pi^{n+1,n}\coloneqq{\arg\min}_{\Pi(\mu,\star)}\mathscr{H}(\cdot|\pi^{n,n})\,,\qquad\pi^{n+1,n+1}\coloneqq{\arg\min}_{\Pi(\star,\nu)}\mathscr{H}(\cdot|\pi^{n+1,n})\,,$

where $\Pi(\mu,\star)$ (resp. $\Pi(\star,\nu)$) is the set of probability
measures $\pi$ on $\mathbb{R}^{2d}$ such that the first marginal is $\mu$,
i.e., $(\mathrm{proj}_{x})_{\sharp}\pi=\mu$ (resp. the second marginal is
$\nu$, i.e., $(\mathrm{proj}_{y})_{\sharp}\pi=\nu$). It is relatively easy to
show (see [47, Section 6]), starting from
$\pi^{0,0}(\mathrm{d}x\mathrm{d}y)\propto\exp(-\psi^{0}(y)-\varphi^{0}(x))\mathrm{d}x\mathrm{d}y$,
that the iterates in (7) are related to Sinkhorn potentials through

(8) $\displaystyle\textstyle\pi^{n+1,n}(\mathrm{d}x\mathrm{d}y)$
$\displaystyle\textstyle\propto\exp(-\nicefrac{{|x-y|^{2}}}{{2T}}-\varphi^{n+1}(x)-\psi^{n}(y))\mathrm{d}x\mathrm{d}y\;,$
(9) $\displaystyle\textstyle\pi^{n+1,n+1}(\mathrm{d}x\mathrm{d}y)$
$\displaystyle\textstyle\propto\exp(-\nicefrac{{|x-y|^{2}}}{{2T}}-\varphi^{n+1}(x)-\psi^{n+1}(y))\mathrm{d}x\mathrm{d}y\;.$

In the sequel, we will refer to the couplings
$(\pi^{n,n},\pi^{n+1,n})_{n\in\mathbb{N}^{*}}$ as Sinkhorn plans.

###  2.1. A first (informal) statement of our results and assumptions

To state the assumptions we impose on marginal distributions, we introduce the
integrated convexity profile $\kappa_{U}:\mathbb{R}^{*}_{+}\to\mathbb{R}$ of a
differentiable function $U:\mathbb{R}^{d}\to\mathbb{R}$ as the function

(10) $\kappa_{U}(r)\coloneqq\inf\biggl{\\{}\frac{\langle\nabla U(x)-\nabla
U(y),\,x-y\rangle}{|x-y|^{2}}\quad\colon\quad|x-y|=r\biggr{\\}}\,.$

Likewise, for a distribution $\zeta$ represented as a Gibbs measure as
$\zeta(\mathrm{d}x)\propto\exp(-U)\mathrm{d}x$, we simply write
$\kappa_{\zeta}$ for $\kappa_{U}$ and referred to this function as to the
integrated log-concavity profile of $\zeta$. The function $\kappa_{U}$ is
often employed to quantify ergodicity of stochastic differential equations
whose drift field is $-\nabla U$, see [26]. The term integrated convexity
profile we coin here is motivated by the observation that
$\kappa_{U}(r)\geq\alpha$ if and only if for any $x,v\in\mathbb{R}^{d}$,
$\left|v\right|=1$,
$\int_{0}^{r}\langle\nabla^{2}U(x+hv)v,v\rangle\mathrm{d}h\geq\alpha r$. The
integrated concavity profile of $U$ is defined in a similar way as
$\ell_{U}=-\kappa_{-U}$, and for $\zeta$ of the form
$\zeta(\mathrm{d}x)\propto\exp(-U)\mathrm{d}x$ we set $\ell_{\zeta}=\ell_{U}$.
Our main results apply when the marginals $\mu,\nu$ satisfy the following
property for $\zeta\in\\{\mu,\nu\\}$,

(11)
$\begin{split}\liminf_{r\rightarrow+\infty}\kappa_{\zeta}(r)>0\;,\quad\liminf_{r\rightarrow
0}r\kappa_{\zeta}(r)=0\;.\end{split}$

We are ready to state our main contributions informally.

######  Theorem 1 (Informal).

Assume that (11) holds for the marginals $\mu,\nu$ specified by (2). In
addition, suppose that $\varphi^{0}\equiv 0$. Then, there exist $A,B,C\geq 0$,
$T_{1}>T_{0}$ and $\gamma\in\left[0,1\right)$ such that for $n\in\mathbb{N}$
and $T>T_{0}$,

(12)
$\displaystyle\textstyle\int|\nabla\varphi^{n}-\nabla\varphi^{\star}|\mathrm{d}\mu+\int|\nabla\psi^{n}-\nabla\psi^{\star}|\mathrm{d}\nu\textstyle$
$\displaystyle\leq C\gamma^{n}$ (13)
$\displaystyle\mathbf{W}_{1}(\pi^{n,n},\,\pi^{\star})+\mathbf{W}_{1}(\pi^{n+1,n},\,\pi^{\star})$
$\displaystyle\leq C\gamma^{n}$

and for $n\in\mathbb{N}$, $x\in\mathbb{R}^{d}$ and $T>T_{1}$

(14)
$\begin{split}|\nabla\varphi^{n}-\nabla\varphi^{\star}|(x)+|\nabla\psi^{n}-\nabla\psi^{\star}|(x)&\leq
C\gamma^{n}(A|x|+B)\\\
\textstyle\mathscr{H}(\pi^{n,n}|\pi^{\star})+\mathscr{H}(\pi^{\star}|\pi^{n,n})+\mathscr{H}(\pi^{n+1,n}|\pi^{\star})+\mathscr{H}(\pi^{\star}|\pi^{n+1,n})&\leq
C\gamma^{n}\;,\end{split}$

where $\\{(\varphi^{n},\psi^{n})\\}_{n\in\mathbb{N}}$ and
$\\{(\pi^{n+1,n},\pi^{n,n})\\}_{n\in\mathbb{N}}$ are defined in (4) and
(8)-(9) respectively.

To complement this informal statement, let us also mention that we deduce from
the exponential convergence of the gradient
$\\{(\nabla\varphi^{n},\nabla\psi^{n})\\}_{n\in\mathbb{N}^{\star}}$ the
exponential convergence of the Sinkhorn iterates
$\\{(\varphi^{n},\psi^{n})\\}_{n\in\mathbb{N}^{\star}}$ to the Schrödinger
potentials $(\varphi^{\star},\psi^{\star})$. Finally, we also establish
pointwise exponential convergence of the Hessian matrices
$\\{(\nabla^{2}\varphi^{n},\nabla^{2}\psi^{n})\\}_{n\in\mathbb{N}^{\star}}$ to
$(\nabla^{2}\varphi^{\star},\nabla^{2}\psi^{\star})$.

Furthermore, while we give in Theorem 1 a qualitative statement only, i.e.,
with non-explicit constants $C$ and $\gamma$, we will give in the next section
quantitative versions of the exponential convergence bounds provided in
Theorem 1.

Convergence rates can be improved if we have additional information on the
integrated concavity profile such as for $\zeta\in\\{\mu,\nu\\}$,

(15) $\limsup_{r\rightarrow+\infty}\ell_{\zeta}(r)<\infty\;,$

or some weaker versions of the above (see the precise statement of A 1 below).
However, we stress that (15) and variants thereof are not necessary
conditions. To obtain the sharpest possible convergence rates, we consider in
the next section two set of functions denoted by $\mathcal{G}$ and
$\tilde{\mathcal{G}}$, that will appear in our assumptions on
$\kappa_{\zeta},\ell_{\zeta}$.

###  2.2. Main results and proof strategy

#### Assumptions

In this section we provide a rigorous statement of the main assumptions we
impose on the marginals $\mu,\nu$. In view of this, it is convenient to
introduce some notation and terminology. We begin by introducing two sets of
functions that will appear in our conditions:
$\mathcal{G}\coloneqq\\{g\in\mathrm{C}^{2}((0,+\infty),\,\mathbb{R}_{+})\,:\,\text{$g$
satisfies \eqref{eq:hyp_mcg}}\\}$ with

($\mathbf{H}_{\mathcal{G}}$) $(r\mapsto r^{1/2}\,g(r^{1/2}))\text{ is non-
decreasing and concave,}\quad\lim_{r\downarrow 0}r\,g(r)=0\;$

and its subset

(16) $\tilde{\mathcal{G}}\coloneqq\biggl{\\{}g\in\mathcal{G}\;\text{ bounded
and s.t.}\quad\lim_{r\downarrow 0}g(r)=0\;,\quad g{{}^{\prime}}\geq
0\quad\text{ and }\quad 2g{{}^{\prime\prime}}+g\,g{{}^{\prime}}\leq
0\biggr{\\}}\subset\mathcal{G}\;.$

The above classes of functions are non-empty and in particular
$\tilde{\mathcal{G}}$ contains $r\mapsto 2\tanh(r/2)$. Though it may not
appear as the most natural at first sight, these choices will become clear in
light of proofs. Indeed, the sets $\mathcal{G}$ and $\tilde{\mathcal{G}}$
enjoy special invariance properties (see Theorem 7 and Section 3) under the
mapping

(17) $g\mapsto-\log P_{T}\exp(-g)\;,$

upon which the proof of the lower bounds on the integrated convexity profiles
of potentials (see Section 3) are built.

Note that the properties prescribed in the definition of $\tilde{\mathcal{G}}$
in particular imply that its elements are sublinear (i.e.,
$\sup_{r>0}g(r)/r<+\infty$) concave functions on $(0,+\infty)$. Indeed, any
$g\in\tilde{\mathcal{G}}$ is clearly non-negative and non-decreasing, which
combined with the differential inequality, implies its concavity. Finally,
since $g(0^{+})=0$, the sublinearity of any $g\in\tilde{\mathcal{G}}$ follows.
We say that a potential $U:\mathbb{R}^{d}\to\mathbb{R}$ is asymptotically
strongly convex if there exist $\alpha_{U}\in\mathbb{R}_{+}^{*}$ and
$\tilde{g}_{U}\in\tilde{\mathcal{G}}$ such that

(18) $\kappa_{U}(r)\geq\alpha_{U}-r^{-1}\,\tilde{g}_{U}(r)$

holds for all $r\geq 0$. We consider the set of asymptotically strongly log-
concave probability measures

(19)
$\mathcal{P}_{\mathrm{alc}}(\mathbb{R}^{d})=\\{\zeta(\mathrm{d}x)=\exp(-U(x))\mathrm{d}x:\,U\in\mathrm{C}^{2}(\mathbb{R}^{d}),\,\text{$U$
asymptotically strongly convex}\\}$

Note that as soon as $\zeta(\mathrm{d}x)=\exp(-U(x))\mathrm{d}x$ there exist
$\beta_{U}\in(0,+\infty]$ and $g\in\mathcal{G}$ such that

(20) $\ell_{U}(r)\leq\beta_{U}+r^{-1}\,{g_{U}}(r)\;$

holds for all $r\geq 0$. This is trivially true as we can choose
$\beta_{U}=+\infty$. However, if the above holds for some $\beta_{U}<+\infty$,
we obtain better lower bounds on the integrated convexity profile of Sinkhorn
potentials and all contraction rates appearing in our main results are better
than those obtained for $\beta_{U}=+\infty$.

######  Remark .

Let us remark that the class of probability measures
$\mathcal{P}_{\mathrm{alc}}(\mathbb{R}^{d})$ contains in particular
probability measures $\zeta$ associated with potentials $U$ satisfying

(21) $\kappa_{U}(r)\geq\begin{cases}\alpha\quad&\text{ if }r>R\\\ \alpha-
C_{U}\quad&\text{ if }r\leq R\,,\end{cases}$

for $\alpha,C_{U},R>0$. Indeed, [18, Proposition 5.1] implies that
$\kappa_{\zeta}(r)\geq\alpha-r^{-1}\;{\tilde{g}}_{L}(r)$ where
${\tilde{g}}_{L}\in\tilde{\mathcal{G}}$ is given by

(22) ${\tilde{g}}_{L}(r)\coloneqq 2\,(L)^{1/2}\,\tanh(r\,L^{1/2}/2)\;\text{
with }\quad
L\coloneqq\inf\\{\bar{L}\,\colon\,R^{-1}\,{\tilde{g}}_{\bar{L}}(R)\geq
C_{U}\\}\;.$

Finally, it is worth mentioning that (21) holds if $U$ can be expressed as the
sum of a strongly convex function and a Lipschitz function with second
derivative bounded from below.

We have now all the concepts and notations to introduce our assumptions.

######  A 1.

The marginals $\mu,\nu$ specified by (2) belong to the set
$\mathcal{P}_{\mathrm{alc}}(\mathbb{R}^{d})$ defined at (19) and have finite
relative entropy with respect to the Lebesgue measure $\mathrm{Leb}$,
$\mathrm{Ent}(\mu),\,\mathrm{Ent}(\nu)<+\infty$.

Under A 1, we denote by $\alpha_{\mu},\beta_{\mu}$, $\tilde{g}_{\mu}$ and
$g_{\mu}$ (resp. $\alpha_{\nu},\beta_{\nu}$, $\tilde{g}_{\nu}$ and $g_{\nu}$)
the constants and functions associated with $\mu$ (resp. $\nu$) such that (18)
and (20) hold for $U_{\mu}$ (resp. $U_{\nu}$).

We prove in Section 3 that the previous set of assumptions imply lower bounds
on the integrated convexity profile of Sinkhron and Schrödinger potentials.
This result is at the core of our proof strategy and generalize the findings
[18], that are valid for Schrödinger but not Sinkhorn potentials and
restricted to the case when $g_{\mu}=g_{\nu}=0$ and
$\tilde{g}_{\mu},\tilde{g}_{\nu}$ are of the form (22), see Remark 3 for more
explanations.

#### Proof strategy

The starting point of our discussion is the following well known result, see
[54, Proposition 2] or [15] for instance. Define
$(x,A)\mapsto\pi_{T}^{x,h}(A)$ as the Markov kernel on
$\mathbb{R}^{d}\times\mathcal{B}(\mathbb{R}^{d})$ whose transition density
w.r.t. Lebesgue measure is proportional to
$(x,y)\mapsto\exp(-h(y)-|x-y|^{2}/(2T))$, i.e., for any $x\in\mathbb{R}^{d}$,
$\pi^{x,h}_{T}$ is defined through

(23)
$\pi_{T}^{x,h}(\mathrm{d}y)\propto\exp\biggl{(}-\frac{|y-x|^{2}}{2T}-h(y)\biggr{)}\mathrm{d}y\;.$

######  Proposition .

Assume A 1 holds. Then for any $n\in\mathbb{N}^{*}$,
$h\in\\{\psi^{\star},\,\psi^{n},\varphi^{\star},\varphi^{n}\\}$ it holds

(24) $\nabla\log
P_{T}\mathrm{e}^{-h}(x)=\,T^{-1}\int(y-x)\,\pi_{T}^{x,h}(\mathrm{d}y)\;,\quad\nabla^{2}\log
P_{T}\mathrm{e}^{-h}(x)=\,-T^{-1}\operatorname{Id}+T^{-2}\,\mathrm{Cov}(\pi_{T}^{x,h})\;.$

The proof of this technical result is given for completeness in Section C.1.
To give a concise sketch of our proof strategy, we record here two
observations. The first is that for any $x\in\mathbb{R}^{d}$ the conditional
distribution, $\pi_{T}^{x,h}$ defined in (23) is the invariant probability
measure for the SDE

(25) $\mathrm{d}Y_{t}=-\biggl{(}\frac{Y_{t}-x}{2T}+\frac{1}{2}\nabla
h(Y_{t})\biggr{)}\mathrm{d}t+\mathrm{d}B_{t}\,.$

The second one is that, defining the adjusted marginals produced along
Sinkhorn as the probability measures

(26) $\text{ $\mu^{n}\coloneqq(\mathrm{proj}_{x})_{\sharp}\pi^{n,n}$ and
$\nu^{n}\coloneqq(\mathrm{proj}_{y})_{\sharp}\pi^{n+1,n}$}\;,$

using (1)-(8)-(9), $\pi^{\star}\in\Pi(\mu,\nu)$ and
$\pi^{n,n}\in\Pi(\star,\nu)$, the Schrödinger and Sinkhorn plans can be
written as

(27)
$\pi^{\star}(\mathrm{d}x\mathrm{d}y)=\mu(\mathrm{d}x)\otimes\pi^{x,\psi^{\star}}_{T}(\mathrm{d}y)\quad\text{
and
}\quad\pi^{n,n}(\mathrm{d}x\mathrm{d}y)=\mu^{n}(\mathrm{d}x)\otimes\pi^{x,\psi^{n}}_{T}(\mathrm{d}y)\;,$

and a direct computation shows that

(28)
$\displaystyle\int\pi_{T}^{x,\psi^{\star}}(\mathrm{d}y)\;\mu(\mathrm{d}x)=\nu(\mathrm{d}y)\;,\quad\quad\int\pi_{T}^{x,\psi^{n}}(\mathrm{d}y)\;\mu^{n}(\mathrm{d}x)=\nu(\mathrm{d}y)\;,$
$\displaystyle\int\pi_{T}^{y,\varphi^{\star}}(\mathrm{d}x)\;\nu(\mathrm{d}y)=\mu(\mathrm{d}x)\;,\quad\quad\int\pi_{T}^{y,\varphi^{n}}(\mathrm{d}x)\;\nu^{n}(\mathrm{d}y)=\mu(\mathrm{d}x)\;.$

Let us now proceed to give a minimal description of our proof strategy. From
(4) and the formula (24) with choices $h=\psi^{n},\,\psi^{\star}$ we deduce
the upper bound

(29) $\int|\nabla\varphi^{n+1}-\nabla\varphi^{\star}|(x)\;\mu(\mathrm{d}x)\leq
T^{-1}\,\int\mathbf{W}_{1}(\pi_{T}^{x,\psi^{n}},\pi_{T}^{x,\psi^{\star}})\;\mu(\mathrm{d}x)\;,$

where $\mathbf{W}_{1}(\cdot,\cdot)$ denotes the Wasserstein distance of order
one. Combining together the lower bounds on the integrated convexity profiles
$\kappa_{\psi^{n}}$ (see Section 3), the representation of
$\pi_{T}^{x,\psi^{n}},\pi_{T}^{x,\psi^{\star}}$ as invariant measures for the
SDE (25) and coupling techniques (see Section 4), we obtain at Section 4 that
the key estimate

(30)
$\mathbf{W}_{1}(\pi_{T}^{x,\psi^{n}},\pi_{T}^{x,\psi^{\star}})\leq\gamma^{\nu}_{n}\,\int|\nabla\psi^{n}-\nabla\psi^{\star}|(y)\;\pi^{x,\psi^{\star}}_{T}(\mathrm{d}y)\,$

holds uniformly on $x\in\mathbb{R}^{d}$ for some $\gamma^{\nu}_{n}>0$.
Integrating the key estimate (30) w.r.t. $\mu$ and invoking (28) yields

(31)
$\int|\nabla\varphi^{n+1}-\nabla\varphi^{\star}|(x)\;\mu(\mathrm{d}x)\leq\frac{\gamma^{\nu}_{n}}{T}\int|\nabla\psi^{n}-\nabla\psi^{\star}|(y)\nu(\mathrm{d}y)\;.$

Repeating the same argument but exchanging the roles of
$\psi^{n},\psi^{\star}$ and $\varphi^{n},\varphi^{\star}$ we obtain

(32)
$\int|\nabla\psi^{n}-\nabla\psi^{\star}|(y)\;\nu(\mathrm{d}y)\leq\frac{\gamma^{\mu}_{n-1}}{T}\int|\nabla\varphi^{n}-\nabla\varphi^{\star}|(x)\mu(\mathrm{d}x)$

for some $\gamma^{\mu}_{n-1}>0$. Combining (32) with (31) allows to establish
exponential convergence provided $T^{-2}\gamma^{\nu}_{n}\gamma^{\mu}_{n-1}<1$
for $n$ large enough.

#### Statement of the main results

Our first contribution is the $\mathrm{L}^{1}$-convergence of Sinkhorn
iterates and $\mathbf{W}_{1}$-convergence of Sinkhorn plans.

######  Theorem 2.

Assume that A 1 holds. Then for any $n\geq 1$

(33)
$\displaystyle\int|\nabla\varphi^{n}-\nabla\varphi^{\star}|\mathrm{d}\mu\leq$
$\displaystyle\;\frac{T}{\gamma_{n-1}^{\mu}}\;\prod_{k=0}^{n-1}\frac{\gamma_{k}^{\mu}\;\gamma_{k}^{\nu}}{T^{2}}\;\int|\nabla\psi^{0}-\nabla\psi^{\star}|\mathrm{d}\nu\;,$
$\displaystyle\int|\nabla\psi^{n}-\nabla\psi^{\star}|\mathrm{d}\nu\leq$
$\displaystyle\;\prod_{k=0}^{n-1}\frac{\gamma_{k}^{\mu}\;\gamma_{k}^{\nu}}{T^{2}}\;\int|\nabla\psi^{0}-\nabla\psi^{\star}|\mathrm{d}\nu\;,$

and

(34) $\displaystyle\mathbf{W}_{1}(\pi^{n,n},\,\pi^{\star})\leq$
$\displaystyle\;T\;\prod_{k=0}^{n-1}\frac{\gamma_{k}^{\mu}\;\gamma_{k}^{\nu}}{T^{2}}\;\int|\nabla\psi^{0}-\nabla\psi^{\star}|\mathrm{d}\nu\;,$
$\displaystyle\mathbf{W}_{1}(\pi^{n+1,n},\,\pi^{\star})\leq$
$\displaystyle\;\gamma_{n}^{\nu}\;\prod_{k=0}^{n-1}\frac{\gamma_{k}^{\mu}\;\gamma_{k}^{\nu}}{T^{2}}\;\int|\nabla\psi^{0}-\nabla\psi^{\star}|\mathrm{d}\nu\;,$

where $(\gamma_{k}^{\mu})_{k\in\mathbb{N}}$ and
$(\gamma_{k}^{\nu})_{k\in\mathbb{N}}$ are non-negative non-increasing
sequences given in Section 2.2, depending on
$\alpha_{\mu},\,\beta_{\nu},\,{\tilde{g}_{\mu}},\,{g}_{\nu},\,T$ and on
$\alpha_{\nu},\,\beta_{\mu},\,{\tilde{g}_{\nu}},\,{g}_{\mu},\,T$ respectively.
Denoting
$\gamma_{\infty}^{\zeta}\coloneqq\lim_{k\to+\infty}\gamma_{k}^{\zeta}$ for
$\zeta\in\\{\mu,\nu\\}$, as a corollary, as soon as
$T^{2}>\gamma_{\infty}^{\mu}\,\gamma_{\infty}^{\nu}$, the asymptotic rate is
strictly less than one, and for any
$T^{-2}\gamma_{\infty}^{\mu}\,\gamma_{\infty}^{\nu}<\lambda<1$, there exists
$C\geq 0$ such that for any $n\in\mathbb{N}^{*}$,

(35)
$\displaystyle\int|\nabla\varphi^{n}-\nabla\varphi^{\star}|\mathrm{d}\mu+\int|\nabla\psi^{n}-\nabla\psi^{\star}|\mathrm{d}\nu\leq
C\lambda^{n}\int|\nabla\psi^{0}-\nabla\psi^{\star}|\mathrm{d}\nu\;.$

######  Remark  (Explicit rates).

For any $\mathfrak{p}\in\\{\mu,\nu\\}$, and for any $k\in\mathbb{N}$ it holds

(36)
$\gamma^{\mathfrak{p}}_{k}\coloneqq\frac{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)^{2}}{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}\biggl{(}\left\|\tilde{g}_{\mathfrak{p}}\right\|_{\infty}\biggl{(}\frac{1}{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)}+\frac{2}{\alpha_{\mathfrak{p},k}+T^{-1}}\biggr{)}\biggr{)}^{2}}\frac{1}{\alpha_{\mathfrak{p},k}+T^{-1}+\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)}\;$

where
$(\alpha_{\mu,n})_{n\in\mathbb{N}}\subseteq(\alpha_{\mu}-T^{-1},\alpha_{\mu}-T^{-1}+(\beta_{\nu}\,T^{2})^{-1}]$
and
$(\alpha_{\nu,n})_{n\in\mathbb{N}}\subseteq(\alpha_{\nu}-T^{-1},\alpha_{\nu}-T^{-1}+(\beta_{\mu}\,T^{2})^{-1}]$
are monotone increasing sequences, built in Theorem 7. If
$\alpha_{\varphi}^{\star}$ and $\alpha_{\psi}^{\star}$ denote their respective
limits, then the asymptotic rates of convergence read as

(37) $\displaystyle\gamma^{\mu}_{\infty}\coloneqq$
$\displaystyle\;\frac{\tilde{g}_{\mu}{{}^{\prime}}(0)^{2}}{\tilde{g}_{\mu}{{}^{\prime}}\biggl{(}\left\|\tilde{g}_{\mu}\right\|_{\infty}\biggl{(}\frac{1}{\tilde{g}_{\mu}{{}^{\prime}}(0)}+\frac{2}{\alpha_{\varphi^{\star}}+T^{-1}}\biggr{)}\biggr{)}^{2}}\frac{1}{\alpha_{\varphi^{\star}}+T^{-1}+\tilde{g}_{\mu}{{}^{\prime}}(0)}\;,$
$\displaystyle\gamma^{\nu}_{\infty}\coloneqq$
$\displaystyle\;\frac{\tilde{g}_{\nu}{{}^{\prime}}(0)^{2}}{\tilde{g}_{\nu}{{}^{\prime}}\biggl{(}\left\|\tilde{g}_{\nu}\right\|_{\infty}\biggl{(}\frac{1}{\tilde{g}_{\nu}{{}^{\prime}}(0)}+\frac{2}{\alpha_{\psi^{\star}}+T^{-1}}\biggr{)}\biggr{)}^{2}}\frac{1}{\alpha_{\psi^{\star}}+T^{-1}+\tilde{g}_{\nu}{{}^{\prime}}(0)}\;.$

From the above expressions it can be deduced that the condition for
exponential convergence $T^{-2}\gamma^{\mu}_{\infty}\gamma^{\nu}_{\infty}<1$
is always satisfied for large enough values of $T$. To be more precise, a
sufficient condition for the exponential convergence of Sinkhorn algorithm is
the following

(38)
$T^{2}>\frac{\tilde{g}_{\mu}{{}^{\prime}}(0)\,(\alpha_{\mu}+\tilde{g}_{\mu}{{}^{\prime}}(0))^{-1}}{\tilde{g}_{\mu}{{}^{\prime}}\biggl{(}\left\|\tilde{g}_{\mu}\right\|_{\infty}\biggl{(}\frac{1}{\tilde{g}_{\mu}{{}^{\prime}}(0)}+\frac{2}{\alpha_{\mu}}\biggr{)}\biggr{)}^{2}}\;\frac{\tilde{g}_{\nu}{{}^{\prime}}(0)\,(\alpha_{\nu}+\tilde{g}_{\nu}{{}^{\prime}}(0))^{-1}}{\tilde{g}_{\nu}{{}^{\prime}}\biggl{(}\left\|\tilde{g}_{\nu}\right\|_{\infty}\biggl{(}\frac{1}{\tilde{g}_{\nu}{{}^{\prime}}(0)}+\frac{2}{\alpha_{\nu}}\biggr{)}\biggr{)}^{2}}\;.$

Expressions simplify when considering $\tilde{g}_{\mathfrak{p}}$ as in (22) in
Section 2.2, i.e.

(39)
$\tilde{g}_{\mathfrak{p}}(r)=2\,(L_{\mathfrak{p}})^{1/2}\,\tanh(r\,L_{\mathfrak{p}}^{1/2}/2)$

for some $L_{\mathfrak{p}}\geq 0$, since the previous asymptotic rates would
read as

(40) $\displaystyle\gamma_{\infty}^{\mu}=$
$\displaystyle\;\frac{1}{\alpha_{\varphi^{\star}}+T^{-1}+L_{\mu}}\;\cosh^{4}\biggl{(}2(L_{\mu})^{1/2}\biggl{(}\frac{1}{L_{\mu}}+\frac{2}{\alpha_{\varphi^{\star}}+T^{-1}}\biggr{)}\biggr{)}\;,$
$\displaystyle\gamma_{\infty}^{\nu}=$
$\displaystyle\;\frac{1}{\alpha_{\psi^{\star}}+T^{-1}+L_{\nu}}\;\cosh^{4}\biggl{(}2(L_{\nu})^{1/2}\biggl{(}\frac{1}{L_{\nu}}+\frac{2}{\alpha_{\psi^{\star}}+T^{-1}}\biggr{)}\biggr{)}\;,$

whereas (38) would read as

(41)
$T^{2}>\frac{\cosh^{4}\biggl{(}\frac{2}{L_{\mu}^{1/2}}+\frac{4\,L_{\mu}^{1/2}}{\alpha_{\mu}}\biggr{)}}{\alpha_{\mu}+L_{\mu}}\;\frac{\cosh^{4}\biggl{(}\frac{2}{L_{\nu}^{1/2}}+\frac{4\,L_{\nu}^{1/2}}{\alpha_{\nu}}\biggr{)}}{\alpha_{\nu}+L_{\nu}}\;.$

Finally, explicit expressions for $(\gamma^{\mu}_{k})_{k\in\mathbb{N}}$,
$(\gamma^{\nu}_{k})_{k\in\mathbb{N}}$ and their limits fully simplify when
$\mu$ and $\nu$ are strongly log-concave, i.e., when
${\tilde{g}}_{\mu},{\tilde{g}}_{\nu},{g}_{\mu},{g}_{\nu}$ all equal to the
null function. We focus on this particular scenario in Appendix A where we
give simple expression for
$(\gamma^{\mu}_{k})_{k\in\mathbb{N}},\,(\gamma^{\nu}_{k})_{k\in\mathbb{N}},\,\gamma_{\infty}^{\mu},\,\gamma_{\infty}^{\nu}$.
In particular we show that Sinkhorn’s algorithm converges as soon as

(42)
$T>\frac{\beta_{\mu}\beta_{\nu}-\alpha_{\mu}\alpha_{\nu}}{\sqrt{\alpha_{\mu}\,\beta_{\mu}\,\alpha_{\nu}\,\beta_{\nu}\,(\alpha_{\mu}+\beta_{\mu})(\alpha_{\nu}+\beta_{\nu})}}\;.$

As $\beta_{\mu}=\beta_{\nu}=+\infty$, this expression simply reads as
$T>(\alpha_{\mu}\,\alpha_{\nu})^{-1/2}$.

By relying on the previous result, we are also able to prove exponential
convergence of the $\mathrm{L}^{1}(\mu^{n})$-norms and
$\mathrm{L}^{1}(\nu^{n})$-norms of the difference of the gradients, along the
adjusted marginals $\mu^{n},\,\nu^{n}$. We have postponed this result to
Section 5.1.

Our coupling approach can also be employed in order to prove pointwise
convergence result. In order to establish such a result the assumptions we
impose on the regularization parameter $T$ are more stringent than the ones we
need for Theorem 2. In addition we must impose an additional assumption (43)
on $\psi^{0}$. As we explain below there is a natural choice for $\psi^{0}$
that guarantees that (43) holds.

######  Theorem 3.

Assume A 1 holds, $\psi^{0}\in\mathrm{C}^{1}(\mathbb{R}^{d})$ and that there
exist two positive constants $A,\,B>0$ such that for any $x\in\mathbb{R}^{d}$

(43) $|\nabla\psi^{0}-\nabla\psi^{\star}|(x)\leq A\;|x|+B\;.$

Then for any $n\in\mathbb{N}^{*}$ and $x\in\mathbb{R}^{d}$

(44) $\displaystyle|\nabla\varphi^{n}-\nabla\varphi^{\star}|(x)\leq$
$\displaystyle\,\frac{T}{\hat{\gamma}^{\mu}_{n-1}}\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\;\hat{\gamma}_{k}^{\nu}}{T^{2}}(A|x|+B)$
$\displaystyle|\nabla\psi^{n}-\nabla\psi^{\star}|(x)\leq$
$\displaystyle\,\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\;\hat{\gamma}_{k}^{\nu}}{T^{2}}\;(A|x|+B)\;,$

where

(45)
${\hat{\gamma}_{k}^{\mu}}\coloneqq\gamma_{k}^{\mu}\;\max\biggl{\\{}(T\,\alpha_{\varphi^{\star}}+1)^{-1}\;,\biggl{(}1+\frac{A}{B}\frac{1+\left\|\tilde{g}_{\mu}\right\|_{\infty}+|\nabla\varphi^{\star}(0)|}{\alpha_{\varphi^{\star}}+T^{-1}}\biggr{)}\biggr{\\}}\;,$

and

(46)
${\hat{\gamma}_{k}^{\nu}}\coloneqq\gamma_{k}^{\nu}\;\max\biggl{\\{}(T\,\alpha_{\psi^{\star}}+1)^{-1}\;,\biggl{(}1+\frac{A}{B}\frac{1+\left\|\tilde{g}_{\nu}\right\|_{\infty}+|\nabla\psi^{\star}(0)|}{\alpha_{\psi^{\star}}+T^{-1}}\biggr{)}\biggr{\\}}\;,$

where
$(\gamma_{k}^{\mu})_{k\in\mathbb{N}},\,(\gamma_{k}^{\nu})_{k\in\mathbb{N}}$
are the non-decreasing non-negative sequences defined in Section 2.2,
$\alpha_{\varphi^{\star}}\in(\alpha_{\mu}-T^{-1},\alpha_{\mu}-T^{-1}+(\beta_{\nu}\,T^{2})^{-1}]$
and
$\alpha_{\psi^{\star}}\in(\alpha_{\nu}-T^{-1},\alpha_{\nu}-T^{-1}+(\beta_{\mu}\,T^{2})^{-1}]$
are given in Theorem 7. In particular, denoting
$\gamma_{\infty}^{\zeta}\coloneqq\lim_{k\to+\infty}\gamma_{k}^{\zeta}$ and
$\hat{\gamma}_{\infty}^{\zeta}\coloneqq\lim_{k\to+\infty}\hat{\gamma}_{k}^{\zeta}$
for $\zeta\in\\{\mu,\nu\\}$, if $T$ is large enough, e.g., if

(47)
$T>\max\biggl{\\{}\alpha_{\mu}^{-1}\,,\;\alpha_{\nu}^{-1}\,,\;\gamma_{\infty}^{\mu}+\frac{\gamma_{\infty}^{\mu}A}{\alpha_{\mu}B}(1+\|{\tilde{g}_{\mu}}\|_{\infty}+|\nabla\varphi^{\star}(0)|)\,,\;\gamma_{\infty}^{\nu}+\frac{\gamma_{\infty}^{\nu}A}{\alpha_{\nu}B}(1+\|{\tilde{g}_{\nu}}\|_{\infty}+|\nabla\psi^{\star}(0)|)\biggr{\\}}\;,$

then $T^{2}>\hat{\gamma}_{\infty}^{\mu}\,\hat{\gamma}_{\infty}^{\nu}$ and as a
result for any
$\lambda\in(T^{-2}\,\hat{\gamma}_{\infty}^{\mu}\hat{\gamma}_{\infty}^{\nu},\,1)$,
there exists $C\geq 0$ such that for any $x\in\mathbb{R}^{d}$ and
$n\in\mathbb{N}^{*}$,

(48)
$\displaystyle|\nabla\varphi^{n}-\nabla\varphi^{\star}|(x)+|\nabla\psi^{n}-\nabla\psi^{\star}|(x)\leq\,C\,\lambda^{n}\;(A|x|+B)\;.$

As for Theorem 2, expressions of
$(\hat{\gamma}_{k}^{\mu})_{k\in\mathbb{N}},\,(\hat{\gamma}_{k}^{\nu})_{k\in\mathbb{N}}$
simplify as $U_{\mu}$ and $U_{\nu}$ are strongly convex. These expressions are
given in Appendix A.

We stress that the previous theorem holds for any smooth initialization
$\psi^{0}\in\mathrm{C}^{1}(\mathbb{R}^{d})$ satisfying (43). A common choice
would be starting at $\psi^{0}=U_{\nu}$, which corresponds in terms of EOT and
entropic potentials as having a null-initialization for Sinkhorn’s
algorithm.111Indeed let us stress here that if $\varphi^{\star}$ and
$\psi^{\star}$ denote the Schrödinger potentials, then the functions $-T\log
P_{T}\exp(-\psi^{\star})$ and $-T\log P_{T}\exp(-\varphi^{\star})$ correspond
to the couple of entropic potentials commonly considered in EOT literature.
This choice agrees with the usual normalization imposed to Sinkhorn iterates
when studying its convergence [24, 11, 23, 10]. Let us also point out that if
one starts Sinkhorn’s algorithm one step before with the null initialization
$\varphi^{0}\coloneqq 0$, then at the first iteration we immediately get
$\psi^{0}=U_{\nu}$. Under this choice, from (3) we deduce that

(49) $\psi^{0}-\psi^{\star}=-\log P_{T}\exp(-\varphi^{\star})$

hence

(50)
$|\nabla\psi^{0}-\nabla\psi^{\star}|(x)=\frac{1}{T}\;\biggl{|}\int(y-x)\,\pi_{T}^{x,\varphi^{\star}}(\mathrm{d}y)\biggr{|}\leq\frac{1}{T}\int|y-x|\,\pi_{T}^{x,\varphi^{\star}}(\mathrm{d}y)\leq\frac{|x|}{T}+\frac{1}{T}\int|y|\,\mathrm{d}\pi_{T}^{x,\varphi^{\star}}(y)\;.$

The latter combined with the bound we are going to give later in (174) shows
that the initialization $\psi^{0}=U_{\nu}$ automatically satisfies the linear
growth condition of Theorem 3. At the same time, integrating (50) w.r.t. $\nu$
and using (27) (exchanging the roles between first and second marginal) yields
to

(51) $\displaystyle\int|\nabla\psi^{0}-\nabla\psi^{\star}|\;\mathrm{d}\nu\leq
T^{-1}\int|x-y|\;\pi_{T}^{x,\varphi^{\star}}(\mathrm{d}y)\;\nu(\mathrm{d}x)=T^{-1}\int|x-y|\;\mathrm{d}\pi^{\star}$

which allows us to state (33) in terms of the first moments of the marginals
$\mu$ and $\nu$. Summarizing the above discussion, if we start from
$\psi^{0}=U_{\nu}$, the previous results read as

######  Corollary .

Assume A 1. If we set the initial Sinkhorn iterate equal to $\psi^{0}=U_{\nu}$
(or equivalently if we start with $\varphi^{0}=0$), then the linear-growth
condition (43) is satisfied with

(52)
$A=T^{-1}\,\frac{T\,\alpha_{\varphi^{\star}}+2}{T\,\alpha_{\varphi^{\star}}+1}\quad\text{
and }\quad
B=\frac{1+\left\|\tilde{g}_{\mu}\right\|_{\infty}+|\nabla\varphi^{\star}(0)|}{T\,\alpha_{\varphi^{\star}}+1}\;,$

where
$\alpha_{\varphi^{\star}}\in(\alpha_{\mu}-T^{-1},\alpha_{\mu}-T^{-1}+(\beta_{\nu}T^{2})^{-1}]$
is given in Theorem 7.

Moreover, the integrated bounds (33) read as

(53)
$\displaystyle\int|\nabla\varphi^{n}-\nabla\varphi^{\star}|(x)\mu(dx)\leq$
$\displaystyle\;\frac{1}{\gamma_{n-1}^{\mu}}\;\prod_{k=0}^{n-1}\frac{\gamma_{k}^{\mu}\;\gamma_{k}^{\nu}}{T^{2}}\;(M_{1}(\mu)+M_{1}(\nu))\;,$
$\displaystyle\int|\nabla\psi^{n}-\nabla\psi^{\star}|(y)\nu(dy)\leq$
$\displaystyle\;\prod_{k=0}^{n-1}\frac{\gamma_{k}^{\mu}\;\gamma_{k}^{\nu}}{T^{2}}\;(M_{1}(\mu)+M_{1}(\nu))\;.$

As a consequence of the previous corollary, starting from $\psi^{0}=U_{\nu}$,
we have exponential pointwise convergence of the gradients as soon as $T$ is
large enough, as mentioned in Theorem 3.

It is possible to infer convergence of Sinkhorn iterates
$(\varphi^{n})_{n\in\mathbb{N}}$ and $(\psi^{n})_{n\in\mathbb{N}}$. from the
convergence of gradients. Since Schrödinger potentials are unique up to a
trivial additive shift, we have to impose a normalization. In particular, we
consider the symmetric normalization

(54)
$\int\varphi^{\star}\,\mathrm{d}\mu+\mathrm{Ent}(\mu)=\int\psi^{\star}\,\mathrm{d}\nu+\mathrm{Ent}(\nu)=\frac{1}{2}\biggl{(}\mathrm{Ent}(\mu)+\mathrm{Ent}(\nu)-\mathscr{H}(\pi^{\star}|\mathrm{R}_{0,T})\biggr{)}\;.$

This normalization has already been used while showing convergence of rescaled
Schrödinger potentials and thier gradients to Kantorovich potentials and the
Brenier map respectively [48, 16]. In what concerns iterates, we work with the
following normalization

(55) $\varphi^{\diamond
n}=\varphi^{n}-\biggl{(}\int\varphi^{n}\mathrm{d}\mu-\int\varphi^{\star}\mathrm{d}\mu\biggr{)}\;,\qquad\psi^{\diamond
n}=\psi^{n}-\biggl{(}\int\psi^{n}\mathrm{d}\nu-\int\psi^{\star}\mathrm{d}\nu\biggr{)}\,,$

This choice guarantees that

(56) $\int\varphi^{\diamond
n}\mathrm{d}\mu+\mathrm{Ent}(\mu)=\int\varphi^{\star}\mathrm{d}\mu+\mathrm{Ent}(\mu)=\int\psi^{\diamond
n}\mathrm{d}\nu+\mathrm{Ent}(\nu)=\int\psi^{\star}\mathrm{d}\nu+\mathrm{Ent}(\nu)\;.$

######  Theorem 4.

Assume that A 1 and (43) hold. Then for any $n\geq 1$ and $x\in\mathbb{R}^{d}$
it holds

(57) $\displaystyle|\varphi^{\diamond n}-\varphi^{\star}|(x)\leq$
$\displaystyle\;\frac{T}{\hat{\gamma}^{\mu}_{n-1}}\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\;\hat{\gamma}_{k}^{\nu}}{T^{2}}\;\biggl{[}A\;|x|^{2}+(A\;M_{1}(\mu)+B)\;|x|+B\;M_{1}(\mu)+2\,A\;M_{2}(\mu)\biggr{]}\;,$
$\displaystyle|\psi^{\diamond n}-\psi^{\star}|(x)\leq$
$\displaystyle\;\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\;\hat{\gamma}_{k}^{\nu}}{T^{2}}\;\biggl{[}A\;|x|^{2}+(A\;M_{1}(\nu)+B)\;|x|+B\;M_{1}(\nu)+2\,A\;M_{2}(\nu)\biggr{]}\;.$

Hence, for $T$ large enough (e.g., (47)), pointwise exponential convergence of
Sinkhorn’s iterates holds. Finally, if the initial iteration is set equal to
$\psi^{0}=U_{\nu}$ (i.e., $\varphi^{0}=0$), then the above bounds hold true
with $A$ and $B$ given at (52).

Let us mention here that a straightforward adaptation of the proof of Theorem
4 implies that this result also holds true under a pointwise normalization
(e.g., $\psi^{\star}(0)=\psi^{\diamond n}(0)=U_{\nu}(0)$) or for the symmetric
zero-mean option considered in [24, 11, 23, 10].

We now move on to consider exponential entropic convergence of Sinkhorn’s
algorithm on the primal side, i.e., for Sinkhorn plans
$(\pi^{n,n})_{n\in\mathbb{N}}$ and $(\pi^{n+1,n})_{n\in\mathbb{N}}$ defined in
(7) and for the adjusted marginals $(\mu^{n})_{n\in\mathbb{N}}$ and
$(\nu^{n})_{n\in\mathbb{N}}$, using a symmetrized version of the relative
entropy. As observed in [16, 33], measuring distances between plans with this
divergence leads to tractable expressions.

######  Theorem 5 (Exponential convergence of Sinkhorn on the primal side).

Assume that A 1 and (43) hold. Then, for any $n\geq 1$ it holds

(58)
$\begin{split}\mathscr{H}(\pi^{n,n}|\pi^{\star})+\mathscr{H}(\pi^{\star}|\pi^{n,n})\leq&\;D(A,B,\mu)\frac{T}{\hat{\gamma}^{\mu}_{n-1}}\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\hat{\gamma}_{k}^{\nu}}{T^{2}}\;,\\\
\mathscr{H}(\pi^{n+1,n}|\pi^{\star})+\mathscr{H}(\pi^{\star}|\pi^{n+1,n})\leq&\;D(A,B,\nu)\;\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\hat{\gamma}_{k}^{\nu}}{T^{2}}\;,\end{split}$

and

(59)
$\begin{split}\mathscr{H}(\pi^{n,n}|\pi^{n,n-1})+\mathscr{H}(\pi^{n,n-1}|\pi^{n,n})=\mathscr{H}(\mu^{n}|\mu)+\mathscr{H}(\mu|\mu^{n})\leq&\;D(A,B,\mu)\frac{T}{\hat{\gamma}^{\mu}_{n-1}}\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\hat{\gamma}_{k}^{\nu}}{T^{2}}\;,\\\
\mathscr{H}(\pi^{n+1,n}|\pi^{n,n})+\mathscr{H}(\pi^{n,n}|\pi^{n+1,n})=\mathscr{H}(\nu^{n}|\nu)+\mathscr{H}(\nu|\nu^{n})\leq&\;D(A,B,\nu)\;\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\hat{\gamma}_{k}^{\nu}}{T^{2}}\;,\end{split}$

where the multiplicative constants $D(A,B,\mu),\,D(A,B,\nu)$ are explicitly
given in (187) and (191).

As a consequence, for $T$ large enough (e.g. (47)), entropic exponential
convergence of Sinkhorn’s plans and adjusted marginals holds. Finally, if the
initial iteration is set equal to $\psi^{0}=U_{\nu}$ (i.e., $\varphi^{0}=0$),
then the above bounds hold true with $A$ and $B$ given at (52).

In Section 5.2 we show how the multiplicative constants $D(A,B,\cdot)$ can be
further improved. The benefit of considering symmetric relative entropies in
Theorem 5 is twofold: not only it allows us to bound these relative entropies
in terms of $\varphi^{n}-\varphi^{\star}$ and $\psi^{n}-\psi^{\star}$, but
also allows us to translate it in terms of $\varphi^{\diamond
n}-\varphi^{\star}$ and $\psi^{\diamond n}-\psi^{\star}$ and therefore apply
the results of Theorem 4.We conclude the presentation of our main results by
showing that the pointwise convergence of
$(\nabla\varphi^{n})_{n\in\mathbb{N}}$ and $(\nabla\psi^{n})_{n\in\mathbb{N}}$
implies the pointwise convergence of the corresponding Hessian matrices
$(\nabla^{2}\varphi^{n})_{n\in\mathbb{N}}$ and
$(\nabla^{2}\psi^{n})_{n\in\mathbb{N}}$ , with the exact same exponential
convergence rate. We will measure this convergence through the Frobenius norm,
i.e., for any $\mathbf{A}\in\mathbb{R}^{d\times d}$

(60)
$\|\mathbf{A}\|_{\mathrm{F}}\coloneqq\sqrt{\operatorname{Tr}(\mathbf{A}\mathbf{A}^{{\sf
T}})}\;.$

The proof of convergence for the Hessians follows a similar scheme than the
one for gradients replacing the representation formula for the derivative of
Sinkhorn potentials with its second order counterpart that follows from the
general relation

(61) $\nabla^{2}\log
P_{T}\exp(-h)(x)=-T^{-1}\operatorname{Id}+T^{-2}\,\mathrm{Cov}(\pi_{T}^{x,h})\;,$

The above is proven as a part of Section 2.2.

######  Theorem 6.

Assume that A 1 and (43) hold. Then, for any $n\geq 1$ it holds

(62)
$\begin{split}\|\nabla^{2}\varphi^{n+1}-\nabla^{2}\varphi^{\star}\|_{\mathrm{F}}(x)\leq
C(x,T,\nu,A,B)\,\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\hat{\gamma}_{k}^{\nu}}{T^{2}}\;,\\\
\|\nabla^{2}\psi^{n+1}-\nabla^{2}\psi^{\star}\|_{\mathrm{F}}(x)\leq
C(x,T,\mu,A,B)\frac{T}{\hat{\gamma}^{\mu}_{n}}\prod_{k=0}^{n}\frac{\hat{\gamma}_{k}^{\mu}\hat{\gamma}_{k}^{\nu}}{T^{2}}\;,\end{split}$

where $C(x,T,\nu,A,B)$ and $C(x,T,\mu,A,B)$ are given in (228). As a
consequence, for $T$ large enough (e.g., (47)), pointwise exponential
convergence of Sinkhorn’s Hessians holds in Frobenius norm. Finally, if the
initial potential is set equal to $\psi^{0}=U_{\nu}$ (i.e., $\varphi^{0}=0$),
then the above bounds hold true with $A$ and $B$ given at (52).

###  2.3. Literature review

As already emphasized in the introduction, when the marginal distributions
$\mu,\nu$ are discrete probability measures, convergence of Sinkhorn’s
algorithm is well-known [30, 9]. In addition, it has also been established
when dealing with compact state spaces or bounded cost functions [12, 23, 7].

When considering non compact spaces with possibly unbounded costs and
marginals, the techniques developed to handle the discrete setting cannot
apply as such and new ideas have emerged. When it comes to results that allow
for unbounded costs, [57] shows qualitative convergences of iterates in
relative entropy and total variation for Sinkhorn plans. More recently, [49]
establishes qualitative convergence both on the primal and dual sides under
mild assumptions. Concerning convergence rates, [41] gives an interpretation
of Sinkhorn’s algorithm as a block coordinate descent on the dual problem and
obtains convergence of marginal distributions in relative entropy at a linear
rate $n^{-1}$ under minimal assumptions. [28] derives polynomial rates of
convergence in Wasserstein distance assuming, among other things, that
marginals admit exponential moments. Lastly, [33] improves existent polynomial
convergence rates for optimal plans with respect to a symmetric relative
entropy.

In contrast to these existing works, the main contribution of this paper is to
establish exponential convergence bounds for the gradients and Hessians of
Sinkhorn iterates as well as for the optimal plans. To the best of our
knowledge these findings are new both in their dual and primal formulation in
that they represent the first exponential convergence results that holds for
unbounded costs and marginals. They also are among the very few results that
yield convergence of derivatives of potentials. On this subject, let us
mention the recent work by the authors of the present paper and collaborators
[36] where the state space is the $d$-dimensional torus and therefore deals
with bounded cost functions. As [36], our proofs are mainly probabilistic, and
differ from other proposed methodologies in that they rely on one-sided
integrated semiconvexity estimates for potentials along Sinkhorn iterates.
These estimates are by themselves a new result, that has potentially several
further implications. Though the current approach and the one we proposed in
[36] are both inspired by coupling methods and stochastic control, there is a
fundamental difference. In [36] exponential convergence is achieved through
Lipschitz estimates on potentials. In the current setup, we make assumptions
on the integrated log-concavity profile of the marginals; these assumption are
of geometric nature and not perturbative.

We conclude this section reporting on results available for multimarginal
entropic optimal transport problem and the natural extension of Sinkhorn
algorthm in this setting. For bounded costs and marginals, (or equivalently
compact spaces) [11] shows well-posedness of the Schrödinger system and smooth
dependence of Schrödinger potentials on the marginal inputs. [24] managed to
show qualitative convergence of Sinkhorn iterates to Schrödinger potentials in
$\mathrm{L}^{p}$-norms using tools from calculus of variations: their results
require bounded costs but apply to multimarginal problems. These results have
been subsequently improved by Carlier [10] who establishes exponential
convergence.

##  3\. Integrated semiconvexity propagation along Sinkhorn

In this section we establish lower bounds on the integrated convexity profile
of Sinkhorn potentials. Before proceeding further, let us point out here that
the results of this section hold under a weaker assumption than A 1.

######  A’ 1.

The two distributions $\mu,\nu\in\mathcal{P}(\mathbb{R}^{d})$ specified by (2)
have finite relative entropy with respect to the Lebesgue measure
$\mathrm{Leb}$, $\mathrm{Ent}(\mu),\,\mathrm{Ent}(\nu)<+\infty$;

######  A’ 2.

  1. (i)

There exist $\alpha_{\nu}\in(0,+\infty)$ and $\beta_{\mu}\in(0,+\infty]$ such
that

(63) $\kappa_{U_{\nu}}(r)\geq\alpha_{\nu}-r^{-1}\,\hat{g}_{\nu}(r)\quad\text{
and }\quad\ell_{U_{\mu}}(r)\leq\beta_{\mu}+r^{-1}\,g_{\mu}(r)\;,$

with $g_{\mu}\in\mathcal{G}$ and $\hat{g}_{\nu}\in\hat{\mathcal{G}}$;

  2. (ii)

There exist $\alpha_{\mu}\in(0,+\infty)$ and $\beta_{\nu}\in(0,+\infty]$ such
that

(64) $\kappa_{U_{\mu}}(r)\geq\alpha_{\mu}-r^{-1}\,\hat{g}_{\mu}(r)\quad\text{
and }\quad\ell_{U_{\nu}}(r)\leq\beta_{\nu}+r^{-1}\,g_{\nu}(r)\;,$

with $g_{\nu}\in\mathcal{G}$ and $\hat{g}_{\mu}\in\hat{\mathcal{G}}$,

where

(65) $\hat{\mathcal{G}}\coloneqq\biggl{\\{}g\in\mathcal{G}\quad\text{ bounded
and s.t.}\quad 2g{{}^{\prime\prime}}+g\,g{{}^{\prime}}\leq
0\biggr{\\}}\subseteq\mathcal{G}\;.$

Notice that the above assumption allows for concave functions having non-null
limit value in the origin (whereas the elements of $\tilde{\mathcal{G}}$ are
sublinear in a neighborhood of the origin), that
$\tilde{\mathcal{G}}\subseteq\hat{\mathcal{G}}$ and hence that A 1 implies A’
1 and A’ 2.

Let us introduce for any fixed $\beta>0$ and any $g\in\mathcal{G}$ and
$\hat{g}\in\hat{\mathcal{G}}$, the following functions for $\alpha,\,s,\,u\geq
0$

(66) $\displaystyle F_{\beta}^{g,\hat{g}}(\alpha,s)=$
$\displaystyle\,\beta\,s+\frac{s}{T(1+T\alpha)}+s^{1/2}\,g(s^{1/2})+\frac{s^{1/2}\,\hat{g}(s^{1/2})}{(1+T\alpha)^{2}}\;,$
$\displaystyle G_{\beta}^{g,\hat{g}}(\alpha,u)=$ $\displaystyle\,\inf\\{s\geq
0\;\colon\;F_{\beta}^{g,\hat{g}}(\alpha,s)\geq u\\}\;,$

with the convention $G_{\beta}^{g,\hat{g}}(\alpha,u)\equiv 0$ whenever
$\beta=+\infty$.

Then, the main result of this section can be stated as follows.

######  Theorem 7.

Assume A’ 1. If A’ 2-$(i)$ holds and if

(67) $\kappa_{\psi^{0}}(r)\geq\alpha_{\nu}-T^{-1}-r^{-1}\hat{g}_{\nu}(r)\;,$

then there exists a monotone increasing sequence
$(\alpha_{\nu,n})_{n\in\mathbb{N}}\subseteq(\alpha_{\nu}-T^{-1},\alpha_{\nu}-T^{-1}+(\beta_{\mu}\,T^{2})^{-1}]$
such that for any $n\geq 1$ and $r>0$ it holds

(68) $\ell_{\varphi^{n}}(r)\leq
r^{-2}\,F_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}(\alpha_{\nu,n},r^{2})-T^{-1}\quad\text{
and }\quad\kappa_{\psi^{n}}(r)\geq\alpha_{\nu,n}-r^{-1}\hat{g}_{\nu}(r)\;,$

with $F_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}$ defined as in (66). Moreover,
the sequence can be explicitly built by setting

(69) $\begin{cases}\alpha_{\nu,0}\coloneqq\alpha_{\nu}-T^{-1}\;,\\\
\alpha_{\nu,n+1}\coloneqq\alpha_{\nu}-T^{-1}+\frac{G_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}(\alpha_{\nu,n},2)}{2\;T^{2}}\;,\quad
n\in\mathbb{N}\;,\end{cases}$

$G_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}$ given in (66). Finally,
$(\alpha_{\nu,n})_{n\in\mathbb{N}}$ converges to
$\alpha_{\psi^{\star}}\in(\alpha_{\nu}-T^{-1},\alpha_{\nu}-T^{-1}+(\beta_{\mu}\,T^{2})^{-1}]$,
fixed point solutions of (69) and for any $r>0$,

(70) $\ell_{\varphi^{\star}}(r)\leq
r^{-2}\,F_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}(\alpha_{\psi^{\star}},r^{2})-T^{-1}\quad\text{
and
}\quad\kappa_{\psi^{\star}}(r)>\alpha_{\psi^{\star}}-r^{-1}\hat{g}_{\nu}(r)\;,$

where $\varphi^{\star}$ and $\psi^{\star}$ are the Schrödinger potentials
introduced in (1).

Similarly, under A’ 2-$(ii)$ there exists a monotone increasing sequence
$(\alpha_{\mu,n})_{n\in\mathbb{N}}\subseteq(\alpha_{\mu}-T^{-1},\alpha_{\mu}-T^{-1}+(\beta_{\nu}\,T^{2})^{-1}]$
such that for any $n\geq 1$ and $r>0$ it holds

(71) $\ell_{\psi^{n}}(r)\leq
r^{-2}\,F_{\beta_{\nu}}^{g_{\nu},\hat{g}_{\mu}}(\alpha_{\mu,n},r^{2})-T^{-1}\quad\text{
and }\quad\kappa_{\varphi^{n}}(r)\geq\alpha_{\mu,n}-r^{-1}\hat{g}_{\mu}(r)\;,$

with $F_{\beta_{\nu}}^{g_{\nu},\hat{g}_{\mu}}$ defined as in (66) and

(72) $\begin{cases}\alpha_{\mu,1}\coloneqq\alpha_{\mu}-T^{-1}\;,\\\
\alpha_{\mu,n+1}\coloneqq\alpha_{\mu}-T^{-1}+\frac{G_{\beta_{\nu}}^{g_{\nu},\hat{g}_{\mu}}(\alpha_{\mu,n},2)}{2\;T^{2}}\;,\quad
n\in\mathbb{N}\;,\end{cases}$

with $G_{\beta_{\nu}}^{g_{\nu},\hat{g}_{\mu}}$ defined as in (66). Finally,
$(\alpha_{\mu,n})_{n\in\mathbb{N}}$ converges to
$\alpha_{\varphi^{\star}}\in(\alpha_{\mu}-T^{-1},\alpha_{\mu}-T^{-1}+(\beta_{\nu}\,T^{2})^{-1}]$,
fixed point solutions of (72) and for any $r>0$,

(73) $\ell_{\psi^{\star}}(r)\leq
r^{-2}\,F_{\beta_{\nu}}^{g_{\nu},\hat{g}_{\mu}}(\alpha_{\varphi^{\star}},r^{2})-T^{-1}\quad\text{
and
}\quad\kappa_{\varphi^{\star}}(r)>\alpha_{\varphi^{\star}}-r^{-1}\hat{g}_{\mu}(r)\;.$

######  Remark .

The above result is an extension of [18, Theorem 1.2] where the author just
provides the limit-bounds (70) and (73) in the case when
$g^{\mu}=g^{\nu}\equiv 0$ and $\hat{g}_{\mu},\tilde{g}_{\nu}$ take the form
(22). In the above result we show that the iterative proof given there can be
actually employed when proving the estimates (68) and (71) along Sinkhorn’s
algorithm.

Let us also mention that our result encompasses [15, Theorem 4] when
considering $\hat{g}_{\nu}\equiv 0$ and $g_{\mu}\equiv 0$ in A’ 2-$(i)$.

We provide a proof of the above theorem at the end of this section. Let us
mention here that we will show in Section 3 that Assumption (67) on $\psi^{0}$
can be essentially dropped; here we simply observe that it is met for a
regular enough initial condition, e.g., for $\psi^{0}=U_{\nu}$.

Let us also point out that the above theorem guarantees existence and
uniqueness for the strong solution of the SDE

(74) $\mathrm{d}Y_{t}=-\biggl{(}\frac{Y_{t}-x}{2T}+\frac{1}{2}\nabla
h(Y_{t})\biggr{)}\mathrm{d}t+\mathrm{d}B_{t}\,.$

for any choice of $h=\psi^{\star},\,\psi^{n},\,\varphi^{\star},\,\varphi^{n}$.
Indeed from Theorem 7 we immediately deduce

######  Corollary .

Under the assumptions of Theorem 7 it holds for any $y\in\mathbb{R}^{d}$

(75)
$\begin{split}\left\langle\nabla\psi^{\star}(y),\,\frac{y}{|y|}\right\rangle\geq\alpha_{\psi^{\star}}\,|y|-\hat{g}_{\nu}(|y|)-|\nabla\psi^{\star}(0)|\;,\quad\left\langle\nabla\psi^{n}(y),\,\frac{y}{|y|}\right\rangle\geq\alpha_{\nu,n}\,|y|-\hat{g}_{\nu}(|y|)-|\nabla\psi^{n}(0)|\;,\\\
\left\langle\nabla\varphi^{\star}(y),\,\frac{y}{|y|}\right\rangle\geq\alpha_{\varphi^{\star}}\,|y|-\hat{g}_{\mu}(|y|)-|\nabla\varphi^{\star}(0)|\;,\quad\left\langle\nabla\varphi^{n}(y),\,\frac{y}{|y|}\right\rangle\geq\alpha_{\mu,n}\,|y|-\hat{g}_{\mu}(|y|)-|\nabla\varphi^{n}(0)|\;.\end{split}$

As a consequence for any even $p\geq 2$ the potential $V_{p}(y)=1+|y|^{p}$ is
a Lyapunov function for (74) with
$h\in\\{\psi^{\star},\,\psi^{n},\,\varphi^{\star},\,\varphi^{n}\\}$. As a
consequence, existence and uniqueness of strong solutions holds for these
SDEs.

###### Proof.

The lower-bounds displayed above are a direct consequence of Theorem 7. Let
$x\in\mathbb{R}^{d}$ be fixed. We will only consider now the case
$h=\psi^{\star}$. The other cases follow the same lines. Since
$\hat{g}_{\nu}\in\hat{\mathcal{G}}$ is bounded, it holds for any
$y\in\mathbb{R}^{d}$

(76) $\begin{split}-\frac{1}{2}\,\left\langle
T^{-1}(y-x)+\nabla\psi^{\star}(y),\,y\right\rangle\leq-\frac{\alpha_{\psi^{\star}}+T^{-1}}{2}\,|y|^{2}+\frac{\|\hat{g}_{\nu}\|_{\infty}+|\nabla
h(0)|+T^{-1}|x|}{2}|y|\;,\end{split}$

and hence there exist $\gamma>0$ and $R>0$ such that

(77) $-\frac{1}{2}\,\left\langle
T^{-1}(y-x)+\nabla\psi^{\star}(y),\,y\right\rangle\leq-\gamma\,|y|^{2}\quad\forall\,|y|\geq
R\;.$

At this stage, [43, Lemma 4.2] guarantees that for any even $p\geq 2$ the
potential $V_{p}(y)=1+|y|^{p}$ is a Lyapunov function for the diffusion (74).

More precisely it holds a geometric drift condition, i.e., for any
$A_{\psi^{\star}}\in(0,p\gamma)$ there exists a finite constant
$B_{\psi^{\star}}=B_{\psi^{\star}}(A_{\psi^{\star}},p)$ such that for any
$y\in\mathbb{R}^{d}$

(78) $\mathcal{L}_{\psi^{\star}}V_{p}(y)\leq-
A_{\psi^{\star}}V_{p}(y)+B_{\psi^{\star}}\;,$

where above $\mathcal{L}_{\psi^{\star}}\coloneqq\Delta/2-\frac{1}{2}\langle
T^{-1}(y-x)+\nabla\psi^{\star}(y),\,\nabla\rangle$ denotes the generator
associated to the SDE (74). Finally, existence and uniqueness of strong
solutions of (74) follows from [56, Theorem 2.1] (see also [44, Section 2]). ∎

The proof of Theorem 7 will be based on a propagation of integrated-convexity
along Hamilton-Jacobi-Bellman (HJB) equations observed in [18], based on
coupling by reflection techniques, which reads as follows

######  Theorem 8 (Theorem 2.1 in [18]).

For any fixed function $\hat{g}\in\hat{\mathcal{G}}$, consider the class of
functions

(79)
$\mathcal{F}_{\hat{g}}\coloneqq\\{h\in\mathrm{C}^{1}(\mathbb{R}^{d})\;\colon\;\kappa_{h}(r)\geq-r^{-1}\,\hat{g}(r)\quad\forall\,r>0\\}\;.$

Then, the class $\mathcal{F}_{\hat{g}}$ is stable under the action of the HJB
flow, i.e.,

(80) $h\in\mathcal{F}_{\hat{g}}\;\Rightarrow\;-\log
P_{T-t}\exp(-h)\in\mathcal{F}_{\hat{g}}\quad\forall\;0\leq t\leq T\;.$

We omit the proof of the above result since it runs exactly as stated in [18].
There it is proven when $\hat{g}$ is of the form (22). However same proofs
allows to reach the conclusion for any function $\hat{g}\in\hat{\mathcal{G}}$
since it only requires $\hat{g}$ to satisfy the differential inequality

(81) $2(\hat{g}){{}^{\prime\prime}}+\hat{g}\,(\hat{g}){{}^{\prime}}\leq 0\;,$

which is an equality in the special case considered there
$\hat{g}(r)=\tanh(r)$.

As a first consequence of the previous theorem we may immediately deduce the
following integrated propagation _convexity-to-concavity_ result.

######  Lemma .

Assume A’ 1. If A’ 2-$(i)$ holds true, $\alpha_{\nu,n}>-T^{-1}$ and if for any
$r>0$

(82) $\kappa_{\psi^{n}}(r)\geq\alpha_{\nu,n}-r^{-1}\,\hat{g}_{\nu}(r)\;,$

then

(83)
$\ell_{\varphi^{n+1}}(r)\leq\beta_{\mu}+g_{\mu}(r)-\frac{\alpha_{\nu,n}}{1+T\alpha_{\nu,n}}+\frac{r^{-1}\,\hat{g}_{\nu}(r)}{(1+T\alpha_{\nu,n})^{2}}=-T^{-1}+r^{-2}\,F_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}(\alpha_{\nu,n},r^{2})\;.$

Similarly if A’ 2-$(ii)$ holds, $\alpha_{\mu,n}>-T^{-1}$ and if for any $r>0$

(84) $\kappa_{\varphi^{n}}(r)\geq\alpha_{\mu,n}-r^{-1}\,\hat{g}_{\mu}(r)\;,$

then

(85)
$\ell_{\psi^{n}}(r)\leq\beta_{\nu}+r^{-1}\,g_{\nu}(r)-\frac{\alpha_{\mu,n}}{1+T\alpha_{\mu,n}}+\frac{r^{-1}\,\hat{g}_{\mu}(r)}{(1+T\alpha_{\mu,n})^{2}}=-T^{-1}+r^{-2}\,F_{\beta_{\nu}}^{g_{\nu},\hat{g}_{\mu}}(\alpha_{\mu,n},r^{2})\;.$

###### Proof.

Let us firstly notice that our assumption (82) is equivalent to stating that

(86)
$\bar{\psi}^{n}\coloneqq\psi^{n}-\frac{\alpha_{\nu,n}}{2}\,|\cdot|^{2}\,\in\mathcal{F}_{\hat{g}_{\nu}}\;,$

and therefore Theorem 8 implies that

(87) $-\log P_{T}\exp(-\bar{\psi}^{n})\in\mathcal{F}_{\hat{g}_{\nu}}\;.$

By recalling that $\varphi^{n+1}$ is defined via (4), in order to conclude it
is enough noticing now that

(88) $\begin{split}-\log P_{T}\exp(-&\psi^{n})(x)-\frac{d}{2}\log(2\pi
T)=-\log\int\exp\biggl{(}-\frac{|x-y|^{2}}{2T}-\frac{\alpha_{\nu,n}}{2}|y|^{2}-\bar{\psi}^{n}(y)\biggr{)}\;\mathrm{d}y\\\
=&\frac{\alpha_{\nu,n}\,|x|^{2}}{2(1+T\alpha_{\nu,n})}-\log\int\exp\biggl{(}-\frac{1+T\alpha_{\nu,n}}{2T}\,|y-(1+T\alpha_{\nu,n})^{-1}x|^{2}-\bar{\psi}^{n}(y)\biggr{)}\;\mathrm{d}y\\\
=&\frac{\alpha_{\nu,n}\,|x|^{2}}{2(1+T\alpha_{\nu,n})}-\log
P_{T/(1+T\alpha_{\nu,n})}\exp(-\bar{\psi}^{n})((1+T\alpha_{\nu,n})^{-1}x)-\frac{d}{2}\log\frac{2\pi
T}{1+T\alpha_{\nu,n}}\;,\end{split}$

and combining it with (87) and A’ 2-$(i)$.

The second part of the statement follows the same lines. ∎

######  Lemma .

Fix $\beta\in(0,+\infty]$ and two functions $g\in\mathcal{G}$ and
$\hat{g}\in\hat{\mathcal{G}}$. Then the following properties hold true.

  1. (1)

For any $\alpha>-T^{-1}$ the function $s\mapsto
F_{\beta}^{g,\hat{g}}(\alpha,s)$ is concave and increasing on $[0,+\infty)$.

  2. (2)

The function $\alpha\mapsto G_{\beta}^{g,\hat{g}}(\alpha,2)$ is positive and
non-decreasing over $(-T^{-1},+\infty)$.

  3. (3)

For any given $a_{0}>0$, the fixed-point problem

(89) $\alpha=a_{0}-T^{-1}+\frac{G_{\beta}^{g,\hat{g}}(\alpha,2)}{2\;T^{2}}$

admits at least one solution on
$(a_{0}-T^{-1},\,a_{0}-T^{-1}+(\beta\,T^{2})^{-1}]$ and, as soon as
$\beta<+\infty$, $a_{0}-T^{-1}$ does not belong to the closure of the set of
solutions of (89).

###### Proof.

  1. (1)

Since $r\mapsto r\,\hat{g}(r)$ and $r\mapsto r\,g(r)$ are non-decreasing and
$\alpha>-T^{-1}$, an explicit computation shows that $s\mapsto
F_{\beta}^{g,\hat{g}}(\alpha,s)$ is an increasing function. The concavity of
the latter function follows from the properties of
$g,\,\hat{g}\in\mathcal{G}$, since for $h=g,\,\hat{g}$ it holds

(90)
$\left.\frac{\mathrm{d}^{2}}{\mathrm{d}u^{2}}\biggl{(}u^{1/2}\,h(u^{1/2})\biggr{)}\right|_{u=s}=\frac{s^{-1/2}}{4}\;\biggl{[}h{{}^{\prime\prime}}(s^{1/2})+s^{-1/2}\,h{{}^{\prime}}(s^{1/2})-s^{-1}\,h(s^{1/2})\biggr{]}\leq
0\;.$

  2. (2)

The proof is by contradiction. Notice that $G_{\beta}^{g,\hat{g}}(\cdot,2)$ is
continuous on $(-T^{-1},+\infty)$ and assume that is not a positive function,
which implies that there exists some $\alpha>-T^{-1}$ such that
$G_{\beta}^{g,\hat{g}}(\alpha,2)=0$ and hence by definition that there exists
a sequence $(s_{n})_{n\in\mathbb{N}}$ converging to zero and such that
$F_{\beta}^{g,\hat{g}}(\alpha,s)\geq 2$, which is clearly impossible since
$\lim_{s\downarrow 0}F_{\beta}^{g,\hat{g}}(\alpha,s)=0$. Hence
$G_{\beta}^{g,\hat{g}}(\cdot,2)$ is a positive function. The monotonicity of
$G^{g_{\mu},\hat{g}_{\nu}}_{\beta_{\mu}}(\cdot,2)$ follows from the fact that
$F_{\beta}^{g,\hat{g}}(\alpha,s)$ is increasing in $s$ and decreasing in
$\alpha\in(-T^{-1},+\infty)$, which implies for any
$\alpha^{\prime}\geq\alpha$ and $u\geq 0$

(91) $\\{s\;\colon\;F_{\beta}^{g,\hat{g}}(\alpha^{\prime},s)\geq
u\\}\subseteq\\{s\;\colon\;F_{\beta}^{g,\hat{g}}(\alpha,s)\geq u\\}\;.$

  3. (3)

Consider the map associated to the fixed-point problem (69), i.e., the
continuous function $H\colon[a_{0}-T^{-1},+\infty)\to\mathbb{R}$ defined for
$\alpha\in(-T^{-1},+\infty)$ as

(92) $H(\alpha)\coloneqq\alpha-
a_{0}+T^{-1}-\frac{G_{\beta}^{g,\hat{g}}(\alpha,2)}{2\;T^{2}}\;.$

Let us now prove that

(93) $H(a_{0}-T^{-1})<0\quad\text{ and
}\quad\lim_{\alpha\to+\infty}H(\alpha)=+\infty\;.$

The first inequality follows from a direct computation, showing that
$G_{\beta}^{g,\hat{g}}(a_{0}-T^{-1},2)>0$. In order to prove the second
statement it is enough noticing that $G_{\beta}^{g,\hat{g}}(\cdot,2)$ is
bounded, which is immediate since $g,\,\hat{g}\geq 0$ implies for any
$\alpha>-T^{-1}$ and $s>0$ that $F_{\beta}^{g,\hat{g}}(\alpha,s)\geq\beta\,s$
and hence that

(94) $G_{\beta}^{g,\hat{g}}(\alpha,2)\leq 2/\beta\;.$

From (93) and the continuity of $G_{\beta}^{g,\hat{g}}(\cdot,2)$ we finally
deduce the existence of some $\bar{\alpha}\in(a_{0}-T^{-1},+\infty)$ such that
$H(\bar{\alpha})=0$, i.e., a fixed point for (89). As a consequence (94)
further implies $\bar{\alpha}\leq a_{0}-T^{-1}+(\beta\,T^{2})^{-1}$. Finally
$a_{0}-T^{-1}$ does not belong to the closure of the set of fixed-points
solutions, because if this was the case then the continuity of
$G_{\beta}^{g,\hat{g}}(\cdot,2)$ would have implied $H(a_{0}-T^{-1})=0$,
clearly in contrast with (93).

∎

As a corollary of the previous lemma we have already proven the following

######  Corollary .

Assume A’ 1. If A’ 2-$(i)$ holds true, then there exists at least one solution
$\alpha^{\star}_{\nu}$ on
$(\alpha_{\nu}-T^{-1},\alpha_{\nu}-T^{-1}+(\beta_{\mu}\,T^{2})^{-1}]$ to the
fixed point associated to (69). Moreover, if $\beta_{\mu}$ is finite then
$\alpha_{\nu}-T^{-1}$ is not an accumulation point for the set of solutions.

Similarly if A’ 2-$(ii)$ holds true, there exists at least one solution
$\alpha^{\star}_{\mu}$ on
$(\alpha_{\mu}-T^{-1},\alpha_{\mu}-T^{-1}+(\beta_{\nu}\,T^{2})^{-1}]$ to the
fixed point associated to (72). Moreover, if $\beta_{\nu}$ is finite then
$\alpha_{\mu}-T^{-1}$ is not an accumulation point for the set of solutions.

The next result is the counterpart to Section 3 and it shows that we do also
have an integrated propagation _concavity-to-convexity_.

######  Lemma .

Assume A’ 1. If A’ 2-$(i)$ holds, $\alpha_{\nu,n}>-T^{-1}$ and if

(95)
$\ell_{\varphi^{n+1}}(r)\leq-T^{-1}+r^{-2}\,F_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}(\alpha_{\nu,n},r^{2})\;,$

then $\alpha_{\nu,n+1}>-T^{-1}$ and for any $r>0$

(96) $\kappa_{\psi^{n+1}}(r)\geq\alpha_{\nu,n+1}-r^{-1}\,\hat{g}_{\nu}(r)\;.$

Similarly if A’ 2-$(ii)$ holds, $\alpha_{\mu,n}>-T^{-1}$ and if

(97)
$\ell_{\psi^{n}}(r)\leq-T^{-1}+r^{-2}\,F_{\beta_{\nu}}^{g_{\nu},\hat{g}_{\mu}}(\alpha_{\mu,n},r^{2})\;,$

then $\alpha_{\mu,n+1}>-T^{-1}$ and for any $r>0$

(98)
$\kappa_{\varphi^{n+1}}(r)\geq\alpha_{\mu,n+1}-r^{-1}\,\hat{g}_{\mu}(r)\;.$

###### Proof.

We are going to prove only the first part of the lemma since the second can be
proven by an analogous reasoning. Firstly, let us consider the function

(99) $\hat{\psi}^{n+1}(y)\coloneqq
T\psi^{n+1}(y)-TU_{\nu}(y)+\frac{|y|^{2}}{2}\;.$

Then (61) implies that its Hessian is given for any $y\in\mathbb{R}^{d}$ by

(100)
$\nabla^{2}\hat{\psi}^{n+1}(y)=\frac{1}{T}\,\mathrm{Cov}(\pi_{T}^{y,\varphi^{n+1}})\;,$

where we recall from the definition of $\pi^{y,\varphi^{n+1}}_{T}$ (23) that

(101)
$\pi_{T}^{y,\varphi^{n+1}}(\mathrm{d}x)\propto\exp\biggl{(}-\frac{|x-y|^{2}}{2T}-\varphi^{n+1}(x)\biggr{)}\,\mathrm{d}x\;.$

Moreover, if we set for notations’ sake
$V^{y,n+1}\coloneqq-\log(\mathrm{d}\pi_{T}^{y,\varphi^{n+1}}/\mathrm{d}\mathrm{Leb})$,
then our assumption implies

(102) $\ell_{V^{y,n+1}}(r)\leq
r^{-2}\,F_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}(\alpha_{\nu,n},r^{2})\quad\forall\,r>0\;.$

In order to prove the desired bound for $\kappa_{\psi^{n+1}}$, we will first
establish a lower bound for the Hessian $\nabla^{2}\hat{\psi}^{n+1}$, i.e., a
lower bound for the covariance matrix (100). In view of that, let us consider
for any fixed $y\in\mathbb{R}^{d}$, the variance
$\mathrm{Var}_{X\sim\pi_{T}^{y,\varphi^{n+1}}}(X_{1})$ where $X_{i}$ denotes
the $i^{th}$ scalar component of the random vector
$X\sim\pi_{T}^{y,\varphi^{n+1}}$. Next, observe that

(103)
$\mathrm{Var}_{X\sim\pi_{T}^{y,\varphi^{n+1}}}(X_{1})\geq\mathbb{E}_{X\sim\pi_{T}^{y,\varphi^{n+1}}}[\mathrm{Var}_{X\sim\pi_{T}^{y,\varphi^{n+1}}}(X_{1}|X_{2},\,\dots,\,X_{d})]\;,$

and notice that for any given $z=(z_{2},\,\dots,\,z_{d})\in\mathbb{R}^{d-1}$
it holds

(104)
$\mathrm{Var}_{X\sim\pi_{T}^{y,\varphi^{n+1}}}(X_{1}|X_{2}=z_{d},\,\dots,\,X_{d}=z_{d})=\frac{1}{2}\,\int_{\mathbb{R}^{2}}|x-\hat{x}|^{2}\;\pi_{T}^{y,\varphi^{n+1}}(\mathrm{d}x|z)\,\pi_{T}^{y,\varphi^{n+1}}(\mathrm{d}\hat{x}|z)$

where $(y,z,A)\mapsto\pi_{T}^{y,h}(A|z)$ is the Markov kernel on
$\mathbb{R}^{d}\times\mathbb{R}^{d-1}\times\mathcal{B}(\mathbb{R})$ whose
transition density w.r.t. Lebesgue measure is proportional to
$\exp(-V^{y,n+1}(x,z))$. If we set $V^{y,z}(x)\coloneqq V^{y,n+1}(x,z)$ we
then have, uniformly in $z\in\mathbb{R}^{d-1}$,

(105)
$\begin{split}1=&\frac{1}{2}\int(\partial_{x}V^{y,z}(x)-\partial_{x}V^{y,z}(\hat{x}))(x-\hat{x})\;\pi_{T}^{y,\varphi^{n+1}}(\mathrm{d}x|z)\,\pi_{T}^{y,\varphi^{n+1}}(\mathrm{d}\hat{x}|z)\\\
=&\frac{1}{2}\int\langle\nabla V^{y,n+1}(x,z)-\nabla
V^{y,n+1}(\hat{x}),\;(x,z)-(\hat{x},z)\rangle\;\pi_{T}^{y,\varphi^{n+1}}(\mathrm{d}x|z)\,\pi_{T}^{y,\varphi^{n+1}}(\mathrm{d}\hat{x}|z)\\\
\overset{\eqref{eq:appo:ell:conditional}}{\leq}&\,\frac{1}{2}\int
F_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}(\alpha_{\nu,n},|x-\hat{x}|^{2})\;\pi_{T}^{y,\varphi^{n+1}}(\mathrm{d}x|z)\,\pi_{T}^{y,\varphi^{n+1}}(\mathrm{d}\hat{x}|z)\\\
\leq&\,F_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}\biggl{(}\alpha_{\nu,n},2\,\mathrm{Var}_{X\sim\pi_{T}^{y,\varphi^{n+1}}}(X_{1}|X_{2}=z_{d},\,\dots,\,X_{d}=z_{d})\biggr{)}\;,\end{split}$

where the last step follows from the concavity of $s\mapsto
F_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}(\alpha_{\nu,n},s)$ (cf. Section 3) and
Jensen’s inequality. By combining the above with (103), since
$\alpha_{\nu,n}>-T^{-1}$ and
$F_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}(\alpha_{\nu,n},\cdot)$ is increasing
(cf. Section 3), we deduce by definition of
$G_{\beta_{m}u}^{g_{\mu},\hat{g}_{\nu}}$ that

(106)
$\mathrm{Var}_{X\sim\pi_{T}^{y,\varphi^{n+1}}}(X_{1})\geq\frac{1}{2}\,G_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}(\alpha_{\nu,n},2)\;.$

Since the definition $\ell_{U}$ is invariant under orthonormal transformation,
for any orthonormal matrix $O$ the functions $\varphi^{n+1}(O\,\cdot)$ satisfy
the condition (95) too. The previous bound and this observation leads to

(107) $\mathrm{Var}_{X\sim\pi_{T}^{y,\varphi^{n+1}}}(\langle
v,\,X\rangle)\geq\frac{1}{2}\,G_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}(\alpha_{\nu,n},2)\quad\forall\,y,\,v\in\mathbb{R}^{d}\;\text{s.t.
}|v|=1\;,$

and hence from (100) we finally deduce

(108) $\langle v,\,\nabla^{2}\hat{\psi}^{n+1}(y)\,v\rangle\geq
G_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}(\alpha_{\nu,n},2)\,\frac{|v|^{2}}{2T}\quad\forall\,y,\,v\in\mathbb{R}^{d}\;.$

By recalling (99) and (69), the above bound concludes our proof. ∎

######  Remark  (A first trivial lower bound).

Under A’ 1 and A’ 2, the above discussion already provides a first trivial
lower bound for $\kappa_{\psi^{n}}$. Indeed (100) tells us that
$\hat{\psi}^{n+1}$ is a convex function for any $n\geq 0$, which combined with
(99) yields to for $n\geq 0$

(109)
$\kappa_{\psi^{n+1}}(r)\geq\alpha_{\nu}-T^{-1}-r^{-1}\hat{g}_{\nu}(r)\;.$

Therefore in Theorem 7 we could consider any initialization $\psi^{0}$ without
prescriptions on its behaviour and run an iteration of Sinkhorn in order to
get $\alpha_{\nu,1}\geq\alpha_{\nu}-T^{-1}$. At this point once can proceed
again with the proof of Theorem 7 with the same (but shifted by $-1$) sequence
of parameters $(\alpha_{\nu,n})_{n\in\mathbb{N}}$.

Let us notice that the same discussion holds for the sequence of
$\varphi^{n}$, which yields to for $n\geq 0$

(110)
$\kappa_{\varphi^{n+1}}\geq\alpha_{\mu}-T^{-1}-r^{-1}\hat{g}_{\mu}(r)\;,$

which for $n=0$ gives for free the base step of the iteration (72).

Finally, let us remark here that since the potentials couple
$(\varphi^{\star},\psi^{\star})$ can be thought as a constant sequence of
Sinkhorn iterates, the above discussion proves also that

(111)
$\kappa_{\psi^{\star}}(r)\geq\alpha_{\nu}-T^{-1}-r^{-1}\hat{g}_{\nu}(r)\quad\text{
and
}\quad\kappa_{\varphi^{\star}}\geq\alpha_{\mu}-T^{-1}-r^{-1}\hat{g}_{\mu}(r)\;.$

Given the above lemmata, we are finally able to prove how lower-bounds for
integrated convexity profiles propagate along Sinkhorn.

###### Proof of Theorem 7.

Let us start showing the first statement Consider the sequence
$(\alpha_{\nu,n})_{n\in\mathbb{N}}$ defined in (69). We will prove our
statement by induction. The case $n=0$ is met under the assumption
$\kappa_{\psi^{0}}(r)\geq\alpha_{\nu}-T^{1}-r^{-1}\hat{g}_{\nu}(r)$ The
inductive step follows by applying consecutively Section 3 and Section 3. As a
direct consequence of item $(ii)$ in Section 3 we deduce that the sequence
$(\alpha_{\nu,n})_{n\in\mathbb{N}}$ is non-decreasing and hence
$\alpha_{\nu,n}\geq\alpha_{\nu,0}=\alpha_{\nu}-T^{-1}$. Since
$G_{\beta_{\mu}}^{g_{\mu},\hat{g}_{\nu}}$ is continuous and
$\alpha_{\nu}-T^{-1}$ is not an accumulation point for the set of solutions of
(89) (cf. item $(iii)$ in Section 3), we deduce that
$\alpha_{\nu,n}>\alpha_{\nu,0}=\alpha_{\nu}-T^{-1}$ for $n\geq 1$ and that the
same holds for its limit $\alpha_{\psi^{\star}}$. The upper bound on
$\alpha_{\nu,n}$ comes for free from (69) and the upper bound (94). The proof
of (70) is obtained in the same way by considering the (constant) Sinkhorn
iterates $(\varphi^{\star},\psi^{\star})$ with the same sequence of
$(\alpha_{\nu,n})_{n\in\mathbb{N}}$.

The proof of the second statement is completely analogous and for this reason
we omit it. The only difference here relies in proving that the base case
$n=1$ holds true, but this has been already proven in the discussion of
Section 3. ∎

##  4\. Wasserstein distance w.r.t a measure with log-concave profile

In this section we consider two probability measures
$\mathfrak{p},\,\mathfrak{q}\in\mathcal{P}(\mathbb{R}^{d})$ that can be again
written with log-densities as

(112)
$\mathfrak{p}(\mathrm{d}x)=\exp(-U_{\mathfrak{p}}(x))\mathrm{d}x\,,\qquad\mathfrak{q}(\mathrm{d}x)=\exp(-U_{\mathfrak{q}}(x))\mathrm{d}x\;.$

######  AO 1.

Assume that
$U_{\mathfrak{p}},\,U_{\mathfrak{q}}\in\mathrm{C}^{1}(\mathbb{R}^{d})$ and
that

  1. (1)

$U_{\mathfrak{q}}$ is coercive, i.e., there exist $\gamma_{\mathfrak{q}}>0$
and $R_{\mathfrak{q}}\geq 0$ such that

(113) $-\frac{1}{2}\langle\nabla
U_{\mathfrak{q}}(x),x\rangle\leq-\gamma_{\mathfrak{q}}\,|x|^{2}\quad\forall\,|x|\geq
R_{\mathfrak{q}}\;.$

  2. (2)

$U_{\mathfrak{p}}$ has an integrated convex profile, i.e., there exist some
$\alpha_{\mathfrak{p}}>0$ and
${\tilde{g}}_{\mathfrak{p}}\in\tilde{\mathcal{G}}$ such that

(114)
$\kappa_{U_{\mathfrak{p}}}(r)\geq\alpha_{\mathfrak{p}}-r^{-1}\,{\tilde{g}}_{\mathfrak{p}}(r)\quad\forall
r>0\;.$

Let us also emphasize here that the convexity of integrated profile assumption
is stronger than the coercivity, since the former implies

(115) $-\frac{1}{2}\langle\nabla
U_{\mathfrak{p}}(x),x\rangle\leq-\frac{\alpha_{\mathfrak{p}}}{2}\,|x|^{2}+\frac{|\nabla
U_{\mathfrak{p}}(0)|+\|{\tilde{g}}_{\mathfrak{p}}\|_{\infty}}{2}|x|\;,$

and hence that the coercive condition holds

(116) $-\frac{1}{2}\langle\nabla
U_{\mathfrak{p}}(x),x\rangle\leq-\gamma_{\mathfrak{p}}|x|^{2}\quad\forall\,|x|\geq
R_{\mathfrak{p}}$

for some $\gamma_{\mathfrak{p}}>0$ and $R_{\mathfrak{p}}>0$. Notice that
$\mathfrak{p}$ and $\mathfrak{q}$ can be seen as invariant measures of the
corresponding SDEs

(117) $\begin{cases}\mathrm{d}X_{t}=-\frac{1}{2}\nabla
U_{\mathfrak{p}}(X_{t})\,\mathrm{d}t+\mathrm{d}B_{t}\;,\\\
\mathrm{d}Y_{t}=-\frac{1}{2}\nabla
U_{\mathfrak{q}}(Y_{t})\,\mathrm{d}t+\mathrm{d}B_{t}\;,\end{cases}$

which admit unique strong solutions, in view of the coercivity of the
corresponding drifts and owing to [56, Theorem 2.1]. Finally let
$(P_{t}^{\mathfrak{p}})_{t\geq 0}$ and $(P_{t}^{\mathfrak{q}})_{t\geq 0}$
denote the corresponding Markov semigroups associated to the above SDEs. Since
$\mathfrak{p}$ and $\mathfrak{q}$ are invariant measures we clearly have
$\mathfrak{p}P^{\mathfrak{p}}_{t}=\mathfrak{p}$ and
$\mathfrak{q}P_{t}^{\mathfrak{q}}=\mathfrak{q}$ for any $t\geq 0$.

The main result of this section is showing how the Wasserstein distance
between $\mathfrak{p}$ and $\mathfrak{q}$ can be bounded w.r.t. the integrated
difference of the drifts appearing in (117).

######  Theorem 9.

Assume AO 1. Then it holds

(118)
$\mathbf{W}_{1}(\mathfrak{p},\;\mathfrak{q})\leq\biggl{(}\frac{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)}{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)}\biggr{)}^{2}\frac{1}{\alpha_{\mathfrak{p}}+\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)}\;\int|\nabla
U_{\mathfrak{p}}-\nabla
U_{\mathfrak{q}}|\mathrm{d}\mathfrak{q}=\biggl{(}\frac{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)}{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)}\biggr{)}^{2}\frac{1}{\alpha_{\mathfrak{p}}+\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)}\;\int\left|\nabla\log\frac{\mathrm{d}\mathfrak{p}}{\mathrm{d}\mathfrak{q}}\right|\mathrm{d}\mathfrak{q}\;,$

with
$R\coloneqq\left\|\tilde{g}_{\mathfrak{p}}\right\|_{\infty}((\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0))^{-1}+2/\alpha_{\mathfrak{p}})$.

###### Proof.

Firstly, let us consider the function

(119)
$f_{\mathfrak{p}}(r)\coloneqq\begin{cases}\tilde{g}_{\mathfrak{p}}(r)\quad&\text{
if }r\leq R\;,\\\
\tilde{g}_{\mathfrak{p}}(R)+\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)\,(r-R)\quad&\text{
otherwise.}\end{cases}$

Notice that $f(0)=0$ and that
$f_{\mathfrak{p}}\in\mathcal{C}^{1}((0,+\infty),\mathbb{R}_{+})\cap\mathcal{C}^{2}((0,R)\cup(R,+\infty),\mathbb{R}_{+})$
with $f_{\mathfrak{p}}{{}^{\prime\prime}}$ having a jump discontinuity in
$r=R$. Moreover, $f_{\mathfrak{p}}$ is non-decreasing and concave and
equivalent to the identity i.e., for any $r>0$ it holds

(120) $\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)\,r\leq
f_{\mathfrak{p}}(r)\leq\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)\,r\;,\quad\text{
since }\;\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)\leq
f_{\mathfrak{p}}{{}^{\prime}}(r)\leq\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)\;.$

Notice that, since $\tilde{g}_{\mathfrak{p}}\in\tilde{\mathcal{G}}$, for any
$r\in(0,R)$ it holds

(121) $\displaystyle
2\;f_{\mathfrak{p}}^{{}^{\prime\prime}}(r)-\frac{r_{t}\;f_{\mathfrak{p}}^{{}^{\prime}}(r)}{2}\kappa_{U_{\mathfrak{p}}}(r)\leq$
$\displaystyle\,2\;\tilde{g}_{\mathfrak{p}}^{{}^{\prime\prime}}(r)-\frac{\alpha_{\mathfrak{p}}}{2}\,\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(r)\,r+\frac{1}{2}{\tilde{g}}_{\mathfrak{p}}(r)\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(r)\leq-\frac{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(r)}{2}\,(\alpha_{\mathfrak{p}}\,r+\tilde{g}_{\mathfrak{p}}(r))$
$\displaystyle\leq$
$\displaystyle\,-\frac{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)}{2}\,\biggl{(}\frac{\alpha_{\mathfrak{p}}}{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)}+1\biggr{)}\;\tilde{g}_{\mathfrak{p}}(r)=-\frac{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)}{2}\,\biggl{(}\frac{\alpha_{\mathfrak{p}}}{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)}+1\biggr{)}\;f_{\mathfrak{p}}(r)\;;$

whereas for any $r>R$ it holds

(122) $\displaystyle
2\;f_{\mathfrak{p}}^{{}^{\prime\prime}}(r)-\frac{r_{t}\;f_{\mathfrak{p}}^{{}^{\prime}}(r)}{2}\kappa_{U_{\mathfrak{p}}}(r)=$
$\displaystyle\,-\frac{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)}{2}\;r\,\kappa_{U_{\mathfrak{p}}}(r)\leq-\frac{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)}{2}\,(\alpha_{\mathfrak{p}}\,r-\tilde{g}_{\mathfrak{p}}(r))$
$\displaystyle\leq$
$\displaystyle\,-\frac{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)}{2}\,\alpha_{\mathfrak{p}}\,\frac{R}{\tilde{g}_{\mathfrak{p}}(R)}\,f_{\mathfrak{p}}(r)+\frac{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)}{2}\,f_{\mathfrak{p}}(r)\;,$

where the last step follows from the concavity of $\tilde{g}_{\mathfrak{p}}$,
which implies for any $r\geq R$ that $\tilde{g}_{\mathfrak{p}}(r)\leq
f_{\mathfrak{p}}(r)$ and the monotonicity of $F(r)\coloneqq r/f(r)$, since for
any $r\geq R$ it holds

(123)
$F{{}^{\prime}}(r)=\frac{f_{\mathfrak{p}}(r)-r\,f_{\mathfrak{p}}{{}^{\prime}}(r)}{f_{\mathfrak{p}}(r)^{2}}=\frac{\tilde{g}_{\mathfrak{p}}(R)-R\,\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)}{f_{\mathfrak{p}}(r)^{2}}\geq
0\;.$

Therefore, by recalling the definition of $R$, for any
$r\in(0,R)\cup(R,+\infty)$ we have shown that

(124)
$2\;f_{\mathfrak{p}}^{{}^{\prime\prime}}(r)-\frac{r_{t}\;f_{\mathfrak{p}}^{{}^{\prime}}(r)}{2}\kappa_{U_{\mathfrak{p}}}(r)\leq-\lambda_{\mathfrak{p}}\,f_{\mathfrak{p}}(r)\quad\text{
with
}\;\lambda_{\mathfrak{p}}\coloneqq\frac{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)}{2}\,\biggl{(}\frac{\alpha_{\mathfrak{p}}}{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)}+1\biggr{)}\;.$

By relying on the above construction, we will prove our thesis considering the
Wasserstein distance

(125)
$\mathbf{W}_{f_{\mathfrak{p}}}(\mathfrak{p},\,\mathfrak{q})\coloneqq\inf_{\pi\in\Pi(\mathfrak{p},\mathfrak{q})}\mathbb{E}_{(X,Y)\sim\pi}\bigl{[}f_{\mathfrak{p}}(|X-Y|)\bigr{]}$

induced by the concave function $f_{\mathfrak{p}}$. The above is indeed a
distance since $f_{\mathfrak{p}}(0)=0$, $f_{\mathfrak{p}}$ is strictly
increasing, concave and hence also subadditive (which implies the triangular
inequality). Moreover, from (120) it follows the equivalence between
$\mathbf{W}_{f_{\mathfrak{p}}}$ and the usual $\mathbf{W}_{1}$, namely

(126)
$\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(R)\;\mathbf{W}_{1}(\mathfrak{p},\,\mathfrak{q})\leq\mathbf{W}_{f_{\mathfrak{p}}}(\mathfrak{p},\,\mathfrak{q})\leq\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)\;\mathbf{W}_{1}(\mathfrak{p},\,\mathfrak{q})\;,$

For any $t\geq 0$ notice that

(127)
$\begin{split}\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{p},\,\mathfrak{q})\leq\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{p},\,\mathfrak{q}\,P_{t}^{\mathfrak{p}})+\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{q}\,P_{t}^{\mathfrak{p}},\,\mathfrak{q})=\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{p}\,P_{t}^{\mathfrak{p}},\,\mathfrak{q}\,P_{t}^{\mathfrak{p}})+\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{q}\,P_{t}^{\mathfrak{p}},\,\mathfrak{q})\;.\end{split}$

In order to bound the first Wasserstein distance appearing in the upper bound
(127), namely
$\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{p}\,P_{t}^{\mathfrak{p}},\,\mathfrak{q}\,P_{t}^{\mathfrak{p}})$,
fix an initial coupling $(X_{0},X^{\mathfrak{q}}_{0})\sim\pi^{\star}$
distributed according to the optimal coupling for
$\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{p},\,\mathfrak{q})$ and consider
the reflection coupling diffusion processes

(128) $\begin{cases}\mathrm{d}X_{t}=-\frac{1}{2}\nabla
U_{\mathfrak{p}}(X_{t})\,\mathrm{d}t+\mathrm{d}B_{t}\\\
\mathrm{d}X^{\mathfrak{q}}_{t}=-\frac{1}{2}\nabla
U_{\mathfrak{p}}(X^{\mathfrak{q}}_{t})\,\mathrm{d}t+\mathrm{d}\hat{B}_{t}\quad\forall\;t\in[0,\tau)\;\text{
and }X_{t}=X^{\mathfrak{q}}_{t}\quad\forall\;t\geq\tau\\\
(X_{0},X^{\mathfrak{q}}_{0})\sim\pi^{\star}\;,\end{cases}$

where $\tau\coloneqq\inf\\{s\geq 0:X^{\mathfrak{q}}_{s}=X_{s}\\}$, and
$(\hat{B}_{t})_{t\geq 0}$ is defined as

(129) $\mathrm{d}\hat{B}_{t}\coloneqq(\operatorname{Id}-2\,e_{t}\,e_{t}^{{\sf
T}}\,\mathbf{1}_{\\{t<\tau\\}})\,\mathrm{d}B_{t}\qquad\text{where}\quad
e_{t}\coloneqq\begin{cases}\frac{Z_{t}}{|Z_{t}|}\quad&\text{ when
}r_{t}>0\,,\\\ u\quad&\text{ when }r_{t}=0\,.\end{cases}$

where $Z_{t}\coloneqq X_{t}-X^{\mathfrak{q}}_{t}$, $r_{t}\coloneqq|Z_{t}|$ and
$u\in\mathbb{R}^{d}$ is a fixed (arbitrary) unit-vector. By Lévy’s
characterization, $(\hat{B}_{t})_{t\geq 0}$ is a $d$-dimensional Brownian
motion. As a result, $X_{t}\sim\mathfrak{p}$ and
$X^{\mathfrak{q}}_{t}\sim\mathfrak{q}P^{\mathfrak{p}}_{t}$ for any $t\geq 0$.
In addition $\mathrm{d}W_{t}\coloneqq e_{t}^{{\sf T}}\mathrm{d}B_{t}$ is a
one-dimensional Brownian motion. Let us notice that for any $t<\tau$ it holds

(130) $\displaystyle\mathrm{d}Z_{t}=$ $\displaystyle\,-2^{-1}\,(\nabla
U_{\mathfrak{p}}(X_{t})-\nabla
U_{\mathfrak{p}}(X^{\mathfrak{q}}_{t}))\,\mathrm{d}t+2\;e_{t}\;\mathrm{d}W_{t}\;,$
$\displaystyle\mathrm{d}r_{t}^{2}=$ $\displaystyle\,-\langle Z_{t},\;\nabla
U_{\mathfrak{p}}(X_{t})-\nabla
U_{\mathfrak{p}}(X^{\mathfrak{q}}_{t})\rangle\;\mathrm{d}t+4\;\mathrm{d}t+4\;\langle
Z_{t},\;e_{t}\rangle\;\mathrm{d}W_{t}\;,$ $\displaystyle\mathrm{d}r_{t}=$
$\displaystyle\,-2^{-1}\,\langle e_{t},\;\nabla U_{\mathfrak{p}}(X_{t})-\nabla
U_{\mathfrak{p}}(X^{\mathfrak{q}}_{t})\rangle\;\mathrm{d}t+2\;\mathrm{d}W_{t}\;.$

Now, an application of Ito-Tanaka formula [55, Chapter VI, Theorem 1.5] to the
concave function
$f_{\mathfrak{p}}\in\mathcal{C}^{1}((0,+\infty),\mathbb{R}_{+})\cap\mathcal{C}^{2}((0,R)\cup(R,+\infty),\mathbb{R}_{+})$,
gives for any $t<\tau$

(131)
$f_{\mathfrak{p}}(r_{t})=f_{\mathfrak{p}}(r_{0})+\int_{0}^{t}f_{\mathfrak{p}}{{}^{\prime}}(r_{s})\,\mathrm{d}r_{s}+\frac{1}{2}\int_{\mathbb{R}}L_{t}^{a}\,\mu_{f_{\mathfrak{p}}}(\mathrm{d}a)\;,$

where $(L_{t}^{a})_{t}$ denotes the right-continuous local time of the
semimartingale $(r_{t})_{t}$, whereas $\mu_{f_{\mathfrak{p}}}$ is the non-
positive measure representing $f_{\mathfrak{p}}{{}^{\prime\prime}}$ in the
sense of distributions, i.e.,
$\mu_{f_{\mathfrak{p}}}([a,b])=f_{\mathfrak{p}}{{}^{\prime}}(b)-f_{\mathfrak{p}}{{}^{\prime}}(a)$
for any $a\leq b$. Let us further notice that the Meyer Wang occupation times
formula [38, Theorem 29.5], which for any measurable function
$H\colon\mathbb{R}\to[0,+\infty)$ reads as

(132)
$\int_{0}^{t}H(r_{s})\,\mathrm{d}[r]_{s}=\int_{\mathbb{R}}H(a)L_{t}^{a}\,\mathrm{d}a\;,$

implies that the random set $\\{s\in[0,\tau]\colon r_{s}=R\\}$ has almost
surely zero Lebesgue measure. Particularly, since $\mu_{f_{\mathfrak{p}}}$ is
non-positive, this combined with the above formula implies

(133)
$\frac{1}{2}\int_{\mathbb{R}}L_{t}^{a}\,\mu_{f_{\mathfrak{p}}}(\mathrm{d}a)\leq\frac{1}{2}\int_{\mathbb{R}}\mathbf{1}_{\\{a\neq
R\\}}\,L_{t}^{a}\,f_{\mathfrak{p}}{{}^{\prime\prime}}(a)\,\mathrm{d}a=2\int_{0}^{t}\mathbf{1}_{\\{r_{s}\neq
R\\}}\,f_{\mathfrak{p}}{{}^{\prime\prime}}(r_{s})\,\mathrm{d}s=2\int_{0}^{t}\,f_{\mathfrak{p}}{{}^{\prime\prime}}(r_{s})\,\mathrm{d}s\;.$

As a byproduct of (131) and (133) we may finally state that for any $t<\tau$
it almost surely holds

(134)
$\begin{split}\mathrm{d}f_{\mathfrak{p}}(r_{t})\leq&-\frac{f_{\mathfrak{p}}^{{}^{\prime}}(r_{t})}{2}\,\langle
e_{t},\;\nabla U_{\mathfrak{p}}(X_{t})-\nabla
U_{\mathfrak{p}}(X^{\mathfrak{q}}_{t})\rangle\;\mathrm{d}t+2\;f_{\mathfrak{p}}^{{}^{\prime\prime}}(r_{t})\;\mathrm{d}t+2\;f_{\mathfrak{p}}^{{}^{\prime}}(r_{t})\;\mathrm{d}W_{t}\\\
\leq&\biggl{(}2\;f_{\mathfrak{p}}^{{}^{\prime\prime}}(r_{t})-\frac{r_{t}\;f_{\mathfrak{p}}^{{}^{\prime}}(r_{t})}{2}\kappa_{U_{\mathfrak{p}}}(r_{t})\biggr{)}\;\mathrm{d}t+2\;f_{\mathfrak{p}}^{{}^{\prime}}(r_{t})\;\mathrm{d}W_{t}\overset{\eqref{eq:appo:diff:general}}{\leq}\;-\lambda_{\mathfrak{p}}\;f_{\mathfrak{p}}(r_{t})\;\mathrm{d}t+2\;f_{\mathfrak{p}}^{{}^{\prime}}(r_{t})\;\mathrm{d}W_{t}\;.\end{split}$

By recalling that $f_{\mathfrak{p}}(r_{t})=f_{\mathfrak{p}}(0)=0$ as soon as
$t\geq\tau$, by taking expectation, integrating over time and Gronwall Lemma
we have finally proven that for any $t\geq 0$

(135)
$\begin{split}\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{p}\,P_{t}^{\mathfrak{p}},\,\mathfrak{q}\,P_{t}^{\mathfrak{p}})\leq\mathbb{E}[{f_{\mathfrak{p}}}(|X_{t}-X^{\mathfrak{q}}_{t}|)]\leq
e^{-\lambda_{\mathfrak{p}}\,t}\;\mathbb{E}[f_{\mathfrak{p}}(|X_{0}-X^{\mathfrak{q}}_{0}|)]=e^{-\lambda_{\mathfrak{p}}\,t}\;\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{p},\,\mathfrak{q})\,.\end{split}$

Let us now provide a bound for the second Wasserstein distance appearing in
the upper bound (127), namely
$\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{q}\,P_{t}^{\mathfrak{p}},\,\mathfrak{q})$,
by relying on synchronous coupling technique. Therefore fix an initial random
variable $Y_{0}\sim\mathfrak{q}$ and a $d$-dimensional Brownian motion
$(B_{t})_{t\geq 0}$, and consider now the diffusion processes

(136) $\begin{cases}\mathrm{d}X^{\mathfrak{q}}_{t}=-\frac{1}{2}\nabla
U_{\mathfrak{p}}(X^{\mathfrak{q}}_{t})\,\mathrm{d}t+\mathrm{d}B_{t}\\\
\mathrm{d}Y_{t}=-\frac{1}{2}\nabla
U_{\mathfrak{q}}(Y_{t})\,\mathrm{d}t+\mathrm{d}B_{t}\\\
X^{\mathfrak{q}}_{0}=Y_{0}=Y_{0}\sim\mathfrak{q}\;.\end{cases}$

Notice that for any $t\geq 0$, $Y_{t}\sim\mathfrak{q}$ whereas
$X^{\mathfrak{q}}_{t}\sim\mathfrak{q}P^{\mathfrak{p}}_{t}$. If we set now
$\bar{r}_{t}\coloneqq|X^{\mathfrak{q}}_{t}-Y_{t}|$ we then have

(137) $\displaystyle\mathrm{d}(X^{\mathfrak{q}}_{t}-Y_{t})=$
$\displaystyle\;-2^{-1}(\nabla U_{\mathfrak{p}}(X^{\mathfrak{q}}_{t})-\nabla
U_{\mathfrak{q}}(Y_{t}))\;\mathrm{d}t\;,$
$\displaystyle\mathrm{d}\bar{r}_{t}^{2}=$ $\displaystyle\;-\langle
X^{\mathfrak{q}}_{t}-Y_{t},\;\nabla
U_{\mathfrak{p}}(X^{\mathfrak{q}}_{t})-\nabla
U_{\mathfrak{q}}(Y_{t})\rangle\;\mathrm{d}t\;.$

At this point we would like to apply the square-root function, however the
latter fails to be $\mathrm{C}^{2}$ in the origin whereas $\bar{r}_{t}$ may be
equal to zero (e.g., we already start with $\bar{r}_{0}=0$). For this reason
we are going to perform an approximation argument. Fix $\delta>0$ and consider
the function $\rho_{\delta}(r)\coloneqq\sqrt{r+\delta}$. Then it holds

(138) $\displaystyle\mathrm{d}\rho_{\delta}(\bar{r}_{t}^{2})=$
$\displaystyle\;-(2\;\rho_{\delta}(\bar{r}_{t}^{2}))^{-1}\langle
X^{\mathfrak{q}}_{t}-Y_{t},\;\nabla
U_{\mathfrak{p}}(X^{\mathfrak{q}}_{t})-\nabla
U_{\mathfrak{q}}(Y_{t})\rangle\;\mathrm{d}t$ $\displaystyle=$
$\displaystyle\;-(2\;\rho_{\delta}(\bar{r}_{t}^{2}))^{-1}\langle
X^{\mathfrak{q}}_{t}-Y_{t},\;\nabla
U_{\mathfrak{p}}(X^{\mathfrak{q}}_{t})-\nabla
U_{\mathfrak{p}}(Y_{t})\rangle\;\mathrm{d}t$
$\displaystyle\qquad\qquad\qquad\qquad-(2\;\rho_{\delta}(\bar{r}_{t}^{2}))^{-1}\langle
X^{\mathfrak{q}}_{t}-Y_{t},\;\nabla U_{\mathfrak{p}}(Y_{t})-\nabla
U_{\mathfrak{q}}(Y_{t})\rangle\;\mathrm{d}t$ $\displaystyle\leq$
$\displaystyle\;-2^{-1}\;\frac{\bar{r}_{t}^{2}}{\rho_{\delta}(\bar{r}_{t}^{2})}\;\kappa_{\mathfrak{p}}(\bar{r}_{t})\;\mathrm{d}t+2^{-1}\;\frac{\bar{r}_{t}}{\rho_{\delta}(\bar{r}_{t}^{2})}\;|\nabla
U_{\mathfrak{p}}-\nabla U_{\mathfrak{q}}|(Y_{t})\mathrm{d}t$
$\displaystyle\leq$
$\displaystyle\;-2^{-1}\;\frac{\bar{r}_{t}^{2}}{\rho_{\delta}(\bar{r}_{t}^{2})}\;(\alpha_{\mathfrak{p}}-\tilde{G}_{\mathfrak{p}})\;\mathrm{d}t+2^{-1}\;\frac{\bar{r}_{t}}{\rho_{\delta}(\bar{r}_{t}^{2})}\;|\nabla
U_{\mathfrak{p}}-\nabla U_{\mathfrak{q}}|(Y_{t})\mathrm{d}t\;,$

where in the last step we have relied on the sublinearity of
${\tilde{g}}_{\mathfrak{p}}\in\tilde{\mathcal{G}}$ (namely that
${\tilde{g}}_{\mathfrak{p}}(r)\leq\tilde{G}_{\mathfrak{p}}\,r$ for some
positive constant ${\tilde{G}}_{\mathfrak{p}}>0$). Therefore it holds

(139)
$\displaystyle\mathrm{d}\rho_{\delta}(\bar{r}_{t}^{2})\leq\;\frac{(\alpha_{\mathfrak{p}}-\tilde{G}_{\mathfrak{p}})^{+}}{2}\;\rho_{\delta}(\bar{r}_{t}^{2})\mathrm{d}t\;+2^{-1}|\nabla
U_{\mathfrak{p}}-\nabla U_{\mathfrak{q}}|(Y_{t})\mathrm{d}t\;.$

By taking expectation and integrating over time the above bound gives

(140)
$\begin{split}\mathbb{E}[\rho_{\delta}(\bar{r}_{t}^{2})]\leq\mathbb{E}[\rho_{\delta}(\bar{r}_{0}^{2})]+\frac{(\alpha_{\mathfrak{p}}-\tilde{G}_{\mathfrak{p}})^{+}}{2}\int_{0}^{t}\mathbb{E}[\rho_{\delta}(\bar{r}_{s}^{2})]\;\mathrm{d}s+2^{-1}\int_{0}^{t}\mathbb{E}[|\nabla
U_{\mathfrak{p}}-\nabla U_{\mathfrak{q}}|(Y_{s})]\;\mathrm{d}s\\\
=\sqrt{\delta}+\frac{(\alpha_{\mathfrak{p}}-\tilde{G}_{\mathfrak{p}})^{+}}{2}\int_{0}^{t}\mathbb{E}[\rho_{\delta}(\bar{r}_{s}^{2})]\;\mathrm{d}s+\frac{t}{2}\,\int|\nabla
U_{\mathfrak{p}}-\nabla
U_{\mathfrak{q}}|\;\mathrm{d}\mathfrak{q}\;,\end{split}$

where in the last step we have relied on the fact that for any $t\geq 0$
$Y_{t}\sim\mathfrak{q}$ and that $\bar{r}_{0}=0$. Therefore Gronwall Lemma
yields to

(141)
$\displaystyle\mathbb{E}[\bar{r}_{t}]\leq\mathbb{E}[\rho_{\delta}(\bar{r}_{t}^{2})]\leq\frac{t}{2}\;\exp\biggl{(}\frac{t}{2}\;(\alpha_{\mathfrak{p}}-\tilde{G}_{\mathfrak{p}})^{+}\biggr{)}\;\int|\nabla
U_{\mathfrak{p}}-\nabla
U_{\mathfrak{q}}|\;\mathrm{d}\mathfrak{q}+\sqrt{\delta}\;\exp\biggl{(}\frac{t}{2}\;(\alpha_{\mathfrak{p}}-\tilde{G}_{\mathfrak{p}})^{+}\biggr{)}\;.$

By letting $\delta$ to zero in the above right-hand-side, we obtain the
desired upper bound

(142)
$\displaystyle\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{q}\,P_{t}^{\mathfrak{p}},\,\mathfrak{q})\overset{\eqref{eq:appo:equivalence:wasser:general}}{\leq}$
$\displaystyle\;\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)\;\mathbf{W}_{1}(\mathfrak{q}\,P_{t}^{\mathfrak{p}},\,\mathfrak{q})\leq\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)\;\mathbb{E}[|X^{\mathfrak{q}}_{t}-Y_{t}|]=\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)\;\mathbb{E}[\bar{r}_{t}]$
$\displaystyle\leq$
$\displaystyle\;\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)\;\frac{t}{2}\;\exp\biggl{(}\frac{t}{2}\;(\alpha_{\mathfrak{p}}-\tilde{G}_{\mathfrak{p}})^{+}\biggr{)}\;\int|\nabla
U_{\mathfrak{p}}-\nabla U_{\mathfrak{q}}|\;\mathrm{d}\mathfrak{q}\;.$

By putting together the last estimate with (127) and (135) we have proven

(143)
$\begin{split}\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{p},\,\mathfrak{q})\leq&\,e^{-\lambda_{\mathfrak{p}}\,t}\;\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{p},\,\mathfrak{q})+\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)\;\frac{t}{2}\;\exp\biggl{(}\frac{t}{2}\;(\alpha_{\mathfrak{p}}-\tilde{G}_{\mathfrak{p}})^{+}\biggr{)}\;\int|\nabla
U_{\mathfrak{p}}-\nabla
U_{\mathfrak{q}}|\;\mathrm{d}\mathfrak{q}\;.\end{split}$

or equivalently

(144)
$\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{p},\,\mathfrak{q})\leq\;\frac{t/2}{1-e^{-\lambda_{\mathfrak{p}}\,t}}\;\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)\;\exp\biggl{(}\frac{t}{2}\;(\alpha_{\mathfrak{p}}-\tilde{G}_{\mathfrak{p}})^{+}\biggr{)}\;\int|\nabla
U_{\mathfrak{p}}-\nabla U_{\mathfrak{q}}|\;\mathrm{d}\mathfrak{q}\;,$

which in the $t$ vanishing limit reads as

(145)
$\mathbf{W}_{{f_{\mathfrak{p}}}}(\mathfrak{p},\,\mathfrak{q})\leq\frac{\tilde{g}_{\mathfrak{p}}{{}^{\prime}}(0)}{2\;\lambda_{\mathfrak{p}}}\;\int|\nabla
U_{\mathfrak{p}}-\nabla U_{\mathfrak{q}}|\;\mathrm{d}\mathfrak{q}\;.$

Combining the above bound with the equivalence (126) concludes the proof. ∎

######  Remark  (Explicit constants for $\tilde{g}_{\mathfrak{p}}$ as in
(22)).

Particularly, when $\tilde{g}_{\mathfrak{p}}$ is as in (22) in Section 2.2,
i.e.

(146) $\tilde{g}_{\mathfrak{p}}(r)=2\,(L)^{1/2}\,\tanh(r\,L^{1/2}/2)$

for some $L$, then in the previous proof we have
$R=2(L)^{1/2}(L^{-1}+2/\alpha_{\mathfrak{p}})$ and therefore

(147)
$\mathbf{W}_{1}(\mathfrak{p},\;\mathfrak{q})\leq\frac{\cosh^{4}(2(L)^{1/2}(L^{-1}+2/\alpha_{\mathfrak{p}}))}{\alpha_{\mathfrak{p}}+L}\;\int|\nabla
U_{\mathfrak{p}}-\nabla U_{\mathfrak{q}}|\mathrm{d}\mathfrak{q}\;.$

Let us conclude this section by showing how the previous result and the
propagation of lower-bounds for integrated convexity profiles (studied in the
previous section) yield to the key estimates that will be employed in the
proof of the main results.

######  Corollary .

Assume A 1. Then, for any $x\in\mathbb{R}^{d}$, it holds

(148)
$\mathbf{W}_{1}(\pi_{T}^{x,\psi^{n}},\pi_{T}^{x,\psi^{\star}})\leq\gamma_{n}^{\nu}\int\left|\nabla\psi^{n}-\nabla\psi^{\star}\right|\mathrm{d}\pi_{T}^{x,\psi^{\star}},$

and similarly

(149)
$\mathbf{W}_{1}(\pi_{T}^{x,\varphi^{n+1}},\pi_{T}^{x,\varphi^{\star}})\leq\gamma_{n}^{\mu}\int\left|\nabla\varphi^{n+1}-\nabla\varphi^{\star}\right|\mathrm{d}\pi_{T}^{x,\varphi^{\star}},$

where $\gamma_{n}^{\mu}$ and $\gamma_{n}^{\nu}$ are given in Section 2.2.

###### Proof.

Inequality (148) follows from the previous theorem when considering
$\mathfrak{p}=\pi_{T}^{x,\psi^{n}}$ and
$\mathfrak{q}=\pi^{x,\psi^{\star}}_{T}$. Indeed these two probabilities are
the invariant measures associated to (117) with
$U_{\mathfrak{p}}(y)=(2T)^{-1}\,|y-x|^{2}+\psi^{n}(y)$ and
$U_{\mathfrak{q}}(y)=(2T)^{-1}\,|y-x|^{2}+\psi^{\star}(y)$ respectively.
Theorem 7 guarantees that there exists
$\alpha_{\nu,n},\,\alpha_{\psi^{\star}}>-T^{-1}$ such that

(150) $\kappa_{U_{\mathfrak{p}}}(r)\geq
T^{-1}+\alpha_{\nu,n}-r^{-1}{\tilde{g}}_{\nu}(r)\quad\text{ and
}\quad\kappa_{U_{\mathfrak{q}}}(r)\geq
T^{-1}+\alpha_{\psi^{\star}}-r^{-1}{\tilde{g}}_{\nu}(r)\;.$

Therefore $\mathfrak{p}$ and $\mathfrak{q}$ satisfy Assumption AO 1 and
Theorem 9 gives (148).

We omit the details for the proof of (149) since it can be obtained in the
same way, this time considering $\mathfrak{p}=\pi_{T}^{x,\varphi^{n+1}}$ and
$\mathfrak{q}=\pi_{T}^{x,\varphi^{\star}}$. ∎

##  5\. Proof of the main results

Let us simply remark here that the sequences
$(\gamma^{\mu}_{n})_{n\in\mathbb{N}}$ and
$(\gamma^{\nu}_{n})_{n\in\mathbb{N}}$ given in Section 2.2 appearing in the
following results are non-increasing, i.e., they provide faster convergence as
the index $n$ increases. It will be clear from our exposition that we will
deduce from the convergence of the gradients, both the convergence of the
iterates on primal and dual formulation as well as the convergence of the
Hessians.

###  5.1. Exponential convergence of the gradients and $\mathbf{W}_{1}$ along
Sinkhorn’s algorithm

The convergence of the gradients will follow by iterating the following
result.

######  Proposition .

Assume A 1 holds. Then for any $n\geq 0$ it holds

(151) $\int|\nabla\varphi^{n+1}-\nabla\varphi^{\star}|\mathrm{d}\mu\leq
T^{-1}\,\gamma_{n}^{\nu}\int|\nabla\psi^{n}-\nabla\psi^{\star}|\mathrm{d}\nu\;,$

and similarly

(152) $\int|\nabla\psi^{n+1}-\nabla\psi^{\star}|\mathrm{d}\nu\leq
T^{-1}\,\gamma_{n}^{\mu}\int|\nabla\varphi^{n+1}-\nabla\varphi^{\star}|\mathrm{d}\mu\;,$

where $\gamma^{\mu}_{n}$ and $\gamma^{\nu}_{n}$ are given in Section 2.2.

###### Proof.

Let us start by showing (151). From (1) we immediatly have

(153) $\begin{cases}\varphi^{n+1}-{\varphi^{\star}}=\log
P_{T}\exp(-\psi^{n})-\log P_{T}\exp(-\psi^{\star})\\\
\psi^{n+1}-\psi^{\star}=\log P_{T}\exp(-\varphi^{n+1})-\log
P_{T}\exp(-\varphi^{\star})\end{cases}$

and since for $h=\varphi^{n+1},\,\varphi^{\star}$ the gradient along the
semigroup has the explicit formulation (cf. Section 2.2)

(154) $\nabla\log
P_{T}\exp(-h)(x)=\frac{1}{T}\int(y-x)\,\pi_{T}^{x,h}(\mathrm{d}y)\,,$

we deduce

(155) $|\nabla\varphi^{n+1}-\nabla\varphi^{\star}|(x)=T^{-1}\biggl{|}\int
y\,\pi_{T}^{x,\psi^{n}}(\mathrm{d}y)-\int
y\,\pi_{T}^{x,\psi^{\star}}(\mathrm{d}y)\biggr{|}\leq
T^{-1}\mathbf{W}_{1}(\pi_{T}^{x,\psi^{n}},\pi_{T}^{x,\psi^{\star}})\,.$

Combining the above observation with Section 4 we end up with

(156)
$|\nabla\varphi^{n+1}-\nabla\varphi^{\star}|(x)\leq\;T^{-1}\;\gamma_{n}^{\nu}\int\left|\nabla\psi^{n}-\nabla\psi^{\star}\right|\mathrm{d}\pi_{T}^{x,\psi^{\star}}\;,$

with $\gamma_{n}^{\nu}$ as introduced in Section 2.2. Combining (28) with
(156) gives (151). The contractive estimate (152) can be proven in the same
manner by relying on (149) by noticing that

(157) $|\nabla\psi^{n+1}-\nabla\psi^{\star}|(x)\leq
T^{-1}\mathbf{W}_{1}(\pi_{T}^{x,\varphi^{n+1}},\pi_{T}^{x,\varphi^{\star}})\;\text{
and
}\;\int\pi_{T}^{x,\varphi^{\star}}(\mathrm{d}y)\;\nu(\mathrm{d}x)=\mu(\mathrm{d}y)\;.$

∎

As a corollary we immediately deduce Theorem 2.

###### Proof of Theorem 2.

As concerns (33), it is enough concatenating the bounds proven in Section 5.1,
as already sketched at (31) and (32).

The $\mathbf{W}_{1}$-convergence bound (34) can be deduce from (33) and
Section 4 since $\pi^{n,n}\in\Pi(\star,\nu)$, $\pi^{n+1,n}\in\Pi(\mu,\star)$
and hence

(158)
$\mathbf{W}_{1}(\pi^{n,n},\pi^{\star})\leq\int\mathbf{W}_{1}(\pi^{x,\varphi^{n}}_{T},\pi^{x,\varphi^{\star}}_{T})\;\mathrm{d}\nu\quad\text{
and
}\mathbf{W}_{1}(\pi^{n+1,n},\pi^{\star})\leq\int\mathbf{W}_{1}(\pi^{x,\psi^{n}}_{T},\pi^{x,\psi^{\star}}_{T})\;\mathrm{d}\mu\;.$

∎

######  Remark  (Convergence along adjusted marginals).

Sinkhorn iterates may be considered as potentials of appropriate Schrödinger
problems. Indeed, the decomposition given in (8) and (9) implies that

  * •

the couple $(\varphi^{n+1},\psi^{n})$ corresponds to a couple of Schrödinger
potentials (as defined in (1)) associated to the Schrödinger problem with
reference measure $\mathrm{R}_{0,T}$ and with marginals $\mu$ and
$\nu^{n}\coloneqq(\mathrm{proj}_{y})_{\sharp}\pi^{n+1,n}$;

  * •

the couple $(\varphi^{n+1},\psi^{n+1})$ corresponds to a couple of Schrödinger
potentials (as defined in (1)) associated to the Schrödinger problem with
reference measure $\mathrm{R}_{0,T}$ and with marginals
$\mu^{n+1}\coloneqq(\mathrm{proj}_{x})_{\sharp}\pi^{n+1,n+1}$ and $\nu$.

This simple observation, the bound (155), the conditional property of
$\pi_{T}^{x,\psi^{n}}$222Let us recall that $\pi^{n,n}$ is the optimal
coupling for the Schrödinger problem with marginals $\mu^{n}$ and $\nu$, that
the corresponding potentials are given by $\varphi^{n},\,\psi^{n}$ and
therefore
$\pi_{T}^{x,\psi^{n}}(\mathrm{d}y)\mu^{n}(\mathrm{d}x)=\pi^{n,n}(\mathrm{d}x\,\mathrm{d}y)$.,
and arguing as in Section 4 (this time with
$\mathfrak{p}=\pi^{x,\psi^{\star}}_{T}$ and
$\mathfrak{q}=\pi^{x,\psi^{n}}_{T}$) prove that

(159)
$\displaystyle\int|\nabla\varphi^{n}-\nabla\varphi^{\star}|\mathrm{d}\mu^{n}\leq$
$\displaystyle\;\frac{\gamma_{\infty}^{\nu}}{T}\;\prod_{k=0}^{n-2}\frac{\gamma_{k}^{\mu}\;\gamma_{k}^{\nu}}{T^{2}}\;\int|\nabla\psi^{0}-\nabla\psi^{\star}|\mathrm{d}\nu\;,$
$\displaystyle\int|\nabla\psi^{n}-\nabla\psi^{\star}|\mathrm{d}\nu^{n}\leq$
$\displaystyle\;\frac{\gamma_{\infty}^{\mu}}{\gamma_{n-1}^{\mu}}\;\prod_{k=0}^{n-1}\frac{\gamma_{k}^{\mu}\;\gamma_{k}^{\nu}}{T^{2}}\;\int|\nabla\psi^{0}-\nabla\psi^{\star}|\mathrm{d}\nu\;.$

We can now move to the pointwise convergence of the gradients stated in
Theorem 3, by relying on contractive techniques similar to the previous ones.

######  Lemma .

Assume A 1 holds. If there are positive constants $A,B>0$ such that for any
$x\in\mathbb{R}^{d}$, $|\nabla\psi^{n}-\nabla\psi^{\star}|(x)\leq A|x|+B$,
then it holds for any $x\in\mathbb{R}^{d}$,

(160)
$|\nabla\varphi^{n+1}-\nabla\varphi^{\star}|(x)\leq{T^{-1}\hat{\gamma}_{n}^{\nu}}\;(A|x|+B)\;,$

with

(161)
${\hat{\gamma}_{n}^{\nu}}\coloneqq\,\gamma_{n}^{\nu}\;\max\biggl{\\{}(T\,\alpha_{\psi^{\star}}+1)^{-1}\;,\biggl{(}1+\frac{A}{B}\frac{1+\|{\tilde{g}_{\nu}}\|_{\infty}+|\nabla\psi^{\star}(0)|}{\alpha_{\psi^{\star}}+T^{-1}}\biggr{)}\biggr{\\}}\;,$

and $\gamma_{n}^{\nu}$ is given in Section 2.2.

Similarly if for any $x\in\mathbb{R}^{d}$,
$|\nabla\varphi^{n+1}-\nabla\varphi^{\star}|(x)\leq A|x|+B$, then it holds for
any $x\in\mathbb{R}^{d}$,

(162)
$|\nabla\psi^{n+1}-\nabla\psi^{\star}|(x)\leq{T^{-1}\hat{\gamma}_{n}^{\mu}}\;(A|x|+B)\;,$

with

(163)
${\hat{\gamma}_{n}^{\mu}}\coloneqq\,\gamma_{n}^{\mu}\;\max\biggl{\\{}(T\,\alpha_{\varphi^{\star}}+1)^{-1}\;,\biggl{(}1+\frac{A}{B}\frac{1+\|{\tilde{g}_{\mu}}\|_{\infty}+|\nabla\varphi^{\star}(0)|}{\alpha_{\varphi^{\star}}+T^{-1}}\biggr{)}\biggr{\\}}\;,$

and $\gamma_{n}^{\mu}$ is given in Section 2.2.

###### Proof.

We will prove only the first inequality since the proof of the second one can
be achieved by following the same argument. Owing to the computations
performed in Section 5.1 and Section 4, let us consider (156) as starting
point of our proof here. Therefore we have

(164)
$|\nabla\varphi^{n+1}-\nabla\varphi^{\star}|(x)\leq\;T^{-1}\;\gamma_{n}^{\nu}\int\left|\nabla\psi^{n}-\nabla\psi^{\star}\right|\mathrm{d}\pi_{T}^{x,\psi^{\star}}\;,$

which combined with our assumption yields to

(165) $|\nabla\varphi^{n+1}-\nabla\varphi^{\star}|(x)\leq
T^{-1}\;\gamma_{n}^{\nu}\biggl{(}A\,\mathbb{E}_{\pi_{T}^{x,\psi^{\star}}}[|Y|]+B\biggr{)}\;.$

Therefore our proof follows once we provide a bound on the above right-hand-
side. In order to do that, let us denote by $Y^{\star}$ the strong solution of

(166)
$\begin{cases}\mathrm{d}Y^{\star}_{t}=-\biggl{(}\frac{Y_{t}^{\star}-x}{2T}+\frac{1}{2}\,\nabla\psi^{\star}(Y^{\star}_{t})\biggr{)}\,\mathrm{d}t+\mathrm{d}B_{t}\\\
Y^{\star}_{0}\sim\pi^{x,\psi^{\star}}_{T}\;.\end{cases}$

Then Ito formula, Section 3 and the boundedness of
${\tilde{g}_{\nu}}\in\tilde{\mathcal{G}}$ imply that

(167)
$\begin{split}\mathrm{d}|Y^{\star}_{t}|^{2}=&-T^{-1}|Y^{\star}_{t}|^{2}\,\mathrm{d}t+T^{-1}\langle
Y^{\star}_{t},\,x\rangle\,\mathrm{d}t-\langle
Y^{\star}_{t},\,\nabla\psi^{\star}(Y^{\star}_{t})\rangle\,\mathrm{d}t+1\,\mathrm{d}t+2\,Y^{\star}_{t}\,\cdot\mathrm{d}B_{t}\\\
\leq&-(\alpha_{\psi^{\star}}+T^{-1})|Y^{\star}_{t}|^{2}\,\mathrm{d}t+(1+T^{-1}|x|+\|{\tilde{g}_{\nu}}\|_{\infty}+|\nabla\psi^{\star}(0)|)|Y_{t}^{\star}|\,\mathrm{d}t+2\,Y^{\star}_{t}\,\cdot\mathrm{d}B_{t}\;,\end{split}$

and therefore for any $\varepsilon\in(0,\alpha_{\psi^{\star}}+T^{-1})$ we have

(168)
$\begin{split}\mathrm{d}|Y^{\star}_{t}|^{2}\leq-(\alpha_{\psi^{\star}}+T^{-1}-\varepsilon)|Y^{\star}_{t}|^{2}\,\mathrm{d}t+(4\varepsilon)^{-1}(1+T^{-1}|x|+\|{\tilde{g}_{\nu}}\|_{\infty}+|\nabla\psi^{\star}(0)|)^{2}\mathrm{d}t+2\,Y^{\star}_{t}\cdot\mathrm{d}B_{t}\;.\end{split}$

If we consider the stopping time $\tau_{M}\coloneqq\inf\\{t\geq
0\;\colon\;|Y_{t}^{\star}|>M\\}$, where we set
$\inf(\emptyset)\coloneqq+\infty$, by integrating over time on
$[0,t\wedge\tau_{M}]$, taking expectation and owing to the Optional Stopping
Theorem we deduce that

(169)
$\begin{split}\mathbb{E}[|Y_{t\wedge\tau_{M}}^{\star}|^{2}]\leq\mathbb{E}[|Y_{0}^{\star}|^{2}]-(\alpha_{\psi^{\star}}+T^{-1}-\varepsilon)\int_{0}^{t}\mathbb{E}[\mathbf{1}_{s\leq\tau_{M}}\,|Y_{s}^{\star}|^{2}]\,\mathrm{d}s\\\
+\frac{t}{4\varepsilon}(1+T^{-1}|x|+\|{\tilde{g}_{\nu}}\|_{\infty}+|\nabla\psi^{\star}(0)|)^{2}\;.\end{split}$

Since $\tau_{M}\uparrow+\infty$ almost surely as $M\uparrow+\infty$ (as a
consequence of Section 3), in the asymptotic regime Fatou Lemma implies for
any $t\geq 0$ that

(170)
$\begin{split}\mathbb{E}[|Y_{t}^{\star}|^{2}]+(\alpha_{\psi^{\star}}+T^{-1}-\varepsilon)\int_{0}^{t}\mathbb{E}[|Y_{s}^{\star}|^{2}]\,\mathrm{d}s\leq\mathbb{E}[|Y_{0}^{\star}|^{2}]+\frac{t}{4\varepsilon}(1+T^{-1}|x|+\|{\tilde{g}_{\nu}}\|_{\infty}+|\nabla\psi^{\star}(0)|)^{2}\;,\end{split}$

which combined with the stationarity of the process
$Y^{\star}_{s}\sim\pi^{x,\psi^{\star}}_{T}$, gives

(171)
$\mathbb{E}_{\pi_{T}^{x,\psi^{\star}}}[|Y|]\leq\mathbb{E}_{\pi_{T}^{x,\psi^{\star}}}[|Y|^{2}]^{1/2}\leq\frac{1+T^{-1}|x|+\|{\tilde{g}_{\nu}}\|_{\infty}+|\nabla\psi^{\star}(0)|}{\sqrt{4\varepsilon(\alpha_{\psi^{\star}}+T^{-1}-\varepsilon)}}\;.$

By minimizing over $\varepsilon\in(0,\alpha_{\psi^{\star}}+T^{-1})$ we finally
get the desired upper-bound

(172)
$\mathbb{E}_{\pi_{T}^{x,\psi^{\star}}}[|Y|]\leq\frac{1+T^{-1}|x|+\|{\tilde{g}_{\nu}}\|_{\infty}+|\nabla\psi^{\star}(0)|}{\alpha_{\psi^{\star}}+T^{-1}}=\frac{|x|}{T\,\alpha_{\psi^{\star}}+1}+\frac{1+\|{\tilde{g}_{\nu}}\|_{\infty}+|\nabla\psi^{\star}(0)|}{\alpha_{\psi^{\star}}+T^{-1}}\;.$

By combining the above bound with (165) we finally conclude that

(173) $\displaystyle|\nabla\varphi^{n+1}-\nabla\varphi^{\star}|(x)\leq
T^{-1}\;\gamma_{n}^{\nu}\biggl{(}\frac{A}{T\,\alpha_{\psi^{\star}}+1}\,|x|+A\frac{1+\|{\tilde{g}_{\nu}}\|_{\infty}+|\nabla\psi^{\star}(0)|}{\alpha_{\psi^{\star}}+T^{-1}}+B\biggr{)}\leq{T^{-1}\hat{\gamma}^{\nu}_{n}}(A\,|x|+B)$

where the last step holds true for the chosen rate
${T^{-1}\hat{\gamma}_{n}^{\nu}}$. This concludes the proof of the first part
of the statement. The proof of the second one is similar and for this reason
we omit it. Let us just mention here that the same reasoning yields to the
moment bound

(174)
$\mathbb{E}_{\pi_{T}^{x,\varphi^{\star}}}[|Y|]<\frac{|x|}{T\,\alpha_{\varphi^{\star}}+1}+\frac{1+\|{\tilde{g}_{\mu}}\|_{\infty}+|\nabla\varphi^{\star}(0)|}{\alpha_{\varphi^{\star}}+T^{-1}}\;.$

∎

###### Proof of Theorem 3.

It is enough concatenating the bounds proven in Section 5.1 and observing that
at each step $\psi^{k}$ and $\varphi^{k}$ satisfy (43) with constants
$A_{k},\,B_{k}$ such that the ratio $A^{k}/B^{k}\equiv A/B$. ∎

###  5.2. Exponential pointwise and entropic convergence of Sinkhorn’s
algorithm

Let us firstly show how the pointwise convergence of the gradients implies the
(pointwise and integrated) convergence of Sinkhorn iterates, i.e., Theorem 4.

###### Proof of Theorem 4.

Owing to the normalizations (54), (55) and to Theorem 3, we immediately deduce
that

(175) $\begin{split}|\varphi^{\diamond
n}-&\varphi^{\star}|(x)=\biggl{|}\varphi^{\diamond n}(x)-\int\varphi^{\diamond
n}\mathrm{d}\mu-\varphi^{\star}(x)+\int\varphi^{\star}\mathrm{d}\mu\biggr{|}=\biggl{|}\int(\varphi^{n}-\varphi^{\star})(x)-(\varphi^{n}-\varphi^{\star})(y)\;\mathrm{d}\mu(y)\biggr{|}\\\
\leq&\int\biggl{|}(\varphi^{n}-\varphi^{\star})(x)-(\varphi^{n}-\varphi^{\star})(y)\biggr{|}\;\mathrm{d}\mu(y)\leq\int\int_{0}^{1}|\nabla(\varphi^{n}-\varphi^{\star})|(y+t(x-y))\;|x-y|\;\mathrm{d}t\;\mathrm{d}\mu(y)\\\
\overset{\eqref{eq:pointwisebound}}{\leq}&\frac{T}{\hat{\gamma}^{\mu}_{n-1}}\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\;\hat{\gamma}_{k}^{\nu}}{T^{2}}\int\int_{0}^{1}(A|y+t(x-y)|+B)|x-y|\;\mathrm{d}t\;\mathrm{d}\mu(y)\\\
\leq&\frac{T}{\hat{\gamma}^{\mu}_{n-1}}\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\;{\hat{\gamma}_{k}^{\nu}}}{T^{2}}\;\biggl{[}A\;|x|^{2}+(A\;M_{1}(\mu)+B)\;|x|+B\;M_{1}(\mu)+2\,A\;M_{2}(\mu)\biggr{]}\;.\end{split}$

The second pointwise bound can be proven in the same manner. ∎

As a consequence of Theorem 4 we may also deduce the convergence of the
$\mathrm{L}^{1}$-norms along the adjusted marginals and along the real
marginals.

######  Corollary .

Assume A 1 holds and (43) for some positive constants $A,\,B>0$. Then for any
$n\geq 1$ it holds

(176) $\begin{split}\|\varphi^{\diamond
n}-\varphi^{\star}\|_{\mathrm{L}^{1}(\mu)}\leq&\;\frac{T}{\hat{\gamma}^{\mu}_{n-1}}\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\;\hat{\gamma}_{k}^{\nu}}{T^{2}}\;\biggl{[}3\,A\;M_{2}(\mu)+(A\;M_{1}(\mu)+B)\,M_{1}(\mu)+B\;M_{1}(\mu)\biggr{]}\;,\\\
\|\psi^{\diamond
n}-\psi^{\star}\|_{\mathrm{L}^{1}(\nu)}\leq&\;\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\;\hat{\gamma}_{k}^{\nu}}{T^{2}}\;\biggl{[}3\,A\;M_{2}(\nu)+(A\;M_{1}(\nu)+B)\,M_{1}(\nu)+B\;M_{1}(\nu)\biggr{]}\;,\end{split}$

and

(177) $\|\varphi^{\diamond
n}-\varphi^{\star}\|_{\mathrm{L}^{1}(\mu^{n})}\leq\;\;C(A,B,\mu)\frac{T}{\hat{\gamma}^{\mu}_{n-1}}\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\;\hat{\gamma}_{k}^{\nu}}{T^{2}}\;,\quad\|\psi^{\diamond
n}-\psi^{\star}\|_{\mathrm{L}^{1}(\nu^{n})}\leq\;C(A,B,\nu)\;\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\;\hat{\gamma}_{k}^{\nu}}{T^{2}}\;,$

where

(178)
$\begin{split}C(A,B,\mu)\coloneqq&\,\biggl{[}3\,A\;M_{2}(\mu)+(A\;M_{1}(\mu)+B)\,M_{1}(\mu)+B\;M_{1}(\mu)\biggr{]}\\\
+&A\,C_{2}(\mu)\biggl{(}\sqrt{\mathcal{H}(\mu^{1}|\mu)}+\frac{\mathcal{H}(\mu^{1}|\mu)}{2}\biggr{)}+(A\,M_{1}(\mu)+B)\,C_{1}(\mu)\,\sqrt{\mathcal{H}(\mu^{1}|\mu)}\end{split}$

and

(179)
$\begin{split}C(A,B,\nu)\coloneqq&\,\biggl{[}3\,A\;M_{2}(\nu)+(A\;M_{1}(\nu)+B)\,M_{1}(\nu)+B\;M_{1}(\nu)\biggr{]}\\\
+&A\,C_{2}(\nu)\biggl{(}\sqrt{\mathcal{H}(\nu^{0}|\nu)}+\frac{\mathcal{H}(\nu^{0}|\nu)}{2}\biggr{)}+(A\,M_{1}(\nu)+B)\,C_{1}(\nu)\,\sqrt{\mathcal{H}(\nu^{0}|\nu)}\;,\end{split}$

with $C_{1}(\zeta)$ and $C_{2}(\zeta)$ are positive constants defined in
Section 5.2 for $\zeta\in\\{\mu,\nu\\}$.

###### Proof.

The proof of the integrated bounds along the marginals (176) is a
straightforward consequence of the pointwise convergence, whereas the bounds
along the adjusted marginals is a consequence of the weighted Csiszár-
Kullback-Pinsker inequalities [8, Theorem 2.1] which imply

(180) $\begin{split}M_{1}&(\mu^{n})\leq
M_{1}(\mu)+C_{1}(\mu)\,\sqrt{\mathcal{H}(\mu^{n}|\mu)}\quad\text{ and }\quad
M_{2}(\mu^{n})\leq
M_{2}(\mu)+C_{2}(\mu)\biggl{(}\sqrt{\mathcal{H}(\mu^{n}|\mu)}+\frac{\mathcal{H}(\mu^{n}|\mu)}{2}\biggr{)}\\\
M_{1}&(\nu^{n})\leq
M_{1}(\nu)+C_{1}(\nu)\,\sqrt{\mathcal{H}(\nu^{n}|\nu)}\quad\text{ and }\quad
M_{2}(\nu^{n})\leq
M_{2}(\nu)+C_{2}(\nu)\biggl{(}\sqrt{\mathcal{H}(\nu^{n}|\nu)}+\frac{\mathcal{H}(\nu^{n}|\nu)}{2}\biggr{)},\end{split}$

where $C_{1}(\mu),\,C_{1}(\nu),\,C_{2}(\mu),\,C_{2}(\nu)$ are positive
constants (independent from $n\in\mathbb{N}$). For sake of clarity we
postponed the proof of the above moment bounds to Section 5.2, below. The
proof of (177) then follows from the fact that the sequences
$(\mathscr{H}(\mu^{n}|\mu))_{n\in\mathbb{N}}$ and
$(\mathscr{H}(\nu^{n}|\nu))_{n\in\mathbb{N}}$ are monotone decreasing along
Sinkhorn [47, Proposition 6.10]. ∎

In the next lemma we show how the weighted Csiszár-Kullback-Pinsker
inequalities imply the above moments inequalities along the adjusted
marginals.

######  Lemma .

Assume A 1 and let $\zeta\in\\{\mu,\nu\\}$. Then for any probability measure
$\mathfrak{p}\in\mathcal{P}(\mathbb{R}^{d})$ such that
$\mathcal{H}(\mathfrak{p}|\zeta)<+\infty$ it holds

(181) $\begin{split}&M_{1}(\mathfrak{p})\leq
M_{1}(\zeta)+C_{1}(\zeta)\,\sqrt{\mathcal{H}(\mathfrak{p}|\zeta)}\quad\text{
and }\quad M_{2}(\mathfrak{p})\leq
M_{2}(\zeta)+C_{2}(\zeta)\biggl{(}\sqrt{\mathcal{H}(\mathfrak{p}|\zeta)}+\frac{\mathcal{H}(\mathfrak{p}|\zeta)}{2}\biggr{)}\\\
&C_{1}(\zeta)\coloneqq\inf_{\sigma_{\zeta}\in(0,\alpha_{\zeta}/2)}\biggl{(}\frac{2}{\sigma_{\zeta}}+\frac{2}{\sigma_{\zeta}}\log\int
e^{\sigma_{\zeta}|x|^{2}}\,\mathrm{d}\zeta\biggr{)}^{1/2}\;,\;C_{2}(\zeta)\coloneqq\inf_{\sigma_{\zeta}\in(0,\alpha_{\zeta}/2)}\biggl{(}\frac{3}{\sigma_{\zeta}}+\frac{2}{\sigma_{\zeta}}\int
e^{\sigma_{\zeta}|x|^{2}}\mathrm{d}\zeta\biggr{)}\;.\end{split}$

###### Proof.

For any $\sigma_{\nu}\in(0,\alpha_{\nu}/2)$, from the weighted Csiszár-
Kullback-Pinsker inequalities [8, Theorem 2.1] applied to the measurable
functions $F_{1}(x)=\sigma_{\nu}^{1/2}|x|$ and
$F_{2}(x)=\sigma_{\nu}|x|^{2}/2$ we immediately deduce

(182)
$\begin{split}M_{1}(\mathfrak{p})&=M_{1}(\nu)+\sigma_{\nu}^{-1/2}\,\int\sigma_{\nu}^{1/2}\,|x|\,\mathrm{d}(\mathfrak{p}-\nu)\leq
M_{1}(\nu)+\sqrt{\mathcal{H}(\mathfrak{p}|\nu)}\,\biggl{(}\frac{2}{\sigma_{\nu}}+\frac{2}{\sigma_{\nu}}\log\int
e^{\sigma_{\nu}|x|^{2}}\,\mathrm{d}\nu\biggr{)}^{1/2}\\\
M_{2}(\mathfrak{p})&=M_{2}(\nu)+\frac{2}{\sigma_{\nu}}\,\int\frac{\sigma_{\nu}}{2}\,|x|^{2}\,\mathrm{d}(\mathfrak{p}-\nu)\leq
M_{2}(\nu)+\biggl{(}\sqrt{\mathcal{H}(\mathfrak{p}|\nu)}+\frac{\mathcal{H}(\mathfrak{p}|\nu)}{2}\biggr{)}\biggl{(}\frac{3}{\sigma_{\nu}}+\frac{2}{\sigma_{\nu}}\int
e^{\sigma_{\nu}|x|^{2}}\mathrm{d}\nu\biggr{)}\end{split}$

which are finite by Appendix C and A 1. Minimizing over
$\sigma_{\nu}\in(0,\alpha_{\nu}/2)$ concludes the proof of the first claim.

Lastly, the moment bounds corresponding to the choice of reference $\mu$, can
be proven in the same way. ∎

Finally, let us conclude the section with the proof of the exponential
convergence of Sinkhorn’s algorithm on the primal side.

###### Proof of Theorem 5.

Let us preliminary point out that as a first consequence of Section 5.2, for
any $n\geq 1$ it holds

(183) $\varphi^{\diamond
n}-\varphi^{\star}\in\mathrm{L}^{1}(\mu)\cap\mathrm{L}^{1}(\mu^{n})\quad\text{
and }\quad\psi^{\diamond
n}-\psi^{\star}\in\mathrm{L}^{1}(\nu)\cap\mathrm{L}^{1}(\nu^{n})\;,$

which will guarantee that the following integrals (and corresponding
summations) are all well-defined.

Now, (1) and (9) imply that

(184)
$\begin{split}\log\frac{\mathrm{d}\pi^{\star}}{\mathrm{d}\pi^{n,n}}(x,y)=&\,\varphi^{n}(x)-\varphi^{\star}(x)+\psi^{n}(y)-\psi^{\star}(y)\;,\end{split}$

and hence the symmetric relative entropies can be rewritten as

(185)
$\begin{split}\mathscr{H}(\pi^{n,n}|\pi^{\star})+\mathscr{H}(\pi^{\star}|\pi^{n,n})=&\int(\varphi^{n}-\varphi^{\star})\oplus(\psi^{n}-\psi^{\star})\;\mathrm{d}\pi^{\star}-\int(\varphi^{n}-\varphi^{\star})\oplus(\psi^{n}-\psi^{\star})\;\mathrm{d}\pi^{n,n}\\\
=&\int(\varphi^{\diamond n}-\varphi^{\star})\oplus(\psi^{\diamond
n}-\psi^{\star})\;\mathrm{d}\pi^{\star}-\int(\varphi^{\diamond
n}-\varphi^{\star})\oplus(\psi^{\diamond
n}-\psi^{\star})\;\mathrm{d}\pi^{n,n}\\\ =&\int(\varphi^{\diamond
n}-\varphi^{\star})\;\mathrm{d}\mu^{n}-\int(\varphi^{\diamond
n}-\varphi^{\star})\;\mathrm{d}\mu\;.\end{split}$

By combining the above with Section 5.2, we then deduce that

(186)
$\begin{split}\mathscr{H}(&\pi^{n,n}|\pi^{\star})+\mathscr{H}(\pi^{\star}|\pi^{n,n})\leq\;\frac{D(A,B,\mu)}{{T^{-1}\hat{\gamma}^{\mu}_{n-1}}}\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\;{\hat{\gamma}_{k}^{\nu}}}{T^{2}}\;,\end{split}$

with

(187)
$\begin{split}D(A,B,\mu)\coloneqq&\,C(A,B,\mu)+3\,A\;M_{2}(\mu)+(A\;M_{1}(\mu)+B)\,M_{1}(\mu)+B\;M_{1}(\mu)\\\
=&\,2\biggl{[}3\,A\;M_{2}(\mu)+(A\;M_{1}(\mu)+B)\,M_{1}(\mu)+B\;M_{1}(\mu)\biggr{]}\\\
+&A\,C_{2}(\mu)\biggl{(}\sqrt{\mathcal{H}(\mu^{1}|\mu)}+\frac{\mathcal{H}(\mu^{1}|\mu)}{2}\biggr{)}+(A\,M_{1}(\mu)+B)\,C_{1}(\mu)\,\sqrt{\mathcal{H}(\mu^{1}|\mu)}\;,\end{split}$

with $C_{1}(\mu),\,C_{2}(\mu)$ being the constants introduced in Section 5.2.

Similarly (1) and (8) imply

(188)
$\log\frac{\mathrm{d}\pi^{\star}}{\mathrm{d}\pi^{n+1,n}}(x,y)=\varphi^{n+1}(x)-\varphi^{\star}(x)+\psi^{n}(y)-\psi^{\star}(y)\;,$

hence

(189)
$\mathscr{H}(\pi^{n+1,n}|\pi^{\star})+\mathscr{H}(\pi^{\star}|\pi^{n+1,n})=\int(\psi^{\diamond
n}-\psi^{\star})\;\mathrm{d}\nu^{n}-\int(\psi^{\diamond
n}-\psi^{\star})\;\mathrm{d}\nu\;.$

The latter combined with Section 5.2 yields to

(190)
$\mathscr{H}(\pi^{n+1,n}|\pi^{\star})+\mathscr{H}(\pi^{\star}|\pi^{n+1,n})\leq
D(A,B,\nu)\;\prod_{k=0}^{n-1}\frac{\hat{\gamma}_{k}^{\mu}\;{\hat{\gamma}_{k}^{\nu}}}{T^{2}}\;,$

with

(191)
$\begin{split}D(A,B,\nu)\coloneqq&\,C(A,B,\nu)+3\,A\;M_{2}(\nu)+(A\;M_{1}(\nu)+B)\,M_{1}(\nu)+B\;M_{1}(\nu)\\\
=&\,2\biggl{[}3\,A\;M_{2}(\nu)+(A\;M_{1}(\nu)+B)\,M_{1}(\nu)+B\;M_{1}(\nu)\biggr{]}\\\
+&A\,C_{2}(\nu)\biggl{(}\sqrt{\mathcal{H}(\nu^{0}|\nu)}+\frac{\mathcal{H}(\nu^{0}|\nu)}{2}\biggr{)}+(A\,M_{1}(\nu)+B)\,C_{1}(\nu)\,\sqrt{\mathcal{H}(\nu^{0}|\nu)}\;.\end{split}$

Finally let us notice that from (8) and (9) we may also deduce that

(192)
$\begin{split}\mathscr{H}(\pi^{n,n}|\pi^{n,n-1})=\mathscr{H}(\mu^{n}|\mu)\quad\text{
and }\quad\mathscr{H}(\pi^{n,n-1}|\pi^{n,n})=\mathscr{H}(\mu|\mu^{n})\;,\\\
\mathscr{H}(\pi^{n+1,n}|\pi^{n,n})=\mathscr{H}(\nu^{n}|\nu)\quad\text{ and
}\quad\mathscr{H}(\pi^{n,n}|\pi^{n+1,n})=\mathscr{H}(\nu|\nu^{n})\;.\end{split}$

and hence the bound for the adjusted marginals is a consequence of the
previous ones, and the trivial inequalities

(193)
$\begin{split}&\mathscr{H}(\pi^{n,n}|\pi^{n,n-1})=\mathscr{H}(\mu^{n}|\mu)\leq\mathscr{H}(\pi^{n,n}|\pi^{\star})\;,\quad\mathscr{H}(\pi^{n,n-1}|\pi^{n,n})=\mathscr{H}(\mu|\mu^{n})\leq\mathscr{H}(\pi^{\star}|\pi^{n,n})\;,\\\
&\mathscr{H}(\pi^{n+1,n}|\pi^{n,n})=\mathscr{H}(\nu^{n}|\nu)\leq\mathscr{H}(\pi^{n+1,n}|\pi^{\star})\;,\quad\mathscr{H}(\pi^{n,n}|\pi^{n+1,n})=\mathscr{H}(\nu|\nu^{n})\leq\mathscr{H}(\pi^{\star}|\pi^{n+1,n})\,.\end{split}$

∎

######  Remark .

Let us remark here that from Theorem 5 we can consider a sharper
multiplicative constant $C_{S}(A,B,\mu)$ in Section 5.2. Indeed, instead of
relying on the monotonicity of relative entropies along Sinkhorn’s algorithm,
in (178) we could define $C_{S}(A,B,\mu)$ as

(194)
$C_{S}(A,B,\mu)\coloneqq\biggl{[}3\,A\;M_{2}(\mu)+(A\;M_{1}(\mu)+B)\,M_{1}(\mu)+B\;M_{1}(\mu)\biggr{]}+\varepsilon_{\mu}(n)\;,$

where

(195) $\varepsilon_{\mu}(n)\coloneqq
A\,C_{2}(\mu)\biggl{(}\sqrt{\mathcal{H}(\mu^{n}|\mu)}+\frac{\mathcal{H}(\mu^{n}|\mu)}{2}\biggr{)}+(A\,M_{1}(\mu)+B)\,C_{1}(\mu)\,\sqrt{\mathcal{H}(\mu^{n}|\mu)}$

is a positive constant exponentially small as $n\uparrow+\infty$ thanks to
Theorem 5. Then, in Theorem 5 instead of (187), we can consider the sharper
multiplicative constant

(196) $D_{S}(A,B,\mu)\coloneqq
2\,\biggl{[}3\,A\;M_{2}(\mu)+(A\;M_{1}(\mu)+B)\,M_{1}(\mu)+B\;M_{1}(\mu)\biggr{]}+\varepsilon_{\mu}(n)\;.$

The same reasoning applies for $C_{S}(A,B,\nu)$ and $D_{S}(A,B,\nu)$.

###  5.3. Exponential convergence of the Hessians along Sinkhorn’s algorithm

For notations sake for any couple of vectors $v,\,w\in\mathbb{R}^{d}$ we will
denote by $v\otimes w=v\,w^{\sf T}$ the matrix given by their tensor product.
Then, for any $v,\,w\in\mathbb{R}^{d}$ the Frobenius norm of their product
reads as

(197) $\|v\otimes w\|_{\mathrm{F}}=|v|\,|w|\;.$

Let us also recall here that for
$h\in\\{\varphi^{n},\,\psi^{n},\,\varphi^{\star},\,\psi^{\star}\\}$ the
gradient and the Hessian along the semigroup has the explicit formulation (cf.
(24) and (61) respectively)

(198) $\begin{split}\nabla\log
P_{T}\exp(-h)(x)=&\,\frac{1}{T}\int(y-x)\,\pi_{T}^{x,h}(\mathrm{d}y)\\\
\nabla^{2}\log
P_{T}\exp(-h)(x)=&\,-T^{-1}\operatorname{Id}+T^{-2}\,\mathrm{Cov}(\pi_{T}^{x,h})\;.\end{split}$

where we recall the conditional distribution $\pi_{T}^{x,h}$ being defined at
(23)

(199)
$\pi_{T}^{x,h}(\mathrm{d}y)\propto\exp\biggl{(}-\frac{|y-x|^{2}}{2T}-h(y)\biggr{)}\mathrm{d}y\;.$

Throught this section we will always assume the validity of the hypothesis of
Theorem 3, i.e., A 1, and (43) are met.

######  Lemma .

There exist two positive constants $C_{1},\,C_{2}>0$ independent of $x$ and
$T$, such that for any coupling
$\pi\in\Pi(\pi^{x,\psi^{\star}}_{T},\,\pi^{x,\psi^{n}}_{T})$ it holds

(200)
$\|\nabla^{2}\varphi^{n+1}-\nabla^{2}\varphi^{\star}\|_{\mathrm{F}}(x)\leq
T^{-2}\,\mathbb{E}_{(Y,Z)\sim\pi}\biggl{[}|Y-Z|\,(|Y|+|Z|+C_{1}\,T^{-1}\,|x|+C_{2})\biggr{]}\;.$

###### Proof.

As a byproduct of (4) and (61) we immediately have for any
$x\in\mathbb{R}^{d}$ and for any coupling
$\pi\in\Pi(\pi_{T}^{x,\psi^{\star}},\pi_{T}^{x,\psi^{n}})$

(201)
$\begin{split}T^{2}\,(\nabla^{2}\varphi^{n+1}-\nabla^{2}\varphi^{\star})(x)=&\;\mathrm{Cov}(\pi_{T}^{x,\psi^{\star}})-\mathrm{Cov}(\pi_{T}^{x,\psi^{n}})\\\
=&\,\mathbb{E}_{(Y,Z)\sim\pi}\biggl{[}Y^{\otimes 2}-Z^{\otimes
2}\biggr{]}+\mathbb{E}_{Z\sim\pi_{T}^{x,\psi^{n}}}[Z]^{\otimes
2}-\mathbb{E}_{Y\sim\pi_{T}^{x,\psi^{\star}}}[Y]^{\otimes 2}\\\
=&\,\mathbb{E}_{(Y,Z)\sim\pi}\biggl{[}(Y-Z)\otimes
Y\biggr{]}+\mathbb{E}_{(Y,Z)\sim\pi}\biggl{[}Z\otimes(Y-Z)\biggr{]}\\\
-&\,\mathbb{E}_{(Y,Z)\sim\pi}\biggl{[}Y-Z\biggr{]}\;\mathbb{E}_{Z\sim\pi_{T}^{x,\psi^{n}}}[Z]-\mathbb{E}_{Y\sim\pi_{T}^{x,\psi^{\star}}}[Y]\;\mathbb{E}_{(Y,Z)\sim\pi}\biggl{[}Y-Z\biggr{]}\end{split}$

By applying the Frobenius norm, and recalling (197), we then deduce that

(202)
$T^{2}\,\|\nabla^{2}\varphi^{n+1}-\nabla^{2}\varphi^{\star}\|_{\mathrm{F}}(x)\leq\mathbb{E}_{(Y,Z)\sim\pi}\biggl{[}|Y-Z|\;\biggl{(}|Y|+|Z|+\mathbb{E}_{Z\sim\pi_{T}^{x,\psi^{n}}}[|Z|]+\mathbb{E}_{Y\sim\pi_{T}^{x,\psi^{\star}}}[|Y|]\biggr{)}\biggr{]}$

In order to bound the last two expected values in the right hand side we will
proceed as in the proof of Section 5.1. Particularly, (172) already proves
that

(203)
$\mathbb{E}_{Y\sim\pi_{T}^{x,\psi^{\star}}}[|Y|]<\frac{|x|}{T\,\alpha_{\psi^{\star}}+1}+\frac{1+\|{\tilde{g}_{\nu}}\|_{\infty}+|\nabla\psi^{\star}(0)|}{\alpha_{\psi^{\star}}+T^{-1}}\;.$

By reasoning in the same way we can prove that

(204)
$\mathbb{E}_{Z\sim\pi_{T}^{x,\psi^{n}}}[|Z|]<\frac{|x|}{T\,\alpha_{\nu,n}+1}+\frac{1+\|{\tilde{g}_{\nu}}\|_{\infty}+|\nabla\psi^{n}(0)|}{\alpha_{\psi^{n}}+T^{-1}}\;,$

Particularly, the pointwise convergence of the gradients of Theorem 3 and the
convergences $\alpha_{\mu,n}\uparrow\alpha_{\varphi^{\star}}$,
$\alpha_{\nu,n}\uparrow\alpha_{\psi^{\star}}$ stated in Theorem 7 yield to the
uniform bound

(205)
$\mathbb{E}_{Y\sim\pi_{T}^{x,\psi^{\star}}}[|Y|]\;\vee\;\sup_{n\in\mathbb{N}}\mathbb{E}_{Z\sim\pi_{T}^{x,\psi^{n}}}[|Z|]\;\leq
C_{1}\,T^{-1}\,|x|+C_{2}\;,$

for some positive constants $C_{1},\,C_{2}>0$ independent of $x$ and $T$. This
concludes our proof. ∎