How to Learn Sinkhorn–Knopp: Matrix Scaling, Doubly Stochastic Matrices, Convergence and Entropic Optimal Transport

Wait, What?

A matrix can be “fixed” by repeatedly normalizing rows and columns—even though every normalization seems to undo part of the previous one.

The Sinkhorn–Knopp algorithm is a deceptively simple iterative method for scaling a nonnegative matrix so that its row and column sums match desired marginals. In the classical square case, the goal is often a doubly stochastic matrix: every row sums to 1 and every column sums to 1. The same scaling idea later became central to computational optimal transport, where entropic regularization turns a large transport problem into alternating multiplicative updates that are highly parallelizable.

Quick Answer

Learn Sinkhorn–Knopp through nonnegative matrices → row/column marginals → diagonal scaling → alternating normalization → convergence conditions → prescribed marginals → entropic optimal transport → numerical stabilization → stopping criteria → production diagnostics. The key is to see the iteration as repeated correction toward two coupled constraints, not as arbitrary renormalization.

1. Start With the Classical Matrix-Scaling Problem

Given a nonnegative matrix A, we want diagonal matrices Dᵣ and D꜀ such that:

B = D_r A D_c

has the desired row and column sums. In the simplest doubly stochastic target, each row and each column should sum to 1. Scaling a row changes the column sums; scaling a column changes the row sums. Sinkhorn–Knopp alternates these corrections until both constraints are approximately satisfied.

2. Work a Tiny Matrix by Hand

Take:

A = [[1, 2],
     [3, 4]]

First divide each row by its row sum. Then compute the new column sums and divide each column by its column sum. The row sums are no longer exactly 1, so repeat. Track the maximum row-sum and column-sum error after each complete cycle. The errors typically shrink toward zero for matrices satisfying the necessary support conditions.

This hand trace makes the algorithm feel less paradoxical: each step enforces one constraint exactly while preserving nonnegativity, and the coupled system converges under appropriate conditions.

3. The Scaling-Vector View Is More Useful

Instead of materializing diagonal matrices, store scaling vectors u and v. For prescribed positive marginals a and b and a nonnegative kernel K, the scaled matrix is:

P = diag(u) K diag(v)

The standard updates are:

u = a / (K v)
v = b / (K^T u)

where division is element-wise. This form is the bridge from classical matrix balancing to modern optimal-transport implementations.

4. What the Invariant Looks Like

After the u update, the row marginals of P match a. After the v update, the column marginals match b. The opposite marginal may drift slightly, but the iteration moves toward simultaneous satisfaction of both constraints.

A learner should be able to explain that invariant before writing vectorized code. If the algorithm is treated as “repeat these two formulas,” numerical bugs become very hard to diagnose.

5. Convergence Is Not Automatic for Every Zero Pattern

The classical theory places structural conditions on the nonzero pattern of the matrix. Positive matrices behave especially well; more general nonnegative matrices require suitable support conditions for a doubly stochastic scaling to exist. Philip Knight’s SIAM analysis surveys convergence and gives explicit rates for important matrix classes.

Professional code should therefore detect impossible or numerically degenerate cases rather than loop forever and call a large iteration count “convergence.”

6. The Optimal-Transport Connection

In discrete optimal transport, we often want a coupling P that moves mass from distribution a to distribution b while minimizing transport cost. The unregularized problem is a linear program. Adding an entropy term produces an entropically regularized problem whose solution can be written using a Gibbs-like kernel:

K_ij = exp(-C_ij / epsilon)

where C is the cost matrix and ε controls regularization. Sinkhorn iterations then scale K so the coupling has the desired marginals.

This is why a classical matrix-balancing algorithm became a modern computational workhorse: the core operations are matrix-vector products and element-wise scaling, which map well to vectorized CPUs, GPUs and accelerators.

7. Understand What ε Changes

Larger ε makes the solution smoother and generally easier to compute, but farther from the sharp unregularized transport optimum. Smaller ε better approximates the original optimal-transport objective but can create numerical underflow and slower or more delicate convergence.

Therefore ε is not just a “speed knob.” It changes the optimization problem. A professional explanation must separate regularization bias from iteration error.

8. Numerical Stability: Work in the Log Domain When Needed

If C/ε is large, exp(-C/ε) can underflow toward zero. Direct scaling vectors may overflow or underflow after many iterations. Stable implementations often use log-domain potentials and log-sum-exp transformations, or periodically absorb scaling factors into dual variables.

This is a recurring numerical lesson: mathematically equivalent formulations can have very different floating-point behaviour.

9. Stopping Criteria Should Measure the Contract

Do not stop merely because u and v change slowly. Measure whether the resulting coupling satisfies its required marginals. Useful diagnostics include:

• maximum absolute row-marginal error;
• maximum absolute column-marginal error;
• L1 marginal violation;
• duality gap or objective diagnostics when available;
• iteration count and numerical warnings.

The stopping rule should match the downstream use. A visualization may tolerate looser marginal error than a scientific inference pipeline.

10. Complexity Comes From Matrix–Vector Work

For a dense n×m kernel, each Sinkhorn half-step performs a matrix-vector product and costs O(nm). The total runtime is that cost multiplied by the number of iterations. Structure changes the picture: sparse kernels, convolutional structure, low-rank approximations or specialized geometry can reduce the effective cost dramatically.

11. How to Learn It Efficiently

Start with a 2×2 or 3×3 matrix and predict the direction of each normalization before calculating it. Run a reference implementation, inspect marginal errors after every iteration, then modify the target marginals. Only after the classical case is secure should you introduce the optimal-transport kernel and ε. This Predict–Run–Investigate–Modify–Make sequence prevents the modern application from hiding the simple matrix-scaling core.

12. Professional Extensions

• ε-scaling: solve a sequence of easier regularized problems as ε decreases.
• Unbalanced transport: relax exact marginal conservation when mass may appear or disappear.
• Batching: solve many transport problems efficiently on accelerators.
• Differentiation: use Sinkhorn layers inside learning systems, while tracking stability and gradient cost.
• Large-scale kernels: exploit sparsity, geometry or low-rank structure instead of materializing every pairwise cost.

Common Failure States

• Normalizing rows and columns without understanding the target marginals.
• Assuming every nonnegative matrix can be scaled to any requested marginals.
• Using a tiny ε and ignoring underflow.
• Stopping because the scaling vectors barely change while marginal violations remain large.
• Comparing regularized and unregularized transport objectives as if they were the same problem.
• Building a full dense cost matrix when the application has exploitable structure.

Practice Ladder

• Beginner: alternately normalize a 2×2 positive matrix by hand.
• Foundation: implement prescribed row/column marginals with scaling vectors.
• Intermediate: add convergence diagnostics and impossible-support tests.
• Advanced: solve a small entropically regularized transport problem and sweep ε.
• Professional: implement a log-stabilized batched version and compare accuracy, marginal error, memory and throughput across dense, sparse and structured kernels.

Learning Hall Boundary

This article owns Sinkhorn–Knopp as alternating diagonal matrix scaling and its entropic optimal-transport use. It does not replace the existing assignment, Hungarian, auction, PageRank, matrix-factorization or Frank–Wolfe jobs.

Evidence Boundary

Sinkhorn and Knopp established the classical matrix-scaling framework; Philip A. Knight’s “The Sinkhorn–Knopp Algorithm: Convergence and Applications,” SIAM Journal on Matrix Analysis and Applications 30(1), 2008, develops convergence results and applications. Gabriel Peyré and Marco Cuturi’s Computational Optimal Transport provides the modern algorithmic bridge from entropy-regularized transport to Sinkhorn scaling, including convergence, stabilization and computational variants.

Professional rule: you understand Sinkhorn when you can state the two marginal constraints, derive the scaling updates, diagnose failure from the matrix support, and distinguish regularization error from numerical and stopping error.