Entropy turns the linear transport problem into a strictly convex smooth problem. Its optimizer has a multiplicative Gibbs form, and enforcing the two marginals reduces to alternating diagonal scalings. The regularization is therefore simultaneously variational, probabilistic, and algorithmic.
The entropic transport problem
Definition. For \(\varepsilon>0\), define \[ \operatorname{OT}_\varepsilon(\mu,\nu) = \inf_{\pi\in\Pi(\mu,\nu)} \left\{ \int c\,\mathrm{d}\pi +\varepsilon\operatorname{KL}(\pi\,\|\,\mu\otimes\nu) \right\}. \] The reference \(\mu\otimes\nu\) makes the entropy term invariant under changes of coordinates and equal to the mutual information of a coupling.
Existence and uniqueness. If the feasible set contains a finite-objective plan and \(c\) is bounded below with the usual lower-semicontinuity assumptions, the entropic problem has a unique minimizer \(\pi_\varepsilon\).
Proof. Compactness of \(\Pi(\mu,\nu)\), lower semicontinuity of the cost integral, and lower semicontinuity of relative entropy give existence by the direct method. Relative entropy is strictly convex on densities with respect to \(\mu\otimes\nu\), while the cost is linear. Two distinct finite-objective minimizers would have a midpoint with strictly smaller entropy and the same averaged cost, contradicting minimality.
Dual formulation
Entropic duality. Under integrability conditions, \[ \operatorname{OT}_\varepsilon(\mu,\nu) = \sup_{\varphi,\psi} \left[ \int\varphi\,\mathrm{d}\mu+\int\psi\,\mathrm{d}\nu -\varepsilon\int \left\{ e^{(\varphi(x)+\psi(y)-c(x,y))/\varepsilon}-1 \right\} \,\mathrm{d}\mu(x)\,\mathrm{d}\nu(y) \right]. \] At optimality, \[ \frac{\mathrm{d}\pi_\varepsilon}{\mathrm{d}(\mu\otimes\nu)}(x,y) = \exp\!\left(\frac{\varphi_\varepsilon(x)+\psi_\varepsilon(y)-c(x,y)}{\varepsilon}\right). \]
Proof. Write \(\pi=r\,\mu\otimes\nu\) and introduce multipliers \(\varphi,\psi\) for the marginals. Pointwise minimization in \(r\ge0\) uses the convex conjugate \[ \sup_{r\ge0}\{sr-\varepsilon(r\log r-r+1)\} =\varepsilon(e^{s/\varepsilon}-1). \] Substituting \(s=\varphi+\psi-c\) gives the dual. Equality in the conjugacy relation gives the exponential density of the optimizer.
The Schrödinger system
Let \(K_\varepsilon(x,y)=e^{-c(x,y)/\varepsilon}\), \(u(x)=e^{\varphi_\varepsilon(x)/\varepsilon}\), and \(v(y)=e^{\psi_\varepsilon(y)/\varepsilon}\). Then \[ \pi_\varepsilon(\mathrm{d}x,\mathrm{d}y) =u(x)K_\varepsilon(x,y)v(y)\,\mu(\mathrm{d}x)\nu(\mathrm{d}y). \] The marginal constraints become \[ u(x)\int K_\varepsilon(x,y)v(y)\,\mathrm{d}\nu(y)=1, \qquad v(y)\int K_\varepsilon(x,y)u(x)\,\mathrm{d}\mu(x)=1. \] This pair of nonlinear integral equations is the Schrödinger system.
Discrete matrix scaling
For discrete marginals \(a\in\Delta_m\), \(b\in\Delta_n\), let \(K_{ij}=e^{-C_{ij}/\varepsilon}>0\). Absorbing the marginal weights into scaling vectors, the unique optimizer has form \[ P_\varepsilon=\operatorname{diag}(u)K\operatorname{diag}(v), \] with \[ u\odot(Kv)=a,\qquad v\odot(K^{\mathsf T}u)=b. \]
Sinkhorn iterations. Starting from \(v^{(0)}>0\), alternate \[ u^{(k+1)}=a\oslash(Kv^{(k)}), \qquad v^{(k+1)}=b\oslash(K^{\mathsf T}u^{(k+1)}), \] where \(\oslash\) is componentwise division. Every row update enforces the first marginal exactly, and every column update enforces the second.
Why Sinkhorn converges
Sinkhorn scaling theorem. If \(K\) has strictly positive entries and \(a,b\) are strictly positive probability vectors, there exist positive \(u,v\), unique up to \((u,v)\mapsto(\lambda u,\lambda^{-1}v)\), such that \(\operatorname{diag}(u)K\operatorname{diag}(v)\in U(a,b)\). Sinkhorn iterates converge to this scaled matrix.
Proof perspectives. In dual coordinates, each update exactly maximizes the strictly concave dual objective in one block of potentials. The objective increases and is bounded above; strict concavity modulo additive constants identifies the limit. Equivalently, each update is a Kullback--Leibler projection onto one affine marginal set. Alternating Bregman projections converge to the unique intersection minimizer. A quantitative proof uses contraction of the positive kernel in Hilbert's projective metric.
Small-regularization limit
Zero-temperature limit. As \(\varepsilon\downarrow0\), \[ \operatorname{OT}_\varepsilon(\mu,\nu)\longrightarrow\operatorname{OT}_0(\mu,\nu). \] Every weak limit point of \(\pi_\varepsilon\) is an unregularized optimal plan. If several unregularized optimizers exist, \(\pi_\varepsilon\) selects the one minimizing \(\operatorname{KL}(\pi\,\|\,\mu\otimes\nu)\), when this selection problem is finite.
Proof. Since the entropy term is nonnegative, \(\operatorname{OT}_\varepsilon\ge\operatorname{OT}_0\). Evaluate the regularized objective at a finite-entropy optimal plan, or approximate an optimal plan by finite-entropy plans, to obtain the limsup inequality. Compactness gives limit points of \(\pi_\varepsilon\); lower semicontinuity of the cost makes them optimal. Comparing two unregularized optimizers at order \(\varepsilon\) gives the entropy selection rule.
Large-regularization limit
As \(\varepsilon\to\infty\), minimizing mutual information dominates the cost. The unique zero-information coupling is \(\mu\otimes\nu\), so \(\pi_\varepsilon\Rightarrow\mu\otimes\nu\) under mild conditions. Entropic transport therefore interpolates between optimal geometric dependence and independence.
Smooth dependence on the marginals
For fixed \(\varepsilon>0\) and positive kernels, dual potentials are unique modulo constants and vary smoothly with positive discrete marginals. The envelope theorem gives \[ \nabla_a\operatorname{OT}_\varepsilon=\varphi_\varepsilon, \qquad \nabla_b\operatorname{OT}_\varepsilon=\psi_\varepsilon \] in normalized coordinates. This smoothness is the main statistical and optimization advantage of regularization, but it becomes poorly conditioned as \(\varepsilon\downarrow0\).
Log-domain stabilization
When \(\varepsilon\) is small, entries of \(K=e^{-C/\varepsilon}\) may underflow. Writing \(u=e^{f/\varepsilon}\), \(v=e^{g/\varepsilon}\) transforms scaling into log-sum-exp updates: \[ f_i=\varepsilon\log a_i-\varepsilon\log\sum_j \exp\!\left(\frac{g_j-C_{ij}}{\varepsilon}\right). \] Subtracting additive constants from \(f,g\) preserves the plan and prevents numerical overflow. This is the stable realization of the dual gauge invariance.
Regularization bias
\(\operatorname{OT}_\varepsilon(\mu,\mu)\) is generally positive, so the regularized cost is not a distance and is biased toward diffuse couplings. Smaller \(\varepsilon\) reduces geometric bias but slows scaling and increases statistical and numerical sensitivity. The next page studies debiasing and separates the fully discrete and semi-discrete computational geometries.