Sinkhorn’s algorithm is wonderfully simple: alternately correct the two marginals of a coupling. Yet its speed can be surprisingly hard to predict. As the entropic temperature \(\varepsilon\) decreases, the transport plan becomes sharper and convergence usually slows down. How much slower: exponentially in \(1/\varepsilon\), or only proportionally to \(1/\varepsilon\)?
Gaussians offer a small laboratory in which this question can be answered exactly. For a quadratic cost, the entropic optimal coupling is Gaussian, and Sinkhorn stays within Gaussian distributions. This closed-form structure is studied in our work with Hicham Janati, Boris Muzellec, and Marco Cuturi [1], and in more detail by Akyildiz, Del Moral, and Miguez [2]. In one dimension, the entire iteration becomes a fractional-linear map.
The conclusion is encouraging: for centered Gaussian iterates, the number of full Sinkhorn cycles needed for a fixed relative reduction of a suitable error is of order \(1/\varepsilon\) at low temperature. A recent paper of Lénaïc Chizat, Alex Delalande, and Tomas Vaškevičius [3] explains why favorable temperature dependence extends well beyond Gaussians.
The problem and the convention
Let \(\alpha=\mathcal N(0,A)\) and \(\beta=\mathcal N(0,B)\), where \(A,B\) are positive-definite \(d\times d\) covariance matrices. We use the cost and regularized objective
\[ c(x,y)=\tfrac12|x-y|^2, \qquad J_\varepsilon(\pi)=\int c\,d\pi +\varepsilon\operatorname{KL}(\pi\mid\alpha\otimes\beta), \] \[ \min_{\pi\in\Pi(\alpha,\beta)}J_\varepsilon(\pi). \]Here \(\Pi(\alpha,\beta)\) is the set of joint laws with marginals \(\alpha,\beta\), and \(\operatorname{KL}(\pi\mid\eta)=\int\log(d\pi/d\eta)\,d\pi\), with value \(+\infty\) when absolute continuity fails. The factor \(1/2\) in the cost matters: changing it changes the temperature scale.
Sinkhorn alternates the two potential updates
\[ \begin{aligned} f_{k+1}(x)&=-\varepsilon\log\int e^{(g_k(y)-c(x,y))/\varepsilon}\,d\beta(y),\\ g_{k+1}(y)&=-\varepsilon\log\int e^{(f_{k+1}(x)-c(x,y))/\varepsilon}\,d\alpha(x). \end{aligned} \]The associated density relative to \(\alpha\otimes\beta\) is \(e^{(f+g-c)/\varepsilon}\). The first update enforces the \(X\)-marginal; the second enforces the \(Y\)-marginal. Throughout this post, one cycle means both updates. After a full cycle, the second marginal is exactly \(\beta\), but the first generally is not yet \(\alpha\).
Gaussian closure: iterate conditional covariances
Quadratic potentials whose Gibbs density is integrable produce Gaussian laws, and integrating a Gaussian produces another quadratic potential. Starting with an admissible centered Gaussian initialization therefore preserves centered Gaussian joint laws. The target covariances \(A,B\) do not evolve; the useful evolving quantities are the conditional covariances
\[ S_k=\operatorname{Cov}(X\mid Y), \qquad R_k=\operatorname{Cov}(Y\mid X). \]We observe the iteration just after each \(Y\)-normalization. It can be initialized at any \(S_0\succ0\) by taking \(Y\sim\beta\) and
\[ X\mid Y=y\ \sim\ \mathcal N(\varepsilon^{-1}S_0y,S_0). \]The usual zero-potential initialization belongs to this family after its first full cycle. This indexing lets us state the recursion without a special first step.
Theorem 1 (the Gaussian Sinkhorn recursion). A full cycle is
\[ \begin{aligned} R_k&=(B^{-1}+\varepsilon^{-2}S_k)^{-1},\\ S_{k+1}&=(A^{-1}+\varepsilon^{-2}R_k)^{-1}. \end{aligned} \]The joint covariance at the end of cycle \(k\) is
\[ \Sigma_k= \begin{pmatrix} S_k+\varepsilon^{-2}S_kBS_k & \varepsilon^{-1}S_kB\\ \varepsilon^{-1}BS_k & B \end{pmatrix}. \]These iterates converge to the unique positive-definite solution of
\[ S_*+\varepsilon^{-2}S_*BS_*=A. \]The limiting joint law therefore has both prescribed marginals and is the entropic optimizer. The proof below gives an explicit matrix-square-root formula for \(S_*\).
Proof: Gaussian conditioning and a matrix quadratic equation
The updates. A joint density in this family has the form
\[ p(x,y)\propto\exp\left( -\tfrac12x^\top P x+\varepsilon^{-1}x^\top y -\tfrac12y^\top Q y\right). \]Here \(P,Q\) are diagonal blocks of the joint precision, not marginal precisions. Completing the square gives \(S=P^{-1}\), \(R=Q^{-1}\). To make the \(X\)-marginal have covariance \(A\), its Schur complement must be \(A^{-1}\): thus \(P_{\rm new}=A^{-1}+\varepsilon^{-2}Q^{-1}\). The following \(Y\)-normalization similarly gives \(Q_{\rm new}=B^{-1}+\varepsilon^{-2}P_{\rm new}^{-1}\). This proves the recursion. Finally, total covariance applied to \(X=\varepsilon^{-1}S_kY+Z\), with independent \(Z\sim\mathcal N(0,S_k)\), gives \(\Sigma_k\).
The fixed point. Substituting the first update into the second shows that a fixed point satisfies \(S+\varepsilon^{-2}SBS=A\). Set
\[ D=\varepsilon^{-2}B^{1/2}SB^{1/2}, \qquad M=\varepsilon^{-2}B^{1/2}AB^{1/2}. \]The equation becomes \(D+D^2=M\). Taking its unique positive-definite solution and undoing the change of variables gives the closed form
\[ D_*=\tfrac12(\sqrt{I+4M}-I), \] \[ S_*=\varepsilon^2 B^{-1/2}D_*B^{-1/2}. \]No commutativity of \(A\) and \(B\) is assumed.
Why all positive-definite initializations converge
In the same variables, the iteration is
\[ D_{k+1}=\Phi(D_k),\qquad \Phi(D)=\big[M^{-1}+(I+D)^{-1}\big]^{-1}. \]Matrix inversion reverses the positive-semidefinite order twice, so \(\Phi\) is order preserving, and \(0\prec\Phi(D)\prec M\) for \(D\succeq0\). The sequence started at \(0\) increases; the one started at \(M\) decreases. Both converge, and continuity forces both limits to solve \(D+D^2=M\). Uniqueness makes their limits equal to \(D_*\). For any positive-definite initialization, \(0\prec D_1\prec M\), hence
\[ \Phi^j(0)\preceq D_{j+1}\preceq\Phi^j(M) \quad(j\ge0). \]The bounding sequences converge to the same limit, so \(D_k\to D_*\).
Optimality. The limiting density has the form \(d\pi_*/d(\alpha\otimes\beta)=e^{(f_*+g_*-c)/\varepsilon}\) and the required marginals. For any feasible \(\pi\) of finite objective, cancellation of the two potential integrals yields
\[ J_\varepsilon(\pi)-J_\varepsilon(\pi_*) =\varepsilon\operatorname{KL}(\pi\mid\pi_*)\ge0. \]Equality forces \(\pi=\pi_*\). \(\square\)
One dimension: a fractional-linear iteration
Now let \(\alpha=\mathcal N(0,a)\), \(\beta=\mathcal N(0,b)\), where \(a,b>0\) are variances, not standard deviations. Write \(s_k>0\) for the conditional variance. The two inversions collapse to
\[ s_{k+1}=F(s_k),\qquad F(s)=\frac{a(bs+\varepsilon^2)}{bs+ab+\varepsilon^2}. \]The current first marginal has variance \(\widetilde a_k=s_k+bs_k^2/\varepsilon^2\), and the cross-covariance is \(bs_k/\varepsilon\). Thus following a single positive number follows the entire coupling.
The two fixed points of the fractional-linear map are
\[ s_+=\frac{\varepsilon}{2b} \big(\sqrt{\varepsilon^2+4ab}-\varepsilon\big), \qquad s_-=-\frac{\varepsilon}{2b} \big(\sqrt{\varepsilon^2+4ab}+\varepsilon\big). \]Only \(s_*=s_+\) is a variance. The negative root is nevertheless useful: measuring the distance to one fixed point relative to the other linearizes the dynamics.
Theorem 2 (an exact scalar convergence law). Define the optimal correlation and the full-cycle factor
\[ \rho=\frac{\sqrt{\varepsilon^2+4ab}-\varepsilon}{2\sqrt{ab}}, \qquad q=\rho^4\in(0,1). \]For every \(s_0>0\), the cross-ratio error satisfies
\[ z_k:=\frac{s_k-s_*}{s_k-s_-}=q^kz_0. \]Consequently the entire trajectory is explicit:
\[ s_k=\frac{s_*-q^kz_0s_-}{1-q^kz_0}. \]Unless \(s_0=s_*\), the asymptotic ratio of conditional-variance errors is
\[ \lim_{k\to\infty} \frac{|s_{k+1}-s_*|}{|s_k-s_*|}=q. \]Proof: move the two fixed points to zero and infinity
A fractional-linear map fixing \(s_*\) and \(s_-\) becomes a multiplication in the coordinate \(z=(s-s_*)/(s-s_-)\). Here direct substitution gives
\[ \frac{F(s)-s_*}{F(s)-s_-} =\left( \frac{\sqrt{\varepsilon^2+4ab}-\varepsilon} {\sqrt{\varepsilon^2+4ab}+\varepsilon} \right)^2\frac{s-s_*}{s-s_-}. \]The squared prefactor is \(\rho^4\), since the product of the two terms \(\sqrt{\varepsilon^2+4ab}\pm\varepsilon\) is \(4ab\). Iteration and solving for \(s_k\) prove the formulas. For \(s_0>0\), \(-\rho^2<z_0<1\), so the denominator never vanishes. Finally,
\[ s_k-s_*=(s_*-s_-)\frac{q^kz_0}{1-q^kz_0}, \]which gives the limiting error ratio. \(\square\)
The same asymptotic factor applies to the cross-covariance error and to the first-marginal variance error: the former is linear in \(s_k\), and the latter has a nonzero derivative at \(s_*\).
What happens at low temperature?
The natural dimensionless temperature is \(\tau=\varepsilon/\sqrt{ab}\). Indeed, the optimal joint covariance is
\[ \Sigma_*= \begin{pmatrix}a&\sqrt{ab}\,\rho\\\sqrt{ab}\,\rho&b\end{pmatrix}, \qquad \rho=e^{-\operatorname{arsinh}(\tau/2)}. \]Low temperature makes \(\rho\) close to one: the coupling concentrates near the deterministic map \(Y=\sqrt{b/a}\,X\). The two marginal corrections then undo much of each other’s progress. Quantitatively,
\[ -\log q=4\operatorname{arsinh}(\tau/2), \qquad q=1-2\tau+O(\tau^2). \]To reduce \(|z_k|\) by a factor \(\eta\in(0,1)\), a nonstationary initialization therefore needs exactly
\[ k_\eta= \left\lceil\frac{\log(1/\eta)}{4\operatorname{arsinh}(\tau/2)}\right\rceil \sim\frac{\sqrt{ab}}{2\varepsilon}\log\frac1\eta \quad(\varepsilon\to0). \]This is a statement about relative reduction of the cross-ratio error. Other errors have the same asymptotic geometric factor, but their initial sizes and the requested absolute tolerance also affect iteration counts. At high temperature, \(q\sim\tau^{-4}\), so only a few cycles are needed.
Why the dual objective gap has the squared rate
For a normalized Gibbs coupling \(\pi_k\), define the dual value
\[ \mathcal D(f,g)=\int f\,d\alpha+\int g\,d\beta -\varepsilon\int\big(e^{(f+g-c)/\varepsilon}-1\big) \,d\alpha\,d\beta. \]Since \(\pi_k\) has total mass one, the last integral is zero. Using the marginals of \(\pi_*\) in the log-density ratio gives the exact identity
\[ \delta_k:=\mathcal D(f_*,g_*)-\mathcal D(f_k,g_k) =\varepsilon\operatorname{KL}(\pi_*\mid\pi_k). \]Both couplings have the same \(Y\)-marginal. Integrating the conditional Gaussian relative entropy, and putting \(u=s_k/s_*\), gives
\[ \frac{\delta_k}{\varepsilon} =\frac12\left[ \log u+\frac1u-1 +\frac{\rho^2}{1-\rho^2}\frac{(u-1)^2}{u}\right]. \]This is a positive constant times \((u-1)^2\), to leading order near \(1\). Thus \(\delta_{k+1}/\delta_k\to q^2\) for a nonstationary trajectory. Notice that \(\pi_k\) generally has only one correct marginal: plugging it into the primal objective would not give a feasible primal optimality gap.
For comparison, if the potentials contain linear terms with coefficients \(\ell,m\), Gaussian integration gives \(\ell_{\rm new}=-R m/\varepsilon\), then \(m_{\rm new}=-S_{\rm new}\ell_{\rm new}/\varepsilon\). At the scalar covariance fixed point their full-cycle multiplier is \(S_*R_*/\varepsilon^2=\rho^2\). Our centered initialization sets these coefficients to zero throughout.
Numerical experimentation
Use the temperature slider to follow this transition. The experiment fixes \(a=b=1\) and starts from \(s_0=2s_*\), an admissible centered, \(Y\)-normalized Gaussian. The left panel shows joint covariance ellipses, not two-dimensional marginals: each axis represents one of the two scalar random variables. The middle panel compares conditional-variance error with the exactly geometric cross-ratio error. The right panel shows the exact number of cycles for a millionfold cross-ratio reduction, together with its low-temperature asymptotic.
The curves evaluate the exact recurrence in cross-ratio coordinates, with analytic interpolation between integer cycles; there is no discretization of the Gaussian tails and no fitted convergence rate. Both ellipses show the same Mahalanobis radius. The initialization depends on temperature so that the initial relative conditional-variance error is always one.
Why a worst-case bound looks much worse
Classical Sinkhorn analysis uses the Hilbert projective metric [4], [5]. For a bounded cost with oscillation \(C=\sup c-\inf c\), the kernel \(K=e^{-c/\varepsilon}\) has projective diameter at most \(2C/\varepsilon\). Birkhoff’s theorem then gives the full-cycle contraction bound
\[ \overline q_{\mathrm H}=\tanh^2\!\left(\frac{C}{2\varepsilon}\right), \qquad 1-\overline q_{\mathrm H}\sim4e^{-C/\varepsilon}. \]The exponentially small quantity is the improvement \(1-\overline q_{\mathrm H}\), not the contraction factor itself: that factor approaches one. This elementary guarantee can lead to a cycle count of order \(e^{C/\varepsilon}\log(1/\eta)\). For half-squared distance on a common ball of radius \(R\), one can take \(C\le2R^2\).
This is a worst-case bound, not a prediction that every instance is exponentially slow. It uses extreme kernel ratios and ignores the geometry of the marginals. Moreover, Gaussian supports are unbounded, so the bounded-cost estimate is actually vacuous for them. The exact calculation above is a different, structure-sensitive analysis: its decay exponent is proportional to \(\varepsilon\), not exponentially small in \(1/\varepsilon\).
Beyond Gaussians
Replacing \(\tfrac12|x-y|^2\) by \(-\langle x,y\rangle\) subtracts fixed marginal integrals, so it preserves the optimizer and the temperature. For this bilinear cost, Chizat, Delalande, and Vaškevičius [3, Theorems 1.1–1.2] prove substantially sharper bounds. With compact supports, a convex source domain, and a log-concave source density, their dual-gap bound has a full-cycle factor of the form \(1-c_0\varepsilon\) at low temperature, where \(c_0>0\) depends on the fixed problem geometry. The target marginal need not have a density.
For a source density bounded above and below by positive constants on a compact convex domain, they obtain a global \(1-c_1\varepsilon^2\) factor, improving eventually to \(1-c_2\varepsilon\). These are structural assumptions, not a claim that smoothness alone always suffices. Their Gaussian example separately establishes the sharp linear dependence on temperature for the asymptotic dual-gap rate; it is not an application of their compact-support theorem.
The Gaussian calculation thus captures a useful lesson: the notorious low-temperature slowdown is real, but its severity depends on more than the smallest entry of the Gibbs kernel. Once the distributional structure enters the argument, a much more informative rate becomes visible.
References
- Hicham Janati, Boris Muzellec, Gabriel Peyré, and Marco Cuturi. Entropic Optimal Transport between Unbalanced Gaussian Measures has a Closed Form. NeurIPS, 2020.
- Ömer Deniz Akyildiz, Pierre Del Moral, and Joaquín Miguez. Gaussian entropic optimal transport: Schrödinger bridges and the Sinkhorn algorithm. Foundations of Data Science 11, 50–120, 2026. Open-access preprint.
- Lénaïc Chizat, Alex Delalande, and Tomas Vaškevičius. Sharper exponential convergence rates for Sinkhorn’s algorithm in continuous settings. Mathematical Programming 215, 809–858, 2026. Open-access preprint.
- Joel Franklin and Jens Lorenz. On the scaling of multidimensional matrices. Linear Algebra and its Applications 114–115, 717–735, 1989.
- Gabriel Peyré and Marco Cuturi. Computational Optimal Transport. Foundations and Trends in Machine Learning 11(5–6), 355–607, 2019.