Intermediate Sinkhorn Iterations

• Intermediate Sinkhorn iterations are a sequence of partial algorithmic updates that progress toward the entropy-regularized optimal transport solution while revealing continuous flow dynamics.
• They bridge discrete matrix scaling steps with continuum limits, exemplified by parabolic Monge-Ampère equations and mirror gradient flows.
• The iterations offer statistically optimal and computationally efficient approximations that underpin adaptive stopping criteria and error bounds in stochastic control tasks.

Intermediate Sinkhorn iterations are the sequence of partial updates generated by the Sinkhorn (or iterative proportional fitting) algorithm before convergence to its regularized optimal transport solution. These iterates play a vital role in both the analysis and practical computation of entropy-regularized optimal transport (OT), matrix scaling, and related stochastic control problems. Theoretical and computational studies have revealed that, under various regimes, intermediate Sinkhorn iterates exhibit structured dynamical properties, fast convergence, statistical optimality, and even admit continuum limits that connect to parabolic Monge-Ampère equations and mirror gradient flows in Wasserstein space.

1. Dynamical Trajectories of Sinkhorn Iterates

At each iteration, the Sinkhorn algorithm alternately updates scalings to approach prescribed marginals—yielding a sequence of couplings or scaling vectors whose marginals or potentials evolve discretely toward OT equilibrium. In the canonical entropic OT problem,

$$\min_{\pi \in \Pi(\mu, \nu)} \int \tfrac12\|x-y\|^2\,d\pi(x,y) + \varepsilon\,\mathrm{KL}(\pi\mid \mu\otimes\nu)$$

the sequence of $X$-marginals (or projections) from iterated Sinkhorn updates forms a discrete trajectory $(\rho_k^\varepsilon)$. As the entropic regularization parameter $\varepsilon \to 0$ and iteration count scales as $k \sim t/\varepsilon$, this sequence fills out an absolutely continuous curve in Wasserstein space $\mathcal{P}_2(\mathbb{R}^d)$ termed the Sinkhorn flow (Deb et al., 2023). The continuous limit is characterized by a continuity equation

$$\frac{\partial}{\partial t}\rho_t + \nabla\cdot(\rho_t v_t) = 0,$$

where the velocity field $v_t$ is determined by the dual potential structure (the mirror coordinate, related to the Brenier map): $$v_t(x) = -(\nabla^2 u_t(x))^{-1} \nabla_x(f(x) + \log \rho_t(x)),$$ with $u_t$ a convex potential such that $(\nabla u_t)_\#\rho_t = \nu$.

2. Parabolic Monge-Ampère Equations and Continuum Limits

Intermediate Sinkhorn iterates, in suitable scaling limits, approximate the time-discretized solution to nonlinear parabolic equations of Monge-Ampère type. For quadratic cost and sufficiently regular data, the evolution of the potential is governed by

$$\frac{\partial}{\partial t} u_t(x) = f(x) - g(\nabla u_t(x)) + \log\det \nabla^2 u_t(x),$$

where $f$ and $g$ are the negative log-densities of source and target. This equation, identified in various works (Berman, 2017, Deb et al., 2023), links the dynamics of Sinkhorn’s dual updates to geometric flows (e.g., Ricci-Kähler flow for the torus case). Practically, after $m \sim k t$ iterations, the interpolated Sinkhorn potentials approximate the PDE solution $u_t(x)$ to within $O(1/k)$, exhibiting uniform and explicit error bounds (Berman, 2017).

3. Statistical and Computational Properties of Intermediate Iterates

Intermediate iterates are not merely theoretical constructs: they are statistically and algorithmically optimal for a wide class of downstream tasks. In Schrödinger bridge estimation from finite samples, the squared total variation error between the law defined by $k$ Sinkhorn steps and the ideal bridge decomposes as $O(1/m + 1/n + r^{4k})$, where $m$, $n$ are the sample sizes and $r \in (0,1)$ a function of regularization and domain (Maeda et al., 26 Oct 2025). The first two terms are statistical estimation error floors; only the third is optimized by further Sinkhorn updates, and it decays exponentially. This provides principled guidance for early stopping: computational effort in increasing $k$ beyond

$k \gtrsim \frac{\log(1/m + 1/n)}{4\log r}$

is wasted, since further error is limited by data availability.

For matrix scaling, intermediate Sinkhorn steps converge in error metrics such as $\ell_1$, $\ell_2$ and KL-divergence at rates that are explicitly characterized (Chakrabarty et al., 2018, He, 13 Jul 2025). In dense matrices ($\gamma > 1/2$), error decreases logarithmically in $n$ and $1/\varepsilon$, with upper complexity $O(\log n - \log \varepsilon)$ iterations; for sparse matrices ($\gamma < 1/2$), convergence is polynomially slower, demonstrating a sharp phase transition (He, 13 Jul 2025). These rates allow certification of the accuracy of intermediate solutions without full convergence.

4. Analytical Structure and Flow Interpretations

Under the vanishing regularization limit ($\varepsilon \to 0$), intermediate Sinkhorn iterates fill out a continuous "mirror flow" in Wasserstein space, minimizing the functional $F(\rho) = \mathrm{KL}(\rho\,||\,e^{-f})$ with mirror $\tfrac12 W_2^2(\rho, \nu)$. The velocity norm along this flow corresponds to the metric derivative in the linearized optimal transport (LOT) metric, not in the canonical $W_2$ distance. The norm is given as

$$\|v_t\|_{L^2(\rho_t)} = \lim_{\delta \to 0} \frac{1}{\delta}\, \mathrm{LOT}_{e^{-g}}(\rho_{t+\delta}, \rho_t),$$

where $$\mathrm{LOT}^2_{e^{-g}}(\rho_{t+\delta}, \rho_t) = \int |T_{t+\delta}(y) - T_t(y)|^2\, d\nu(y)$$ and $T_t$ is the optimal transport map from $e^{-g}$ to $\rho_t$ (Deb et al., 2023). This positions the path of Sinkhorn iterates as geodesic in LOT geometry.

The flow admits an associated McKean-Vlasov stochastic differential equation

$$dX_t = \left( -\frac{\partial f}{\partial u_t}(X_t) - \frac{\partial g}{\partial u_t}(X_t^{u_t}) + \frac{\partial h_t}{\partial u_t}(X_t) \right) dt + \sqrt{2 (\nabla^2 u_t(X_t))^{-1}}\, dB_t,$$

whose marginals evolve according to Sinkhorn flow, revealing links to nonlinear diffusions and stochastic control.

5. Exponential Convergence and Regularity Conditions

Exponential convergence of intermediate iterates to the exact regularized OT or Schrödinger bridge solution holds under regularity and convexity assumptions. Essential conditions include:

• A logarithmic Sobolev inequality for the target measure $e^{-f}$, i.e.,

$$\mathrm{KL}(\rho\|e^{-f}) \le \frac{1}{2\alpha} I(\rho|e^{-f}),$$

where $I$ is the relative Fisher information.

• Uniform lower bounds on the determinant of Hessians of the intermediate potentials, $\inf_{x,t} \det \nabla^2 u_t(x) \geq c_0 > 0$.

These guarantee exponential decay in KL divergence,

$$\mathrm{KL}(\rho_t\|e^{-f}) \leq \mathrm{KL}(\rho_0\|e^{-f}) e^{-ct},$$

and thus in Wasserstein distance $W_2(\rho_t,e^{-f}) = O(e^{-ct/2})$, for some $c>0$ (Deb et al., 2023). Pointwise exponential rates also hold for iterates and their gradients on the torus, with constants tied to geometric parameters (Greco et al., 2023).

6. Applications and Broader Implications

Intermediate Sinkhorn iterates are practically significant as computationally efficient, provably accurate, and statistically optimal approximations. They inform solver design for high-dimensional OT, dynamic stochastic control, Schrödinger bridge problems, and generative modeling. The structural convergence results justify the adoption of early stopping or adaptive iteration schedules, especially in large-scale and sample-based scenarios, and underpin the analysis of variance reduction, numerical stability, and computational cost. By connecting discrete iterates to known geometric PDEs, mirror descent flows, and stochastic dynamics, this body of work reveals the foundational role of intermediate Sinkhorn iterates in both the analysis and practical implementation of entropic regularization and OT-based algorithms.

---

Conceptual Focus	Main Result or Formula
Sinkhorn flow (scaling limit)	$$\frac{\partial \rho_t}{\partial t} + \nabla\cdot(\rho_t v_t) = 0$$, $v_t = -(\nabla^2 u_t)^{-1} \nabla(f + \log \rho_t)$, $(\nabla u_t)_\#\rho_t = \nu$
Parabolic Monge-Ampère evolution	$$\frac{\partial u_t}{\partial t}(x) = f(x) - g(\nabla u_t(x)) + \log\det \nabla^2 u_t(x)$$
LOT metric derivative	$$\|v_t\|_{L^2(\rho_t)} = \lim_{\delta\to 0} \frac{1}{\delta} \mathrm{LOT}_{e^{-g}}(\rho_{t+\delta},\rho_t)$$
Statistical error bound (Schr. bridge)	$$\mathbb{E}\left[\mathrm{TV}^2(P_{\infty,\infty}^{*,[0,\tau]}, P_{m,n}^{k,[0,\tau]})\right] = O\left(\frac{1}{m}+\frac{1}{n}+r^{4k}\right)$$, $r<1$, $k$= #Sinkhorn
Exponential convergence (KL)	$$\mathrm{KL}(\rho_t \| e^{-f}) \le \mathrm{KL}(\rho_0 \| e^{-f}) e^{-ct}$$
McKean-Vlasov (Sinkhorn diffusion)	$$dX_t = \left(-\frac{\partial f}{\partial u_t}(X_t) - \frac{\partial g}{\partial u_t}(X_t^{u_t}) + \frac{\partial h_t}{\partial u_t}(X_t)\right)dt + \sqrt{2(\nabla^2 u_t(X_t))^{-1}}\,dB_t$$

---

In summary, intermediate Sinkhorn iterations constitute a rich mathematical object, governing the transient dynamics between arbitrary initialization and regularized OT equilibrium, and their scaling and continuum limits elucidate deep geometric, analytic, and statistical properties of entropic optimal transport algorithms (Deb et al., 2023, Maeda et al., 26 Oct 2025, Berman, 2017, Greco et al., 2023, He, 13 Jul 2025).