Matrices with 0-1 Entries and Prescribed Row and Column Sums

ON THE NUMBER OF MATRICES AND A RANDOM MATRIX WITH PRESCRIBED ROW AND COLUMN SUMS AND 0-1 ENTRIES 2009 Abstract. We consider the set Σ(R, C) of all m × n matrices having 0-1 entries and prescribed row sums R = (r1 , . . . , rm ) and column sums C = (c1 , . . . , cn ). We prove an asymptotic estimate for the cardinality |Σ(R, C)| via the solution to a convex optimization problem. We show that if Σ(R, C) is sufficiently large, then a random matrix D ∈ Σ(R, C) sampled from the uniform probability measure in Σ(R, C) with high probability is close to a particular matrix Z = Z(R, C) that maximizes the sum of entropies of entries among all matrices with row sums R, column sums C and entries between 0 and 1. Similar results are obtained for 0-1 matrices with prescribed row and column sums and assigned zeros in some positions. 1. Introduction and main results Matrices with 0-1 entries and prescribed row and column sums is a classical object which appears in many branches of pure and applied mathematics. In combinatorics, such matrices encode hypergraphs with prescribed degrees of vertices and related structures, see, for example, [LW01]. In algebra, certain structural constants in the ring of symmetric functions and, consequently, in the representation theory of the symmetric and general linear groups are expressed as numbers of 0-1 matrices with prescribed row and column sums, see Chapter 1 of [Ma95]. In statistics, 0-1 matrices with prescribed row and column sums are known as binary contingency tables, see [C+05]. 1991 Mathematics Subject Classification. 05A16, 05C30, 60C05, 15A15, 15A52. Key words and phrases. 0-1 matrices, entropy, asymptotic estimates, integer flows. This research was partially supported by NSF Grant DMS 0400617 and a United States - Israel BSF grant 2006377. Typeset by AMS-TEX 1 Let R = (r1 , . . . , rm ) be a positive integer m-vector and let C = (c1 , . . . , cn ) be a positive integer n-vector such that m X i=1 0 < ri < n ri = n X cj = N and j=1 for i = 1, . . . , m and 0 < cj < m for j = 1, . . . , n. Let Σ(R, C) be the set of all m × n matrices (binary contingency tables) D = (dij ) such that n X dij = ri for i = 1, . . . , m, j=1 m X dij = cj for j = 1, . . . , n and i=1 dij ∈ {0, 1}. In words: Σ(R, C) is the set of 0-1 matrices with row sums R and column sums C. Vectors R and C are called margins of a matrix D ∈ Σ(R, C). Our first main result provides an estimate of the cardinality of Σ(R, C). (1.1) Theorem. Let us define the function   ! n m Y Y Y −c i  F (x, y) = x−r yj j   (1 + xi yj ) i for i=1 j=1 x = (x1 , . . . , xn ) and ij y = (y1 , . . . , yn ) and let α(R, C) = inf x1 ,... ,xm >0 y1 ,... ,yn >0 F (x, y). Then for the number |Σ(R, C)| of m × n zero-one matrices with row sums R and column sums C we have !  n cj  m n−ri Y Y cj (mn)! (n − ri )   α(R, C). α(R, C) ≥ |Σ(R, C)| ≥ (mn)mn (n − ri )! cj ! i=1 j=1 Let us estimate the ratio between the lower and the upper bounds for |Σ(R, C)| using Stirling’s formula √
s!s−s = e−s 2πs 1 + O(s−1 ) . Since e−mn m Y i=1  ! n Y en−ri  ecj  = 1, j=1 2 the “ e−s ” contributions from Stirling’s formula cancel each other out and we obtain α(R, C) ≥ |Σ(R, C)| ≥ (mn)−γ(m+n) α(R, C) for some absolute constant γ > 0. We note that in many interesting cases we have |Σ(R, C)| = 2Ω(mn) , see also Section 3.1, in which case the estimate of Theorem 1.1 captures the logarithmic order of |Σ(R, C)|. Let us substitute xi = esi , y = eti in F (x, y). Then ln F (x, y) = G(s, t), where m n X X X
G(s, t) = − ri s i − cj tj + ln 1 + esi +tj i=1 for s = (s1 , . . . , sm ) j=1 ij and t = (t1 , . . . , tn ) . One can observe that G(s, t) is a convex function on Rm ×Rn , hence to compute the infimum of G(s, t) one can use any of the efficient convex optimization algorithms, see, for example, [NN94]. Suppose that margins R, C are such that the set Σ(R, C) is not empty and let us consider Σ(R, C) as a finite probability space with the uniform measure. Let us pick a random matrix D ∈ Σ(R, C). What is D likely to look like? This question is of some interest to statistics: a binary contingency table D = (dij ) may represent certain statistical data (for example, dij may be equal to 1 or 0 depending on whether or not Darwin finches of the i-th species can be found on the j-th Galapagos island, as in [C+05]). One can condition on the row and column sums and ask what is special about a particular table D ∈ Σ(R, C), considering all tables in Σ(R, C) as equiprobable, see [C+05]. To answer this question we need to know what a random table D ∈ Σ(R, C) looks like. We prove that with high probability D is close to a particular matrix Z with row sums R and column sums C and entries between 0 and 1, which we call the maximum entropy matrix. (1.2) The maximum entropy matrix. For 0 ≤ x ≤ 1 let us consider the entropy function 1 1 . H(x) = x ln + (1 − x) ln x 1−x As is known, H is a strictly concave function with H(0) = H(1) = 0. For an m × n matrix X = (xij ) such that 0 ≤ xij ≤ 1 for all i, j, we define X H(X) = H (xij ) . ij Assume that Σ(R, C) is non-empty. Let us consider the polytope P(R, C) of matrices X = (xij ) such that m n X X xij = cj for j = 1, . . . , n and xij = ri for i = 1, . . . , m, j=1 i=1 0 ≤ xij ≤ 1 for all i, j. 3 Since H(X) is strictly concave, it attains a unique maximum Z = Z(R, C) on P(R, C), which we call the maximum entropy matrix with margins (R, C). For example, if all ri are equal, then by the symmetry argument we must have Z = (zij ) where zij = cj /m for all i, j. The following observation characterizes the maximum entropy matrix as the solution to the problem that is convex dual to the optimization problem of Theorem 1.1. (1.3) Lemma. Suppose that the polytope P(R, C) has a non-empty interior, that is, contains a matrix Y = (yij ) such that 0 < yij < 1 for all i, j. Then the infimum α(R, C) in Theorem 1.1 is attained at a particular point x∗ = (ξ1 , . . . , ξm ) and y∗ = (η1 , . . . , ηn ). For the maximum entropy matrix Z = (zij ) we have (1.3.1) zij = ξi ηj 1 + ξi ηj for all i, j and, moreover, (1.3.2) α(R, C) = eH(Z) . Conversely, if the infimum α(R, C) in Theorem 1.1 is attained at a certain point x∗ = (ξ1 , . . . , ξm ) and y∗ = (η1 , . . . , ηn ) then for the maximum entropy matrix Z = (zij ) equations (1.3.1) and (1.3.2) hold. The condition that the polytope P(R, C) has a non-empty interior is equivalent to the requirement that for every choice of 1 ≤ k ≤ m and 1 ≤ l ≤ n there is
0 0 0 a matrix D ∈ Σ(R, C), D = dij , such that d0kl = 0 and there is a matrix
D1 ∈ Σ(R, C), D1 = d1ij , such that d1kl = 1. One can take Y to be the average of all matrices D ∈ Σ(R, C). In other words, we require the set Σ(R, C) to be reasonably large. We also observe that if ri cj < N for all i, j (recall that N is the total sum of the matrix entries) one can choose yij = ri cj /N . We prove that with high probability a random matrix D ∈ Σ(R, C) is close to the maximum entropy matrix Z as far as sums over subsets of entries are concerned. For a subset n o S ⊂ (i, j) : i = 1, . . . , m, j = 1, . . . , n and an m × n matrix A = (aij ), let us denote X σS (A) = aij , (i,j)∈S the sum of the entries of A indexed by S. In what follows, we are interested in the case of the density N/mn separated from 0. Without loss of generality, we assume that n ≥ m. 4 (1.4) Theorem. Let us fix numbers κ > 0 and 0 < δ < 1. Then there exists a number q = q(κ, δ) such that the following holds. Let (R, C) be margins such that n ≥ m > q and the polytope P(R, C) has a non-empty interior, and let Z ∈ P(R, C) be the maximum entropy matrix. Let  S ⊂ (i, j) : i = 1, . . . , m, j = 1, . . . , n be a subset such that σS (Z) ≥ δmn and let ln n ǫ = δ√ . m If ǫ ≤ 1 then n Pr D ∈Σ(R, C) : o (1 − ǫ)σS (Z) ≤ σS (D) ≤ (1 + ǫ)σS (Z) ≥ 1 − 2n−κn . Let us associate with a non-negative, non-zero m × n matrix A = (aij ) a finite probability space on the ground set {(i, j) : i = 1, . . . , m, j = 1, . . . , n} with Pr {(i, j)} = aij /N , where N > 0 is the total sum of matrix entries. Theorem 1.4 asserts that the probability space associated with the maximum entropy matrix Z reasonably well approximates the probability space associated with a random binary contingency table D ∈ Σ(R, C) as far as events S whose probability is separated from 0 are concerned. The following interpretation of the maximum entropy matrix was suggested to the author by J.A. Hartigan, see [BH09]. (1.5) Theorem. Let Z = (zij ) be the m×n maximum entropy matrix with margins (R, C) and let us suppose that the polytope P(R, C) has a non-empty interior. Let X = (xij ) be the random m × n matrix of independent Bernoulli random variables such that E X = Z. In other words, Pr {xij = 1} = zij and Pr {xij = 0} = 1 − zij independently for all i, j. Then the probability mass function of X is constant on the set Σ(R, C) of binary contingency tables with margins (R, C), and, moreover,  Pr X = D = e−H(Z) for all D ∈ Σ(R, C). The distribution of the random matrix X in Theorem 1.5 can be characterized as the maximum entropy distribution in the class consisting of all probability distributions on the set {0, 1}m×n of matrices with 0-1 entries whose expectations lie in the affine subspace consisting of the matrices with row sums R and column sums C, see [BH09]. 5 2. Extensions and ramifications Our results hold in a somewhat greater generality. Let us fix an m × n nonnegative matrix W = (wij ), which we call the matrix of weights. Let us consider the following partition function X |Σ(R, C; W )| = Y wij . D∈Σ(R,C) i,j D=(dij ) dij =1 In particular, if wij = 1 for all i, j then |Σ(R, C; W )| = |Σ(R, C)|. If wij ∈ {0, 1} then the partition function counts binary contingency tables with zeros assigned to some positions: the value of |Σ(R, C; W )| is equal to the number of m × n matrices D = (dij ) such that the row sums of D are R, the column sums of D are C, dij ∈ {0, 1} for all i, j, and, additionally, dij = 0 if wij = 0. In combinatorial terms, the set Σ(R, C; W ) can be interpreted as the set of all subgraphs with prescribed degrees of vertices of a given bipartite graph. Binary contingency tables with preassigned zeros are of interest in statistics, see [C+05]. We prove the following result. (2.1) Theorem. Let us define the function F (x, y; W ) = m Y xi−ri  n Y j=1 i=1 for ! x = (x1 , . . . , xn ) and   Y −c yj j   (1 + wij xi yj ) ij y = (y1 , . . . , yn ) and let α(R, C; W ) = inf x1 ,... ,xm >0 y1 ,... ,yn >0 F (x, y; W ). Then for the partition function |Σ(R, C; W )| we have α(R, C; W ) ≥ |Σ(R,C; W )| ≥ (mn)! (mn)mn m Y (n − ri )n−ri i=1 (n − ri )!  cj n Y cj  α(R, C; W ).  c ! j=1 j ! As before, the function obtained as the result of the substitution xi = eti , yj = esj in ln F (x, y; W ), G(s, t; W ) = − for m X i=1 s = (s1 , . . . , sm ) ri s i − n X j=1 cj tj + X ln 1 + wij esi +tj ij and t = (t1 , . . . , tn ) 6
is convex on Rm × Rn , hence computing α(R, C; W ) is a convex optimization problem. Let us assume now that wij ∈ {0, 1} for all (i, j) and let us consider the set Σ(R, C; W ) of all m × n binary contingency tables D = (dij ) with the additional constraint that dij = 0 if wij = 0. Assuming that Σ(R, C; W ) is not empty, we consider this set as a finite probability space with the uniform measure. We call matrix W the pattern. We are interested in what a random table D ∈ Σ(R, C; W ) looks like. We define the maximum entropy matrix as before. (2.2) The maximum entropy matrix. Suppose that the set Σ(R, C; W ) is nonempty. Let us consider the polytope P(R, C; W ) of m × n matrices X = (xij ) such that n m X X xij = ri for i = 1, . . . , m, xij = cj for j = 1, . . . , n, j=1 i=1 0 ≤ xij ≤ 1 for all i, j and xij = 0 whenever wij = 0. Thus P(R, C; W ) is a face of polytope P(R, C) of Section 1.2. Let H(X) be the entropy function of Section 1.2. Since H(X) is strictly concave, it attains a unique maximum Z = Z(R, C; W ) on polytope P(R, C; W ), which we call the maximum entropy matrix with margins (R, C) and pattern W . (2.3) Lemma. Suppose that the polytope P(R, C; W ) contains a matrix Y = (yij ) such that 0 < yij < 1 whenever wij = 1, in which case we say that P(R, C; W ) has a non-empy interior. Then the infimum α(R, C; W ) in Theorem 2.1 is attained at a certain point x∗ = (ξ1 , . . . , ξm ) and y∗ = (η1 , . . . , ηn ). The maximum entropy matrix Z = (zij ) satisfies ξi ηj for all i, j such that wij = 1 (2.3.1) zij = 1 + ξi ηj Moreover, (2.3.2) α(R, C; W ) = eH(Z) . Conversely, if the infimum α(R, C; W ) is attained at a point x∗ = (ξ1 , . . . , ξm ) and y∗ = (η1 , . . . , ηn ), then for the maximum entropy matrix Z = (zij ) the equations (2.3.1) and (2.3.2) hold. For P(R, C; W ) to have a non-empty interior is equivalent to the requirement that for every pair k, l such that wkl = 1 there is a matrix D0 ∈ Σ(R, C; W
),
0 0 D = dij , such that d0kl = 0 and there is a matrix D1 ∈ Σ(R, C; W ), D1 = d1ij , such that d1kl = 1. In other words, we require the set Σ(R, C; W ) to be reasonably large. We prove an analogue of Theorem 1.4. We consider subsets n o S ⊂ (i, j) : wij = 1 . As before, we denote by σS (A) the sum of the entries of a matrix A indexed by the subset S. 7 (2.4) Theorem. Let us fix numbers κ > 0 and 0 < δ < 1. Then there exists a number q = q(κ, δ) such that the following holds. Let (R, C) be margins such that n ≥ m > q and the polytope P(R, C; W ) has a non-empty interior, and let Z ∈ P(R, C; W ) be the maximum entropy matrix. Let  S ⊂ (i, j) : wij = 1 be a subset such that σS (Z) ≥ δmn and let ln n ǫ = δ√ . m If ǫ ≤ 1 then n Pr D ∈Σ(R, C; W ) : o (1 − ǫ)σS (Z) ≤ σS (D) ≤ (1 + ǫ)σS (Z) ≥ 1 − 2n−κn . The statement of the theorem is, of course, vacuous unless pattern W contains Ω(mn) ones. There is an analogue of Theorem 1.5. (2.5) Theorem. Suppose that the polytope P(R, C; W ) has a non-empty interior and let Z ∈ P(R, C; W ) be the maximum entropy matrix. Let X = (xij ) be the random m × n matrix of independent Bernoulli random variables such that E X = Z, that is, Pr {xij = 1} = zij and Pr {xij = 0} = 1 − zij independently for all i, j. Then the probability mass function of X is constant on the set Σ(R, C; W ) and, moreover,  Pr X = D = e−H(Z) for all D ∈ Σ(R, C; W ). 3. Comparisons with the literature There is a vast literature on 0-1 matrices with prescribed row and column sums and with or without zeros in prescribed positions, see for example, Chapter 16 of [LW01], [Ne69], [Be74], [GC77], recent [CM05], [G+06], [C+08], [GM09] and references therein. A simple and efficient criterion for the existence of a 0-1 matrix with prescribed row and column sums is given by the classical Gale-Ryser Theorem; in the case of enforced zeros, the question reduces to the existence of a network flow, see for example, Chapter 16 of [LW01]. Estimating the number of such matrices also attracted a lot of attention. Precise asymptotic formulas for the number of matrices were obtained in sparse cases for which ri ≪ n and cj ≪ m [Ne69], [Be74], [G+06], the regular case of all row sums ri equal and all column sums cj equal [C+08] and cases close to regular [C+08], [GM09]. Formulas of Theorems 1.1 and 2.1 are not so precise but they are applicable to a wide class of margins (R, C) and they uncover some interesting features of the numbers |Σ(R, C)| and |Σ(R, C; W )|. The following construction provides some insight into the combinatorial interpretation of the number α(R, C) from Theorem 1.1. 8 (3.1) Cloning the margins. Let us fix some margins R, C for which the set Σ(R, C) is not empty, and, moreover, the polytope P(R, C) contains an interior point, so the conditions of Lemma 1.3 are satisfied. Let R = (r1 , . . . , rm ) and C = (c1 , . . . , cn ). For a positive integer k, let us define the km-vector   Rk = kr1 , . . . , kr1 , . . . , krm , . . . , krm  {z } {z } | | k times and the kn-vector k times   Ck = kc1 , . . . , kc1 , . . . , kcn , . . . , kcn  . | {z } {z } | k times k times In other words, we obtain margins (Rk , Ck ) if we choose a matrix Y ∈ P(R, C) and then create a new block matrix Yk by arranging k 2 copies of Y into a km × kn matrix. Then Rk is the vector of row sums of Yk and Ck is the vector of column sums of Yk . Clearly, the conditions of Lemma 1.3 are satisfied for (Rk , Ck ). Theorem 1.1 then implies that (3.1.1) lim 1/k2 k−→+∞ |Σ(Rk , Ck )| = α(R, C). Indeed, the infimum α(R, C) is attained at a certain point x∗ = (ξ1 , . . . , ξm ) and y∗ = (η1 , . . . , ηn ) . It is not hard to see that the infimum α (Rk , Ck ) is attained at   x∗k = ξ1 , . . . , ξ1 , . . . , ξm , . . . , ξm  | {z } | {z } k times k times   and yk∗ = η1 , . . . , η1 , . . . , ηn , . . . , ηn  . | {z } | {z } k times k times (3.2) Asymptotic repulsion in the space of matrices. A natural candidate for an approximation of |Σ(R, C)| is the “independence estimate” (3.2.1) I(R, C) =  mn N −1 Y n   m  Y m n , ri j=1 cj i=1 see [GC77], [G+06], and [C+08]. 9 The intuitive meaning of (3.2.1) is as follows. Let us consider the set of all m × n matrices with 0-1 entries and with the total sum of entries equal to N as a finite probability space with the uniform measure. Let us consider the two events in this space: the event R consisting of the matrices with row sums R and the event C consisting of the matrices with column sums C. One can see that  −1 Y  −1 Y n   m   m mn n mn and Pr (C) = Pr (R) = cj N ri N j=1 i=1 and that  mn Pr (R ∩ C). |Σ(R, C)| = N Thus the value of (3.2.1) equals |Σ(R, C)| if the events R and C are independent. It turns out that (3.2.1) indeed approximates |Σ(R, C)| reasonably well in the sparse and near-unform cases, see [G+06] and [C+08]. However, for generic R and C, the independence estimate I(R, C) overestimates |Σ(R, C)| by a 2Ω(mn) factor. To see why, let us fix some margins R = (r1 , . . . , rm ) and C = (c1 , . . . , cn ) such that not all row sums ri are equal and not all column sums cj are equal and the conditions of Lemma 1.3 are satisfied. Let us consider the cloned margins Rk and Ck as in Section 3.1. Applying Stirling’s formula, we get  lim (3.2.2) 1/k2 I (Rk , Ck ) k−→+∞    m r  X N i = exp −mnH +n H mn n i=1 n  c  X j , +m H m j=1 where H is the entropy function, see Section 1.2. To compare (3.2.2) and (3.1.1) we use Lemma 1.3 and the multivariate entropy function H(p1 , . . . , pk ) = k X i=1 pk ln 1 , pk where p1 , . . . , pk are non-negative numbers such that p1 + . . . + pk = 1. Thus H(x) = H(x, 1 − x) for 0 ≤ x ≤ 1 and we rewrite (3.2.2) as 1 lim ln I (Rk , Ck ) k−→+∞ k 2   r n − r1 rm  n − rm 1 + (mn − N )H ,... , ,... , =N H N N mn − N mn − N    c m − c1 cn m − cn 1 + (mn − N )H ,... , ,... , +N H N N mn − N mn − N − N ln N − (mn − N ) ln(mn − N ). 10 On the other hand, applying Lemma 1.3, we can rewrite (3.1.1) as 1 ln |Σ (Rk , Ck )| k−→+∞ k 2   z  1 − zij ij + (mn − N )H − N ln N =N H N mn − N − (mn − N ) ln(mn − N ), lim where Z = (zij ) is the maximum entropy matrix for margins (R, C). We now use some classical entropy inequalities, see, for example, [Kh57]. Namely, by the inequality relating the entropies of two partitions of a probability space and the entropy of their intersection, we have z  r c rm  cn  ij 1 1 H ≤H +H ,... , ,... , N N N N N with the equality if and only if ri cj (3.2.3) zij = for all i, j N and       n − r1 m − c1 n − rm m − cn 1 − zij ≤H +H ,... , ,... , H mn − N mn − N mn − N mn − N mn − N with the equality if and only if (n − ri )(m − cj ) for all i, j. mn − N However, if we have both (3.2.3) and (3.2.4), we must have (ri m−N )(cj n−N ) = 0, so unless all row sums ri are equal or all column sums cj are equal, we have (3.2.4) 1 − zij = lim k−→+∞ 1/k2 |Σ(Rk , Ck )| < lim k−→+∞ I (Rk , Ck ) 1/k2 . Therefore, as k grows, the independence estimate (3.2.1) overestimates the number 2 of 0-1 matrices with row sums Rk and column sums Ck by a factor of 2Ω(k ) . In probabilistic terms, as k grows, the event Rk consisting of the 0-1 matrices with row sums Rk and the event Ck consisting of the 0-1 matrices with column sums Ck repel each other (the events are negatively correlated) instead of being asymptotically independent. The procedure of cloning described in Section 3.1 produces margins of increasing size with the following features: the density remains separated from 0 and 1, and if the margins were non-uniform initially, they stay away from uniform. One can show that for more general sequences of margins that share these two features, we have the asymptotic repulsion of the event consisting of the 0-1 matrices with prescribed row sums and the event consisting of the 0-1 matrices with prescribed column sums. This is in contrast to the case of contingency tables (non-negative integer matrices with prescribed row and column sums), where we have the asymptotic attraction of the events [Ba09]. 11 (3.3) Randomized counting and sampling. Jerrum, Sinclair, and Vigoda [J+04] showed how to apply their algorithm for computing the permanent of a non-negative matrix to construct a fully polynomial randomized approximation scheme (FPRAS) to compute |Σ(R, C)| and, more generally, |Σ(R, C; W )|, where W is a 0-1 pattern, see also [B+07]. Furthermore, they obtained a polynomial time algorithm for sampling a random D ∈ Σ(R, C) and D ∈ Σ(R, C; W ) from a “nearly uniform” distribution. This problem arises naturally in statistics, see, for example, [C+05]. The estimates of Theorem 1.1 and Theorem 2.1 are not nearly as precise as those of [J+04], but they are deterministic, easily computable, and amenable to analysis. Similarly, we do not provide a sampling algorithm but show instead in Theorems 1.4 and 2.4 what a random matrix is likely to look like. (3.4) An open question. Theorem 1.5 allows us to interpret Theorem 1.4 as a law of large numbers for binary contingency tables: with respect to sums σS (D) for sufficiently “heavy” sets S of indices, a random binary contingency table D ∈ Σ(R, C) behaves approximately as the matrix of independent Bernoulli random variables whose expectation is the maximum entropy matrix Z = (zij ). Similar concentration results can be obtained for other well-behaved functions on binary contingency tables. One can ask whether the distribution of a particular entry dij of a random table D ∈ Σ(R, C) converges in distribution to the Bernoulli distribution with expectation zij as the dimensions m and n of the table grow in some regular way, for example, when the margins are cloned as in Section 3.1. Our approach, based on estimating combinatorial quantities via solutions to optimization problems, reminds one of that of Gurvits [Gu08]. The appearance of entropy in combinatorial counting problems reminds one of recent papers of Cuckler and Kahn [CK09a], [CK09b], although methods and results seem to be quite different. In the rest of the paper, we prove the results stated in Sections 1 and 2. 4. Preliminaries: permanents and scaling Let A = (aij ) be an n × n matrix. The permanent of A is defined by the expression n X Y aiσ(i) , per A = σ∈Sn i=1 where Sn is the symmetric group of all permutations of the set {1, . . . , n}. The relevance of permanents to us is that both values of |Σ(R, C)| and |Σ(R, C; W )| can be expressed as permanents of mn × mn matrices. This result is not new, for |Σ(R, C)| it was observed, for example, in [JS90]. For |Σ(R, C; W )|, where W is a 0-1 pattern, a construction is presented in [J+04]. We give a general construction for |Σ(R, C; W )|, where W is an arbitrary matrix, which is slightly different from that of [J+04]. 12 (4.1) Lemma. Let us choose margins R = (r1 , . . . , rm ), C = (c1 , . . . , cn ) and an m × n matrix W = (wij ) of weights. Let us construct an mn × mn matrix A = A(R, C; W ) as follows. The rows of A are split into disjoint m blocks having n − r1 , . . . , n − rm rows respectively (blocks of type I) and n blocks having c1 , . . . , cn rows respectively (blocks of type II). The columns of A are split into m disjoint blocks of n columns in each. For i = 1, . . . , m the entry of A that lies in a row from the i-th block of rows of type I and in a column from the i-th block of columns is equal to 1. For i = 1, . . . , m and j = 1, . . . , n the entry of A that lies in a row from the j-th block of rows of type II and the j-th column from the i-th block of columns is equal to wij . All other entries of A are 0s. Then  ! n m Y 1 Y 1   per A. |Σ(R, C; W )| = (n − r )! c ! i j j=1 i=1 Proof. First, we express |Σ(R, C; W )| as a coefficient in a certain polynomial. Let x1 , . . . , xn be formal variables and let X er (x1 , . . . , xn ) = xi1 · · · xir 1≤i1 <... 0 and a function φ(ǫ), 0 < ǫ < ǫ0 , such that lim φ(ǫ) = 1 ǫ−→0+ and for any n × n non-negative matrix B = (bij ) such that n X bij = 1 for j = 1, . . . , n i=1 and 1−ǫ ≤ n X j=1 bij ≤ 1 + ǫ 15 for i = 1, . . . , n for some 0 ≤ ǫ < ǫ0 , we have per B ≥ n! φ(ǫ). nn  From [L+00], one can choose ǫ0 = 1/n and φ(ǫ) = (1 − ǫn)n . 5. Proofs of Theorems 1.1 and 2.1 We prove Theorem 2.1 only since Theorem 1.1 is a particular case of Theorem 2.1. We start with a straightforward observation. (5.1) Lemma. We have Y (1 + wij xi yj ) = ij X (R,C) |Σ(R,C; W )|xR yC , where xR = xr11 · · · xrmm , yC = y1c1 · · · yncn , and the sum is taken over all margins R, C.  Next, we need a technical lemma. (5.2) Lemma. Let W = (wij ) be an m × n non-negative matrix such that α(R, C; W ) > 0. Then, for any ǫ > 0 there exist points x = x(ǫ) and y = y(ǫ), x = (x1 , . . . , xn ) and y = (y1 , . . . , yn ), such that −ri + −cj + n X j=1 m X i=1 wij xi yj <ǫ 1 + wij xi yj for i = 1, . . . , m wij xi yj <ǫ 1 + wij xi yj for j = 1, . . . , n xi , y j > 0 for all and i, j. Proof. Let us consider the function G(s, t; W ) = − for m X i=1 s = (s1 , . . . , sm ) ri s i − n X j=1 cj tj + X ln 1 + wij esi +tj ij and t = (t1 , . . . , tn ) . 16
Then G(s, t; W ) is convex and inf G(s, t) = ln α(R, C; W ) > −∞. s∈Rm t∈Rn Hence G(s, t) is bounded from below, it is also easy to check that the Hessian of G remains bounded on Rm × Rn . Therefore, the gradient of G(s, t) can get arbitrarily close to 0. That is, for any ǫ > 0 there are points s(ǫ) = (s1 (ǫ), . . . , sm (ǫ)) and t(ǫ) = (t1 (ǫ), . . . , tn (ǫ)) such that ∂ G(s, t) ∂si ∂ G(s, t) ∂tj s=s(ǫ),t=t(ǫ) < ǫ for i = 1, . . . , m and s=s(ǫ),t=t(ǫ) < ǫ for j = 1, . . . , n (it suffices to choose s(ǫ) and t(ǫ) so that the value of G(s(ǫ), t(ǫ)) is sufficiently close to the infimum). In other words, −ri + n X wij esi (ǫ)+tj (ǫ) < ǫ for 1 + wij esi (ǫ)+tj (ǫ) m X wij esi (ǫ)+tj (ǫ) < ǫ for 1 + wij esi (ǫ)+tj (ǫ) j=1 i = 1, . . . , m and −cj + i=1 j = 1, . . . , n. We now let xi = xi (ǫ) = esi (ǫ) for i = 1, . . . , m and yj = yj (ǫ) = etj (ǫ) for j = 1, . . . , n.  (5.3) Proof of Theorem 2.1. The upper bound α(R, C; W ) ≥ |Σ(R, C; W )| follows from Lemma 5.1. Let us prove the lower bound. If α(R, C; W ) = 0 then |Σ(R, C; W )| = 0 and the lower bound follows. Hence we assume that α(R, C; W ) > 0. Let A = A(R, C; W ) be the mn × mn block matrix constructed in Lemma 4.1. Let us consider the mn × mn block matrix B(ǫ) obtained from A as follows. For 17 ǫ > 0, let x(ǫ) = (x1 , . . . , xm ) and y(ǫ) = (y1 , . . . , yn ) be the point constructed in Lemma 5.2. For i = 1, . . . , m we multiply every row of A in the i-th block of type I by 1 ; xi (n − ri ) For j = 1, . . . , n we multiply every row of A in the j-th block of type II by yj cj for j = 1, . . . , n; For i = 1, . . . , m and j = 1, . . . , n we multiply the j-th column in the i-th block of columns of A by xi . 1 + wij xi yj This choice of scaling factors is, basically, a lucky guess made in the hope to match the structure of the function F (x, y; W ). Thus we have  ! n m Y −c c Y per A = yj j cj j  xin−ri (n − ri )n−ri  j=1 i=1  × and hence |Σ(R, C; W )| = (5.3.1) Y ij   x−1 i (1 + wij xi yj ) per B(ǫ) m Y (n − ri )n−ri !  c n Y cj j   cj c j j=1 i=1
× F x(ǫ), y(ǫ); W per B(ǫ) !  n cj  m n−ri Y Y cj (n − ri )   ≥ cj (n − r )! c i i=1 j=1 j (n − ri )! × α(R, C; W ) per B(ǫ). Finally, we claim that B(ǫ) is close to a doubly stochastic matrix. Indeed, For i = 1, . . . , m and j = 1, . . . , n the entry of B(ǫ) that lies in a row from the i-th block of rows of type I and in the j-th column from the i-th block of columns is equal to 1 . (n − ri )(1 + wij xi yj ) 18 For i = 1, . . . , m and j = 1, . . . , n the entry of B(ǫ) that lies in a row from the j-th block of rows of type II and the j-th column from the i-th block of columns is equal to wij xi yj . cj (1 + wij xi yj ) All other entries of B(ǫ) are 0s. Let us compute the row sums of B(ǫ). For a row in the i-th block of rows of type I the sum equals ai = n X j=1 Since n X j=1 1 . (n − ri )(1 + wij xi yj ) n n X 1 + wij xi yj X wij xi yj 1 = − , 1 + wij xi yj 1 + wij xi yj j=1 1 + wij xi yj j=1 by the inequalities of Lemma 5.2, we have |ai − 1| < ǫ ≤ ǫ for n − ri i = 1, . . . , m. For a row in the j-th block of rows of type II the sum equals bj = m X i=1 wij xi yj . cj (1 + wij xi yj ) By the inequalities of Lemma 5.2, we have |bj − 1| < ǫ ≤ ǫ for cj j = 1, . . . , n. Let us compute the column sums of B(ǫ). For the j-th column from the i-th block of columns the sum equals (n − ri ) 1 wij xi yj + cj = 1. (n − ri )(1 + wij xi yj ) cj (1 + wij xi yj ) Clearly, B(ǫ) is non-negative and hence by Lemma 4.3, we have per B(ǫ) ≥ (mn)! φ(ǫ) where (mn)mn The proof now follows by (5.3.1) as ǫ −→ 0+. 19 lim φ(ǫ) = 1. ǫ−→0+  6. Proofs of Lemmas 1.3 and 2.3 We prove Lemma 2.3 only since Lemma 1.3 is a particular case of Lemma 2.3. Proof of Lemma 2.3. Since H ′ (x) = ln(1 − x) − ln x, the value of the derivative at x = 0 is +∞ (we consider the right derivative there), the value of the derivative at x = 1 is −∞ (we consider the left derivative there) and the value of the derivative is finite for any 0 < x < 1. Suppose that for the maximum entropy matrix Z we have zij ∈ {0, 1} for some i, j such that wij = 1. If Y ∈ P(R, C; W ), Y = (yij ), is a matrix such that 0 < yij < 1 whenever wij = 1 then H (ǫY + (1 − ǫ)Z) > H(Z) for a sufficiently small ǫ > 0, which contradicts to the choice of Z. Hence 0 < zij < 1 whenever wij = 1. Therefore, the gradient of H(X) at X = Z is orthogonal to the affine subspace of matrices X = (xij ) having row sums R, column sums C, and such that xij = 0 whenever wij = 0. Hence (6.1) ln 1 − zij = λi + µj zij for all i, j such that wij = 1 and some λ1 , . . . , λm and µ1 , . . . , µn . Hence e−λi e−µj zij = 1 + e−λi e−µj whenever wij = 1. Therefore, X e−λi e−µj = ri 1 + e−λi e−µj for i = 1, . . . , m X e−λi e−µj = cj 1 + e−λi e−µj for j = 1, . . . , n. j: wij =1 i: wij =1 Therefore, s∗ = (−λ1 , . . . , −λm ) and t∗ = (−µ1 , . . . , −µn ) is a critical point of G(s, t; W ) = − m X i=1 ri s i − n X cj tj + j=1 20 X (i,j): wij =1
ln 1 + esi +tj . Since G is convex, (s∗ , t∗ ) is also a minimum point. Therefore, the point x∗ = (ξ1 , . . . , ξm ) and y∗ = (η1 , . . . , ηn ) where ξi = e−λi i = 1, . . . , m and ηj = e−µj for for j = 1, . . . , n is a minimum point of m Y F (x, y; W ) = xi−ri !  n Y −cj  yj j=1 i=1  Y (1 + xi yj ) (i,j): wij =1 and satisfies X ξi ηj = ri 1 + ξi ηj for i = 1, . . . , m X ξi ηj = cj 1 + ξi ηj for j = 1, . . . , n. j: wij =1 (6.2) i: wij =1 Conversely, if x∗ = (ξ1 , . . . , ξm ) and y∗ = (η1 , . . . , ηn ) is a point where the minimum of F (x, y; W ) is attained, then, setting the gradient of ln F to 0, we obtain equations (6.2). Letting zij = ξi ηj 1 + ξi ηj when wij = 1 and zij = 0 when wij = 0, we obtain a matrix Z ∈ P(R, C; W ). Moreover, the gradient of H(X) at X = Z satisfies (6.1) with λi = − ln ξi and µj = − ln ηj , so Z is the maximum entropy matrix. We now check: X X H(Z) = − zij ln zij − (1 − zij ) ln (1 − zij ) (i,j): wij =1 =− =− X (i,j): wij =1 m X i=1 (i,j): wij =1 X ξi ηj 1 1 ξi ηj ln − ln 1 + ξi ηj 1 + ξi ηj 1 + ξi ηj 1 + ξi ηj (i,j): wij =1  X ln ξi    j: wij =1 + X n X  ξi ηj  ξi ηj  X − ln η   j  1 + ξi ηj 1 + ξi ηj i: j=1 wij =1 ln (1 + ξi ηj ) (i,j): wij =1 =− m X i=1 ri ln ξi −  n X cj ln ηj + j=1 X (i,j): wij =1 21 ln (1 + ξi ηj ) by (6.2). Hence H(Z) = ln F (x∗ , y∗ ; W ) and the proof follows.  7. Proofs of Theorems 1.5 and 2.5 We prove Theorem 2.5 only, since Theorem 1.5 is a particular case of Theorem 2.5. From formula (6.1), we have 1 − zij = eλi +µj zij for all i, j such that wij = 1 and some λ1 , . . . , λm and µ1 , . . . , µn . Then, for all i, j such that wij = 1 and any dij ∈ {0, 1}, we have  d 1−d Pr xij = dij =zijij (1 − zij ) ij = (1 − zij ) = (1 − zij ) e−(λi +µj )dij .  1 − zij zij −dij Consequently, for any D ∈ Σ(R, C; W ), D = (dij ), we have  Pr X = D = Y i,j: wij =1  = Y (1 − zij ) e−(λi +µj )dij i,j: wij =1  m Y (1 − zij ) e−λi ri i=1 !  n Y j=1  e−µj cj  . On the other hand, e−H(Z) = Y i,j: wij =1  z zijij (1 − zij ) 1−zij   −zij  1 − zij  zij i,j: wij =1 i,j: wij =1    ! n m Y Y Y = (1 − zij ) e−µj cj  , e−λi ri  = Y (1 − zij )  i,j: wij =1 Y i=1 which completes the proof. j=1  22 8. Proofs of Theorems 1.4 and 2.4 We prove Theorem 2.4 only since Theorem 1.4 is a particular case of Theorem 2.4. We will use standard large deviation inequalities for bounded random variables, see, for example, Corollary 5.2 of [Mc89]. (8.1) Lemma. Let Y1 , . . . , Yk be independent random variables such that 0 ≤ Yi ≤ 1 for i = 1, . . . , k. Let Y = Y1 + . . . + Yk and let a = E Y . Then, for 0 ≤ ǫ ≤ 1 we have  Pr Y ≥ (1 + ǫ)a  Pr Y ≤ (1 − ǫ)a   1 2 and ≤ exp − ǫ a 3   1 2 ≤ exp − ǫ a . 2  (8.2) Proof of Theorem 2.4. Let X = (xij ) be the m × n matrix of independent Bernoulli random variables such that E X = Z, as in Theorem 2.5. By Theorem 2.5, the distribution of X conditioned on Σ(R, C; W ) is uniform and hence  Pr D ∈Σ(R, C; W ) : σS (D) ≤ (1 − ǫ)σS (Z)  Pr X : σS (X) ≤ (1 − ǫ)σS (Z) and X ∈ Σ(R, C; W )  . = Pr X : X ∈ Σ(R, C; W ) Similarly,  Pr D ∈Σ(R, C; W ) : σS (D) ≥ (1 + ǫ)σS (Z)  Pr X : σS (X) ≥ (1 + ǫ)σS (Z) and X ∈ Σ(R, C; W )  = . Pr X : X ∈ Σ(R, C; W ) By Theorem 2.5, Lemma 2.3 and Theorem 2.1, we get  Pr X ∈ Σ(R, C; W ) = e−H(Z) |Σ(R, C; W )| ≥ (mn)! (mn)mn i=1 ≥ (mn)−γ(m+n) for some absolute constant γ > 0. 23 !  n cj  Y cj   (n − ri )! c ! j j=1 m Y (n − ri )n−ri Therefore,  Pr D ∈Σ(R, C; W ) : σ(D) ≤ (1 − ǫ)σS (Z)  ≤ (mn)γ(m+n) Pr X : σS (X) ≤ (1 − ǫ)σS (Z) and similarly  Pr D ∈Σ(R, C; W ) : σ(D) ≥ (1 + ǫ)σS (Z)  ≤ (mn)γ(m+n) Pr X : σS (X) ≥ (1 + ǫ)σS (Z) . (8.2.1) By Lemma 8.1,  Pr X : σS (X) ≤ (1 − ǫ)σS (Z) (8.2.2) 1 ≤ exp − ǫ2 σS (Z) 2 and  Pr X : σS (X) ≥ (1 + ǫ)σS (Z) Hence for  δ ln n ǫ= √ m   1 2 ≤ exp − ǫ σS (Z) . 3  and σS (Z) ≥ δmn we have ǫ2 σS (Z) ≥ δ 3 n ln2 n. (8.2.3) Combining (8.2.1)–(8.2.3), we conclude that for any κ > 0 and all sufficiently large n ≥ m > q(κ, δ) we have  Pr D ∈ Σ(R, C; W ) : σS (D) ≤ (1 − ǫ)σS (Z)  Pr D ∈ Σ(R, C; W ) : σS (D) ≥ (1 + ǫ)σS (Z) as required. ≤ n−κn ≤ n and −κn  Acknowledgment I am grateful to Alex Samorodnitsky for emphasizing the importance of matrix scaling in permanent computations during many enlightening conversations. I benefitted from conversations with Catherine Greenhill about asymptotic enumeration during the 2008 Schloss Dagstuhl meeting on design and analysis of randomized and approximation algorithms. After the first version of this paper was written, John Hartigan pointed out to connections with the maximum entropy principle and suggested Theorem 1.5 to me (cf. [BH09]), which led to a substantial simplification of the original proof and some strengthening of Theorems 1.4 and 2.4. 24 References [Ba96] [Ba07] [Ba09] [BH09] [Be74] [B+07] [C+08] [CM05] [C+05] [CK09a] [CK09b] [GC77] [G+06] [GM09] [Gu04] [Gu08] [JS90] [J+04] [Kh57] [L+00] A. Barvinok, Two algorithmic results for the traveling salesman problem, Mathematics of Operations Research 21 (1996), 65–84. A. Barvinok, Integration and optimization of multivariate polynomials by restriction onto a random subspace, Foundations of Computational Mathematics 7 (2007), 229– 244. A. Barvinok, Asymptotic estimates for the number of contingency tables, integer flows, and volumes of transportation polytopes, International Mathematics Research Notices. IMRN 2009 (2009), 348–385. A. Barvinok and J. Hartigan, Maximum entropy Gaussian approximation for the number of integer points and volumes of polytopes, preprint arXiv:0903.5223 (2009). E.A. Bender, The asymptotic number of non-negative integer matrices with given row and column sums, Discrete Mathematics 10 (1974), 217–223. I. Bezáková, N. Bhatnagar, and E. Vigoda, Sampling binary contingency tables with a greedy start, Random Structures & Algorithms 30 (2007), 168–205. E.R. Canfield, C. Greenhill, and B.D. McKay, Asymptotic enumeration of dense 01 matrices with specified line sums, Journal of Combinatorial Theory. Series A 115 (2008), 32–66. E. R. Canfield and B.D. McKay, Asymptotic enumeration of dense 0-1 matrices with equal row sums and equal column sums, Research Paper 29, 31 pp., Electronic Journal of Combinatorics 12 (2005). Y. Chen, P. Diaconis, S.P. Holmes, and J.S. Liu, Sequential Monte Carlo methods for statistical analysis of tables, Journal of the American Statistical Association 100 (2005), 109–120. W. Cuckler and J. Kahn, Hamiltonian cycles in Dirac graphs, Combinatorica 29 (2009), 299–326. W. Cuckler and J. Kahn, Entropy bounds for perfect matchings and Hamiltonian cycles, Combinatorica 29 (2009), 327–335. I.J. Good and J.F. Crook, The enumeration of arrays and a generalization related to contingency tables, Discrete Mathematics 19 (1977), 23–45. C. Greenhill, B.D. McKay, and X. Wang, Asymptotic enumeration of sparse 0-1 matrices with irregular row and column sums, Journal of Combinatorial Theory. Series A 113 (2006), 291–324. C. Greenhill and B.D. McKay, Random dense bipartite graphs and directed graphs with specified degrees, Random Structures & Algorithms 35 (2009), 222–249. L. Gurvits, Classical complexity and quantum entanglement, Journal of Computer and System Sciences 69 (2004), 448–484. L. Gurvits, Van der Waerden/Schrijver-Valiant like conjectures and stable (aka hyperbolic) homogeneous polynomials: one theorem for all, Research Paper 66, Electronic Journal of Combinatorics 15(1) (2008). M. Jerrum and A. Sinclair, Fast uniform generation of regular graphs, Theoretical Computer Science 73 (1990), 91–100. M. Jerrum, A. Sinclair, and E. Vigoda, A polynomial-time approximation algorithm for the permanent of a matrix with nonnegative entries, Journal of the ACM 51 (2004), 671–697. A.I. Khinchin, Mathematical Foundations of Information Theory, Dover Publications, Inc., New York, N. Y., 1957. N. Linial, A. Samorodnitsky, and A. Wigderson, A deterministic strongly polynomial algorithm for matrix scaling and approximate permanents, Combinatorica 20 (2000), 545–568. 25 [LW01] [Ma95] [Mc89] [Ne69] [NN94] J.H. van Lint and R.M. Wilson, A Course in Combinatorics. Second edition, Cambridge University Press, Cambridge, 2001. I.G. Macdonald, Symmetric Functions and Hall Polynomials. Second edition. With contributions by A. Zelevinsky, Oxford Mathematical Monographs, Oxford Science Publications. The Clarendon Press, Oxford University Press, New York, 1995. C. McDiarmid, On the method of bounded differences, Surveys in combinatorics, 1989 (Norwich, 1989), London Mathematical Society Lecture Note Series, vol. 141, Cambridge Univ. Press, Cambridge, 1989, pp. 148–188. P.E. O’Neil, Asymptotics and random matrices with row-sum and column-sum restrictions, Bulletin of the American Mathematical Society 75 (1969), 1276–1282. Yu. Nesterov and A. Nemirovskii, Interior-Point Polynomial Algorithms in Convex Programming, SIAM Studies in Applied Mathematics, 13, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1994. Department of Mathematics, University of Michigan, Ann Arbor, MI 48109-1043, USA E-mail address: barvinok@umich.edu 26