Fortet–Sinkhorn Iteration in Optimal Transport

• Fortet–Sinkhorn iteration is a scaling algorithm that alternates row and column normalization to produce doubly stochastic matrices or solve entropic optimal transport problems.
• The method applies iterative proportional fitting and dual coordinate ascent, ensuring convergence through metrics like the Hilbert projective metric and KL-potential decreases.
• Recent advances include accelerated schemes and extensions to constrained and structured settings, with applications in imaging, stochastic processes, and semi‐supervised learning.

Fortet–Sinkhorn iteration is the classical matrix-scaling, iterative proportional fitting, or alternate minimization procedure that alternately normalizes rows and columns of a positive matrix to prescribed marginals. In modern optimal transport, the same iteration computes the diagonal scalings of the Gibbs kernel $K=e^{-C/\varepsilon}$ in the entropically regularized transport plan $P=\operatorname{diag}(u)K\operatorname{diag}(v)$, and it is therefore the canonical algorithm behind the Sinkhorn solver for entropic OT. Historically, very similar scaling iterations were already studied by Fortet; hence the literature often speaks of “Fortet–Sinkhorn iterations” or the “Fortet–Sinkhorn algorithm” (Karlsson et al., 2016, Chen et al., 2020).

1. Classical matrix scaling and the Sinkhorn limit

For a positive matrix $A=(a_{i,j})$, the classical scaling operators are defined from the row and column sums

$$\mathrm{row}_i(A)=\sum_{j=1}^n a_{i,j},\qquad \mathrm{col}_j(A)=\sum_{i=1}^m a_{i,j}.$$

If $A$ has positive row sums, the row-scaling matrix is

$$X(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{row}_1(A)},\dots,\frac{1}{\mathrm{row}_m(A)}\right), \qquad R(A)=X(A)A,$$

so $R(A)$ is row stochastic. If $A$ has positive column sums, the column-scaling matrix is

$$Y(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{col}_1(A)},\dots,\frac{1}{\mathrm{col}_n(A)}\right), \qquad C(A)=AY(A),$$

so $C(A)$ is column stochastic. The Fortet–Sinkhorn procedure alternates these two operations: $P=\operatorname{diag}(u)K\operatorname{diag}(v)$0 with $P=\operatorname{diag}(u)K\operatorname{diag}(v)$1 and $P=\operatorname{diag}(u)K\operatorname{diag}(v)$2 (Nathanson, 2019).

For a positive square matrix $P=\operatorname{diag}(u)K\operatorname{diag}(v)$3, there exist positive diagonal matrices $P=\operatorname{diag}(u)K\operatorname{diag}(v)$4 and $P=\operatorname{diag}(u)K\operatorname{diag}(v)$5 such that $P=\operatorname{diag}(u)K\operatorname{diag}(v)$6 is doubly stochastic, and this doubly stochastic matrix is unique; it is called the Sinkhorn limit $P=\operatorname{diag}(u)K\operatorname{diag}(v)$7. The diagonal scalings are unique up to a global scalar: if both $P=\operatorname{diag}(u)K\operatorname{diag}(v)$8 and $P=\operatorname{diag}(u)K\operatorname{diag}(v)$9 are doubly stochastic, then $A=(a_{i,j})$0 and $A=(a_{i,j})$1 for some $A=(a_{i,j})$2. For symmetric positive $A=(a_{i,j})$3, there exists a unique positive diagonal matrix $A=(a_{i,j})$4 such that $A=(a_{i,j})$5 is doubly stochastic. The alternating row/column scaling sequence converges to $A=(a_{i,j})$6, and the corresponding diagonal factors converge as well (Nathanson, 2019).

These facts place Fortet–Sinkhorn iteration within the older matrix-balancing tradition, independently of optimal transport. The optimal transport interpretation did not replace this classical viewpoint; rather, it re-expressed the same scaling mechanism in a Gibbs-kernel geometry.

2. Entropic optimal transport and dual coordinate ascent

In the discrete Kantorovich formulation, two nonnegative histograms $A=(a_{i,j})$7 and $A=(a_{i,j})$8 with equal total mass are coupled by a transport plan $A=(a_{i,j})$9 satisfying

$$\mathrm{row}_i(A)=\sum_{j=1}^n a_{i,j},\qquad \mathrm{col}_j(A)=\sum_{i=1}^m a_{i,j}.$$0

The unregularized OT cost is

$$\mathrm{row}_i(A)=\sum_{j=1}^n a_{i,j},\qquad \mathrm{col}_j(A)=\sum_{i=1}^m a_{i,j}.$$1

Entropic regularization replaces this LP by

$$\mathrm{row}_i(A)=\sum_{j=1}^n a_{i,j},\qquad \mathrm{col}_j(A)=\sum_{i=1}^m a_{i,j}.$$2

with

$$\mathrm{row}_i(A)=\sum_{j=1}^n a_{i,j},\qquad \mathrm{col}_j(A)=\sum_{i=1}^m a_{i,j}.$$3

Because this objective is strictly convex in $$\mathrm{row}_i(A)=\sum_{j=1}^n a_{i,j},\qquad \mathrm{col}_j(A)=\sum_{i=1}^m a_{i,j}.$$4, the solution is unique, and any optimal plan has the factorized form

$$\mathrm{row}_i(A)=\sum_{j=1}^n a_{i,j},\qquad \mathrm{col}_j(A)=\sum_{i=1}^m a_{i,j}.$$5

with $$\mathrm{row}_i(A)=\sum_{j=1}^n a_{i,j},\qquad \mathrm{col}_j(A)=\sum_{i=1}^m a_{i,j}.$$6, $$\mathrm{row}_i(A)=\sum_{j=1}^n a_{i,j},\qquad \mathrm{col}_j(A)=\sum_{i=1}^m a_{i,j}.$$7. The marginal constraints become

$$\mathrm{row}_i(A)=\sum_{j=1}^n a_{i,j},\qquad \mathrm{col}_j(A)=\sum_{i=1}^m a_{i,j}.$$8

The Fortet–Sinkhorn updates are therefore

$$\mathrm{row}_i(A)=\sum_{j=1}^n a_{i,j},\qquad \mathrm{col}_j(A)=\sum_{i=1}^m a_{i,j}.$$9

where $A$0 denotes elementwise division (Karlsson et al., 2016).

The same iteration admits a dual interpretation. For entropic OT, the dual problem can be written as

$A$1

Alternately maximizing this dual in $A$2 and $A$3 yields the scaling updates after exponentiation. In this sense, Fortet–Sinkhorn iteration is block-coordinate ascent on the dual of entropically regularized OT (Karlsson et al., 2016).

A complementary formulation uses the concave Lyapunov potential

$A$4

whose gradient satisfies

$A$5

for

$A$6

Classical Sinkhorn is then alternating maximization in $A$7 and $A$8, with each block update enforcing one set of marginals exactly (Tang et al., 2024).

3. Convergence, complexity, and differentiation

For empirical entropic OT, Fortet–Sinkhorn/IPF is exponentially contractive in the Hilbert projective metric, with contraction factor

$A$9

This is the classical projective-metric picture of Sinkhorn convergence (Maeda et al., 26 Oct 2025).

A different analysis, based on a KL-potential decrease and Pinsker-type control of marginal errors, yields a dimension-free iteration bound: for a positive matrix $$X(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{row}_1(A)},\dots,\frac{1}{\mathrm{row}_m(A)}\right), \qquad R(A)=X(A)A,$$0 with $$X(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{row}_1(A)},\dots,\frac{1}{\mathrm{row}_m(A)}\right), \qquad R(A)=X(A)A,$$1 and smallest entry $$X(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{row}_1(A)},\dots,\frac{1}{\mathrm{row}_m(A)}\right), \qquad R(A)=X(A)A,$$2, the Sinkhorn algorithm reaches

$$X(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{row}_1(A)},\dots,\frac{1}{\mathrm{row}_m(A)}\right), \qquad R(A)=X(A)A,$$3

in

$$X(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{row}_1(A)},\dots,\frac{1}{\mathrm{row}_m(A)}\right), \qquad R(A)=X(A)A,$$4

iterations. This analysis underlies near-linear time approximation algorithms for OT and also motivates the greedy Greenkhorn variant, which updates only the most violated row or column while retaining the same theoretical guarantees (Altschuler et al., 2017).

Differential properties have also been analyzed at the level of the iterates. For parametric entropic OT, if $$X(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{row}_1(A)},\dots,\frac{1}{\mathrm{row}_m(A)}\right), \qquad R(A)=X(A)A,$$5 is the coupling produced by $$X(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{row}_1(A)},\dots,\frac{1}{\mathrm{row}_m(A)}\right), \qquad R(A)=X(A)A,$$6 Sinkhorn–Knopp iterations and $$X(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{row}_1(A)},\dots,\frac{1}{\mathrm{row}_m(A)}\right), \qquad R(A)=X(A)A,$$7 is the exact entropic OT solution, then

$$X(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{row}_1(A)},\dots,\frac{1}{\mathrm{row}_m(A)}\right), \qquad R(A)=X(A)A,$$8

and

$$X(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{row}_1(A)},\dots,\frac{1}{\mathrm{row}_m(A)}\right), \qquad R(A)=X(A)A,$$9

locally uniformly in $R(A)$0, with $R(A)$1 (Pauwels et al., 2022). In parallel, an implicit differentiation framework formulates the Sinkhorn layer through its KKT system and computes analytical gradients jointly with respect to the cost matrix $R(A)$2 and the marginals $R(A)$3, while providing error bounds for approximate forward inputs (Eisenberger et al., 2022).

A common misconception is that Sinkhorn differentiation is merely a numerical convenience. The available results instead show that both unrolled derivatives and implicit derivatives admit rigorous convergence control in the entropic regime.

4. Generalized and constrained Fortet–Sinkhorn schemes

Fortet–Sinkhorn iteration extends beyond the fixed-marginal OT problem. A general variational form is

$R(A)$4

with $R(A)$5 proper, convex, l.s.c. The dual becomes

$R(A)$6

so the first block update remains the classical Sinkhorn scaling, while the second block is replaced by a subgradient condition involving $R(A)$7. For the quadratic choice

$R(A)$8

the generalized update is elementwise and involves the Wright $R(A)$9-function: $A$0 The resulting coordinate ascent computes the proximal operator of $A$1, converges to the unique maximizer, and has locally q-linear convergence. In inverse problems, this proximal map can be embedded in a generalized Douglas–Rachford splitting scheme; in limited-angle computerized tomography, the OMT prior transports mass from a prior image to compensate for missing measurements, avoiding the “fade-in/fade-out” artifacts characteristic of an $A$2 prior, although entropic smoothing still induces blur (Karlsson et al., 2016).

A different constrained extension appears in entropic vector quantile regression. There the dual involves three potentials $A$3, with $A$4 enforcing a mean-independence constraint. One Sinkhorn-type algorithm directly adapts the classical Schrödinger system and requires solving an implicit functional equation for $A$5 at each iteration; a second algorithm replaces this implicit step with a projected gradient update on $A$6. For both algorithms, and for general compactly supported marginals, linear convergence holds in both the dual objective value and the iterates, and the analysis relies on explicit quantitative bounds on the dual potentials and Sinkhorn iterates (Kato et al., 23 Mar 2026).

These constructions show that the defining feature of Fortet–Sinkhorn iteration is not the fixed-marginal setting itself, but the alternating solution of entropic dual blocks with exponential coupling structure.

5. Schrödinger bridges, mirror flows, and continuous-time limits

The Schrödinger bridge problem gives a historical and conceptual antecedent of modern Sinkhorn practice. Fortet’s algorithm was developed for Schrödinger systems, and the discrete Sinkhorn iteration is its matrix-scaling counterpart. In the two-sample Schrödinger bridge setting, intermediate Sinkhorn iterates define a sequence of bridge estimators, and the squared total variation error satisfies the non-asymptotic bound

$A$7

so the algorithmic error decays exponentially in the number of Sinkhorn steps while the statistical error scales with the sample sizes. This perspective also identifies intermediate iterates of DSBM-IMF and BM2 with Sinkhorn-bridge estimators based on the $A$8-th Fortet–Sinkhorn potentials (Maeda et al., 26 Oct 2025).

A complementary line of work reinterprets Sinkhorn as mirror descent on the space of couplings. With

$A$9

the classical step-size-$$Y(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{col}_1(A)},\dots,\frac{1}{\mathrm{col}_n(A)}\right), \qquad C(A)=AY(A),$$0 Sinkhorn update is a mirror descent step. Allowing arbitrary step sizes $$Y(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{col}_1(A)},\dots,\frac{1}{\mathrm{col}_n(A)}\right), \qquad C(A)=AY(A),$$1 yields a $$Y(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{col}_1(A)},\dots,\frac{1}{\mathrm{col}_n(A)}\right), \qquad C(A)=AY(A),$$2-Sinkhorn iteration, and the continuous-time limit is the Sinkhorn$$Y(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{col}_1(A)},\dots,\frac{1}{\mathrm{col}_n(A)}\right), \qquad C(A)=AY(A),$$3 flow / Schr flow, for which

$$Y(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{col}_1(A)},\dots,\frac{1}{\mathrm{col}_n(A)}\right), \qquad C(A)=AY(A),$$4

This framework unifies the Sinkhorn algorithm with Wasserstein mirror flow and the mean-field Schrödinger equation, and extends in parallel to $$Y(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{col}_1(A)},\dots,\frac{1}{\mathrm{col}_n(A)}\right), \qquad C(A)=AY(A),$$5-IPF for dynamical Schrödinger bridges (Karimi et al., 2023).

In the small-regularization regime, with $$Y(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{col}_1(A)},\dots,\frac{1}{\mathrm{col}_n(A)}\right), \qquad C(A)=AY(A),$$6, the marginals of the Fortet–Sinkhorn iteration converge to a continuous Sinkhorn flow in $$Y(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{col}_1(A)},\dots,\frac{1}{\mathrm{col}_n(A)}\right), \qquad C(A)=AY(A),$$7: $$Y(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{col}_1(A)},\dots,\frac{1}{\mathrm{col}_n(A)}\right), \qquad C(A)=AY(A),$$8 The limit is a Wasserstein mirror gradient flow of relative entropy, is equivalent to a parabolic Monge–Ampère-type equation, and can be represented by a McKean–Vlasov diffusion whose marginals follow the Sinkhorn flow (Deb et al., 2023).

6. Accelerations, structured kernels, and representative applications

Recent work has targeted the longstanding trade-off between the low per-iteration cost of Fortet–Sinkhorn scaling and the potentially large number of iterations needed for convergence. One direction is acceleration within the same scaling geometry: Acc-Sinkhorn views row scaling as exact minimization in an inner variable and column scaling as a unit-step dual mirror descent step on a reduced objective. This yields a Hessian-driven Nesterov acceleration with the same scaling form and per-iteration cost, an $$Y(A)=\operatorname{diag}\!\left(\frac{1}{\mathrm{col}_1(A)},\dots,\frac{1}{\mathrm{col}_n(A)}\right), \qquad C(A)=AY(A),$$9 rate under a verifiable stability condition, and complexity

$C(A)$0

for an $C(A)$1-approximation of unregularized OT, improved from

$C(A)$2

for Sinkhorn. On synthetic problems, color transfer, and word alignment, reported speedups are $C(A)$3–$C(A)$4 at small regularization (Xu et al., 28 May 2026).

A second direction augments Fortet–Sinkhorn with Newton steps. Sinkhorn-Newton-Sparse introduces early stopping for the matrix-scaling stage and then switches to a Newton-type method on the dual potential. The key claim is that the Hessian is approximately sparse near the optimum, so sparsification preserves an $C(A)$5 per-iteration cost while yielding orders-of-magnitude faster convergence in practice (Tang et al., 2024). A third direction exploits kernel structure exactly: for one-dimensional $C(A)$6 on a uniform grid, the Gibbs kernel belongs to Lower-/Upper-ColLinear Triangular matrix classes, enabling $C(A)$7 matrix-vector products, closure under Hadamard products and matrix scaling, and hence $C(A)$8 Sinkhorn-type inner loops in IPOT; the two-dimensional block extension gives $C(A)$9 complexity (Liao et al., 2022).

These advances do not eliminate the intrinsic regularization trade-off. Small $P=\operatorname{diag}(u)K\operatorname{diag}(v)$00 better approximates unregularized OT but leads to ill-conditioning, slow convergence, and numerical instability; large $P=\operatorname{diag}(u)K\operatorname{diag}(v)$01 increases smoothing and can blur transport plans and reconstructions. In finite-sample Schrödinger bridge estimation, once the algorithmic term $P=\operatorname{diag}(u)K\operatorname{diag}(v)$02 is well below the statistical term $P=\operatorname{diag}(u)K\operatorname{diag}(v)$03, further Fortet–Sinkhorn iterations are statistically wasted (Karlsson et al., 2016, Maeda et al., 26 Oct 2025).

Representative applications illustrate the breadth of the method:

Setting	Role of Fortet–Sinkhorn iteration	Paper
Limited-angle CT	generalized Sinkhorn computes the OMT proximal block inside Douglas–Rachford splitting; the OMT prior transports mass from prior to target geometry	(Karlsson et al., 2016)
Nested distance of stochastic processes	local entropy-regularized OT subproblems are solved recursively by Sinkhorn; reported speedups range from about $P=\operatorname{diag}(u)K\operatorname{diag}(v)$04 to $P=\operatorname{diag}(u)K\operatorname{diag}(v)$05	(Pichler et al., 2021)
Semi-supervised classification	label allocation is cast as an augmented OT problem with class-capacity and allocation-fraction constraints and solved by log-domain Sinkhorn inside annealed self-training	(Tai et al., 2021)

Fortet–Sinkhorn iteration therefore occupies a dual role. It remains the classical scaling procedure for positive matrices, converging to the Sinkhorn limit, and it is simultaneously the central computational primitive of entropic OT, Schrödinger bridges, generalized proximal transport solvers, and a wide class of constrained or structured transport models.