Better and Simpler Error Analysis of the Sinkhorn-Knopp Algorithm for Matrix Scaling

Given a non-negative $n \times m$ real matrix $A$, the {\em matrix scaling} problem is to determine if it is possible to scale the rows and columns so that each row and each column sums to a specified target value for it. This problem arises in many algorithmic applications, perhaps most notably as a preconditioning step in solving a linear system of equations. One of the most natural and by now classical approach to matrix scaling is the Sinkhorn-Knopp algorithm (also known as the RAS method) where one alternately scales either all rows or all columns to meet the target values. In addition to being extremely simple and natural, another appeal of this procedure is that it easily lends itself to parallelization. A central question is to understand the rate of convergence of the Sinkhorn-Knopp algorithm. In this paper, we present an elementary convergence analysis for the Sinkhorn-Knopp algorithm that improves upon the previous best bound. In a nutshell, our approach is to show a simple bound on the number of iterations needed so that the KL-divergence between the current row-sums and the target row-sums drops below a specified threshold $\delta$, and then connect the KL-divergence with $\ell_1$ and $\ell_2$ distances. For $\ell_1$, we can use Pinsker's inequality. For $\ell_2$, we develop a strengthening of Pinsker's inequality, called (KL vs $\ell_1/\ell_2$) in the paper, which lower bounds the KL-divergence by a combination of $\ell_1$ and $\ell_2$ distance. This inequality may be of independent interest. The idea of studying Sinkhorn-Knopp convergence via KL-divergence is not new and has indeed been previously explored. Our contribution is an elementary, self-contained presentation of this approach and an interesting new inequality that yields a significantly stronger convergence guarantee for the extensively studied $\ell_2$-error.


Introduction
In the matrix scaling problem one is given an n × m non-negative matrix A, and positive integer vectors r ∈ Z n >0 and c ∈ Z m >0 with the same ℓ 1 norm The objective is to determine if there exist diagonal matrices R ∈ R n×n and S ∈ R m×m such that the ith row of the matrix RAS sums to r i for all 1 ≤ i ≤ n and the jth column of RAS sums to c j for all 1 ≤ j ≤ m. Of special importance is the case when n = m and r ≡ c ≡ 1 n , the n-dimensional all-ones vector -the (1, 1)-matrix scaling problem wishes to scale the rows and columns of A to make it doubly stochastic. This problem arises in many different areas ranging from transportation planning [12,26] to quantum mechanics [32,1]; we refer the reader to a recent comprehensive survey by Idel [15] for more examples.
One of the most natural algorithms for the matrix scaling problem is the following Sinkhorn-Knopp algorithm [33,34], which is known by many names including the RAS method [4] and the Iterative Proportional Fitting Procedure [30]. The algorithm starts off by multiplicatively scaling all the columns by the columns-sum times c j to get a matrix A (0) with column-sums c. Subsequently, for t ≥ 0, it obtains the B (t) by scaling each row of A (t) by the respective row-sum times r i , and obtain A (t+1) by scaling each column of B (t) by the respective column sums time c j . More precisely, The above algorithm is simple and easy to implement and each iteration takes O(nnz(A)), the number of non-zero entries of A. Furthermore, it has been known for almost five decades [33,34,13,35] that if A is (r, c)-scalable then the above algorithm asymptotically 1 converges to a right solution. More precisely, given ε > 0, there is some finite t by which one obtains a matrix which is "ε-close to having row-and column-sums r and c".
However, the rate of convergence of this simple algorithm is still not fully understood. Since the rate depends on how we measure "ε-closeness", we look at two natural error definitions. For any t, let r (t) := A (t) 1 m denote the vector of row-sums of A (t) . Similarly, we define c (t) := B (t) ⊤ 1 n to be the vector of the column-sums of B (t) . Note that The error of the matrix A t (the error of matrix B t similarly defined) is In this note, we give simple convergence analysis for both error norms. Our result is the following.
which is (r, c)-scalable, and any ε > 0, the Sinkhorn-Knopp algorithm For the special case of n = m and r ≡ c ≡ 1 n , we get the following as a corollary.

in time
Remark 1.1. To our knowledge, the ℓ 1 -error hasn't been explicitly studied in the literature, although for small ε ∈ (0, 1) the same can be deduced from previous papers on matrix scaling [20,14,19,16]. One of our main motivations to look at ℓ 1 -error arose from the connections to perfect matchings in bipartite graphs as observed by Linial, Samorodnitsky and Wigderson [20]. For the ℓ 2 error, which is the better studied notion in the matrix scaling literature, the best analysis is due to Kalantari et al [18,19]. They give aÕ(ρh 2 /ε 2 ) upper bound on the number of iterations for the general problem, and for the special case when m = n and the square matrix has positive permanent (see [18]), they give aÕ(ρ(h 2 − nh + n)/ε 2 ) upper bound. Thus, for (1, 1)-scaling, they get the same result as in Corollary 2. We get a quadratic improvement on h in the general case, and we think our proof is more explicit and simpler. 1 Computationally, this asymptotic viewpoint is unavoidable in the sense that there are simple examples for which the unique matrix scaling matrices need to have irrational entries. For instance, consider the following example from Rothblum and Schneider [29]. The matrix The unique R and S matrices are Remark 1.2. Both parts of Theorem 1 and Corollary 2 are interesting in certain regimes of error. When the error ε is "small" (say, ≤ 1) so that 1/ε 2 ≥ 1/ε, then statement 2 of Corollary 2 implies statement 1 by Cauchy-Schwarz. However, this breaks down when ε is "large" (say ε = δn for some constant δ > 0). In that case, statement 1 implies that in O(ln n/δ 2 ) iterations, the ℓ 1 -error is ≤ δn, but Statement 2 only implies that in O(ln n/δ 2 ) iterations, the ℓ 2 norm is ≤ δn. This "large ℓ 1 -error regime" is of particular interest for an application to approximate matchings in bipartite graphs discussed below.
Applications to Parallel Algorithms for Bipartite Perfect Matching. As a corollary, we get the following application, first pointed by Linial et al [20], to the existence of perfect matchings in bipartite graphs. Let A be the adjacency matrix of a bipartite graph If G has a perfect matching, then clearly there is a doubly stochastic matrix X in the support of A. This suggests the algorithm of running the Sinkhorn-Knopp algorithm to A, and the following claim suggests when to stop. Note that each iteration can be run in O(1) parallel time with m-processors.
Corollary 3 (Fast Parallel Approximate Matchings). Given a bipartite graph G of max-degree ∆ and an ε ∈ (0, 1), O(ln ∆/ε 2 )-iterations of Sinkhorn-Knopp algorithm suffice to distinguish between the case when G has a perfect matching and the case when the largest matching in G has size at most n(1 − ε).
Thus the approximate perfect matching problem in bipartite graphs is in NC for ε as small as polylogarithmic in n. This is not a new result and can indeed be obtained from the works on parallel algorithms for packing-covering LPs [21,36,3,23], but the Sinkhorn-Knopp algorithm is arguably simpler.

Perspective
As mentioned above, the matrix scaling problem and in particular the Sinkhorn-Knopp algorithm has been extensively studied over the past 50 years. We refer the reader to Idel's survey [15] and the references within for a broader perspective; in this subsection we mention the most relevant works.
We have already discussed the previously best known, in their dependence on h, analysis for the Sinkhorn-Knopp algorithm in Remark 1.1. For the special case of strictly positive matrices, better rates are known. Kalantari and Khachiyan [16] showed that for positive matrices and the (1, 1)-scaling problem, the Sinkhorn-Knopp algorithm obtains ℓ 2 error ≤ ε in O( √ n ln(1/ν)/ε)-iterations; this result was extended to the general matrix scaling problem by Kalantari et al [19]. In a different track, Franklin and Lorenz [13] show that in fact the dependence on ε can be made logarithmic, and thus the algorithm has "linear convergence", however their analysis 2 has a polynomial dependence of (1/ν). All these results use the positivity crucially and seem to break down even with one 0 entry. The Sinkhorn-Knopp algorithm has polynomial dependence on the error parameter and therefore is a "pseudopolynomial" time approximation. We conclude by briefly describing bounds obtained by other algorithms for the matrix scaling problem whose dependence on ε is logarithmic rather than polynomial. Kalantari and Khachiyan [17] describe a method based on the ellipsoid algorithm which runs in time O(n 4 ln(n/ε) ln(1/ν)). Nemirovskii and Rothblum [25] describe a method with running time O(n 4 ln(n/ε) ln ln(1/ν)). The first strongly polynomial time approximation scheme (with no dependence on ν) was due to Linial, Samoridnitsky, and Wigderson [20] who gave ã O(n 7 ln(h/ε)) time algorithm. Rote and Zachariasen [28] reduced the matrix scaling problem to flow problems to give a O(n 4 ln(h/ε)) time algorithms for the matrix scaling problem. To compare, we should recall that Theorem 1 shows that our algorithm runs in time O(nnz(A)h 2 /ε 2 ) time.
Very recently, two independent works obtain vastly improved running times for matrix scaling. Cohen et al [9] givẽ O(nnz(A) 3/2 ) time algorithm, while Allen-Zhu et al [2] give aÕ(n 7/3 + nnz(A) · (n + n 1/3 h 1/2 )) time algorithm; the tildes in both the above running times hide the logarithmic dependence on ε and ν. Both these algorithms look at the matrix scaling problem as a convex optimization problem and perform second order methods.
We briefly describe the entropy minimization viewpoint. Given two non-negative matrices M and N let us define the Kullback-Leibler divergence 3 between M and N as follows with the convention that the summand is zero if both M ij and N ij are 0, and is ∞ if M ij > 0 and N ij = 0. Let Φ r be the set of n × m matrices whose row-sums are r and let Φ c be the set of n × m matrices whose column sums are c.
Given matrix A suppose we wish to find the matrix A * = arg min B∈Φr∩Φc D(B, A). One algorithm for this is to use the method of alternate projections with respect to the KL-divergence [7] (also known as I-projections [10]) which alternately finds the matrices in Φ r and Φ c closest in the KL-divergence sense to the current matrix at hand, and then sets the minimizer to be the current matrix. It is not too hard to see (see Idel [15], Observation 3.17 for a proof) that the above alternate projection algorithm is precisely the Sinkhorn-Knopp algorithm. Therefore, at least in this sense, the right metric to measure the distance to optimality is not the ℓ 1 or the ℓ 2 error as described in the previous section, but the rather the KL-divergence between the normalized vectors as described below. Let π (t) r := r (t) /h be the n-dimensional probability vector whose ith entry is r (t) i /h; similarly define the mdimensional vector π (t) c . Let π r denote the n-dimensional probability vector with the ith entry being r i /h; similarly define π c . Recall that the KL-divergence between two probability distributions p, q is defined as D KL (p||q) := n i=1 p i ln(q i /p i ). The following theorem gives the convergence time for the KL-divergence.
is (r, c)-scalable, then for any δ > 0 there is a t ≤ T = ⌈ ln(1+2nρ/ν) Proof. Let Z := RAS be a matrix with row-sums r and column-sums c for diagonal matrices R, S. Recall A 0 is the matrix obtained by column-scaling A. Note that the minimum non-zero entry of A 0 is ≥ ν/n. Proof. By definition, For a fixed j, the vectors  q min is the smallest non-zero entry of q. For our purpose, we note that the minimum non-zero probability of the A (0) j distribution being ≥ ν/nρ. Therefore, the second summand is at most ln(1 + 2nρ/ν) giving us D(Z, A (0) ) ≤ 1 h m j=1 c j · ln(1 + 2nρ/ν) = ln(1 + 2nρ/ν).
Proof. The LHS of the first equality is simply r ). The other equation follows analogously.
Theorem 1 follows from Theorem 4 using connections between the KL-divergence and the ℓ 1 and ℓ 2 norms. One is the following famous Pinsker's inequality which allows us to easily prove part 1 of Theorem 1. Given any two probability distributions p, q, Proof of Theorem 1, Part 1. Apply (Pinsker) on the vectors π r and π (t) r to get Set δ := ε 2 2h 2 and apply Theorem 4. In O h 2 ln(nρ/ν) ε 2 time we would get a matrix with δ > D KL (π r ||π (t) r ) which from the above inequality would imply ||r (t) − r|| 1 ≤ ε.
To prove Part 2, we need a way to relate the ℓ 2 norm and the KL-divergence. In order to do so, we prove a different lower bound which implies Pinsker's inequality (with a worse constant), but is significantly stronger in certain regimes. To the best of our knowledge this is a new bound which may be of independent interest in other domains. Below we state the version which we need for the proof of Theorem 1, part 2. This is an instantiation of the general inequality Lemma 3.1 whcih we prove in Section 3.
Proof of Theorem 1, Part 2. We apply Lemma 2.3 on the vectors π r and π (t) r . Lemma 2.3 gives us where C = 1 − ln 2. If the second summand in the parenthesis of the RHS is ≥ 1 2 ||r (t) − r|| 2 2 , then we get where we used the weak fact that the sum of some positive numbers is at least the square-root of the sum of their squares. In any case, we get the following To complete the proof of part 2 of Theorem 1, set δ := and apply Theorem 4. In O ρh ln (nρ/ν) · 1 ε + 1 ε 2 time we would get a matrix with δ ≥ D KL (π r ||π (t) r ). If the minimum of the RHS of (2) is the first term, then we get ||r (t) − r|| 2 2 ≤ ε 2 implying the ℓ 2 -error is ≤ ε. If the minimum is the second term, then we get ||r (t) − r|| 2 ≤ ε √ 2ρ < ε since ρ ≥ 1.

New Lower Bound on the KL-Divergence
We now establish a new lower bound on KL-divergence which yields (KL vs ℓ 1 /ℓ 2 ) as a corollary.
Proof of Lemma 3.1: We need the following fact which follows from calculus; we provide a proof later for completeness.
Define η i := qi−pi pi . Note that A θ = {i : η i > θ} and B θ is the rest. We can write the KL-divergence as follows For i ∈ A θ , since η i > θ, we upper bound (1 + η i ) ≤ e a θ ηi using Fact 3.2. For i ∈ B θ , that is η i ≤ θ, we upper bound (1 + η i ) ≤ e ηi−b θ η 2 i using Fact 3.2. Lastly, we note i p i η i = 0 since p, q both sum to 1, implying i∈B θ p i η i = − i∈A θ p i η i . Putting all this in the definition above we get The proof of inequality (3) follows by noting that b θ = 1−a θ θ .
Proof of Lemma 3.2. The proof of both facts follow by proving non-negativity of the relevant function in the relevant interval. Recall a θ = ln(1 + θ)/θ and b θ = 1 θ (1 − a θ ). We start with the following three inequalities about the log-function.
Consider the function g(t) = e t(1−bt) − (1 + t). Note that g(0) = g(θ) = 0. We break the argument in two parts: we argue that g(t) is strictly positive for all t ≤ 0, and that g(t) is strictly positive for t ∈ (0, θ). This will prove the second part of Fact 3.2.
For the second part, we first note that g ′ (θ) < 0 since b > 1 2(1+θ) . That is, g is strictly decreasing at θ. On the other hand g is increasing at θ. To see this, looking at g ′ is not enough since g ′ (0) = 0. However, g ′′ (0) > 0 since b < 1/2. This means that 0 is a strict (local) minimum for g implying g is increasing at 0. In sum, g vanishes at 0 and θ, and is increasing at 0 and decreasing at θ. This means that if g does vanish at some r ∈ (0, θ), then it must vanish once again in [r, θ) for the it to be decreasing at θ. In particular, g ′ must vanish three times in (0, θ) and thus four times in [0, θ) since g ′ (0) = 0. This in turn implies g ′′ vanishes three times in [0, θ) which is a contradiction since g ′′ is a quadratic in t multiplied by a positive term.