Entropic Regularization: Sinkhorn Algorithm

Entropic regularization makes optimal transport smooth, strictly convex and scalable. This chapter first explains the discrete KL-regularized problem, derives Sinkhorn’s alternating matrix scaling algorithm, and then rewrites the same construction as a relative-entropy projection problem. It then records the general continuous formulation, develops the dual soft-transform picture, explains the path-space Schrodinger problem behind the static coupling formulation, and presents the main convex regularization variants and the debiased Sinkhorn divergence. A final section records a less standard viewpoint: after fixing the potential gauge, the finite-dimensional Sinkhorn equations admit a local holomorphic continuation to complex values of the temperature.

The presentation connects the older matrix-scaling literature Sinkhorn, 1964Sinkhorn & Knopp, 1967Sinkhorn, 1967 with modern entropic OT Cuturi, 2013Peyré & Cuturi, 2019.

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))

Entropic Regularization for Discrete Measures¶

Entropy turns a possibly non-unique linear program into a unique smooth problem. The price is bias, but the reward is differentiability and fast scaling algorithms.

Using this entropy as a regularizing function gives the approximate transport value

Equivalently, the regularizer is $\epsilon\sum_{i,j}\P_{i,j}\log \P_{i,j}$. It penalizes concentrated couplings and makes the objective strictly convex on the relative interior of the transport polytope.

Smoothing Effect¶

The entropy acts as a barrier for positivity and makes $\mathcal{L}_{\C}^{\epsilon}(a,b)$ smooth in $a$, $b$, and $\C$ as long as these variables stay in the relative interior. As $\epsilon\to+\infty$, the minimizer converges to the independent coupling $a\otimes b$; as $\epsilon\to0$, it approaches the optimal face of the original transport linear program.

Figure Div visualizes this temperature-dependent path and contrasts it with a generic logarithmic barrier on linear-programming slacks.

show_book_figure("sinkhorn-entropy-lp-geometry")

Entropic regularization and slack barriers. Large $\epsilon$ selects an interior reference point, while small $\epsilon$ moves the minimizer toward a low-cost face of the transport polytope. The second row gives the analogous entropy-on-slacks picture for a generic linear program.

Entropy Barriers Versus Generic LP Barriers¶

For a generic linear program $\min_z \ell^\top z$ with constraints $Az\le b$, one can introduce positive slacks $s=b-Az$ and penalize them by an entropy. This is a useful analogy, but it is not the standard self-concordant interior-point barrier. The canonical barrier is the Burg, or reverse-KL, barrier $-\sum_i\log s_i$, which leads to Newton systems.

Optimal transport is special because entropy is placed on the entries of $\P$, while the constraints are only row and column marginals. This separable structure turns Bregman projections into diagonal rescalings, giving the Sinkhorn iterations.

Sinkhorn’s Algorithm¶

Sinkhorn’s algorithm is alternating normalization of rows and columns. The key point is that the optimizer of the entropic problem has a multiplicative scaling form.

In matrix notation, $\P=\operatorname{diag}(u)K\operatorname{diag}(v)$. The marginal constraints become

Solving each equation in turn gives Sinkhorn’s algorithm:

The division is entrywise. The scaling vectors are not unique: multiplying $u$ by $\lambda>0$ and $v$ by $1/\lambda$ leaves $\P$ unchanged.

Figure Div exposes the alternating feasibility mechanism on a small matrix: each row or column normalization enforces one marginal exactly while generally perturbing the other.

show_book_figure("sinkhorn-marginal-errors")

Marginal constraints during Sinkhorn scaling. Row normalizations align the red source marginal and leave a blue defect; column normalizations align the blue target marginal and leave a red defect.

The interactive demo exposes the alternating row/column normalization directly. Change the half-step count to see the current coupling acquire one marginal, lose the other, and then converge toward both.

Figure Div shows the same alternating projection mechanism on a dense one-dimensional discretization, where the marginal defects appear as continuous side curves.

show_book_figure("sinkhorn-continuous-marginal-scaling")

Dense Sinkhorn scaling for one-dimensional Gaussian-mixture marginals. The violet side curves are the current row and column sums; the red and blue curves are the prescribed marginals.

After convergence, the regularization strength controls how much of the Gibbs kernel remains visible in the optimal plan. Small $\epsilon$ produces a concentrated transport band, while larger $\epsilon$ spreads the same marginals into a smoother coupling.

Figure Div compares these converged plans at four temperatures while keeping both marginals fixed.

show_book_figure("sinkhorn-coupling-iterations")

Final Sinkhorn couplings for the same one-dimensional marginals and four regularization strengths. Decreasing $\epsilon$ sharpens the plan toward an optimal-transport graph; increasing $\epsilon$ keeps more of the product structure.

Before that, Figure Div tracks the same Sinkhorn scaling in dual variables.

show_book_figure("sinkhorn-potentials-iterations")

KL-normalized dual potentials along the scaling iteration. The logarithmic scaling potentials stabilize as the row/column normalizations converge.

The next interactive demo keeps the iteration count high and varies the temperature. It is the quickest way to see the geometry-bias tradeoff: low temperature is geometric and sharp, high temperature is smooth and closer to independence.

Complexity bounds for Sinkhorn and comparisons with accelerated first-order methods are discussed in Altschuler et al., 2017Dvurechensky et al., 2018Knight, 2008. For a dense $n\times m$ problem, each iteration costs one multiplication by $K$ and one by $K^\top$, so the cost scales like $Cnm$ for $\C$ iterations. For fixed positive $\epsilon$, the marginal error eventually has a linear regime, but small $\epsilon$ makes the Gibbs kernel more peaked and scaling harder.

Figure Div complements the complexity discussion by plotting the marginal defect across half-steps and showing how smaller temperatures slow the observed linear regime.

show_book_figure("sinkhorn-linear-rate-epsilon")

Marginal violation along Sinkhorn half-steps for several values of $\epsilon$. Smaller $\epsilon$ gives sharper transport geometry but slower scaling.

Reformulation Using Relative Entropy¶

The KL formulation identifies Sinkhorn as a projection method. It also prepares the continuous and unbalanced settings, where a reference measure is essential.

Relative Entropy¶

A convenient tool to reformulate and normalize discrete entropy is relative entropy. It turns entropy regularization into a finite-dimensional projection problem and admits a direct measure-theoretic extension.

For matrices with the same total mass, the affine terms cancel and

On fixed-mass couplings, taking $Q=\mathbf 1_{n\times m}$ is equivalent to subtracting the Shannon--Boltzmann entropy.

KL Reformulation of Regularized OT¶

Choosing the tensor product $\a\otimes\b=(\a_i\b_j)_{i,j}$ as reference measure leads to the normalized problem

For every $\P\in\CouplingsD(\a,\b)$,

Hence this problem has exactly the same minimizer as the original entropic OT problem, while its optimal value is $\MKD_\C^\epsilon(\a,\b)+\epsilon(\HD(\a)+\HD(\b))$. The normalization becomes substantive in unbalanced OT, where changing the reference measure is no longer merely an additive shift.

The tensor-product reference is nevertheless useful when supports vary. It makes explicit which entries may vanish and passes cleanly to the continuous formulation.

Figure Div shows how the corresponding KL-normalized dual potentials deform with temperature, from nearly hard Kantorovich potentials to smoother log-sum-exp profiles.

show_book_figure("sinkhorn-dual-potentials-epsilon")

KL-normalized Sinkhorn dual potentials for one-dimensional Gaussian-mixture histograms. For $\epsilon=0.010$ the curves are already close to the unregularized one-dimensional Kantorovich potentials; increasing $\epsilon$ turns this hard $c$-transform geometry into smoother log-sum-exp potentials.

Figure Div illustrates the two limiting regimes established above: the plan approaches a sparse optimal coupling as $\epsilon\downarrow0$ and the product coupling as $\epsilon$ grows.

show_book_figure("sinkhorn-plan-epsilon")

Entropically regularized couplings between the red disk and blue annulus point clouds. The plans are strictly positive for every $\epsilon>0$, but the visible mass pattern evolves from nearly radial and sparse to diffuse as $\epsilon$ increases.

General Formulation¶

The continuous formulation replaces matrices by measures and discrete KL by relative entropy. This section records the measure-theoretic problem, explains how the temperature $\epsilon$ connects exact transport to the independent product coupling, and states the two asymptotic regimes that are useful later: a large-temperature expansion around independence and a small-temperature expansion around quadratic optimal transport.

Measure Formulation¶

The only structural change from the discrete problem is that matrix entries are replaced by densities with respect to a product reference measure. For probability measures $\alpha$ and $\beta$, define

For fixed balanced marginals, the specific product reference only matters up to additive constants, provided the reference marginals are mutually absolutely continuous with $\alpha$ and $\beta$. Its support still matters: it determines which couplings have finite entropy.

Probabilistic Interpretation¶

With this terminology, the entropic problem is

Large $\epsilon$ favors nearly independent endpoints, while small $\epsilon$ suppresses endpoint randomness and recovers an optimal Monge--Kantorovich coupling in the limit. When the unregularized quadratic problem has a Brenier map, this limiting coupling is deterministic.

Sinkhorn for General Measures¶

The multiplicative scaling structure extends from matrices to positive functions. Define the Gibbs kernel and its two integral operators by

Fubini’s theorem gives the adjoint identity $\int_\Xx u\mathcal K_\epsilon v\d\al =\int_\Yy v\mathcal K_\epsilon^*u\d\be$. For compact marginal supports and continuous $c$, Propositions Proposition: Continuous Entropic Duality and Proposition: Existence and Uniqueness of Entropic Dual Potentials provide optimal dual potentials $(f_\epsilon,g_\epsilon)$. Set $u_\epsilon=e^{f_\epsilon/\epsilon}$ and $v_\epsilon=e^{g_\epsilon/\epsilon}$. The continuous density law (54) then becomes

Because $\al\otimes\be$ already contains the prescribed marginals, their target densities with respect to $\al$ and $\be$ are both one. Thus

Starting, for instance, from $v^{(0)}=1$, continuous Sinkhorn alternately enforces these two identities:

Equivalently, pointwise,

Given $v^{(\ell)}$, the first update produces an intermediate coupling with $\Xx$-marginal $\al$; the second produces one with $\Yy$-marginal $\be$, while generally perturbing the first marginal again. The scalings retain the gauge $(u,v)\mapsto(\lambda u,v/\lambda)$.

There is one useful situation in which no iteration is required: choose the target by applying the normalized Gibbs kernel itself to the source.

If $k_\epsilon(x,\cdot)$ is already normalized with respect to $\beta_0$, then $Z_\epsilon=1$ and (30) reduces to

For example, let $\X=\Y=\RR^d$, let $\beta_0$ be Lebesgue measure, and take $c(x,y)=\norm{x-y}^2$. Then $Z_\epsilon=(\pi\epsilon)^{d/2}$ and

Equivalently, the closed-form coupling is the law of $(X,X+\sqrt{\epsilon/2}\,G)$ for independent $X\sim\alpha$ and $G\sim\Gaussian(0,\Id)$. Time-indexed Gaussian blurrings are the forward noising mechanism behind diffusion models, with an additional deterministic rescaling for variance-preserving Ornstein--Uhlenbeck schedules; see Connection with diffusion models..

For discrete measures, setting $\uD_i=\a_i u(x_i)$ and $\vD_j=\b_jv(y_j)$ gives $\P_{i,j}=\a_i\b_j u(x_i)\K_{i,j}v(y_j)=\uD_i\K_{i,j}\vD_j$, so the functional iteration reduces exactly to the matrix Sinkhorn iteration (8). Its logarithmic interpretation as continuous dual block ascent is derived in Dual Sinkhorn for General Measures. Continuous convergence is revisited through a generalized Fortet-type monotonicity argument in Section Sinkhorn Convergence: Monotone Point of View; the finite-dimensional linear rate is studied through Hilbert’s metric in Section Sinkhorn Convergence: Linear Hilbert Metric Rate.

Convergence with $\epsilon$¶

The continuous problem has the same qualitative temperature limits as the finite-dimensional problem, but the zero-temperature selection is subtler. For quadratic transport between smooth densities, the limiting OT plan is typically supported on a graph and is therefore singular with respect to $\alpha\otimes\beta$. Thus the robust statement is weak convergence of minimizers. This is the standard $\Gamma$-convergence mechanism for entropic OT Léonard, 2012Carlier et al., 2017; the density hypothesis below isolates the only approximation point needed in the proof.

The proof is the standard $\Gamma$-convergence argument: the entropy is nonnegative, finite-entropy couplings provide recovery sequences, and Pinsker’s inequality from Theorem Theorem: Pinsker Inequality turns the large-$\epsilon$ entropy bound into total-variation convergence.

Large-Temperature Expansion¶

When $\epsilon$ is large, the entropy dominates and the optimal plan is a small perturbation of the product coupling. The useful object is the part of the cost that cannot be absorbed into row and column potentials. Let $r=\alpha\otimes\beta$ and define

and

The coefficient $c_0$ has zero conditional means. The second-order term follows by expanding the constrained exponential tilt and using this conditional orthogonality.

Small-Temperature Expansion for Smooth Densities¶

At small temperature, entropic transport is a viscous perturbation of quadratic optimal transport. The expansion contains an $\epsilon\log\epsilon$ term from the Gaussian normalization of Brownian bridges, an endpoint entropy correction, and a Fisher-information term along the McCann interpolation. The formula below translates the small-noise Schrödinger expansion to the convention $\|x-y\|^2+\epsilon\operatorname{KL}(\cdot|\alpha\otimes\beta)$ Conforti & Tamanini, 2021Chizat et al., 2020.

The Brownian entropy used in the proof is relative to the sigma-finite endpoint measure $p_T(x,y)\,\d x\d y$, and therefore uses $\mathscr H(\pi|\xi)=\int\log(\d\pi/\d\xi)\,\d\pi$ rather than the finite-measure generalized KL above. This distinction is what fixes the Gaussian normalization and the $\epsilon\log\epsilon$ coefficient.

Dual of Sinkhorn¶

The dual point of view replaces couplings by potentials and soft $c$-transforms. It is the right formulation for stabilized implementations and differentiation.

Discrete Dual¶

The KL-normalized problem has the dual

The optimal potentials are linked to the scaling variables through

Discrete Soft $c$-Transforms¶

For fixed $g$, maximizing the dual with respect to $f$ gives

This is a smoothed minimum.

Exponentiating the alternating soft-transform iterations recovers Sinkhorn’s algorithm. For small $\epsilon$, one must compute the log-sum-exp terms with the usual stabilization trick: subtract the minimum before exponentiating and add it back afterward.

Figure Div visualizes the corresponding soft minimum: decreasing $\epsilon$ sharpens the smooth best response toward the hard $c$-transform envelope.

show_book_figure("sinkhorn-soft-c-transform-epsilon")

Soft $c$-transforms for decreasing temperatures. A positive $\epsilon$ replaces the hard lower envelope by a log-sum-exp soft minimum.

Continuous Dual and Soft-Transforms¶

The continuous formula follows from the same entropy conjugacy as its matrix counterpart; it is not a discretization heuristic.

This is the smooth counterpart of the hard feasibility constraint $f\oplus g\le c$ from the Kantorovich dual.

Dual Sinkhorn for General Measures¶

The soft transforms are not only regularized analogues of hard $c$-transforms: they are the exact block-maximization steps of the continuous dual objective (53). Indeed, for fixed $g$ and $h\in\Cc(\X)$,

Hence exact maximization first over $f$ and then over $g$ gives

The transforms are those of Definition Definition: Continuous Soft $c$-Transforms. Their decorations record their domains: the $\bar c$-transform sends a potential on $\Y$ to one on $\X$, whereas the $c$-transform sends a potential on $\X$ to one on $\Y$.

This dual iteration is exactly the logarithmic form of the continuous scaling iteration (27). Set $u^{(\ell)}=e^{f^{(\ell)}/\epsilon}$ and $v^{(\ell)}=e^{g^{(\ell)}/\epsilon}$. Exponentiating the dual updates and using the kernel operators (24) gives

The coupling density reconstructed from the current potentials is therefore $u^{(\ell)}(x)k_\epsilon(x,y)v^{(\ell)}(y)$ with respect to $\alpha\otimes\beta$, as in (25). At a fixed point, its two marginals are $\alpha$ and $\beta$; equivalently, the potentials jointly maximize the continuous dual.

Neural Dual Solvers¶

The convex-potential structure above suggests a sample-based alternative to evaluating soft transforms on a grid or on all pairs of samples. For the bilinear cost $c(x,y)=-\dotp{x}{y}$, the signs of the dual potentials are convex: writing $\Phi=-f$ and $\Psi=-g$, the zero-temperature constraint $f(x)+g(y)\leq-\dotp{x}{y}$ becomes

For the quadratic cost this is the same statement after subtracting the quadratic terms. One can therefore maximize the continuous dual over parameterized convex potentials, estimating the integrals by stochastic samples.

A useful parameterization is given by input-convex neural networks (ICNNs) Amos et al., 2017Makkuva et al., 2020. The construction mirrors elementary closure rules for convex functions: nonnegative linear combinations preserve convexity, composition with a convex nondecreasing scalar nonlinearity preserves convexity, and the ReLU $r\mapsto\max(r,0)$ is both convex and nondecreasing. Thus a feed-forward network with nonnegative hidden-to-hidden weights and affine skip connections from the input defines a convex function of its input. This gives a flexible cone of convex trial potentials, although the finite-dimensional optimization over the network weights is not a convex optimization problem. Universal-approximation statements must be read with this distinction in mind: max-affine functions are dense among continuous convex functions on compact convex sets, and ICNN-type architectures are designed to inherit this approximation principle. General ReLU universal approximation results, such as the width $d+1$ theorem on compact subsets of $\RR^d$ Hanin, 2019, provide useful background but do not by themselves enforce convexity. In practice, neural dual solvers trade exact Sinkhorn scaling for amortized stochastic optimization of the dual potentials.

The path-space meaning of the static Schrodinger problem now follows the dual construction. A noisy reference dynamics first defines a probability law on trajectories; after optimizing out the conditional law of the path given its endpoints, only the endpoint coupling remains.

Path-Space Schrodinger Problem¶

Schrodinger’s reciprocal problem is naturally posed on paths rather than on endpoint pairs. The Sinkhorn problem appears after the path law is reduced to its two endpoint marginals.

Unregularized Path-Space Transport¶

Let $\Omega=C([0,1];\X)$ be a path space and let $e_t(\omega)=\omega_t$ be the evaluation maps. Given a path action $\mathcal A:\Omega\to[0,+\infty]$, the unregularized path-space problem is

For quadratic Wasserstein geometry on $\RR^d$,

The endpoint cost induced by the action is

For the quadratic action, the minimizing path is the straight segment and $c_{\mathcal A}(x,y)=\norm{x-y}^2$.

Entropic Path-Space Problem¶

Let $R^\epsilon\in\mathcal P(\Omega)$ be a reference path law, for instance a Brownian or Langevin dynamics at noise level $\epsilon$. The dynamic Schrodinger bridge problem is the entropy projection

It asks for the most likely path law, relative to the prior dynamics, among all path laws matching the observed endpoints Schrödinger, 1931Léonard, 2012Léonard, 2014Chen et al., 2016.

Stochastic-Control Interpretation¶

For Brownian references, the entropy projection has a direct control meaning. To fix constants, suppose that $R^\epsilon$ is the law of

so that the generator is $(\epsilon/2)\Delta$. Under the standard finite-entropy assumptions, a path law $M$ with the same initial marginal and $M\ll R^\epsilon$ can be represented, in the weak sense, by an adapted drift $u_t$ of finite energy,

Conversely, such a drift defines an absolutely continuous change of law when the usual Girsanov integrability conditions hold. Girsanov’s formula then gives

If the reference initial law is not $\alpha$, an additional term $\epsilon\operatorname{KL}(\alpha|R^\epsilon_0)$ appears, but it is fixed by the endpoint constraint. Hence, up to this fixed endpoint cost, the Schrodinger bridge can be read as

Thus the bridge is the least energetic change of drift that steers the Brownian prior from $\alpha$ to $\beta$. The optimizer may be chosen Markovian. If $h_t$ denotes the positive space-time harmonic Schrodinger factor associated with the terminal constraint, then the optimal feedback drift has the Doob-transform form

Equivalently, with the value potential $\phi_t=\epsilon\log h_t$, one has $u_t^\star=\nabla\phi_t$ and $\phi_t$ solves the viscous Hamilton--Jacobi equation $\partial_t\phi_t+\frac12\norm{\nabla\phi_t}^2+\frac{\epsilon}{2}\Delta\phi_t=0$. This control viewpoint and the endpoint-coupling reduction below are two faces of the same object: the controlled diffusion describes the whole path law, whereas the static Sinkhorn coupling records only its two endpoints.

Viscous Benamou--Brenier Formulations¶

The dynamic problem also has viscous Benamou--Brenier formulations. Here $v_t$ denotes the forward drift called $u_t$ in the control formulation, whereas $u_t$ denotes the associated current velocity. If the uncontrolled noise is $\sqrt{\sigma}\,\d B_t$, its generator is $(\sigma/2)\Delta$ and

and one minimizes a kinetic action. After absorbing the diffusion into the velocity $u_t=v_t-\frac{\sigma}{2}\nabla\log\rho_t$, the same minimizers are obtained from

up to endpoint entropy terms. The extra term is a Fisher-information penalty.

Brownian Bridges and Sinkhorn Couplings¶

For the convention $\d X_t=\sqrt{\epsilon}\,\d B_t$, the endpoint kernel is $\exp(-\norm{x-y}^2/(2\epsilon))$ and the corresponding static cost is $\norm{x-y}^2/2$. Renaming the static temperature gives the equivalent kernel

After rewriting this prior with respect to $\alpha\otimes\beta$, the endpoint problem is exactly the continuous Sinkhorn problem up to an additive constant. Sinkhorn computes which endpoints should be paired; the path-space Schrodinger bridge then connects each pair by a Brownian bridge.

Figure Div illustrates this endpoint-to-path lifting on a small discrete example.

show_book_figure("sinkhorn-path-space-bridges")

Endpoint couplings lifted to Brownian bridges. Increasing $\epsilon$ both softens the endpoint coupling and amplifies the Brownian fluctuations between paired endpoints.

Marginal-Dependent Problems¶

The balanced Sinkhorn problem fixes both marginals exactly. Many nearby models instead optimize the transported marginals, but only through penalties or constraints applied separately to each marginal. The useful point, emphasized by the generalized scaling algorithms of Chizat, Peyré, Schmitzer and Vialard Chizat et al., 2018, is that entropic OT remains a diagonal-scaling problem whenever these marginal terms admit simple KL-proximal maps.

Let $\mathcal F$ and $\mathcal G$ be proper convex lower semicontinuous functionals on finite nonnegative measures on $\Xx$ and $\Yy$. The unregularized marginal-dependent transport problem is

where $\pi_1=(\mathrm p_{\Xx})_\sharp\pi$ and $\pi_2=(\mathrm p_{\Yy})_\sharp\pi$ are the two marginals of $\pi$. Entropic regularization turns this problem into the scaling-friendly form

where $\alpha$ and $\beta$ are reference measures. In this subsection, when the total mass of $\pi$ is not fixed, the KL term is understood in the generalized sense associated with $\varphi(s)=s\log s-s+1$. Thus, if $\lambda=\alpha\otimes\beta$ and $\pi=\rho\lambda+\pi^\perp$ is the Lebesgue decomposition of $\pi$ with respect to $\lambda$, then

with value $+\infty$ when $\pi^\perp\ne0$. On probability couplings this coincides with the usual relative entropy.

Balanced OT is recovered by taking $\mathcal F=\iota_{\{\alpha\}}$ and $\mathcal G=\iota_{\{\beta\}}$. Unbalanced OT replaces these hard indicators by marginal divergences, as developed later in Section Unbalanced OT. An entropic JKO step fixes the first marginal to the previous iterate and puts the energy on the second marginal, for instance $\mathcal F=\iota_{\{\alpha_t\}}$ and $\mathcal G=E$, with cost $c/(2\tau)$; this is the static counterpart of the minimizing-movement schemes of Chapter Paragraph. Barycenters are the multi-coupling extension: several couplings share one unknown marginal and are treated by the generalized Sinkhorn updates of Section OT Barycenters.

In finite dimension, with reference weights satisfying $a_i,b_j>0$ and proper convex functions $\mathsf F:\RR_+^n\to\RR\cup\{+\infty\}$, $\mathsf G:\RR_+^m\to\RR\cup\{+\infty\}$, the entropic version becomes

Equivalently, if $\K_{ij}=a_i b_j e^{-\C_{ij}/\epsilon}$, the terms involving $\C$ can be absorbed into the Gibbs reference and the problem is, up to the additive constant $\epsilon\sum_{i,j}(a_i b_j-\K_{ij})$,

When $\mathsf F=\iota_{\{\a\}}$ and $\mathsf G=\iota_{\{\b\}}$, the first two dual terms are $\langle\mathbf f,\a\rangle+\langle\mathbf g,\b\rangle$, and one recovers the usual entropic dual. The classical Sinkhorn update is the special case in which the KL proximal maps return the prescribed marginals $\a$ and $\b$.

The usefulness of this formulation is that many KL-proximal maps are explicit.

• Hard marginal constraint. If $\mathsf F=\iota_{\{\a\}}$, then $\operatorname{prox}_{\mathsf F/\epsilon}^{\KLD}(z)=\a$.

• KL marginal relaxation. If $\mathsf F(r)=\tau\KLD(r|\a)$ with $\tau>0$, then, coordinatewise,

so $\mathbf u\leftarrow(\a\oslash z)^{\tau/(\tau+\epsilon)}$, the damped scaling of unbalanced Sinkhorn.

• Pointwise bounds. If $\mathsf F=\iota_{\{\ell\le r\le u\}}$, then the proximal map is the coordinatewise clipping $\operatorname{prox}_{\mathsf F/\epsilon}^{\KLD}(z)_i=\min\{u_i,\max\{\ell_i,z_i\}\}$.

• Total-variation marginal relaxation. If $\mathsf F(r)=\tau\norm{r-\a}_1$ and $\lambda=\tau/\epsilon$, then

• Fixed total mass. If $\mathsf F=\iota_{\{\langle r,\ones\rangle=m\}}$, then $\operatorname{prox}_{\mathsf F/\epsilon}^{\KLD}(z)=m z/\langle z,\ones\rangle$.

These examples explain why generalized Sinkhorn algorithms remain practical: the expensive operation is still multiplication by $\K$ or $\K^\top$, while the model-specific part is a low-dimensional KL-proximal update on a marginal.

Heat Kernels and Hopf--Cole Transforms¶

The Gaussian kernel used by Sinkhorn is also the Euclidean heat kernel. This viewpoint clarifies when entropic OT admits fast grid and surface implementations, and it places soft minima and Hopf--Cole transforms in the same heat-kernel calculus.

Geodesics in Heat¶

On $\mathbb R^d$, the heat kernel for $\partial_t u=\Delta u$ is