Sinkhorn Normalization in Matrix Scaling and OT

• Sinkhorn Normalization is a method that scales nonnegative matrices via alternating row and column normalizations to achieve prescribed marginals and form doubly stochastic matrices.
• It computes regularized transport plans through diagonal scaling of a Gibbs kernel, with convergence analyzed under the Hilbert projective metric and geometric contraction rates.
• The technique underpins several applications in machine learning, retrieval, and networked data by serving as a differentiable surrogate for permutations, matchings, and constrained assignments.

Sinkhorn Normalization (SN), also called the Sinkhorn–Knopp algorithm, matrix scaling, or iterative proportional fitting, is the alternating row/column normalization procedure that rescales a nonnegative matrix to satisfy prescribed marginals. In the square uniform case, the limit is a doubly stochastic matrix in the Birkhoff polytope; in entropic optimal transport, the same iteration computes the regularized transport plan through diagonal scalings of a Gibbs kernel; and in modern machine learning it serves as a differentiable surrogate for permutations, matchings, and constrained attention or assignment matrices (Mazzilli et al., 2022, Tang et al., 2024, Adams et al., 2011).

1. Classical matrix scaling problem

Classically, SN solves the $(\mathbf r,\mathbf c)$-scaling problem. Given a nonnegative matrix $A\in\mathbb R_+^{n\times m}$ and positive target marginals $\mathbf r\in\mathbb R_+^n$, $\mathbf c\in\mathbb R_+^m$ with $\sum_i r_i=\sum_j c_j$, the goal is to find diagonal scaling matrices $D_r$ and $D_c$ such that

$B=D_r A D_c$

has row sums $r_i$ and column sums $c_j$. Writing $A\in\mathbb R_+^{n\times m}$0 and $A\in\mathbb R_+^{n\times m}$1, the standard alternating updates are

$A\in\mathbb R_+^{n\times m}$2

Equivalently,

$A\in\mathbb R_+^{n\times m}$3

with entry-wise division denoted by $A\in\mathbb R_+^{n\times m}$4. Under the total support condition, the scaling factors remain strictly positive and the scaled matrix converges to the unique doubly stochastic solution (Mazzilli et al., 2022).

The same procedure can be written as alternate row- and column-scaling operators $A\in\mathbb R_+^{n\times m}$5 and $A\in\mathbb R_+^{n\times m}$6. For a nonnegative matrix $A\in\mathbb R_+^{n\times m}$7, $A\in\mathbb R_+^{n\times m}$8 rescales each row to sum to $A\in\mathbb R_+^{n\times m}$9, and $\mathbf r\in\mathbb R_+^n$0 rescales each column to sum to $\mathbf r\in\mathbb R_+^n$1. Starting from $\mathbf r\in\mathbb R_+^n$2, one alternates $\mathbf r\in\mathbb R_+^n$3 and $\mathbf r\in\mathbb R_+^n$4, obtaining the familiar row/column normalization sequence (Cohen et al., 2019).

A recurring misconception is that exact finite termination can happen after an arbitrary number of alternations. The sharp finite-step result is stricter: if an alternating sequence reaches an $\mathbf r\in\mathbb R_+^n$5-doubly stochastic matrix in finitely many steps, then it does so in at most two scalings (Cohen et al., 2019).

2. Entropic optimal transport formulation

In optimal transport, SN appears as the solver for the entropically regularized transport problem. Given histograms $\mathbf r\in\mathbb R_+^n$6 with $\mathbf r\in\mathbb R_+^n$7, a cost matrix $\mathbf r\in\mathbb R_+^n$8, and regularization strength $\mathbf r\in\mathbb R_+^n$9, the problem is

$\mathbf c\in\mathbb R_+^m$0

The unique minimizer has the scaling form

$\mathbf c\in\mathbb R_+^m$1

and enforcing the marginal constraints yields the Sinkhorn updates

$\mathbf c\in\mathbb R_+^m$2

For finite $\mathbf c\in\mathbb R_+^m$3, the entropy term makes the problem strictly convex and enables fast matrix-scaling methods; as $\mathbf c\in\mathbb R_+^m$4, the regularized plan converges to the unregularized optimal transport plan (Tang et al., 2024).

The dual-variable viewpoint makes the dynamics explicit. Defining

$\mathbf c\in\mathbb R_+^m$5

the updates become

$\mathbf c\in\mathbb R_+^m$6

These iterations alternate maximization of the concave Lyapunov function

$\mathbf c\in\mathbb R_+^m$7

whose gradients are $\mathbf c\in\mathbb R_+^m$8 and $\mathbf c\in\mathbb R_+^m$9 (Tang et al., 2024).

A complementary geometric interpretation views SN as alternating projection in Kullback–Leibler divergence onto the row-sum and column-sum affine constraints of the transport polytope. In the continuum formulation, it is the time-$\sum_i r_i=\sum_j c_j$0 Euler, or Trotter–Euler, discretization of coupled nonlinear integral equations for Schrödinger potentials; after semi-discretization, these flows reduce exactly to the discrete Sinkhorn updates (Modin, 2023).

3. Convergence theory, rates, and complexity regimes

The classical global convergence statement is projective: if $\sum_i r_i=\sum_j c_j$1 entrywise and $\sum_i r_i=\sum_j c_j$2 are prescribed marginals with equal total mass, then there exist unique positive scalings $\sum_i r_i=\sum_j c_j$3, $\sum_i r_i=\sum_j c_j$4 such that $\sum_i r_i=\sum_j c_j$5 has row sums $\sum_i r_i=\sum_j c_j$6 and column sums $\sum_i r_i=\sum_j c_j$7. In the Hilbert projective metric,

$\sum_i r_i=\sum_j c_j$8

the Sinkhorn map is a contraction, and one has

$\sum_i r_i=\sum_j c_j$9

with

$D_r$0

This is the standard global geometric convergence mechanism (Modin, 2023).

For entropic OT, the same linear behavior can be slow in practice when the regularization is small. One explicit rate statement is

$D_r$1

so driving the marginal KL divergence below $D_r$2 may require

$D_r$3

iterations in worst-case theory. The same source notes that empirical dependence on the tolerance $D_r$4 is often polynomial, such as $D_r$5 or $D_r$6 (Tang et al., 2024). From the geometric viewpoint, as $D_r$7, $D_r$8 becomes increasingly ill-conditioned, the contraction factor approaches $D_r$9, and the number of iterations to reach a fixed accuracy grows like $D_c$0 (Modin, 2023).

A more recent refinement isolates a sharp density threshold. For an $D_c$1 normalized matrix whose density parameter satisfies $D_c$2, Sinkhorn–Knopp reaches $D_c$3-error at most $D_c$4 in

$D_c$5

iterations, yielding $D_c$6 total time. By contrast, for every $D_c$7 there exists a $D_c$8-dense matrix requiring $D_c$9 iterations for $B=D_r A D_c$0-error and $B=D_r A D_c$1 iterations for $B=D_r A D_c$2-error. This establishes a phase transition at $B=D_r A D_c$3: logarithmic iteration complexity above the threshold and polynomial dependence below it (He, 13 Jul 2025).

4. Accelerated and generalized algorithms

The most direct acceleration in the provided literature is "Accelerating Sinkhorn Algorithm with Sparse Newton Iterations" (Tang et al., 2024). The proposed Sinkhorn-Newton-Sparse (SNS) algorithm is explicitly two-stage. First, $B=D_r A D_c$4 standard Sinkhorn steps provide a warm start close enough to the optimum and produce a transport plan $B=D_r A D_c$5 that admits a small-error sparse approximation. Second, a Newton-type stage uses the Hessian of the Lyapunov potential,

$B=D_r A D_c$6

which is dense in principle but approximately sparse in many practical instances because mass concentrates on a small subset of pairs $B=D_r A D_c$7. Sparsifying the Hessian by thresholding or retaining only the top $B=D_r A D_c$8 entries yields $B=D_r A D_c$9 per-iteration complexity, the same order as Sinkhorn. The paper reports that SNS converges orders of magnitude faster across empirical transport tasks, and once in the Newton regime the convergence in $r_i$0 becomes super-linear, indeed super-exponential in practice (Tang et al., 2024).

A different generalization arises in unbalanced optimal transport. There, classical Sinkhorn may be slow because the dual potentials drift: the coupling

$r_i$1

is invariant under $r_i$2, while the unbalanced dual objective is not. Translation invariant Sinkhorn introduces an over-parameterized dual

$r_i$3

and the invariant objective

$r_i$4

The resulting $r_i$5-Sinkhorn removes the drift and yields a provably accelerated algorithm; in the KL-penalty case, the composite update is contractive at rate $r_i$6, whereas the standard $r_i$7-Sinkhorn rate is described as $r_i$8 (Séjourné et al., 2022).

These developments clarify an important methodological point. SN is not a single fixed algorithmic object but a scaling framework whose core invariants are preserved while the update rule is modified to exploit problem structure: sparse Hessians in regularized OT, translation invariance in unbalanced OT, or other domain-specific symmetries. This suggests why recent work frequently treats SN both as a classical iterative procedure and as a template for designing faster or more stable normalization dynamics.

5. Applications in machine learning, retrieval, and networked data

In learning to rank, SN is used as a differentiable route into the Birkhoff polytope. "Ranking via Sinkhorn Propagation" observes that expectations of rank-linear gains depend only on the marginals of a distribution over permutation matrices, hence only on a doubly stochastic matrix $r_i$9. The method applies incomplete Sinkhorn normalization $c_j$0 to a nonnegative pre-Sinkhorn matrix, computes the expected gain

$c_j$1

and backpropagates through the row and column Jacobians in $c_j$2 time per query (Adams et al., 2011). Closely related, Gumbel–Sinkhorn networks define

$c_j$3

treating SN as a continuous analogue of the softmax operator for latent matchings. As $c_j$4, $c_j$5 approaches a hard permutation, and the method is demonstrated on sorting numbers, solving jigsaw puzzles, and identifying neural signals in worms (Mena et al., 2018).

In structured matching, GCNNMatch embeds a Sinkhorn block directly into online multi-object tracking. After feature extraction and affinity scoring, the method forms a raw score matrix $c_j$6 with a slack row and column to model births and deaths, then alternates row and column softmax normalizations with prescribed slack margins. The reported hyperparameters are entropic weight $c_j$7, slack initialization $c_j$8, and $c_j$9 Sinkhorn iterations (Papakis et al., 2020).

In cross-modal retrieval, SN is used to balance retrieval probabilities and reduce hubness. The formulation replaces a similarity matrix $A\in\mathbb R_+^{n\times m}$00 by the entropy-regularized transport plan $A\in\mathbb R_+^{n\times m}$01 with uniform marginals, equivalently balancing both queries and targets. The same work proves that Inverted Softmax is the single-marginal version obtained by dropping the row constraint. On Flickr30k, the skewness of the $A\in\mathbb R_+^{n\times m}$02-occurrence distribution is reported as baseline $A\in\mathbb R_+^{n\times m}$03, IS $A\in\mathbb R_+^{n\times m}$04, and SN $A\in\mathbb R_+^{n\times m}$05; the Dual Bank Sinkhorn Normalization (DBSN) extension augments the query bank with a target bank in query-agnostic settings (Pan et al., 4 Aug 2025).

In economic complexity, the Fitness–Complexity algorithm is shown to be algebraically equivalent to Sinkhorn–Knopp under $A\in\mathbb R_+^{n\times m}$06, $A\in\mathbb R_+^{n\times m}$07, and $A\in\mathbb R_+^{n\times m}$08. The associated potential

$A\in\mathbb R_+^{n\times m}$09

yields an energy interpretation in which the relative energy cost

$A\in\mathbb R_+^{n\times m}$10

acts as a barrier separating feasible from infeasible exports. The equivalence also clarifies the scale invariance of the algorithm and the role of normalization choices in time-series comparisons (Mazzilli et al., 2022).

In Transformer models, SN has been proposed as a replacement for row-stochastic softmax attention. Alternating row and column normalizers produces a doubly stochastic attention matrix $A\in\mathbb R_+^{n\times m}$11 satisfying $A\in\mathbb R_+^{n\times m}$12 and $A\in\mathbb R_+^{n\times m}$13. The reported analysis shows that SN preserves rank more effectively than standard Softmax attention, although pure self-attention without skips still exhibits doubly-exponential rank decay to one with depth. The same source emphasizes that skip connections are crucial and notes that, in practice, $A\in\mathbb R_+^{n\times m}$14–$A\in\mathbb R_+^{n\times m}$15 Sinkhorn iterations suffice for the normalization step (Lapenna et al., 9 Apr 2026).

6. Structural extensions beyond the standard setting

One line of work modifies the target structure rather than the application domain. "Normalizing Diffusion Kernels with Optimal Transport" introduces symmetric Sinkhorn normalization for a positive symmetric smoothing operator $A\in\mathbb R_+^{n\times m}$16, where $A\in\mathbb R_+^{n\times m}$17 is symmetric and $A\in\mathbb R_+^{n\times m}$18 is a mass matrix. The goal is a single positive diagonal scaling $A\in\mathbb R_+^{n\times m}$19 such that

$A\in\mathbb R_+^{n\times m}$20

preserves constant signals, remains symmetric with respect to the $A\in\mathbb R_+^{n\times m}$21-weighted inner product, has nonnegative entries, and has spectrum contained in $A\in\mathbb R_+^{n\times m}$22. The fixed-point iteration computes

$A\in\mathbb R_+^{n\times m}$23

and the paper states that 5–10 iterations typically suffice (Kessler et al., 8 Jul 2025). This contrasts with standard alternating row/column normalization, which generally breaks symmetry during the intermediate steps.

Another extension treats continuous-state Gaussian models exactly. In "Gaussian entropic optimal transport: Schrödinger bridges and the Sinkhorn algorithm," the iterative proportional fitting procedure remains finite-dimensional when $A\in\mathbb R_+^{n\times m}$24, $A\in\mathbb R_+^{n\times m}$25, and the reference kernel are Gaussian. The iterates stay Gaussian, the covariance updates reduce to Riccati-type recursions related to the Kalman filter, and the per-iteration cost is $A\in\mathbb R_+^{n\times m}$26. At the limit, the Schrödinger bridge and the entropic transport map admit closed-form Gaussian expressions (Akyildiz et al., 2024).

Distributed and asynchronous variants further broaden the framework. "From Local Updates to Global Balance" studies arbitrary sequences in which only one row or one column is normalized at each step. If every row and every column is visited infinitely often and the initial matrix has support, the iterates converge to a doubly stochastic matrix. The same paper uses this fact to analyze a Decentralized Random Walk on a directed graph, where a local walker rescales only the row and column of its current vertex; the limiting transition matrix is doubly stochastic, so the uniform distribution is invariant (Aletti et al., 3 Jun 2025).

A more radical structural extension replaces positivity-preserving diagonal scalings by phase-preserving ones. "Sinkhorn normal form for unitary matrices" proves that every unitary $A\in\mathbb R_+^{n\times m}$27 admits diagonal unitaries $A\in\mathbb R_+^{n\times m}$28 such that

$A\in\mathbb R_+^{n\times m}$29

has all row sums and column sums equal to one; equivalently,

$A\in\mathbb R_+^{n\times m}$30

with $A\in\mathbb R_+^{n\times m}$31 and $A\in\mathbb R_+^{n\times m}$32 diagonal unitaries. The proof is non-constructive and uses symplectic topology via the non-displaceability of the Clifford torus. The same work notes that, unlike the classical positive case, no elementary convergence proof of an alternating algorithm is known in the unitary setting (Idel et al., 2014).

Taken together, these variants show that SN is best understood as a family of normalization principles centered on diagonal rescaling, marginal constraints, and fixed-point geometry. The standard doubly stochastic matrix-scaling problem remains the canonical case, but the same logic extends to sparse Newton acceleration, unbalanced transport, diffusion operators, Gaussian Schrödinger bridges, asynchronous local updates, and even unitary normal forms.