---
title: Sinkhorn Normalization in Matrix Scaling and OT
url: https://www.emergentmind.com/topics/sinkhorn-normalization-sn
type: topic
---

# Sinkhorn Normalization in Matrix Scaling and OT

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 [2212.12356] [2401.12253] [1106.1925].

## 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 \(D_r^{(k)}=\mathrm{diag}(u^{(k)})\) and \(D_c^{(k)}=\mathrm{diag}(v^{(k)})\), the standard alternating updates are
\[
u_i^{(k)}=\frac{r_i}{\sum_{j=1}^m A_{ij}v_j^{(k-1)}},\qquad
v_j^{(k)}=\frac{c_j}{\sum_{i=1}^n A_{ij}u_i^{(k)}}.
\]
Equivalently,
\[
D_r^{(k)}=\mathrm{diag}\!\bigl(r\oslash(A\,D_c^{(k-1)}\mathbf 1_m)\bigr),\qquad
D_c^{(k)}=\mathrm{diag}\!\bigl(c\oslash(A^T\,D_r^{(k)}\mathbf 1_n)\bigr),
\]
with entry-wise division denoted by \(\oslash\). Under the total support condition, the scaling factors remain strictly positive and the scaled matrix converges to the unique doubly stochastic solution [2212.12356].

The same procedure can be written as alternate row- and column-scaling operators \(R\) and \(C\). For a nonnegative matrix \(A\), \(R(A)\) rescales each row to sum to \(1\), and \(C(A)\) rescales each column to sum to \(1\). Starting from \(A^{(0)}\), one alternates \(C\) and \(R\), obtaining the familiar row/column normalization sequence [1912.00095].

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 \((r,c)\)-doubly stochastic matrix in finitely many steps, then it does so in at most two scalings [1912.00095].

## 2. Entropic optimal transport formulation

In optimal transport, SN appears as the solver for the entropically regularized transport problem. Given histograms \(r,c\in\mathbb R^n\) with \(\sum_i r_i=\sum_j c_j=1\), a cost matrix \(C\in\mathbb R^{n\times n}\), and regularization strength \(\eta>0\), the problem is
\[
\min_{P\in\mathbb R^{n\times n}}
\sum_{i,j}C_{ij}P_{ij}+\frac1\eta\sum_{i,j}P_{ij}\ln P_{ij},
\qquad
\text{s.t. }P\mathbf 1=r,\;P^T\mathbf 1=c,\;P_{ij}\ge 0.
\]
The unique minimizer has the scaling form
\[
P^\star=\operatorname{diag}(u)\,K\,\operatorname{diag}(v),\qquad
K_{ij}=\exp(-\eta C_{ij}),
\]
and enforcing the marginal constraints yields the Sinkhorn updates
\[
u\leftarrow r\./(Kv),\qquad v\leftarrow c\./(K^T u).
\]
For finite \(\eta\), the entropy term makes the problem strictly convex and enables fast matrix-scaling methods; as \(\eta\to\infty\), the regularized plan converges to the unregularized optimal transport plan [2401.12253].

The dual-variable viewpoint makes the dynamics explicit. Defining
\[
P=\exp\bigl(\eta(-C+x\mathbf 1^T+\mathbf 1 y^T)-1\bigr),
\]
the updates become
\[
x\gets x+\tfrac1\eta\bigl(\ln r-\ln(P\mathbf1)\bigr),\qquad
y\gets y+\tfrac1\eta\bigl(\ln c-\ln(P^T\mathbf1)\bigr).
\]
These iterations alternate maximization of the concave Lyapunov function
\[
f(x,y)=
-\frac1\eta\sum_{i,j}\exp\bigl(\eta(-C_{ij}+x_i+y_j)-1\bigr)
+\sum_i r_i x_i+\sum_j c_j y_j,
\]
whose gradients are \(\nabla_x f=r-P\mathbf1\) and \(\nabla_y f=c-P^T\mathbf1\) [2401.12253].

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-\(1\) 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 [2309.09089].

## 3. Convergence theory, rates, and complexity regimes

The classical global convergence statement is projective: if \(K>0\) entrywise and \(a,b>0\) are prescribed marginals with equal total mass, then there exist unique positive scalings \(U=\mathrm{Diag}(u)\), \(V=\mathrm{Diag}(v)\) such that \(UKV\) has row sums \(a\) and column sums \(b\). In the Hilbert projective metric,
\[
d_H(x,y)=\log\!\bigl(\max_i x_i/y_i\bigr)-\log\!\bigl(\min_i x_i/y_i\bigr),
\]
the Sinkhorn map is a contraction, and one has
\[
d_H\bigl(v^{(k)},v^\ast\bigr)\le \alpha^k d_H\bigl(v^{(0)},v^\ast\bigr),
\]
with
\[
\alpha=\tanh\!\Bigl(\tfrac14\log\!\bigl(\tfrac{\max K}{\min K}\bigr)\Bigr).
\]
This is the standard global geometric convergence mechanism [2309.09089].

For entropic OT, the same linear behavior can be slow in practice when the regularization is small. One explicit rate statement is
\[
\mathrm{KL}(r\|P_t\mathbf1)+\mathrm{KL}(c\|P_t^T\mathbf1)
=
O\bigl((1-e^{-O(\eta\|C\|_\infty)})^t\bigr),
\]
so driving the marginal KL divergence below \(\alpha\) may require
\[
O\bigl(e^{O(\eta\|C\|_\infty)}\ln(1/\alpha)\bigr)
\]
iterations in worst-case theory. The same source notes that empirical dependence on the tolerance \(\epsilon\) is often polynomial, such as \(O(\epsilon^{-1})\) or \(O(\epsilon^{-1/2})\) [2401.12253]. From the geometric viewpoint, as \(\epsilon\to 0\), \(K=\exp(-C/\epsilon)\) becomes increasingly ill-conditioned, the contraction factor approaches \(1\), and the number of iterations to reach a fixed accuracy grows like \(O(1/\epsilon)\) [2309.09089].

A more recent refinement isolates a sharp density threshold. For an \(n\times n\) normalized matrix whose density parameter satisfies \(\gamma>1/2\), Sinkhorn–Knopp reaches \(\ell_1\)-error at most \(\varepsilon\) in
\[
O\!\Bigl((2\gamma-1)^{-5}\,(\log n-\log\varepsilon-\log\rho)\Bigr)
\]
iterations, yielding \(\widetilde O(n^2)\) total time. By contrast, for every \(\gamma<1/2\) there exists a \(\gamma\)-dense matrix requiring \(\Omega(\sqrt n/\varepsilon)\) iterations for \(\ell_2\)-error and \(\Omega(n/\varepsilon)\) iterations for \(\ell_1\)-error. This establishes a phase transition at \(\gamma=1/2\): logarithmic iteration complexity above the threshold and polynomial dependence below it [2507.09711].

## 4. Accelerated and generalized algorithms

The most direct acceleration in the provided literature is "Accelerating Sinkhorn Algorithm with Sparse Newton Iterations" [2401.12253]. The proposed Sinkhorn-Newton-Sparse (SNS) algorithm is explicitly two-stage. First, \(N_1\) standard Sinkhorn steps provide a warm start close enough to the optimum and produce a transport plan \(P\) that admits a small-error sparse approximation. Second, a Newton-type stage uses the Hessian of the Lyapunov potential,
\[
\nabla^2 f(x,y)=
\eta
\begin{pmatrix}
\mathrm{diag}(P\mathbf1) & P \\
P^T & \mathrm{diag}(P^T\mathbf1)
\end{pmatrix},
\]
which is dense in principle but approximately sparse in many practical instances because mass concentrates on a small subset of pairs \((i,j)\). Sparsifying the Hessian by thresholding or retaining only the top \(\lambda n^2\) entries yields \(O(n^2)\) 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 \(f\) becomes super-linear, indeed super-exponential in practice [2401.12253].

A different generalization arises in unbalanced optimal transport. There, classical Sinkhorn may be slow because the dual potentials drift: the coupling
\[
\pi_{ij}\propto \exp[(f_i+g_j-C_{ij})/\epsilon]
\]
is invariant under \((f,g)\mapsto(f+\tau,g-\tau)\), while the unbalanced dual objective is not. Translation invariant Sinkhorn introduces an over-parameterized dual
\[
G(\bar f,\bar g,\lambda)=F(\bar f+\lambda,\bar g-\lambda)
\]
and the invariant objective
\[
H(\bar f,\bar g)=\max_{\lambda\in\mathbb R}G(\bar f,\bar g,\lambda).
\]
The resulting \(H\)-Sinkhorn removes the drift and yields a provably accelerated algorithm; in the KL-penalty case, the composite update is contractive at rate \(\bar\kappa=\kappa_1\kappa_2<1\), whereas the standard \(F\)-Sinkhorn rate is described as \(\simeq (1+\epsilon/\rho)^{-2}\) [2201.00730].

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 \(\Pi\). The method applies incomplete Sinkhorn normalization \(Z^T(A)\) to a nonnegative pre-Sinkhorn matrix, computes the expected gain
\[
U=\sum_{j,k}\Pi_{j,k}\,\ell(r_j,k),
\]
and backpropagates through the row and column Jacobians in \(O(TJ^2)\) time per query [1106.1925]. Closely related, Gumbel–Sinkhorn networks define
\[
S_\tau(X)=\arg\max_{P\in U_n}\langle P,X\rangle+\tau h(P),
\]
treating SN as a continuous analogue of the softmax operator for latent matchings. As \(\tau\to 0^+\), \(S_\tau(X)\) approaches a hard permutation, and the method is demonstrated on sorting numbers, solving jigsaw puzzles, and identifying neural signals in worms [1802.08665].

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 \(\hat S\in\mathbb R^{(M+1)\times(N+1)}\) 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 \(l=5\), slack initialization \(s_{\text{slack}}=0.2\), and \(T=8\) Sinkhorn iterations [2010.00067].

In cross-modal retrieval, SN is used to balance retrieval probabilities and reduce hubness. The formulation replaces a similarity matrix \(S\in\mathbb R^{m\times n}\) by the entropy-regularized transport plan \(P^\ast\) 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 \(k\)-occurrence distribution is reported as baseline \(2.24\), IS \(0.50\), and SN \(0.07\); the Dual Bank Sinkhorn Normalization (DBSN) extension augments the query bank with a target bank in query-agnostic settings [2508.02538].

In economic complexity, the Fitness–Complexity algorithm is shown to be algebraically equivalent to Sinkhorn–Knopp under \(A=M_{cp}\), \(r_c=1\), and \(c_p=1\). The associated potential
\[
g(x,y)=\sum_{c,p}x_c M_{cp} y_p-\sum_c r_c\ln x_c-\sum_p c_p\ln y_p
\]
yields an energy interpretation in which the relative energy cost
\[
E_{cp}=\ln\frac{Q_p}{F_c}
\]
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 [2212.12356].

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 \(P\) satisfying \(P\mathbf1=\mathbf1\) and \(\mathbf1^T P=\mathbf1^T\). 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, \(30\)–\(100\) Sinkhorn iterations suffice for the normalization step [2604.07925].

## 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 \(S=KM\), where \(K\) is symmetric and \(M\) is a mass matrix. The goal is a single positive diagonal scaling \(\Lambda\) such that
\[
Q=\Lambda S\Lambda
\]
preserves constant signals, remains symmetric with respect to the \(M\)-weighted inner product, has nonnegative entries, and has spectrum contained in \([0,1]\). The fixed-point iteration computes
\[
d_i=\sum_j S_{ij}\lambda_j,\qquad
\lambda_i\leftarrow \sqrt{\frac{\lambda_i}{d_i}},
\]
and the paper states that 5–10 iterations typically suffice [2507.06161]. 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 \(\mu\), \(\nu\), 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 \(O(d^3)\). At the limit, the Schrödinger bridge and the entropic transport map admit closed-form Gaussian expressions [2412.18432].

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 [2506.08035].

A more radical structural extension replaces positivity-preserving diagonal scalings by phase-preserving ones. "Sinkhorn normal form for unitary matrices" proves that every unitary \(U\in U(n)\) admits diagonal unitaries \(D_1,D_2\) such that
\[
A=D_1UD_2
\]
has all row sums and column sums equal to one; equivalently,
\[
U=e^{i\phi}LAR
\]
with \(L\) and \(R\) 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 [1408.5728].

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.

Source: https://www.emergentmind.com/topics/sinkhorn-normalization-sn