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

Given a non-negative n × m real matrix A, the 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 δ, and then connect the KL-divergence with ℓ_1 and ℓ_2 distances. For ℓ_1, we can use Pinsker's inequality. For ℓ_2, we develop a strengthening of Pinsker's inequality, called (KL vs ℓ_1/ℓ_2) in the paper, which lower bounds the KL-divergence by a combination of ℓ_1 and ℓ_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 ℓ_2-error.

• 20 publications
• 20 publications
02/15/2022

### A short note on an inequality between KL and TV

The goal of this short note is to discuss the relation between Kullback–...
06/29/2019

### A New Lower Bound for Kullback-Leibler Divergence Based on Hammersley-Chapman-Robbins Bound

In this paper, we derive a useful lower bound for the Kullback-Leibler d...
11/04/2019

### Proximal Langevin Algorithm: Rapid Convergence Under Isoperimetry

We study the Proximal Langevin Algorithm (PLA) for sampling from a proba...
12/11/2020

### Scaling positive random matrices: concentration and asymptotic convergence

It is well known that any positive matrix can be scaled to have prescrib...
10/10/2021

### Convergence of Random Reshuffling Under The Kurdyka-Łojasiewicz Inequality

We study the random reshuffling (RR) method for smooth nonconvex optimiz...
06/01/2022

### Split-kl and PAC-Bayes-split-kl Inequalities

We present a new concentration of measure inequality for sums of indepen...
03/19/2017

### Bernoulli Rank-1 Bandits for Click Feedback

The probability that a user will click a search result depends both on i...

## 1 Introduction

In the matrix scaling problem one is given an non-negative 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 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 [16] 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 [5] and the Iterative Proportional Fitting Procedure [33]. 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,

 A(0)ij:=Aij∑ni=1Aij⋅cj      ∀t≥0,   B(t)ij:=A(t)ij∑mj=1A(t)ij⋅ri,   A(t+1)ij:=B(t)ij∑ni=1B(t)ij⋅cj

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 [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 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

 ℓ1-error:  error1(At):=||r(t)−r||1             ℓ2-error:  error2(At):=||r(t)−r||2

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 Sinkhorn-Knopp algorithm

1. in time returns a matrix or with -error .

2. 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.

###### Corollary 2.

Given a matrix which is -scalable, and any , the Sinkhorn-Knopp algorithm

1. in time returns a matrix or with -error .

2. in time returns a matrix or with -error .

Here is the maximum number of non-zeros 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 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 [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 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.

###### 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 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 .

Thus the approximate perfect matching problem in bipartite graphs is in NC for as small as polylogarithmic in . This is not a new result and can indeed be obtained from the works on parallel algorithms for packing-covering LPs [22, 40, 3, 24], but the Sinkhorn-Knopp algorithm is arguably simpler.

### 1.1 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 [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 Sinkhorn-Knopp 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 Sinkhorn-Knopp 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 analysis222 [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 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 [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 Allen-Zhu 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 Sinkhorn-Knopp Algorithm

There have been many approaches (see Idel [16], 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

 D(M,N):=1h∑1≤i≤n ∑1≤j≤mMijln(MijNij) (1)

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 [8] (also known as -projections [11]) 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 [16], 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 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 KL-divergence between two probability distributions

is defined as . The following theorem gives the convergence time for the KL-divergence.

###### Theorem 4.

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 .

###### Proof.

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 .

###### Proof.

By definition,

 D(Z,A(t))=1hm∑j=1n∑i=1Zijln⎛⎜⎝ZijA(t)ij⎞⎟⎠=1hm∑j=1cjn∑i=1Zijcjln⎛⎜⎝ZijA(t)ij⎞⎟⎠

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, [34]) 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 . ∎

###### Lemma 2.2.
 D(Z,A(t))−D(Z,B(t))=DKL(πr||π(t)r)   and   D(Z,B(t))−D(Z,A(t+1))=DKL(πc||π(t)c)
###### Proof.

The LHS of the first equality is simply

 1hm∑j=1n∑i=1Zijln⎛⎜⎝B(t)ijA(t)ij⎞⎟⎠ = 1hm∑j=1n∑i=1Zijln⎛⎝rir(t)i⎞⎠ = 1hn∑i=1ln⎛⎝rir(t)i⎞⎠m∑j=1Zij = n∑i=1(rih)⋅ln⎛⎝ri/hr(t)i/h⎞⎠

since . The last summand is precisely . The other equation follows analogously. ∎

The above two lemmas easily imply the theorem. If for all , both and , then substituting in Lemma 2.2 and summing we get contradicting Lemma 2.1. ∎

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 ,

 DKL(p||q)≥12⋅||p−q||21 (Pinsker)
###### Proof of Theorem 1, Part 1.

Apply (Pinsker) on the vectors and to get

 DKL(πr||π(t)r)≥12h2||r(t)−r||21

Set and apply Theorem 4. In time we would get a matrix with which from the above inequality would imply . ∎

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.

###### Lemma 2.3.

Given any pair of probability distributions over a finite domain, define and . Then,

 DKL(p||q)≥(1−ln2)⋅(∑i∈A|qi−pi|+∑i∈B(qi−pi)2pi) (KL vs ℓ1/ℓ2)
###### Proof of Theorem 1, Part 2.

We apply Lemma 2.3 on the vectors and . Lemma 2.3 gives us

 DKL(πr||π(t)r) ≥ C⋅⎛⎝1h∑i∈A|r(t)i−ri| + 1h∑i∈B(r(t)i−ri)2ri⎞⎠ ≥ Ch(∑i∈A|r(t)i−ri| + 1ρ∑i∈B(r(t)i−ri)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

 (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 KL-Divergence

We now establish a new lower bound on KL-divergence 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 KL-divergence as follows

 DKL(p||q):=n∑i=1piln(pi/qi)=−n∑i=1piln(1+ηi)

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 non-negativity of the relevant function in the relevant interval. Recall and . We start with the following three inequalities about the log-function.

 For all z>0,    z+z2/2>(1+z)ln(1+z)>z     and     ln(1+z)>z−z2/2 (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 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

 ||p−q||21=(∑i∈A|qi−pi|+∑i∈B|qi−pi|)2≤2(∑i∈A|qi−pi|)2+2(∑i∈Bpi|qi−pi|pi)2

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

 DKL(p||q)≥(1−ln2)2||p−q||21

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.

The following inequality is known (see Reiss [30] p 99, Pollard [28] Chap 3.3, or the webpage [39] for a proof).

 For any two distributions p,q,   DKL(p||q)≥D2Hellinger(p,q) (KL-vs-Hellinger)

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,

 ∑i∈A(qi−pi)=∑i∈A(√qi+√pi)(√qi−√pi)=Θ(∑i∈A(√qi−√pi)2)

For any , let where . Via a Taylor series expansion it is not hard to check in this range of . Observing that

 pi(qi−pipi)2=η2ipi       and       (√pi−√qi)2=2pi(1+ηi2−√1+ηi)

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 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.
• [4] 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.
• [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. I-divergence 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 deflation-inflation 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 e-prints, 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 nonnegative-matrix 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 Twenty-fifth 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 11-15, 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. Springer-Verlag, 1989.
• [31] 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.
• [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 11-15, 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. [Online; accessed 15-January-2018].
• [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.