1 Introduction
In the matrix scaling problem one is given an nonnegative matrix
, and positive integer vectors
and with the same norm . The objective is to determine if there exist diagonal matrices and such that the th row of the matrix sums to for all and the th column of sums to for all . Of special importance is the case when and , the dimensional allones vector – the matrix scaling problem wishes to scale the rows and columns of to make it doubly stochastic. This problem arises in many different areas ranging from transportation planning [13, 27] to quantum mechanics [35, 1]; we refer the reader to a recent comprehensive survey by Idel [16] for more examples.One of the most natural algorithms for the matrix scaling problem is the following SinkhornKnopp algorithm [36, 37], which is known by many names including the RAS method [5] and the Iterative Proportional Fitting Procedure [33]. The algorithm starts off by multiplicatively scaling all the columns by the columnssum times to get a matrix with columnsums . Subsequently, for , it obtains the by scaling each row of by the respective rowsum times , and obtain by scaling each column of by the respective column sums time . More precisely,
The above algorithm is simple and easy to implement and each iteration takes , the number of nonzero entries of . Furthermore, it has been known for almost five decades [36, 37, 14, 38] that if is scalable then the above algorithm asymptotically^{1}^{1}1Computationally, 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 [32]. The matrix is with . The unique and matrices are and , respectively, giving . converges to a right solution. More precisely, given , there is some finite by which one obtains a matrix which is “close to having row and columnsums and ”.
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 , let denote the vector of rowsums of . Similarly, we define to be the vector of the columnsums of . Note that for all . The error of the matrix (the error of matrix similarly defined) is
In this note, we give simple convergence analysis for both error norms. Our result is the following.
Theorem 1.
Given a matrix which is scalable, and any , the SinkhornKnopp algorithm

in time returns a matrix or with error .

in time returns a matrix or with error .
Here , , , and is the maximum number of nonzeros in any column of .
For the special case of and , we get the following as a corollary.
Corollary 2.
Given a matrix which is scalable, and any , the SinkhornKnopp algorithm

in time returns a matrix or with error .

in time returns a matrix or with error .
Here is the maximum number of nonzeros in any column of .
Remark 1.1.
To our knowledge, the error hasn’t been explicitly studied in the literature (but see last paragraph of Section 1.1), although for small the same can be deduced from previous papers on matrix scaling [21, 15, 20, 17]. One of our main motivations to look at error arose from the connections to perfect matchings in bipartite graphs as observed by Linial, Samorodnitsky and Wigderson [21]. For the error, which is the better studied notion in the matrix scaling literature, the best analysis is due to Kalantari et al [19, 20]. They give a upper bound on the number of iterations for the general problem, and for the special case when and the square matrix has positive permanent (see [19]), they give a upper bound. Thus, for scaling, they get the same result as in Corollary 2. We get a quadratic improvement on in the general case, and we think our proof is more explicit and simpler.
Remark 1.2.
Both parts of Theorem 1 and Corollary 2 are interesting in certain regimes of error. When the error is “small” (say, ) so that , then statement 2 of Corollary 2 implies statement 1 by CauchySchwarz. However, this breaks down when is “large” (say for some constant ). In that case, statement 1 implies that in iterations, the error is , but Statement 2 only implies that in iterations, the norm is . This “large 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 [21], to the existence of perfect matchings in bipartite graphs. Let be the adjacency matrix of a bipartite graph with iff . If
has a perfect matching, then clearly there is a doubly stochastic matrix
in the support of . This suggests the algorithm of running the SinkhornKnopp algorithm to , and the following claim suggests when to stop. Note that each iteration can be run in parallel time with processors where is the number of edges.Lemma 1.1.
If we find a column (or row) stochastic matrix in the support of such that , then has a matching of size .
Proof.
Suppose is column stochastic. Given , consider . On the other hand, . Therefore, for every , . The claim follows by approximate Hall’s theorem. ∎
Corollary 3 (Fast Parallel Approximate Matchings).
Given a bipartite graph of maxdegree and an , iterations of SinkhornKnopp algorithm suffice to distinguish between the case when has a perfect matching and the case when the largest matching in has size at most .
1.1 Perspective
As mentioned above, the matrix scaling problem and in particular the SinkhornKnopp algorithm has been extensively studied over the past 50 years. We refer the reader to Idel’s survey [16] 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 , analysis for the SinkhornKnopp algorithm in Remark 1.1. For the special case of strictly positive matrices, better rates are known. Kalantari and Khachiyan [17] showed that for positive matrices and the scaling problem, the SinkhornKnopp algorithm obtains error in iterations; this result was extended to the general matrix scaling problem by Kalantari et al [20]. In a different track, Franklin and Lorenz [14] show that in fact the dependence on can be made logarithmic, and thus the algorithm has “linear convergence”, however their analysis^{2}^{2}2 [14] never make the base of the logarithm explicit, but their proof shows it can be as large as . has a polynomial dependence of . All these results use the positivity crucially and seem to break down even with one entry.
The SinkhornKnopp 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 [18] describe a method based on the ellipsoid algorithm which runs in time . Nemirovskii and Rothblum [26] describe a method with running time . The first strongly polynomial time approximation scheme (with no dependence on ) was due to Linial, Samoridnitsky, and Wigderson [21] who gave a time algorithm. Rote and Zachariasen [31] reduced the matrix scaling problem to flow problems to give a time algorithms for the matrix scaling problem. To compare, we should recall that Theorem 1 shows that our algorithm runs in time time.
Very recently, two independent works obtain vastly improved running times for matrix scaling. Cohen et al [10] give time algorithm, while AllenZhu et al [2] give a 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. After the first version of this paper was made public, we were pointed out another recent paper by Altschuler, Weed and Rigollet [4] who also study the error and obtain the same result as part 1 of our Theorem. Indeed their proof techniques are very similar to what we use to prove part 1.
2 Entropy Minimization Viewpoint of the SinkhornKnopp Algorithm
There have been many approaches (see Idel [16], Section 3 for a discussion) towards analyzing the SinkhornKnopp algorithm including convex optimization and logbarrier methods [17, 20, 23, 6], nonlinear PerronFrobenius theory [25, 38, 14, 9, 17], topological methods [29, 7], connections to the permanent [21, 19], and the entropy minimization method [8, 11, 12, 15] which is what we use for our analysis.
We briefly describe the entropy minimization viewpoint. Given two nonnegative matrices and let us define the KullbackLeibler divergence^{3}^{3}3The KLdivergence is normally stated between two distributions and doesn’t have the factor. Also the logarithms are usually base . between and as follows
(1) 
with the convention that the summand is zero if both and are , and is if and . Let be the set of matrices whose rowsums are and let be the set of matrices whose column sums are . Given matrix suppose we wish to find the matrix . One algorithm for this is to use the method of alternate projections with respect to the KLdivergence [8] (also known as projections [11]) which alternately finds the matrices in and closest in the KLdivergence 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 [16], Observation 3.17 for a proof) that the above alternate projection algorithm is precisely the SinkhornKnopp algorithm. Therefore, at least in this sense, the right metric to measure the distance to optimality is not the or the error as described in the previous section, but the rather the KLdivergence between the normalized vectors as described below.
Let be the
dimensional probability vector whose
th entry is ; similarly define the dimensional vector . Let denote the dimensional probability vector with the th entry being ; similarly define. Recall that the KLdivergence between two probability distributions
is defined as . The following theorem gives the convergence time for the KLdivergence.Theorem 4.
If the matrix is scalable, then for any there is a with either or . Recall, , , and is the maximum number of nonzeros in any column of .
Proof.
Let be a matrix with rowsums and columnsums for diagonal matrices . Recall is the matrix obtained by columnscaling . Note that the minimum nonzero entry of is .
Lemma 2.1.
and for all .
Proof.
By definition,
For a fixed , the vectors and are probability vectors, and therefore the above is a sum of weighted KLdivergences which is always nonnegative. For the upper bound, one can use the fact (Inequality 27, [34]) that for any two distributions and , where is the smallest nonzero entry of . For our purpose, we note that the minimum nonzero probability of the distribution being . Therefore, the second summand is at most giving us . ∎
Lemma 2.2.
Proof.
The LHS of the first equality is simply
since . The last summand is precisely . The other equation follows analogously. ∎
Theorem 1 follows from Theorem 4 using connections between the KLdivergence and the and 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 ,
(Pinsker) 
Proof of Theorem 1, Part 1.
To prove Part 2, we need a way to relate the norm and the KLdivergence. 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. This 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.
Lemma 2.3.
Given any pair of probability distributions over a finite domain, define and . Then,
(KL vs ) 
Proof of Theorem 1, Part 2.
where . If the second summand in the parenthesis of the RHS is , then we get . Otherwise, we have , where we used the weak fact that the sum of some positive numbers is at least the squareroot of the sum of their squares. In any case, we get the following
(2) 
To complete the proof of part 2 of Theorem 1, set and apply Theorem 4. In time we would get a matrix with . If the minimum of the RHS of (2) is the first term, then we get implying the error is . If the minimum is the second term, then we get since . ∎
3 New Lower Bound on the KLDivergence
We now establish a new lower bound on KLdivergence which yields (KL vs ) as a corollary.
Lemma 3.1.
Let and be two distributions over a finite element universe. For any fixed , define the sets and . Then we have the following inequality
(3) 
When , we get (KL vs ).
Proof of Lemma 3.1:.
We need the following fact which follows from calculus; we provide a proof later for completeness.
Lemma 3.2.
Given any , define and . Then,

For ,

For ,
Define . Note that and is the rest. We can write the KLdivergence as follows
For , since , we upper bound using Lemma 3.2. For , that is , we upper bound using Lemma 3.2. Lastly, we note since both sum to , implying . Putting all this in the definition above we get
The proof of inequality (3) follows by noting that . ∎
Proof of Lemma 3.2.
The proof of both facts follow by proving nonnegativity of the relevant function in the relevant interval. Recall and . We start with the following three inequalities about the logfunction.
(4) 
The third inequality in (4) implies and thus, . The first inequality in (4) implies which in turn implies . For brevity, henceforth let us lose the subscript on and .
Consider the function . Note that which is increasing in since . So, for any , we have , by the second inequality in (4). Therefore, is increasing when . The first part of Fact 3.2 follows since by definition of .
Consider the function . Note that . We break the argument in two parts: we argue that is strictly positive for all , and that is strictly positive for . This will prove the second part of Fact 3.2.
The first derivative is and the second derivative is . Since , we have , and thus for , . Therefore, is strictly increasing for . However, , and so for all . This implies is strictly decreasing in the interval . Noting , we get for all . This completes the first part of the argument.
For the second part, we first note that since . That is, is strictly decreasing at . On the other hand is increasing at . To see this, looking at is not enough since . However, since . This means that is a strict (local) minimum for implying is increasing at . In sum, vanishes at and , and is increasing at and decreasing at . This means that if does vanish at some , then it must vanish once again in for the it to be decreasing at . In particular, must vanish three times in and thus four times in since . This in turn implies vanishes three times in which is a contradiction since is a quadratic in multiplied by a positive term.
We end by proving (4). This also follows the same general methodology. Define and . Differentiating, we get for all , and for all . Thus, is increasing, and is decreasing, in . The first two inequalities of (4) follow since . To see the third inequality, define and observe which is if . Thus is strictly increasing, and the third inequality of (4) follows since .
∎
3.1 Comparison with other wellknown inequalities
We connect (KL vs ) with two well known lower bounds on the KLDivergence. First we compare with Pinsker’s inequality (Pinsker). To see that (KL vs ) generalizes (Pinsker) with a weaker constant, note that
The first parenthetical term above, since it is , is at most the first summation in the parenthesis of (KL vs ). The second parenthetical term above, by CauchySchwarz, is at most the second summation in the parenthesis of (KL vs ). Thus (KL vs ) implies
On the other hand, the RHS of (KL vs ) can be much larger than that of (Pinsker). For instance, suppose for all , , and for , . The RHS of (Pinsker) is while that of (KL vs ) is which is the correct order of magnitude for .
The KLdivergence between two distributions is also at least the Hellinger distance between them. Before proceeding, let us define this distance.
The following inequality is known (see Reiss [30] p 99, Pollard [28] Chap 3.3, or the webpage [39] for a proof).
(KLvsHellinger) 
It seems natural to compare the RHS of (KL vs ) and (KLvsHellinger) (we thank Daniel Dadush for bringing this to our attention).
As the subsequent calculation shows, the RHS of (KL vs ) is in fact . In particular, this implies one can obtain (by reverse engineering the argument below) part 2 of Theorem 2 via the application of (KLvsHellinger) as well.
For the set , we know . Therefore,
For any , let where . Via a Taylor series expansion it is not hard to check in this range of . Observing that
we get that the RHS of (KL vs ) is .
Acknowledgements
We thank Daniel Dadush for asking the connection between our inequality and Hellinger distance, and Jonathan Weed for letting us know of [4].
References
 [1] S. Aaronson. Quantum computing and hidden variables. Phys. Rev. A, 71:032325, Mar 2005.
 [2] Z. Allen Zhu, Y. Li, R. Oliveira, and A. Wigderson. Much faster algorithms for matrix scaling. 58th IEEE Annual Symposium on Foundations of Computer Science, FOCS, 2017.

[3]
Z. Allen Zhu and L. Orecchia.
Using optimization to break the epsilon barrier: A faster and simpler widthindependent algorithm for solving positive linear programs in parallel.
In Proceedings of the TwentySixth Annual ACMSIAM Symposium on Discrete Algorithms, SODA, San Diego, CA, USA, pages 1439–1456, 2015.  [4] J. Altschuler, J. Weed, and P. Rigollet. Nearlinear time approximation algorithms for optimal transport via sinkhorn iteration. In I. Guyon, U. von Luxburg, S. Bengio, H. M. Wallach, R. Fergus, S. V. N. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, 49 December 2017, Long Beach, CA, USA, pages 1961–1971, 2017.
 [5] M. Bacharach. Estimating nonnegative matrices from marginal data. International Economic Review, 6(3):294–310, 1965.
 [6] H. Balakrishnan, I. Hwang, and C. J. Tomlin. Polynomial approximation algorithms for belief matrix maintenance in identity management. In 2004 43rd IEEE Conference on Decision and Control (CDC) (IEEE Cat. No.04CH37601), volume 5, pages 4874–4879 Vol.5, Dec 2004.
 [7] R. Bapat and T. Raghavan. An extension of a theorem of Darroch and Ratcliff in loglinear models and its application to scaling multidimensional matrices. Linear Algebra and its Applications, 114:705 – 715, 1989. Special Issue Dedicated to Alan J. Hoffman.
 [8] L. Bregman. The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming. USSR Computational Mathematics and Mathematical Physics, 7(3):200 – 217, 1967.
 [9] R. A. Brualdi, S. V. Parter, and H. Schneider. The diagonal equivalence of a nonnegative matrix to a stochastic matrix. Journal of Mathematical Analysis and Applications, 16(1):31 – 50, 1966.
 [10] M. B. Cohen, A. Madry, D. Tsipras, and A. Vladu. Matrix scaling and balancing via box constrained newton’s method and interior point methods. 58th IEEE Annual Symposium on Foundations of Computer Science, FOCS, 2017.
 [11] I. Csiszar. Idivergence geometry of probability distributions and minimization problems. Ann. Probab., 3(1):146–158, 02 1975.
 [12] I. Csiszar. A geometric interpretation of Darroch and Ratcliff’s generalized iterative scaling. Ann. Statist., 17(3):1409–1413, 09 1989.
 [13] W. E. Deming and F. F. Stephan. On a least squares adjustment of a sampled frequency table when the expected marginal totals are known. Ann. Math. Statist., 11(4):427–444, 12 1940.
 [14] J. Franklin and J. Lorenz. On the scaling of multidimensional matrices. Linear Algebra and its Applications, 114:717 – 735, 1989. Special Issue Dedicated to Alan J. Hoffman.
 [15] L. Gurvits and P. N. Yianilos. The deflationinflation method for certain semidefinite programming and maximum determinant completion problems. Technical report, NEC Research Institute, 4 Independence Way, Princeton, NJ 08540, 1998.
 [16] M. Idel. A review of matrix scaling and Sinkhorn’s normal form for matrices and positive maps. ArXiv eprints, Sept. 2016.
 [17] B. Kalantari and L. Khachiyan. On the rate of convergence of deterministic and randomized RAS matrix scaling algorithms. Oper. Res. Lett., 14(5):237–244, 1993.
 [18] B. Kalantari and L. Khachiyan. On the complexity of nonnegativematrix scaling. Linear Algebra and its Applications, 240:87 – 103, 1996.
 [19] B. Kalantari, I. Lari, F. Ricca, and B. Simeone. On the complexity of general matrix scaling and entropy minimization via the RAS algorithm. Technical Report, n.24, Department of Statistics and Applied probability, La Sapienza University, Rome, 2002.
 [20] B. Kalantari, I. Lari, F. Ricca, and B. Simeone. On the complexity of general matrix scaling and entropy minimization via the RAS algorithm. Math. Program., 112(2):371–401, 2008.
 [21] N. Linial, A. Samorodnitsky, and A. Wigderson. A deterministic strongly polynomial algorithm for matrix scaling and approximate permanents. Combinatorica, 20(4):545–568, 2000.

[22]
M. Luby and N. Nisan.
A parallel approximation algorithm for positive linear programming.
In
Proceedings of the Twentyfifth Annual ACM Symposium on Theory of Computing
, STOC ’93, pages 448–457, New York, NY, USA, 1993. ACM.  [23] S. M. Macgill. Theoretical properties of biproportional matrix adjustments. Environment and Planning A, 9(6):687–701, 1977.
 [24] M. W. Mahoney, S. Rao, D. Wang, and P. Zhang. Approximating the Solution to Mixed Packing and Covering LPs in Parallel Time. In 43rd International Colloquium on Automata, Languages, and Programming, ICALP 2016, July 1115, 2016, Rome, Italy, pages 52:1–52:14, 2016.
 [25] M. Menon. Reduction of a matrix with positive elements to a doubly stochastic matrix. Proceedings of the American Mathematical Society, 18(2):244–247, 1967.
 [26] A. Nemirovskii and U. Rothblum. On complexity of matrix scaling. Linear Algebra and its Applications., 302:435–460, 1999.
 [27] J. d. D. Ortúzar and L. G. Willumsen. Modelling Transport. John Wiley & Sons, Ltd, 2011.
 [28] D. Pollard. A User’s Guide to Measure Theoretic Probability. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2001.
 [29] T. Raghavan. On pairs of multidimensional matrices. Linear Algebra and its Applications, 62:263 – 268, 1984.
 [30] R.D. Reiss. Approximate distributions of order statistics. SpringerVerlag, 1989.
 [31] G. Rote and M. Zachariasen. Matrix scaling by network flow. In Proceedings of the Eighteenth Annual ACMSIAM Symposium on Discrete Algorithms, SODA ’07, pages 848–854, Philadelphia, PA, USA, 2007. Society for Industrial and Applied Mathematics.
 [32] U. Rothblum and H. Schneider. Scalings of matrices which have prespecified row sums and column sums via optimization. Linear Algebra and its Applications, 114:737 – 764, 1989. Special Issue Dedicated to Alan J. Hoffman.
 [33] L. Ruschendorf. Convergence of the iterative proportional fitting procedure. Ann. Statist., 23(4):1160–1174, 1995.
 [34] I. Sason and S. Verdú. Upper bounds on the relative entropy and rényi divergence as a function of total variation distance for finite alphabets. In 2015 IEEE Information Theory Workshop  Fall (ITW), Jeju Island, South Korea, October 1115, 2015, pages 214–218. IEEE, 2015.
 [35] E. Schrödinger. Über die umkehrung der naturgesetze. Preuss. Akad. Wiss., Phys.Math. Kl., 1931.
 [36] R. Sinkhorn. Diagonal equivalence to matrices with prescribed row and column sums. The American Mathematical Monthly, 74(4):402–405, 1967.
 [37] R. Sinkhorn and P. Knopp. Concerning nonnegative matrices and doubly stochastic matrices. Pacific J. Math., 21(2):343–348, 1967.
 [38] G. W. Soules. The rate of convergence of sinkhorn balancing. Linear Algebra and its Applications, 150:3 – 40, 1991.
 [39] Xi’an. Statistics stack exchange. https://stats.stackexchange.com/questions/130432/differencesbetweenbhattacharyyadistanceandkldivergence, Dec 27, 2014. [Online; accessed 15January2018].
 [40] N. E. Young. Sequential and parallel algorithms for mixed packing and covering. In 42nd Annual Symposium on Foundations of Computer Science, FOCS, Las Vegas, Nevada, USA, pages 538–546, 2001.