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, while Fortet’s order argument gives qualitative fixed-point convergence. Robust Bregman estimates then give a non-asymptotic $O(1/\ell)$ dual-gap bound, and Hilbert’s metric gives a global linear contraction when the kernel is uniformly positive. A local spectral analysis sharpens this factor and motivates over-relaxation and variable projection. 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 (103), and their alternating dual-maximization scheme is (106). 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/\log(nm)$, so the temperature decreases with the requested accuracy and with the number of support points.

The assumption $C\geq0$ is harmless for transport costs. For a signed matrix, subtracting $\min_{i,j}C_{ij}$ leaves Sinkhorn iterates and dual gaps unchanged up to gauge; the proposition then applies to the normalized nonnegative cost.

The preceding rate becomes useful when Sinkhorn serves as an approximate solver for exact OT. With the KL-normalized dual, the iterate itself can be compared with the unregularized cost, without any additive entropy correction.

For a dense $n\times n$ problem, one Sinkhorn cycle costs $O(n^2)$ operations. The corollary therefore gives

arithmetic complexity for additive accuracy $\delta$. At fixed accuracy and normalized cost scale, this is nearly quadratic in $n$, compared with the exact $O(n^3)$ Hungarian method of Proposition Proposition: Correctness and Complexity of the Hungarian Primal-Dual Method and the $O(n^3)$ work per scaling phase of auction in Proposition Proposition: Complexity of $\varepsilon$-Scaling. The accuracy comparison is reversed: Sinkhorn’s present guarantee is polynomial in $1/\delta$, whereas Hungarian is exact, scaled auction is exact for integer costs once its tolerance is below $1/n$, and self-concordant interior-point analysis in Section Linear-Programming Algorithms has logarithmic dependence on $1/\delta$ but requires costly Newton solves. These assignment comparisons concern uniform square problems. They show why Sinkhorn is especially effective for large-scale problems requiring low or moderate precision rather than highly accurate exact solutions.

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, so it is a metric on projective classes.

Before applying this to nonlinear Sinkhorn scaling, consider ordinary power iteration. For a positive square kernel $K$, normalize $Ku$ on any fixed positive affine section of the projective cone. The preceding contraction theorem and completeness of this cone give a unique fixed ray $u^\star>0$, hence the Perron--Frobenius eigenvector up to scale, and

Figure Div displays this contraction simultaneously for the entire three-state simplex. It uses $J_3=\mathbf1_3\mathbf1_3^\top/3$ and $K_{\rho,\delta}=\rho I_3+(1-\rho)J_3+\delta A$, where $A$ is skew-cyclic, together with a positive non-normal kernel. These choices respectively show isotropic contraction, contraction with rotation, and anisotropic contraction.

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 the initial simplex to blue near the Perron eigenvector, marked by a star. The panels show an isotropic baseline, a rotating kernel, and an anisotropic non-normal kernel.

Sinkhorn Contraction¶

A Sinkhorn cycle composes two positive kernel maps with entrywise inversion and multiplication. The kernel maps contract Hilbert distance, whereas the latter two operations are projective isometries. A complete row--column cycle thus contracts by the square of the kernel factor. The following theorem, due to Franklin and Lorenz Franklin & Lorenz, 1989, makes this argument precise.

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.

Local Convergence Analysis and Acceleration¶

The global projective rate of Sinkhorn Convergence: Linear Hilbert Metric Rate depends only on the kernel and can be pessimistic near the solution. Linearization reveals a sharper common structure: the dual Hessian, the Jacobian of a complete Sinkhorn cycle, and the curvature of the entropic semi-dual are governed by the same conditional-expectation operator. This gives the exact local factor and motivates blockwise over-relaxation and variable projection.

Conditional Operator and Maximal Correlation¶

Let $L_0^2(\alpha)=\{h\in L^2(\alpha):\int h\,d\alpha=0\}$ and define $L_0^2(\beta)$ similarly.

The two operators are adjoint in the marginal $L^2$ spaces. Conditional Jensen gives $\|T_\pi\|\leq1$, and constants are preserved, so $0\leq\sigma(\pi)\leq1$. In the discrete setting,

with adjointness in the weighted Euclidean inner products. Thus $\sigma_\epsilon$ is the largest nonconstant singular value of the normalized optimal coupling, equivalently the second singular value of $\operatorname{diag}(a)^{-1/2}P_\epsilon \operatorname{diag}(b)^{-1/2}$.

Dual Hessian and Exact Local Rate¶

The next result linearizes the dual objective and the two soft transforms at the same optimum.

The law of total variance gives the useful gap formula

Sinkhorn therefore slows down when the conditional laws $Y\mid X=x$ become nearly deterministic, as happens when $\epsilon\to0$ in a smooth Monge regime.

As expected, a global contraction factor must bound the derivative of the iteration at its fixed point. The next proposition makes this comparison precise, while emphasizing that the kernel-only global estimate can be much more pessimistic than the coupling-dependent local rate.

The global inequality is robust but often pessimistic. In a smooth nondegenerate Monge regime, the conditional laws have width $O(\sqrt\epsilon)$; a coercive Laplace expansion then suggests $1-\sigma_\epsilon^2\asymp\epsilon$, whereas the kernel-only Hilbert factor can approach one exponentially fast.

Blockwise Over-Relaxation¶

A minimal acceleration extrapolates each block toward its soft-transform maximizer. For $\omega>0$,

The case $\omega=1$ is ordinary Sinkhorn; $1<\omega<2$ over-relaxes both half-steps.

Writing $\delta_\epsilon=1-\sigma_\epsilon^2$, the accelerated factor obeys $r_\star=1-2\sqrt{\delta_\epsilon}+2\delta_\epsilon +O(\delta_\epsilon^{3/2})$, compared with $1-\delta_\epsilon$ for ordinary Sinkhorn.

show_book_figure("sinkhorn-overrelaxation")

Blockwise over-relaxation accelerates local Sinkhorn convergence. Left: the small-temperature entropic coupling. Middle: marginal residuals after a common ordinary-Sinkhorn warm start for five relaxation parameters. Right: the continuous curve is the theoretical local gain $1-r_{\rm loc}(\omega)$; filled circles mark the five displayed parameters and crosses are empirical gains fitted from the marginal residuals. The palette switches to teal--blue beyond $\omega_\star$ to distinguish the post-optimal branch.

Entropic Semi-Dual and Variable Projection¶

Eliminating $f$ by its soft $\bar c$-transform gives

This is variable projection (VarPro): one optimizes out a block whenever its conditional optimizer is available. Its Schur-complement Hessian cannot worsen local conditioning Golub & Pereyra, 1973Golub & Pereyra, 2003.

VarPro was introduced for separable nonlinear least-squares problems. A representative example is training a two-layer network through

using the weights $(U,V)$ and activation $\sigma$ of the models studied in Section Training Two-Layer MLPs as Wasserstein Flows. For fixed $V$, the outer weights $U$ are obtained by linear least squares before optimizing the reduced, generally nonconvex objective in $V$ Kim & Lee, 2008. In entropic OT, $(U,V)$ is replaced by the potential pair $(f,g)$, the objective $\mathcal D_\epsilon(f,g)$ is jointly concave, and the eliminated block has the closed-form maximizer $f=g^{\bar c,\epsilon}$.

For a generic VarPro problem, the Schur-complement improvement is only qualitative because the full and reduced spectra depend on several unrelated Hessian blocks. Here both diagonal blocks in (69) are identities, while the off-diagonal blocks are $T_\epsilon$ and $T_\epsilon^*$. Both condition numbers therefore depend only on $\sigma_\epsilon$: the full-dual bound is $2/(1-\sigma_\epsilon)$, whereas the semi-dual bound is $1/(1-\sigma_\epsilon^2)$. Eliminating $f$ thus improves the bound by the explicit factor $2(1+\sigma_\epsilon)$.

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.

Gaussian log densities and soft transforms are quadratic, so every dual Sinkhorn iterate is a quadratic polynomial and the algorithm closes on finitely many coefficients.

Closure of the iterates suggests, but does not by itself prove, Gaussian optimality. Entropy maximization at fixed covariance supplies the variational argument and yields the following closed form.

Debiasing the preceding closed form cancels the separate trace terms and produces a smooth spectral approximation of the Bures covariance term.

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.

The one-dimensional isotropic case permits an exact non-asymptotic analysis along the invariant family of quadratic potentials.