Entropic Regularization: Convergence

This chapter focuses on algorithmic convergence for entropic optimal transport: the marginals and the temperature are fixed, and one asks how Sinkhorn iterates, soft transforms and related monotone scaling procedures approach their regularized fixed point. Statistical convergence is a separate question, because the empirical marginals then change with the sample size; this is the topic of Paragraph.

The chapter revisits Sinkhorn convergence through several complementary lenses. Bregman projections explain the alternating-projection geometry, Fortet’s order argument gives qualitative fixed-point convergence, and an order-theoretic M-function viewpoint explains why Sinkhorn-like scaling can still converge even when the equations are no longer the gradient of a convex potential. Robust Bregman estimates then give a non-asymptotic dual-gap bound, while Hilbert’s metric gives a clean linear contraction when the kernel is uniformly positive. The last sections discuss Gaussian closed forms and continuous $\varepsilon$-Sinkhorn flows as model cases where the fixed-point structure becomes explicit.

from pathlib import Path
import sys

from IPython.display import Image as DisplayImage
from IPython.display import display

here = Path.cwd()
myst_dir = None
for candidate in [here, here.parent, here / "myst", here.parent / "myst", here.parent.parent / "myst"]:
    if (candidate / "ot4ml_web.py").exists():
        myst_dir = candidate.resolve()
        sys.path.insert(0, str(myst_dir))
        break

if myst_dir is None:
    raise RuntimeError("Could not locate myst/ot4ml_web.py")

repo_root = myst_dir.parent
thumbnails = repo_root / "notebooks-figures" / "thumbnails"

def show_book_figure(name, width=760):
    display(DisplayImage(filename=str(thumbnails / f"{name}.png"), width=width))

Sinkhorn Convergence: Bregman Point of View¶

Sinkhorn can be read as alternating Bregman projections. The main geometric message is simple: each row or column rescaling is the KL projection onto one affine marginal constraint. The convergence mechanism then follows from the Pythagorean identity for Bregman divergences.

For simplicity, this section is written for discrete measures. The same ideas carry over to general measures, and the robust-rate section below expresses the constants through cost and potential oscillations rather than through the number of grid points.

Alternating KL Projections¶

The projection viewpoint explains Sinkhorn as repeated enforcement of one marginal constraint at a time. It is not specific to entropy, although KL is the case where the projections reduce to elementary row and column scalings.

The following matrix construction is the finite-dimensional counterpart of the measure-valued Definition: Measure Bregman Divergence.

For the quadratic generator $\Phi(\P)=\frac12\|\P\|_{\mathrm F}^2$ on $\mathbb R^{n\times m}$, one recovers half the squared Euclidean distance, $B_\Phi(\P\mid Q)=\frac12\|\P-Q\|_{\mathrm F}^2$. For the negative entropy $\Phi(\P)=\sum_{i,j}\P_{i,j}\log \P_{i,j}$ on the nonnegative orthant, extended by $+\infty$ outside it, one obtains $B_\Phi(\P\mid Q)=\operatorname{KL}(\P\mid Q)$. For $\P\geq0$ and $Q>0$, this identity is understood through the lower-semicontinuous extension, with $0\log0=0$.

Bregman divergences are useful because their geometry can encode constraints. A Legendre-type generator blows up, or has an infinite derivative, at the boundary of its domain. For negative entropy, positivity is therefore built into the divergence, so one projects onto affine marginal constraints without separately handling non-negativity.

Linear Tilts and Gibbs References¶

Adding a linear cost to a Bregman penalty merely shifts the reference point in dual coordinates. The usual Gibbs--KL reformulation is exactly the entropy specialization.

For the negative entropy, take $Q=a\otimes b$. The tilted reference is

Thus

On the transport polytope, scaling $K_{a,b}^\epsilon$ is equivalent to scaling the Gibbs kernel $K=e^{-C/\epsilon}$ because the factors $a_i$ and $b_j$ can be absorbed into the Sinkhorn scalings. The unique entropic optimizer is the KL projection of this tilted Gibbs reference onto the coupling constraints:

Cyclic Projection Convergence¶

Given two closed convex constraint sets $\mathcal C_1$ and $\mathcal C_2$, the cyclic method projects first onto $\mathcal C_1$, then onto $\mathcal C_2$, and repeats. Full iterates satisfy the second constraint and half-iterates satisfy the first; both are attained together only in the limit.

The following algorithm records the general cyclic iteration. The two nonnegative constraint defects vanish on their respective sets and provide a practical stopping test.

The convergence mechanism is the classical one of Bregman projections Bregman, 1967. General convex constraints determine a feasible limit, while affine constraints give the stronger closest-point characterization.

Row and Column Scalings¶

We now apply the general Bregman-projection framework to entropic OT. The two affine constraints impose the source and target marginals, while negative entropy turns their cyclic projections into explicit row and column rescalings. Denote these constraint sets by

Positivity is a separate constraint, so

For negative entropy, however, the Legendre generator has effective domain $\Omega=\mathbb R_+^{n\times m}$ and is $+\infty$ outside $\Omega$. Thus $\operatorname{Proj}_{\mathcal C_i}^{\operatorname{KL}}$ already minimizes over $\mathcal C_i\cap\Omega$: the generator enforces positivity, which need not be repeated in the projector notation. The cyclic iteration therefore specializes to

These two projectors are explicit: they rescale respectively the rows and the columns.

For positive histograms and a strictly positive Gibbs kernel, iterative proportional fitting keeps the half-steps positive and converges Ruschendorf, 1995Rüschendorf & Thomsen, 1998. Since the two marginal sets are affine, Proposition: Convergence of Cyclic Bregman Projections, applied with $\P^{(0)}=K_{a,b}^\epsilon$, identifies the limit as the entropic optimizer.

Defining

the two projection steps are the usual Sinkhorn updates on the scaling vectors. In practice one stores the vectors and multiplies by the Gibbs kernel, often exploiting separable, sparse, low-rank or geometric structure.

The Bregman proof is geometric, but its direct finite-dimensional linear-rate constants can degrade with dimension and with small $\epsilon$. The robust dual analysis below gives a dimension-free qualitative message: before any local linear regime becomes visible, one can still guarantee an $O(1/\ell)$ dual gap whose constants depend on the cost range and potential oscillation.

Other Divergences¶

The simplicity of the KL construction relies on negative entropy encoding nonnegativity through its effective domain. Suppose instead that $\Phi$ is Legendre on the full space $\mathbb R^{n\times m}$. Positivity must then be included explicitly in the marginal constraints:

Thus $\mathcal U(a,b)=\mathcal C_{a,+}^1\cap\mathcal C_{b,+}^2$. These sets are convex but no longer affine, and their $B_\Phi$-projectors generally have no closed-form scaling formula. The affine clause of Proposition: Convergence of Cyclic Bregman Projections no longer applies: ordinary cyclic projections may converge to a feasible point without computing $\operatorname{Proj}_{\mathcal U(a,b)}^{B_\Phi}(\P^{(0)})$.

For the Euclidean generator $\Phi(\P)=\frac12\|\P\|_{\mathrm F}^2$, projection onto $\mathcal C_{a,+}^1$ is a rowwise simplex projection:

The column formula is analogous. These thresholds are not multiplicative scalings, but sorting or selection computes them efficiently Blondel et al., 2018.

To recover the closest-point projection onto the intersection, Dykstra’s algorithm alternates the same Bregman projectors while carrying correction variables in dual coordinates. Under the standard finite-dimensional assumptions, it converges to $\operatorname{Proj}_{\mathcal C_{a,+}^1\cap\mathcal C_{b,+}^2}^{B_\Phi} (\P^{(0)})$ Dykstra, 1985Censor & Reich, 1998Bauschke & Lewis, 2000. Benamou, Carlier, Cuturi, Nenna, and Peyré developed this construction systematically for regularized transport Benamou et al., 2015. Through Fenchel duality, Bregman--Dykstra is block-coordinate ascent on the marginal dual potentials; the correction variables store the dual memory absent from plain cyclic projection Bregman et al., 1999.

For the quadratic generator and the product reference $\xi=a\otimes b$,

This is closely related, but not identical, to the quadratic $\phi$-divergence generated by $\phi(r)=\frac12(r-1)^2$:

for positive $a_i,b_j$. The difference is the fixed inverse product-marginal weighting in the second quadratic norm. Both models yield positive-part thresholding and potentially sparse plans. The section on Other Convex Regularizers develops their dual laws and alternating dual-maximization algorithms.

Sinkhorn Convergence: Monotone Point of View¶

This section isolates the order structure shared by generalized Sinkhorn updates. It first introduces the abstract language of topical maps, then applies it to the generalized soft transforms of Other Convex Regularizers, and finally derives a monotone fixed-point argument valid beyond KL regularization.

Variation Seminorm and Topical Maps¶

Potentials are defined only up to additive constants, so their natural gauge-invariant size is their oscillation.

The order maps relevant to Sinkhorn commute with this additive gauge.

Generalized Sinkhorn Maps¶

For the $\phi$-divergence regularized transport problem, the one-sided block updates are expressed using the Legendre transform $\phi^*$ defined in (47). They are the generalized soft $c$-transforms of (121), and their alternating dual-maximization scheme is (124). After choosing a consistent extremal minimizer whenever a scalar update is not unique, one full cycle acting on the first potential is

For the KL generator, this is the usual double Sinkhorn map.

Topicality gives nonexpansiveness, not a strict contraction. For KL, positivity of the Gibbs kernel supplies the stronger Hilbert-metric contraction proved in Sinkhorn Convergence: Linear Hilbert Metric Rate; no comparable strict factor follows from the order properties alone for a general $\phi$.

Monotone Convergence for Generalized Sinkhorn¶

The following order argument traces back to Fortet’s proof of the Schrodinger system Fortet, 1940Essid & Pavon, 2019Léonard, 2019. It is not based on an optimization principle: a subsolution generates an increasing orbit, while a fixed point provides an upper barrier.

The subsolution and supersolution conditions are invariant under additive shifts, but a shift cannot turn an arbitrary initialization into either one. The proposition therefore isolates the genuinely order-theoretic convergence mechanism; unrestricted initialization requires an additional compactness, ascent, or contraction argument.

Sinkhorn Convergence: Sublinear Robust Rate¶

The preceding projection and monotonicity arguments establish convergence but do not provide a quantitative rate. The Hilbert-metric analysis in Sinkhorn Convergence: Linear Hilbert Metric Rate will give geometric, or linear, convergence, the canonical asymptotic behavior of a strictly contractive fixed-point iteration; however, the gap between its global contraction factor and one typically becomes exponentially small in the cost oscillation divided by $\epsilon$. We therefore first establish a slower $O(1/\ell)$ dual-gap rate whose constant grows only as $1/\epsilon$ Peyré, 2026Altschuler et al., 2017Dvurechensky et al., 2018. This robust dependence is crucial when Sinkhorn approximates unregularized discrete OT: Corollary Corollary: Approximating Unregularized OT by Regularized Dual Costs chooses $\epsilon=\delta/(2\log(nm))$, so the temperature decreases with the requested accuracy and with the number of support points.

The same identities give computable diagnostics. If $r=P\mathbf1_m$, the KL projection onto the row constraint satisfies $\operatorname{KL}(\operatorname{Proj}(P)\mid P)=\operatorname{KL}(a\mid r)$; the column identity is analogous. Each observed dual ascent is therefore a marginal KL defect, and Theorem Theorem: Pinsker Inequality turns it into an $\ell^1$ residual. The residual is not itself the remaining dual gap: a certificate also needs the quotient-radius estimate used above.

Sinkhorn Convergence: Linear Hilbert Metric Rate¶

Hilbert’s projective metric measures positive scaling vectors modulo global multiplication. A positive kernel contracts this geometry, yielding a global linear rate for Sinkhorn scaling rays.

Projective Contraction¶

Multiplying either vector by a positive scalar leaves $\mathsf H$ unchanged. It is therefore a metric only on projective classes.

Nested Simplex Images¶

The contraction is already visible on the three-state probability simplex. With the column-vector convention, every positive column-stochastic matrix $K$ maps $\simplex_3$ into its interior, and

The first two panels use the explicit family

This family is doubly stochastic and is strictly positive when $|\delta|<(1-\rho)/3$. We use

The control $K_1$ acts as $\rho I$ on the zero-sum tangent plane. Since the restriction of $A$ has eigenvalues $\pm\mathrm i\sqrt3$, the non-Perron eigenvalues of $K_{\rho,\delta}$ are $\rho\pm\mathrm i\sqrt3\delta$. Thus $K_2$ contracts isotropically while turning gradually. The third panel instead uses the positive doubly stochastic, non-normal kernel

Its non-Perron eigenvalues are approximately .9222 and .7978. Their unequal moduli make the images progressively slender, while non-normality adds shear and a gradual transient turn. Writing

gives $r_1=\rho$, $r_2=\sqrt{\rho^2+3\delta^2}\simeq.9003$, and $r_3\simeq.9222$. This Euclidean asymptotic rate is distinct from the global Birkhoff factor $\lambda(K_i)$ in Hilbert’s metric.

Figure Div contrasts isotropic contraction, rotation, and anisotropic non-normal contraction.

# Isotropic, rotating, and anisotropic positive kernels.
show_book_figure("sinkhorn-birkhoff-simplex-contraction")

Positive Markov kernels contract the three-state simplex and can simultaneously rotate its image. Color progresses from red at $\ell=0$ to blue at $\ell=15$, and the star is the common stationary vector $\mathbf1_3/3$. The isotropic $K_1$ keeps parallel edges, $K_2$ rotates an essentially homothetic triangle, and the non-normal $K_3$ both turns and strongly elongates its image because its tangent modes contract at unequal rates. The outer boundary is only a geometric reference, since the Hilbert metric becomes finite after the first positive image.

The theorem applies to positive linear maps between proper cones. Related nonlinear projective results require order preservation and homogeneity; a generic affine map is not covered.

Nonlinear Sinkhorn Images of the Simplex¶

The projective contraction can be visualized simultaneously for every possible left-scaling ray. Using the row and column maps from the proof, define the full-cycle map

It satisfies $F_u(su)=sF_u(u)$ for every $s>0$. It therefore induces the projective self-map

and the normalized complete-cycle iterates obey $\widehat u^{(\ell+1)}=\widehat F_u(\widehat u^{(\ell)})$. The image sets are nested because

After one cycle they lie in the positive cone, and the preceding theorem gives

The figure reuses the three kernels $K_i$ from Div, with the same uniform Sinkhorn marginals $a=b=\mathbf1_3/3$ in every panel. This isolates the passage from the linear action $K_i^\ell$ to the nonlinear balancing map $\widehat F_{u,i}^{\,\ell}$. The symmetric control remains aligned, $K_2$ turns the curved images, and the unequal tangent rates of $K_3$ produce a pronounced anisotropic collapse. The normalized Sinkhorn fixed ray and stationary Markov vector both equal $\mathbf1_3/3$ in this doubly stochastic example, although they need not coincide in general.

Figure Div reuses the three kernels $\K_i$ from Figure Div, with the same uniform Sinkhorn marginals $\a=\b=\ones_3/3$ in every panel.

# The same three kernels under complete Sinkhorn cycles.
show_book_figure("sinkhorn-projective-scaling-simplex")

Complete Sinkhorn cycles curve, turn and contract the simplex of normalized left scalings. Panel $i$ reuses $K_i$ from the preceding linear figure and the common marginals $a=b=\mathbf1_3/3$. Color is the densely sampled boundary of $\widehat F_{u,i}^{\,\ell}(\simplex_3)$, progressing from the red triangle at $\ell=0$ to the blue curve at $\ell=15$; the star is $\widehat u_i^\star=\mathbf1_3/3$. Reciprocal scaling bends all sixteen boundaries; $K_2$ adds a gradual turn, whereas the non-normal $K_3$ combines turning with a much stronger collapse across one tangent direction. For these invertible kernels the curves are the actual boundaries of the nested image sets. The Hilbert estimate begins after the first positive cycle.

Dual-Potential Form¶

The KL-normalized potentials satisfy $f_\ell=\epsilon\log(u^{(\ell)}\oslash a)$ and $g_\ell=\epsilon\log(v^{(\ell)}\oslash b)$. Hence

The temperature factor is essential when passing to coupling densities:

For $K_\epsilon=e^{-C/\epsilon}$ and $R=\max C-\min C$,

Thus the global rate becomes exponentially pessimistic as $\epsilon\downarrow0$. Sharper continuous analyses obtain polynomial rates under additional semiconcavity or log-concavity assumptions Chizat et al., 2026. In practice, monitor the marginal not just normalized: the source residual at $P^{(\ell)}$ and the target residual at $P^{(\ell+1/2)}$ are the two meaningful posterior diagnostics.

Entropic Optimal Transport Between Gaussians¶

Gaussian marginals provide an explicit finite-dimensional model of Sinkhorn’s behavior. The soft $c$-transform preserves quadratic potentials, the optimal entropic coupling is Gaussian, and the value can be written with matrix square roots Janati et al., 2020. This is the entropic counterpart of the Gaussian $\Wass_2$ and Bures formula.

The controls below expose exactly the quantities in the formula: $\epsilon$ sets the singular-value shrinkage, anisotropy changes the eigenvalues, and the angle changes the covariance misalignment.

This scalar calculation illustrates the general Gaussian convergence picture of Chizat, Delalande and Vaskevicius Chizat et al., 2026: the rate improves when $\epsilon$ is large or the covariance scales overlap well, and deteriorates in the small-temperature limit where the entropic coupling approaches a deterministic Brenier map.

Continuous $\varepsilon$-Sinkhorn Flow¶

This section studies a simultaneous high-resolution, many-iteration limit. It is not a continuous-time interpolation of a fixed-temperature algorithm: the grid is refined while the temperature and the fictitious time step both vanish as $1/k$.

Parabolic Monge--Ampere Limit¶

For the quadratic torus cost $c(x,y)=d_{\mathbb T^d}(x,y)^2/2$, Berman’s scaling discretizes both marginals on a grid of mesh $1/k$, sets $\epsilon_k=1/k$, and assigns duration $1/k$ to each Sinkhorn update Berman, 2020. The $m$-th log-potential is observed at time $t=m/k$. In this coupled limit, Sinkhorn becomes a parabolic Monge--Ampere flow.

To make the scaling explicit, let $\alpha^{(k)}$ and $\beta^{(k)}$ be the positive grid discretizations and define

The multiplicative increment is

The normalized update subtracts the spatial mean of $S_ku$.

No explicit $\epsilon$ remains in the normalized limiting PDE: it records the vanishing temperature before rescaling. The Kahler analogue is the parabolic complex Monge--Ampere equation.

Gaussian Closure¶

Gaussian marginals give a finite-dimensional test case for the continuous flow. The theorem above is stated on the flat torus, but the same local Laplace calculation can be read formally on $\RR^d$ for confining Gaussian densities. Write $\alpha=\Gaussian(m_\alpha,\Sigma_\alpha)$ and $\beta=\Gaussian(m_\beta,\Sigma_\beta)$, and restrict the potential to the quadratic ansatz for which

Taking the spatial gradient of the continuous $\varepsilon$-Sinkhorn PDE removes the additive gauge. Since

coefficient matching in the identity $\partial_tT_t=\nabla\partial_tu_t$ gives

Thus the parabolic Monge--Ampere equation reduces, on the Gaussian ansatz, to a Riccati evolution for the linear part of the transport. The image mean is $q_t$, and the image covariance is

At equilibrium,

which is equivalent to $B\Sigma_\alpha B=\Sigma_\beta$. The stationary map is therefore the Gaussian Brenier map, and the endpoint covariance is governed by the same Bures--Wasserstein geometry as in Section Entropic Optimal Transport Between Gaussians. This Gaussian reduction should be viewed as the finite-dimensional covariance shadow of the vanishing-temperature continuous Sinkhorn limit, not as the fixed-temperature Gaussian Sinkhorn formula itself.

In one dimension the flow reduces to

as long as $1+u_t''>0$. This scalar case is useful for visualization because the potential curves can be plotted directly and the positivity condition is exactly the monotonicity of $x\mapsto x+u_t'(x)$.

Figure Div shows this evolution from the zero initialization for two smooth pairs of marginals.

show_book_figure("sinkhorn-continuous-epsilon-flow")

Continuous $\varepsilon$-Sinkhorn flow in one dimension. The curves are snapshots of the gauge-fixed potential $u_t$, initialized at $u_0=0$, under the parabolic Monge--Ampere equation obtained from Berman’s high-resolution, vanishing-temperature Sinkhorn scaling. Time is encoded from red to blue; the faint bottom silhouettes show the source density in red and the target density in blue.

Monotone Clearing Beyond Variational Sinkhorn¶

The preceding convergence mechanisms mostly used variational structure: Sinkhorn is alternating KL projection, coordinate ascent on a dual objective, or a contraction in a projective metric. There is another, more algebraic, convergence mechanism which keeps the scaling form but discards the existence of an objective. In Galichon’s equilibrium-flow viewpoint, and in related work on substitutability and inverse isotonicity, one studies nonlinear market-clearing equations whose Jacobian has a substitute sign structure Galichon & Jacquet, 2024Galichon et al., 2022. This is the nonlinear analogue of the classical theory of nonsingular M-matrices and M-functions Moré & Rheinboldt, 1973Plemmons, 1977. The relevance for OT is that Sinkhorn is the canonical scaling example, but the same monotone clearing proof also covers fixed-point equations that are not first-order conditions of any convex regularized transport problem.

Two-block clearing maps¶

Write the unknowns as two blocks $z=(u,v)\in\RR^n\times\RR^m$, where $u$ and $v$ will be signed log-scalings, and let

The parallel coordinate-clearing map $T$ is defined by solving, for each coordinate,

This is the Jacobi version: all coordinates are cleared against the old values of the other coordinates. In two-block scaling problems, all coordinates in $u$ decouple when $v$ is fixed, and all coordinates in $v$ decouple when $u$ is fixed. The more common alternating Sinkhorn sweep is the Gauss--Seidel composition of these two block clearings; the order argument below is stated for the parallel map to keep the notation short. The same construction extends to any finite number of blocks.

The M-function assumption says that cross-effects have the sign of substitutes, and that own effects dominate them strongly enough to prevent a global reversal of order. Galichon, Samuelson and Vernet formulate this idea through nonreversingness and unified gross substitutes; for single-valued maps this is the inverse-isotone structure used below. A degenerate $M_0$-function keeps the same order structure but allows a null gauge direction, which is exactly what happens for balanced Sinkhorn before a potential normalization is imposed.

Sinkhorn as the canonical $M_0$-system¶

Let $K=\exp(-C/\epsilon)>0$ and write

The signed convention $v=-\log s$ makes the clearing equations

Then $Q=0$ is exactly $\P\mathbf 1=a$ and $\P^\top\mathbf 1=b$, and coordinate clearing gives the usual Sinkhorn scalings