---
title: Sinkhorn–Knopp Algorithm
url: https://www.emergentmind.com/topics/sinkhorn-knopp-algorithm-27605038-92db-4352-9e75-9be3eb1a52de
type: topic
---

# Sinkhorn–Knopp Algorithm

The Sinkhorn–Knopp (SK) algorithm is a foundational iterative procedure for scaling the rows and columns of a nonnegative matrix to achieve prescribed sums. Originally developed to construct doubly stochastic matrices, it has evolved into an indispensable tool in matrix analysis, entropic optimal transport, numerical linear algebra, economic complexity, and computational geometry. The method alternates between rescaling rows and columns and admits elegant connections to convex optimization, information geometry, and dynamical systems. Its convergence guarantees, phase complexity, extensions to structural constraints, and diverse algorithmic accelerations have produced a rich and actively developing body of literature.

## 1. Formal Definition, Iterative Scheme, and Existence

Given a nonnegative matrix $A\in\mathbb{R}_+^{n\times m}$ and target row sum vector $r\in\mathbb{R}_{>0}^n$, column sum vector $c\in\mathbb{R}_{>0}^m$ with $\|r\|_1 = \|c\|_1 = h$, the matrix scaling problem is to find diagonal matrices $R, S$ such that $R A S$ has row sums $r$ and column sums $c$. A matrix is said to be $(r,c)$-scalable if such $R,S$ exist. The classical Sinkhorn–Knopp iteration operates as follows [1801.02790]:

- **Initialization (column scaling):**
  $$
  A^{(0)}_{ij} := A_{ij} \cdot \frac{c_j}{\sum_{i'} A_{i'j}}
  $$
- **Row-scaling (step $t$):**
  $$
  B^{(t)}_{ij} := A^{(t)}_{ij} \cdot \frac{r_i}{\sum_{j'} A^{(t)}_{ij'}}
  $$
- **Column-scaling (step $t$):**
  $$
  A^{(t+1)}_{ij} := B^{(t)}_{ij} \cdot \frac{c_j}{\sum_{i'} B^{(t)}_{i'j}}
  $$
This process alternates until both marginal constraints are met within a specified tolerance.

For strictly positive $A$, convergence to the unique scaling is guaranteed (up to normalization) by the Sinkhorn theorem and follows from Perron–Frobenius theory [2309.09089]. For matrices with zeros, $(r,c)$-scalability must be verified on the support graph.

## 2. Interpretation as Alternating Bregman Projections and Mirror Descent

The SK iteration is a manifestation of alternating Bregman projections (specifically, I-projections for KL-divergence) onto the affine spaces defined by the row and column sum constraints. In the context of entropic optimal transport, the problem may be expressed as [1909.06918, 2309.09089]:

$$
\min_{X \ge 0} \; \langle C, X \rangle + \epsilon \sum_{i,j} X_{ij} \log X_{ij}
\quad \text{s.t.} \quad X\mathbf{1} = r, \, X^\top \mathbf{1} = c
$$

Defining the Gibbs kernel $K_{ij} = \exp(-C_{ij}/\epsilon)$, the unique minimizer admits the factorization $X^* = \operatorname{diag}(u) K \operatorname{diag}(v)$, with $u, v > 0$.

Under the mirror-descent framework, with negative entropy as the mirror map and KL-divergence for the penalty, the SK updates are recovered as stochastic mirror descent steps. Each iteration corresponds to a Bregman projection onto the respected marginal constraints, yielding closed-form multiplicative rescalings [1909.06918]. The process generalizes to multidimensional arrays and admits acceleration by considering full gradients or alternating orderings of projections.

## 3. Convergence Guarantees, Error Bounds, and Complexity Phase Transition

The convergence of the SK algorithm is quantified via KL-divergence between iterates and the target marginals. The key result for $(r,c)$-scalable $A$ is [1801.02790]:

- For any $\delta > 0$, after at most
  $$
  T = \left\lceil \frac{\ln(1 + 2 \Delta \rho / \nu)}{\delta} \right\rceil
  $$
  alternating steps, at least one of the KL-divergences $D_{KL}(\tilde r^* || \tilde r^{(t)})$ or $D_{KL}(\tilde c^* || \tilde c^{(t)})$ falls below $\delta$. Here, $\rho = \max_i r_i$, $\nu = \min_{A_{ij} > 0} A_{ij} / \max_{ij} A_{ij}$, $\Delta = \max_j | \{ i : A_{ij} > 0 \} |$.

- $\ell_1$ and $\ell_2$ error bounds are obtained via Pinsker's inequality and a new KL–$\ell_1/\ell_2$ bound. For $\ell_1$, $T = O(h^2 \ln(\Delta\rho/\nu) / \epsilon_1^2)$ suffices for $\| r^{(t)} - r \|_1 \leq \epsilon_1$. For $\ell_2$, using the stronger bound, $T = O(\rho h \ln(\Delta\rho/\nu) (1/\epsilon_2 + 1/\epsilon_2^2))$ suffices for $\| r^{(t)} - r \|_2 \leq \epsilon_2$.

A fundamental discovery is the **phase transition in complexity** based on the matrix density $\gamma$ [2507.09711]. For normalized $A$ with density $\gamma > 1/2$ (i.e., each row and column has at least $\gamma n$ significant entries), SK achieves
$$
O(\log n - \log \varepsilon)
$$
iterations and $\tilde O(n^2)$ total time to reach marginal error $\varepsilon$. For $\gamma < 1/2$, pathological constructions require
$$
\tilde \Omega(\sqrt{n} / \varepsilon)
$$
iterations for $\ell_2$ error, with a sharp threshold at $\gamma = 1/2$. Thus, on dense matrices SK rapidly converges, while in the sparse regime it is intrinsically slow.

## 4. Algorithmic Variants and Accelerations

### Overrelaxation

Overrelaxed SK algorithms introduce a relaxation parameter $\omega > 1$ to accelerate local convergence [2012.12562, 1711.01851]:
$$
u_{k+1} = u_k^{1 - \omega} \odot (a / (K v_k))^{\omega}, \quad
v_{k+1} = v_k^{1 - \omega} \odot (b / (K^\top u_{k+1}))^{\omega}
$$
with $0 < \omega < 2/(1 + \Lambda(K))$ ensuring global convergence. The optimal value for local acceleration is
$$
\omega^{\text{opt}} = 2/(1 + \sqrt{1 - \mu_2^2})
$$
where $\mu_2$ is the subdominant eigenvalue of the local Jacobian. Online adaptive strategies estimate $\mu_2$ via residual tracking [2012.12562].

A Lyapunov-based global convergence guarantee for overrelaxed SK is established by monitoring the KL-divergence to the true solution and ensuring the monotonic decrease of this functional with a locally chosen $\omega$ [1711.01851].

### Arnoldi and Krylov-Type Methods

For matrices with clustered dominant singular values (nearly reducible structure), the standard fixed-point SK iteration slows dramatically. Arnoldi-type accelerations and "power $\ell$" schemes apply high-order or block Krylov iterations to the associated local nonlinear eigenvalue problem for the scaling factors [1810.09700]. By doing so, superlinear convergence rates can be achieved, particularly effective in the nearly decomposable regime.

## 5. Extensions to Structure and Large-Scale Computation

### Constrained Optimal Transport and Zero Patterns

SK naturally accommodates structural constraints by encoding forbidden (zero) assignments directly into the kernel. Let $\mathcal{Z} \subset \{1, ..., m\} \times \{1, ..., n\}$ enforce $P_{ij} = 0$ for $(i,j) \in \mathcal{Z}$. The constrained kernel $K$ is defined as in the unconstrained case but with $K_{ij} = 0$ for forbidden pairs. The scaling iterates are unchanged in form [2404.00003]:
$$
u \leftarrow a / (K v), \quad v \leftarrow b / (K^\top u)
$$
The induced structure is preserved at every step, and convergence to the unique solution is guaranteed if a strictly positive feasible coupling exists.

### High-Dimensional and Structured Fast Solvers

In high dimensions or where $K$ is a Gaussian or translation-invariant kernel, per-iteration complexity can be reduced to $O(n \log n)$ using the non-equally spaced Fast Fourier Transform (NFFT) or related convolutional accelerations [2201.07524]. For measures supported on generic grids, NFFT approximates matrix-vector products required in SK at nearly linear cost.

In applications with extreme sparsity—such as Word Mover's Distance for massive vocabularies—SK is efficiently parallelized using sparse kernel fusion (sparse-dense-dense matrix multiplication compositions) and shared memory architectures, with demonstrated $O(V^2 t)$ scaling and strong parallel efficiency [2005.06727].

## 6. Analytical Connections and Generalizations

- **Geometric and PDE interpretation:** In the continuous setting, and under entropic regularization, the discrete SK iterates approximate the solution to a fully nonlinear parabolic Monge–Ampère PDE in the large sample limit [1712.03082]. This establishes a deep connection to geometric flows such as (twisted) Ricci flows and modern transport theory.
- **Energy-minimization viewpoint:** Both SK and related iterative scaling approaches are recovered as block coordinate descent algorithms for a convex logarithmic-barrier, or energy, functional in the "potential" variables (log-scalings) [2212.12356].
- **Mirror descent and Bregman geometry:** SK is a concrete instance of mirror descent with the negative entropy regularizer and, more generally, represents alternating Bregman projections in the space of row- and column-sum constraints [1909.06918].
- **Equivalence to Fitness–Complexity:** In economic complexity, the Fitness–Complexity algorithm is algebraically equivalent to SK with uniform marginals; both admit the same convergence and scale invariance properties [2212.12356].

## 7. Derivative Convergence, Sensitivity, and Automatic Differentiation

It is established that not only the primal iterates of the SK algorithm but also their derivatives (w.r.t. all problem data: marginals, costs, regularization) converge linearly to the derivatives of the entropic optimal transport solution [2207.12717]. This justifies automatic differentiation and backpropagation through truncated SK iterations in machine learning pipelines, guaranteeing accurate sensitivity and gradient information once a sufficient iteration depth is reached.

---

**Key References:**

- [1801.02790]: Error analysis and tight KL/p-norm convergence rates
- [2012.12562], [1711.01851]: Overrelaxation and adaptive accelerated SK
- [2404.00003]: Constrained/structured SK with zero-pattern
- [2005.06727], [2201.07524]: Parallelized and NFFT-accelerated Sinkhorn
- [2507.09711]: Complexity phase transition at matrix density $\gamma=1/2$
- [1810.09700]: Krylov/Arnoldi-type accelerations
- [2207.12717]: Derivative convergence of Sinkhorn iterates
- [2309.09089], [1712.03082]: Geometric, dynamical, and PDE perspectives
- [2212.12356]: Equivalence between SK and Fitness–Complexity algorithms

The SK algorithm is thus a central, robust, and theoretically well-understood primitive for matrix scaling and entropic transport, adaptable to a wide variety of constraints, domains, and large-scale computational settings.

Source: https://www.emergentmind.com/topics/sinkhorn-knopp-algorithm-27605038-92db-4352-9e75-9be3eb1a52de