In the matrix scaling problem one is given an non-negative matrix
, and positive integer vectorsand 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 all-ones 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  for more examples.
One of the most natural algorithms for the matrix scaling problem is the following Sinkhorn-Knopp algorithm [36, 37], which is known by many names including the RAS method  and the Iterative Proportional Fitting Procedure . The algorithm starts off by multiplicatively scaling all the columns by the columns-sum times to get a matrix with column-sums . Subsequently, for , it obtains the by scaling each row of by the respective row-sum 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 non-zero entries of . Furthermore, it has been known for almost five decades [36, 37, 14, 38] that if is -scalable then the above algorithm asymptotically111Computationally, 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 . 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 column-sums 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 row-sums of . Similarly, we define to be the vector of the column-sums 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.
Given a matrix which is -scalable, and any , the Sinkhorn-Knopp algorithm
in time returns a matrix or with -error .
in time returns a matrix or with -error .
Here , , , and is the maximum number of non-zeros in any column of .
For the special case of and , we get the following as a corollary.
Given a matrix which is -scalable, and any , the Sinkhorn-Knopp algorithm
in time returns a matrix or with -error .
in time returns a matrix or with -error .
Here is the maximum number of non-zeros in any column of .
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 . 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 ), 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.
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 Cauchy-Schwarz. 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 , 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 matrixin the support of . This suggests the algorithm of running the Sinkhorn-Knopp 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.
If we find a column (or row) stochastic matrix in the support of such that , then has a matching of size .
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 max-degree and an , -iterations of Sinkhorn-Knopp algorithm suffice to distinguish between the case when has a perfect matching and the case when the largest matching in has size at most .
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  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 Sinkhorn-Knopp algorithm in Remark 1.1. For the special case of strictly positive matrices, better rates are known. Kalantari and Khachiyan  showed that for positive matrices and the -scaling problem, the Sinkhorn-Knopp algorithm obtains error in -iterations; this result was extended to the general matrix scaling problem by Kalantari et al . In a different track, Franklin and Lorenz  show that in fact the dependence on can be made logarithmic, and thus the algorithm has “linear convergence”, however their analysis222  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 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  describe a method based on the ellipsoid algorithm which runs in time . Nemirovskii and Rothblum  describe a method with running time . The first strongly polynomial time approximation scheme (with no dependence on ) was due to Linial, Samoridnitsky, and Wigderson  who gave a time algorithm. Rote and Zachariasen  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  give time algorithm, while Allen-Zhu et al  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  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 Sinkhorn-Knopp Algorithm
There have been many approaches (see Idel , Section 3 for a discussion) towards analyzing the Sinkhorn-Knopp algorithm including convex optimization and log-barrier methods [17, 20, 23, 6], non-linear Perron-Frobenius 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 non-negative matrices and let us define the Kullback-Leibler divergence333The KL-divergence is normally stated between two distributions and doesn’t have the factor. Also the logarithms are usually base . between and as follows
with the convention that the summand is zero if both and are , and is if and . Let be the set of matrices whose row-sums 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 KL-divergence  (also known as -projections ) which alternately finds the matrices in and 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 , 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 or the error as described in the previous section, but the rather the KL-divergence between the normalized vectors as described below.
Let be the
-dimensional probability vector whoseth entry is ; similarly define the -dimensional vector . Let denote the -dimensional probability vector with the th entry being ; similarly define
. Recall that the KL-divergence between two probability distributionsis defined as . The following theorem gives the convergence time for the KL-divergence.
If the matrix is -scalable, then for any there is a with either or . Recall, , , and is the maximum number of non-zeros in any column of .
Let be a matrix with row-sums and column-sums for diagonal matrices . Recall is the matrix obtained by column-scaling . Note that the minimum non-zero entry of is .
and for all .
For a fixed , the vectors and are probability vectors, and therefore the above is a sum of -weighted KL-divergences which is always non-negative. For the upper bound, one can use the fact (Inequality 27, ) that for any two distributions and , where is the smallest non-zero entry of . For our purpose, we note that the minimum non-zero probability of the distribution being . Therefore, the second summand is at most giving us . ∎
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 KL-divergence 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 ,
Proof of Theorem 1, Part 1.
To prove Part 2, we need a way to relate the 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. 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.
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 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 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 KL-Divergence
We now establish a new lower bound on KL-divergence which yields (KL vs ) as a corollary.
Let and be two distributions over a finite -element universe. For any fixed , define the sets and . Then we have the following inequality
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.
Given any , define and . Then,
Define . Note that and is the rest. We can write the KL-divergence as follows
The proof of inequality (3) follows by noting that . ∎
Proof of Lemma 3.2.
The proof of both facts follow by proving non-negativity of the relevant function in the relevant interval. Recall and . We start with the following three inequalities about the log-function.
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 well-known inequalities
We connect (KL vs ) with two well known lower bounds on the KL-Divergence. 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 Cauchy-Schwarz, 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 KL-divergence between two distributions is also at least the Hellinger distance between them. Before proceeding, let us define this distance.
It seems natural to compare the RHS of (KL vs ) and (KL-vs-Hellinger) (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 (KL-vs-Hellinger) 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 .
We thank Daniel Dadush for asking the connection between our inequality and Hellinger distance, and Jonathan Weed for letting us know of .
-  S. Aaronson. Quantum computing and hidden variables. Phys. Rev. A, 71:032325, Mar 2005.
-  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.
Z. Allen Zhu and L. Orecchia.
Using optimization to break the epsilon barrier: A faster and simpler width-independent algorithm for solving positive linear programs in parallel.In Proceedings of the Twenty-Sixth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA, San Diego, CA, USA, pages 1439–1456, 2015.
-  J. Altschuler, J. Weed, and P. Rigollet. Near-linear 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, 4-9 December 2017, Long Beach, CA, USA, pages 1961–1971, 2017.
-  M. Bacharach. Estimating nonnegative matrices from marginal data. International Economic Review, 6(3):294–310, 1965.
-  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.
-  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.
-  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.
-  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.
-  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.
-  I. Csiszar. I-divergence geometry of probability distributions and minimization problems. Ann. Probab., 3(1):146–158, 02 1975.
-  I. Csiszar. A geometric interpretation of Darroch and Ratcliff’s generalized iterative scaling. Ann. Statist., 17(3):1409–1413, 09 1989.
-  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.
-  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.
-  L. Gurvits and P. N. Yianilos. The deflation-inflation method for certain semidefinite programming and maximum determinant completion problems. Technical report, NEC Research Institute, 4 Independence Way, Princeton, NJ 08540, 1998.
-  M. Idel. A review of matrix scaling and Sinkhorn’s normal form for matrices and positive maps. ArXiv e-prints, Sept. 2016.
-  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.
-  B. Kalantari and L. Khachiyan. On the complexity of nonnegative-matrix scaling. Linear Algebra and its Applications, 240:87 – 103, 1996.
-  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.
-  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.
-  N. Linial, A. Samorodnitsky, and A. Wigderson. A deterministic strongly polynomial algorithm for matrix scaling and approximate permanents. Combinatorica, 20(4):545–568, 2000.
M. Luby and N. Nisan.
A parallel approximation algorithm for positive linear programming.
Proceedings of the Twenty-fifth Annual ACM Symposium on Theory of Computing, STOC ’93, pages 448–457, New York, NY, USA, 1993. ACM.
-  S. M. Macgill. Theoretical properties of biproportional matrix adjustments. Environment and Planning A, 9(6):687–701, 1977.
-  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 11-15, 2016, Rome, Italy, pages 52:1–52:14, 2016.
-  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.
-  A. Nemirovskii and U. Rothblum. On complexity of matrix scaling. Linear Algebra and its Applications., 302:435–460, 1999.
-  J. d. D. Ortúzar and L. G. Willumsen. Modelling Transport. John Wiley & Sons, Ltd, 2011.
-  D. Pollard. A User’s Guide to Measure Theoretic Probability. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2001.
-  T. Raghavan. On pairs of multidimensional matrices. Linear Algebra and its Applications, 62:263 – 268, 1984.
-  R.-D. Reiss. Approximate distributions of order statistics. Springer-Verlag, 1989.
-  G. Rote and M. Zachariasen. Matrix scaling by network flow. In Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’07, pages 848–854, Philadelphia, PA, USA, 2007. Society for Industrial and Applied Mathematics.
-  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.
-  L. Ruschendorf. Convergence of the iterative proportional fitting procedure. Ann. Statist., 23(4):1160–1174, 1995.
-  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 11-15, 2015, pages 214–218. IEEE, 2015.
-  E. Schrödinger. Über die umkehrung der naturgesetze. Preuss. Akad. Wiss., Phys.-Math. Kl., 1931.
-  R. Sinkhorn. Diagonal equivalence to matrices with prescribed row and column sums. The American Mathematical Monthly, 74(4):402–405, 1967.
-  R. Sinkhorn and P. Knopp. Concerning nonnegative matrices and doubly stochastic matrices. Pacific J. Math., 21(2):343–348, 1967.
-  G. W. Soules. The rate of convergence of sinkhorn balancing. Linear Algebra and its Applications, 150:3 – 40, 1991.
-  Xi’an. Statistics stack exchange. https://stats.stackexchange.com/questions/130432/differences-between-bhattacharyya-distance-and-kl-divergence, Dec 27, 2014. [Online; accessed 15-January-2018].
-  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.