Sampling Contingency Tables

Martin Dyer^{∗}  Ravi Kannan^{†}  John Mount^{‡}

February 3, 1995

1 Introduction

Given positive integers r_{1}, r_{2}, … r_{m} and c_{1}, c_{2}, … c_{n}, let I(r, c) be the set of m×n arrays with nonnegative integer entries and row sums r_{1}, r_{2}, … r_{m} respectively and column sums c_{1}, c_{2}, … c_{n} respectively. Elements of I(r, c) are called contingency tables with these row and column sums.

We consider two related problems on contingency tables. Given r_{1}, r_{2}, … r_{m} and c_{1}, c_{2}, … c_{n},

1) Determine |I(r, c)|.

2) Generate randomly an element of I(r, c), each with probability 1/|I(r, c)|.

The counting problem is of combinatorial interest in many contexts. See for example, the survey by Diaconis and Gangolli [4]. We show that even in the case when m or n is 2, this problem is #P hard.

∗School of Computer Studies, Leeds University, Leeds, England

†School of Computer Science, Carnegie-Mellon University, Pittsburgh, PA 15213, USA. Partially supported by NSF Grant CCR 9208597.

‡School of Computer Science, Carnegie-Mellon University, Pittsburgh, PA 15213, USA. Partially supported by NSF Grant CCR 9208597.

The random sampling problem is of interest in Statistics. In particular, it comes up centrally in an important new method formulated by Diaconis and Efron [3] for testing independence in two-way tables. (See also Mehta and Patel [17] for a more classical test of independence). We show that this problem can be solved in approximately in polynomial time provided the row and column sums are sufficiently large, in particular, provided r_{i} ∈ Ω(n^{2}m) and c_{j} ∈ Ω(m^{2}n). This then implies under the same conditions, a polynomial time algorithm to solve the counting problem (1) approximately.

The algorithm we give relies on the random walk approach, and is closely related to those for volume computation and sampling from log-concave distributions [?, ?, ?, ?, ?, ?, ?].

2 Preliminaries and Notation

The number of rows will always be denoted by m and the number of columns by n. We will let N denote (n − 1)(m − 1). For convenience, we number the coordinates of any vector in R^{nm} by a pair of integers (i, j) or just ij where i runs from 1 through m and j from 1 through n. Suppose u is an m vector of reals (the row sums) and v is an n vector of reals (the column sums). We define

V (u, v) = {x ∈ R^{nm} : ∑_{j} x_{ij} = u_{i} for i = 1, 2, … m; ∑_{i} x_{ij} = v_{j} for j = 1, 2, … n}.

V (u, v) can be thought of as the set of m × n matrices with row and column sums specified by u and v respectively. Let

P (r, c) = V (r, c) ∩ {x : x_{ij} ≥ 0 for i = 1, 2, … m, j = 1, 2, … n}.

We call such a polytope a “contingency polytope”.

Then, I(u, v) denotes the set of vectors in P (u, v) with all integer coordinates.

Note that the above sets are all trivially empty if ∑_{i} u_{i} ≠ ∑_{j} v_{j}. So we will assume throughout that ∑_{i} u_{i} = ∑_{j} v_{j}.

We will need another quantity denoted by α(u, v) defined for u, v satisfying u_{i} > 2n∀i and v_{j} > 2m∀j :

α(u, v) = max_{i,j} (\frac{u_{i} + 2n}{u_{i} − 2n}, \frac{v_{j} + 2m}{v_{j} − 2m})^{(n−1)(m−1)} .

We will abbreviate α(u, v) to α when u, v are clear from the context. An easy calculation shows that if u_{i} is Ω(n^{2}m) ∀i and v_{j} is Ω(nm^{2}) ∀j, then α is O(1).

On input u, v (we assume that u_{i} > 2n∀i and v_{j} > 2m∀j), our algorithm runs for time bounded by α(u, v) times a polynomial in n , m , max_{ij}(log u_{i}, log v_{j}), and (1/ϵ) where ϵ will be an error parameter specifying the desired degree of accuracy. (See section 5, first paragraph for a description of the input/output specifications of the algorithm.)

3 Hardness of counting

We will show that exactly counting contingency tables is hard, even in a very restricted case. We will use the familiar notion of # P-completeness introduced by Valiant [20].

Theorem 1 The problem of determining the exact number of contingency tables with prescribed row and column sums is is #P-complete, even in the 2 × n case.

Proof It is easy to see that this problem is in #P. We simply guess mn integers x_{ij} in the range [0, J] (where J is the sum of the row sums) and check whether they satisfy the constraints. The number of accepting computations is the number of tables.

For proving hardness, we proceed as follows : Given positive integers a_{1}, a_{2}, …, a_{n−1}, b, it is shown in [?] that it is #P-hard to compute the (n − 1)-dimensional volume of the polytope

∑_{j=1}^{n−1} a_{j}y_{j} ≤ b, 0 ≤ y_{j} ≤ 1 (j = 1, 2, …, n − 1).

It follows that it is #P-hard to compute the (n − 1)-dimensional volume of the polytope

∑_{j=1}^{n} a_{j}y_{j} = b, 0 ≤ y_{j} ≤ 1 (j = 1, 2, …, n),

where a_{n} = b. Hence, by substituting x_{1j} = a_{j}y_{j}, x_{2j} = a_{j}(1 − y_{j}), that it follows that it is #P-hard to compute the (n − 1)-dimensional volume of the polytope P (r, c) (with 2 rows and n columns) where

r = (b, ∑_{j=1}^{n−1} a_{j}), c = (a_{1}, …, a_{n−1}, b).

Now the number ν(r, c) of integer points in P (r, c) is clearly the number of contingency tables with the required row and column sums. But, for integral r, c, P (r, c) is a polytope with integer vertices. (It is the polytope of a 2 × n transportation problem [8].) Consider the family of polytopes P (tr, tc) for t = 1, 2, . . .. It is well known [9, Chapter 12] that the number of integer points in P (tr, tc) will be a polynomial in t, the Ehrhart polynomial, of degree (n − 1), the dimension of P (r, c). Moreover the coefficient of t^{n−1} will be the (n − 1)-dimensional volume of P (r, c). It is also straightforward to show that the coefficients of this polynomial are of size polynomial in the length of the description of P (r, c).

Suppose therefore we could count the number of 2 × n contingency tables for arbitrary r, c. Then we could compute ν(r, c) for arbitrary integral r, c. Hence we could compute ν(tr, tc) for t = 1, 2, …, n and thus determine the coefficients of the Ehrhart polynomial. But this would allow us to compute the volume of P (r, c), which we have seen is a #P-hard quantity.

□

Remark If both n, m are fixed, we can apply a recent result of Barvinok’s [?] (that the number of lattice points in a fixed dimensional polytope can be counted exactly in polynomial time) to compute |I(u, v)|. Diaconis and Efron [3] give an explicit formula for |I(u, v)| in case n = 2. Their formula has exponentially many terms (as expected from our hardness result). Mann has a similar formula for n = 3. Sturmfels demonstrates a structural theorem which allows for an effective counting scheme that runs in polynomial time for n, m fixed.[19] Sturmfels has used the technique to count 4 × 4 contingency tables in the literature. We will, in later subsequent drafts of this paper, discuss our work counting 4 × 4 contingency tables.

4 Sampling Contingency tables : Reduction to continuous sampling

This section reduces the problem of sampling from the discrete set of contingency tables to the problem of sampling with near-uniform density from a contingency polytope.

To this end, we first take a natural basis for the lattice of all integer points in V (u, v) and associate with each lattice point a parallelepiped with respect to this basis. Then we produce a convex set P^{⊙} which has two properties : a) each parallelepiped associated with a point in I(u, v) is fully contained in P^{⊙} and b) the volume of P^{⊙} divided by the total volume of the parallelepipeds associated with points in I(u, v) is at most α(r, c). Now the algorithm is simple : pick a random point y from P^{⊙} with near uniform density (this can be done in polynomial time), find the lattice point x in whose parallelepiped y lies and accept x if it is in I(u, v) [acceptance occurs with probability at least 1/α(u, v) ]; otherwise, reject and rerun. The volume of each parallelepiped is the same, so the probability distribution of the x in I(u, v) is near uniform. This general approach has been used also for the simpler context of sampling from the 0-1 solutions to a knapsack problem [10].

[While P^{⊙} is a convex set and by the general methods, we can sample from it in polynomial time with near uniform density, it will turn out that P^{⊙} is also isomorphic to a contingency polytope. In the next section, we show how to exploit this to get a better polynomial time algorithm than the general one.]

We will have to build up some geometric facts before producing the P^{⊙}.

Let U be the lattice

{x ∈ R^{nm} : ∑_{j} x_{ij} = 0 for i = 1, 2, … m; ∑_{i} x_{ij} = 0 for j = 1, 2, … n; x_{ij} ∈ Z}.

For 1 ≤ i ≤ m − 1 and 1 ≤ j ≤ n − 1, let b(ij) be the vector in R^{nm} given by b(ij)_{ij} = 1, b(ij)_{i+1,j} = −1, b(ij)_{i,j+1} = −1, b(ij)_{i+1,j+1} = 1 and b(ij)_{kl} = 0 for kl other than the 4 above.

Any vector x in V (0, 0) can be expressed as a linear combination of the b(ij)’s as follows (the reader may check this by direct calculation)

x = ∑_{k=1}^{m−1} ∑_{l=1}^{n−1} (∑_{i=1}^{k} ∑_{j=1}^{l} x_{ij}) b(kl).  (1)

It is also easy to check that the b(ij) are all linearly independent. This implies that the subspace V (0, 0) has dimension (n − 1)(m − 1) and so does the affine set V (u, v) which is just a translate of V (0, 0). Also, we see that if x is an integer vector, then the above linear combination has integer coeffients; so, the b(ij) form a basis of the lattice U.

It is easy to see that if u, v are positive vectors, then the dimension of P (u, v) is also N. To see this, it suffices to come up with an x ∈ R^{nm} with row and column sums given by u, v and with each entry x_{ij} strictly positive because then we can add any small real multiples of the b(ij) to x and still remain in P (u, v). Such an x is easy to obtain : for example, we can choose x_{11}, x_{12}, … x_{1n} to satisfy 0 < x_{1j} < min(u_{1}, v_{j}) and summing to u_{1}. Then we subtract this amount from the column sums and repeat the process on the second row etc.

We will denote by Vol(P (u, v)) the (N dimensional) volume of P (u, v).

Obviously, if either u or v has a non-integer coordinate, then I(u, v) is empty.

Lemma 1 : If p, q are m vectors of positive reals and s, t are n vectors of positive reals with q ≥ p and t ≥ s (componentwise), then

Vol(P (q, t)) ≥ Vol(P (p, s)).

Proof By induction on the number of coordinates of (q t) that are strictly greater than the corresponding coordinates of (p s). After changing row and column numbers if necessary, assume without loss of generality that q_{1} − p_{1} is the least among all POSITIVE components of (q t) − (p s). After permuting columns if necessary, we may also assume that we have t_{1} − s_{1} > 0. Let P_{1} = P (p, s). Let P_{2} be defined by

P_{2} = {x ∈ R^{nm} : x_{11} ≥ −(q_{1} − p_{1}); x_{ij} ≥ 0 for (ij) ≠ (11)} ∩ V (p, s).

Let p′ be an m−vector defined by p′_{1} = q_{1}; p′_{i} = p_{i} for i = 2, 3, … m. Let s′ be defined by s′_{1} = s_{1} + q_{1} − p_{1}; s′_{j} = s_{j} for j = 2, 3, … n. Denote by P_{3} the set P (p′, s′). Then P_{2} and P_{3} have the same volume as seen by the isomorphism x_{11} → x_{11} + (q_{1} − p_{1}). Clearly, P_{2} contains P_{1}. Also the row sums and column sums in P_{3} are no greater than the corresponding row and column sums in (q t). Applying the inductive assumption to P_{3} and P (q, t), we have the lemma. □

Let u, v be fixed and consider V (u, v). Let

U(u, v) = V (u, v) ∩ Z^{nm}.

We may partition V (u, v) into fundamental parallelepipeds, one corresponding to each point of U(u, v). Namely,

V (u, v) = ∪_{x∈U(u,v)}{x+∑_{i,j} λ_{ij}b(ij) : 0 ≤ λ_{ij} < 1 for i = 1, 2, … m; j = 1, 2, … n}.

We call {x + ∑_{i,j} λ_{ij}b(ij) : 0 ≤ λ_{ij} < 1 for i = 1, 2, … m; j = 1, 2, … n}, “the parallelepiped” associated with x and denote it by F (x).

Claim : For any x ∈ I(u, v), and for any y belonging to F (x) we have x_{ij} − 2 < y_{ij} < x_{ij} + 2.

Proof : This follows from the fact that at most four of the b(ij) ’s have a nonzero entry in each position and two of them have a +1, the other two a -1. □

Lemma 2 : For nonnegative vectors u ∈ R^{m} and v ∈ R^{n} and any l ∈ R^{nm}, let

P (u, v, l) = V (u, v) ∩ {x ∈ R^{nm} : x_{ij} ≥ l_{ij}}.

Then,

(i) for any x in I(u, v), we have F (x) ⊆ P (u, v, −\underline{2}) where \underline{2} is the nm-vector of all 2 ’s.

(ii) P (u, v, +\underline{2}) ⊆ ∪_{x∈I(u,v)}F (x). [P (u, v, −\underline{2}) will play the role of P^{⊙} in the brief description of the algorithm given earlier.]

Proof : (i) Since x is in I(u, v), we have x_{ij} ≥ 0 and so for any y ∈ F (x), we have y_{ij} ≥ −2 proving (i). For (ii), observe that for any y ∈ V (u, v), there is a unique x ∈ U(u, v) such that y ∈ F (x). If y_{ij} ≥ 2, then we must have x_{ij} ≥ 0, so the corresponding x belongs to I(u, v) as claimed.

Lemma 3 : For u ∈ R^{m} and v ∈ R^{n} with u_{i} > 2n for all i and v_{j} > 2m for all j, we have

\frac{Vol(∪_{x∈I(u,v)}F (x))}{Vol(P (u, v, −\underline{2}))} ≥ \frac{1}{α(u, v)}.

Proof : By (ii) of the lemma above, Vol(∪_{x∈I(u,v)}F (x)) is at least Vol(P (u, v, +\underline{2})). But P (u, v, +\underline{2}) is isomorphic to P (u−2n\underline{1}, v−2m\underline{1}) where \underline{1} ’s are vectors of all 1 ’s of suitable dimensions as seen from the substitution x′_{ij} = x_{ij} − 2. Similarly, we have that P (u, v, −\underline{2}) is isomorphic to P (u + 2n\underline{1}, v + 2m\underline{1}). Let ρ = α^{1/N}. Consider the set ρP (u − 2n\underline{1}, v − 2m\underline{1}) = \{ρx : x ∈ P (u − 2n\underline{1}, v − 2m\underline{1})\}. This is precisely the set P (ρ(u − 2n\underline{1}), ρ(v − 2m\underline{1})). By the definition of ρ, we have ρ(u_{i} − 2n) ≥ u_{i} + 2n∀i and ρ(v_{j} − 2m) ≥ v_{j} + 2m∀j. So by lemma [?], the volume of ρP (u − 2n\underline{1}, v − 2m\underline{1}) is at least the volume of P (u + 2n\underline{1}, v + 2m\underline{1}). But ρP (u − 2n\underline{1}, v − 2m\underline{1}) is a dilation of the N dimensional object P (u − 2n\underline{1}, v − 2m\underline{1}) by a factor of ρ and thus has volume precisely equal to ρ^{N} = α times the volume of P (u − 2n\underline{1}, v − 2m\underline{1}) completing the proof.

□

Since P (u, v, −2) is isomorphic to P (u + 2n\underline{1}, v + 2m\underline{1}), the essential problem we have is one of sampling from a contingency polytope with uniform density.

Proposition 1 For every x ∈ U(r, c), we have vol_{N}(F (x)) = 1.

Proof : Note that if z = ∑_{i=1}^{m−1} ∑_{j=1}^{n−1} λ_{ij}b(ij), then we have z_{ij} = λ_{ij} − λ_{i−1,j} − λ_{i,j−1} + λ_{i−1,j−1} (where if one of the subscripts is 0, we take the λ to be zero too. Thus for each i, 1 ≤ i ≤ m − 1 and each j, 1 ≤ j ≤ n − 1, the “height” of the parallelepiped F (x) in the direction of z_{ij} perpendicular to the “previous” ((a, b) < (i, j) if a < i or b < j) coordinates is 1 proving our proposition.

□

5 The Algorithm

We now describe the algorithm which given row sums u_{1}, u_{2}, … u_{m} and column sums v_{1}, v_{2}, … v_{n}, and an error parameter ϵ > 0, draws a set W of samples from I(u, v) (the set of integer contingency tables with the given row and column sums) with the following property :

if g : I(u, v) → R is any function, then the true mean ḡ of g and the sample mean g̃ (defined below) satisfy

E((|g̃ − ḡ|)^{2}) ≤ ϵg_{0}^{2}

where

ḡ = \frac{∑_{w∈I(u,v)} g(w)}{|I(u, v)|} g̃ = \frac{∑_{w∈W} g(w)}{|W |} g_{0} = max_{w,x∈I(u,v)} g(w) − g(x).

As stated in the last section, we will first draw samples from P (u+2n\underline{1}, v+2m\underline{1}). We let r = u + 2n\underline{1} and c = v + 2m\underline{1}.

We assume that u_{i} > 2n∀i and v_{j} > 2m∀j throughout this section.

The function g often takes on values 0 or 1. This is so in the motivating application described by Diaconis and Efron [3] where one wants to estimate the proportion of contingency tables that have (the so-called) χ^{2} measure of deviation from the independent greater than the observed table. So, in this case, g will be 1 for a table if its χ^{2} measure is greater than that of the observed, otherwise, it will be zero. The evaluation of g for a particular table is in general simple as it is in this case.

At the end of the section, we describe how to use the sampling algorithm to estimate |I(u, v)|.

For reasons that will become clear later, we rearrange the rows and columns so that

r_{1}, r_{2}, … r_{m−1} ≤ r_{m} c_{1}, c_{2}, … c_{n−1} ≤ c_{n}.  (2)

We assume the above holds from now on. An m×n table x in P (r, c) is obviously specified completely by its entries in all but the last row and column. From this, it is easy to see that P (r, c) is isomorphic to the polytope Q(r, c) in R^{N} (recall the notation N = (m − 1)(n − 1)) which is defined as the set of x satisfying the following inequalities :

∑_{j=1}^{n−1} x_{ij} ≤ r_{i}∀i ∑_{i=1}^{m−1} x_{ij} ≤ c_{j}∀j ∑_{i=1}^{m−1} ∑_{j=1}^{n−1} x_{ij} ≥ ∑_{i=1}^{m−1} r_{i} − c_{n} x_{ij} ≥ 0.

In the above, as well as in the rest of the section, i will run through 1, 2, … m − 1 and j will run through 1, 2, … n − 1 unless otherwise specified.

The algorithm will actually pick points in Q(r, c) with near uniform density. To do so, the algorithm first scales each coordinate so that Q(r, c) becomes “wellrounded”. Let

ρ_{ij} = Min (\frac{r_{i}}{n − 1}, \frac{c_{j}}{m − 1}) for i = 1, 2, … m − 1 and j = 1, 2, … n − 1. (3)

The scaling is given by

y_{ij} = \frac{N}{ρ_{ij}} x_{ij} for i = 1, 2, … m − 1 and j = 1, 2, … n − 1.  (4)

This transformation transforms the polytope Q(r, c) to a polytope QQ = QQ(r, c) “in y space”. Note that QQ(r, c) is the set of y that satisfy the following constraints :

\frac{1}{N} ∑_{j} ρ_{ij}y_{ij} ≤ r_{i}∀i  \frac{1}{N} ∑_{i} ρ_{ij}y_{ij} ≤ c_{j}∀j  (5)

\frac{1}{N} ∑_{i,j} ρ_{ij}y_{ij} ≥ ∑_{i} r_{i} − c_{n} y_{ij} ≥ 0.

Instead of working in QQ(r, c), the algorithm will work with a larger convex set. To describe the set, first consider the convex set P′(r, c) obtained from P (r, c) by discarding the constraint x_{mn} ≥ 0. The corresponding set Q′(r, c) in N space is given by upper bounds on the row and column sums and nonnegativity constraints. The scaling above transforms Q′(r, c) to the following set QQ′(r, c) :

QQ′(r, c) = \{y : \frac{1}{N} ∑_{j} ρ_{ij}y_{ij} ≤ r_{i}∀i \frac{1}{N} ∑_{i} ρ_{ij}y_{ij} ≤ c_{j}∀j y_{ij} ≥ 0\} (6)

We will impose the following log-concave function F : P′(r, c) → R_{+} :

F (x) = Min (1, e^{Mx_{mn}})  where M = 2(\frac{n − 1}{r_{m}} + \frac{m − 1}{c_{n}})  (7)

which is 1 on P (r, c) and falls off exponentially outside. The corresponding function on QQ′ is given by :

G(y) = Min (1, e^{M(\frac{1}{N} ∑_{i,j} ρ_{ij}y_{ij}−∑_{i} r_{i}+c_{n})}) .  (8)

We call a cube of the form {y : 0.4s_{ij} ≤ y_{ij} < 0.4(s_{ij} + 1)} in y space where s_{ij} are all integers a “lattice cube” (it is a cube of side 0.4); its center p where p_{ij} = 0.4s_{ij} + 0.2 will be called a lattice cube center (lcc). Let L be the set of lattice cubes that intersect QQ′(r, c). We will interchangeably use L to denote the set of lcc’s of such cubes. (Note that the lcc itself may not be in QQ′.) If y is an lcc, we denote by C(y) the lattice cube of which it is the center. Note that for a particular lcc y, it is easy to check whether it is in L - we just round down all coordinates (to an integer multiple of 0.4) and check if the point so obtained is in QQ′. In our algorithm, each step will only modify one coordinate of y; in this case, by keeping running row and column sums, we can check in O(1) time whether the new y is in L.

A simple calculation shows :

Proposition 2 : For any lcc y and any z ∈ C(y), we have e^{−0.8}G(y) ≤ G(z) ≤ e^{0.8}G(y).

Now we are ready to present the algorithm. There are two steps of the algorithm of which the first is the more important. Some details of execution of the algorithm are given after the algorithm. In the second step, we may reject in one of three places. If we do so, then, this trial does not produce a final sample.

The Algorithm

The first paragraph of this section describes the input / output behavior of the algorithm. See Remark following the Theorem for the choice of T_{0}, T_{1}.

I. Let T_{0}, T_{1} ≥ 1. Their values are to determined by the accuracy needed.

Run the random walk Y^{(1)}, Y^{(2)}, … Y^{(T_{0}+T_{1})} with L as state space described below for T_{0} + T_{1} steps starting at any state of L.

II. for i = T_{0} + 1 to i = T_{0} + T_{1} do the following :

a. Pick Z^{(i)} with uniform density from C(Y^{(i)}). Reject with probability 1 − \frac{G(Z^{(i)})}{e^{.8}G(Y^{(i)})}. [Reject means end the execution for this i.]

b. Reject if Z^{(i)} is not in QQ(r, c).

c. Let p^{(i)} be the point in P (r, c) corresponding to Z^{(i)}. If p^{(i)} is in a F (w) for some w ∈ I(u, v), then add w^{(i)} = w to the set W of samples to be output; otherwise reject.

□

Random walk for step I

Step I is a random walk on L. Two lcc ’s in L will be called “adjacent” if they differ in one coordinate by 0.4 and are equal on all the other coordinates. The transition probabilities for the random walk for step I are given by :

Pr(y → y′) = \frac{1}{2N} Min (1, \frac{G(y′)}{G(y)}) for y, y′ adjacent lcc ’s in L

Pr(y → y′) = 0 for y, y′ not adjacent and

Pr(y → y) = 1 − ∑_{y′∈L\{y}} Pr(y → y′).

This is an irreducible time-reversible Metropolis Markov Chain and from the identity

G(y) Pr(y → y′) = G(y′) Pr(y′ → y),

it follows that the steady state probabilities, π, exist and are proportional to G(·).

An execution of the random walk produces Y^{(i)}, i = 1, 2, … where Y^{(i+1)} is chosen from Y^{(i)} according to the transition probababilities.

In the next section, we will use the results of Frieze, Kannan, Polson [[12]] to show that the second largest absolute value of an eigenvalue of the Markov Chain (we denote this quantity by θ) satisfies (see Theorem 3):

(1 − θ)^{−1} ≤ N^{4}(n + m − 2)(35 + o(1)),

where the o(1) term goes to zero as n + m, N → ∞. As pointed out there, for N ≥ 10 and n + m ≥ 11, we have (1 − θ)^{−1} ≤ 36N^{4}(n + m − 2). The quantity 1 − θ is called the “spectral gap”.

The following theorem is proved using a result of Aldous [[?]]. It involves a quantity A_{1} which is the limiting probability of acceptance in the above algorithm. A sequence of lemmas (4-7) is then used to bound A_{1} from below. This lower bound may be just plugged into the Theorem. See Remark below for a more practical way to estimate A_{1}.

Theorem 2 The set of samples W produced by the algorithm satisfies :

if g : I(u, v) → R is any function, with g_{0} = max_{w,x∈I(u,v)} g(w) − g(x), and the mean ḡ = ^{\frac{∑_{w∈I(u,v)} g(w)}{|I(u,v)|}}, we have

E((\frac{1}{|W |} ∑_{w∈W} g(w)−ḡ)^{2}) ≤ \frac{1}{T_{1}A_{1}} \frac{(11.5 − 900α(u, v) log θ)}{− log θ} (1 + θ^{T_{0}}e^{5N}(5N)^{N}e^{4}) g_{0}^{2}.

Further,

A_{1} ≥ e^{−2.8− \frac{1.6}{N}} /α(u, v).

Remark To use the theorem, we must choose T_{0}, T_{1} to be large enough. This can be done in a fairly standard manner : suppose we want

E((\frac{1}{|W |} ∑_{w∈W} g(w) − ḡ)^{2})) ≤ ϵ.

Then, note that using Theorem [?], we can choose

T_{0} = (5N + N log(5N) + \frac{4}{N})(35 + o(1))N^{4}(n + m − 2)

T_{1} = \frac{(840 + o(1))}{A_{1}}g_{0}^{2}N^{4}(n + m − 2)/ϵ.

We may estimate A_{1} from the run of the algorithm which potentially fares much better than the lower bound we have given.

Proof of the Theorem

Let S = ∑_{w∈W} (g(w) − ḡ). Let f : L × I(u, v) → [0 1] be such that f(y, w) is the probability that given Y^{(i)} = y in step I of the algorithm, we pick w^{(i)} = w. Letting τ denote the linear transformation such that for each x ∈ P (r, c), τ(x) = z gives us the corresponding point in QQ, (the exact description of τ is not needed here), we have

f(y, w) = (2.5)^{N} ∫_{z∈τ(F (w))∩C(y)} \frac{G(z)}{e^{.8}G(y)}dz.

Let h(y) = ∑_{w∈I(u,v)} f(y, w)(g(w)−ḡ). [This is the expected value of g(w^{(i)})−ḡ given Y^{(i)} = y where we count a rejected Y^{(i)} as giving us the value 0.]

Let

S′ = ∑_{i=T_{0}+1}^{T_{0}+T_{1}} h(Y^{(i)}).

[Note that the expectation of S given the Y^{(i)} is precisely S′.] Note that

∑_{y} π(y)h(y) = ∑_{y} ∑_{w} π(y)f(y, w)(g(w)−ḡ) = ∑_{y} \frac{G(y)}{∑_{L} G} ∑_{w} ∫_{z∈τ(F (w))∩C(y)} \frac{G(z)}{e^{.8}G(y)}(g(w)−ḡ)dz

= \frac{1}{e^{.8} ∑_{L} G} ∑_{w}(g(w) − ḡ) ∫_{z∈τ(F (w))} G(z)dz = ∑_{w} \frac{vol(τ(F (w)))}{e^{.8} ∑_{L} G}(g(w) − ḡ) = 0

because G(z) is 1 on all of QQ which contains τ(F (w)); also the volume of τ(F (w)) is the same for all w ∈ I.

If we had started the chain in the steady state, the expected value of S′ would be T_{1} ∑_{y} π(y)h(y) = 0. Thus as T_{1} tends to infinity, the limit of S′/T_{1} is 0. Let A(y) be the probability of accepting y. Let A_{1} = ∑_{y} π(y)A(y). We also need the variance of h (wrt π) denoted h_{2} = ∑_{y} π(y)(h(y))^{2}.

h_{2} = ∑_{y} π(y)(h(y))^{2} ≤ g_{0}^{2} ∑_{y} π(y)(A(y))^{2} ≤ g_{0}^{2}A_{1}

It follows from Aldous’s Proposition 4.2 (using the fact that his τ_{e}/T_{1} will be less than 1 here which implies that his α(τ_{e}/T_{1}) is at most 2(τ_{e}/T_{1})(1 + e^{−1}))

E((S′)^{2}) ≤ T_{1} (1 + \frac{θ^{T_{0}}}{π_{∗}} ) \frac{2(1 + e^{−1})}{(− log θ)}g_{0}^{2}A_{1},

where, π_{∗} is the minimum of all π(y), y ∈ L. [We note that Aldous lets T_{0} actually be a random variable to avoid the dependence on negative eigenvalues of the chain. In our case, we will be able to get an upper bound on θ, the second largest absolute value of an eigenvalue, so this is not necessary. A simple modification of Aldous’s argument which we do not give here implies the above inequality even though our T_{0} is a deterministic quantity.]

It also follows by similar arguments that

E((|W | − T_{1}A_{1})^{2}) ≤ T_{1} (1 + \frac{θ^{T_{0}}}{π_{∗}} ) \frac{2(1 + e^{−1})}{(− log θ)}A_{1}.  (9)

Consider now S. Recall that

E(S|Y^{(i)}, i = 1, 2, … T_{0} + T_{1}) = S′.

Given Y^{(i)}, i = 1, 2, … T_{0}+T_{1}, the process of producing each w^{(i)} is independent. So, we have

Var(S|Y^{(i)}, i = 1, 2, … T_{0}+T_{1}) = ∑_{i} Var(g(w^{(i)})−ḡ|Y^{(i)}) ≤ E((g(w^{(i)})−ḡ))^{2} ≤ T_{1}g_{o}^{2}.

So we have, (using the inequality (A + B)^{2} ≤ 11A^{2} + 1.1B^{2})

E(S^{2}) = E(((S − S′) + S′)^{2}) ≤ 11E((S − S′)^{2}) + 1.1E((S′)^{2})
≤ 1.1(1 − 135α(u, v) log θ)T_{1}A_{1} (1 + \frac{θ^{T_{0}}}{π_{∗}} ) \frac{2(1 + e^{−1})}{(− log θ)}g_{0}^{2},  (10)

where we have used a lower bound on A_{1} which we derive now. First note that

A_{1} = \frac{vol(τ(F (w)))}{e^{.8} ∑_{L} G}|I(u, v)|.

We have ∑_{L} G ≤ e^{.2} ∫_{T} G, where T is the union of cubes in L.

In a sequence of lemmas below, we show that ∫_{T} G ≤ e^{2.8+ \frac{1.6}{N}} vol(QQ). So, using Lemma 3,

A_{1} ≥ e^{−2.8−1.6/N} \frac{vol(∪_{I(u,v)}τ(F (w)))}{vol(QQ)} = e^{−2.8−1.6/N} \frac{vol(∪_{I(u,v)}F (w))}{vol(P (r, c))} ≥ \frac{e^{−2.8−1.6/N}}{α(u, v)} .

We will now use equations 9 and 10 to argue the conclusion of the Theorem. To this end, let

E ((\frac{1}{|W |} ∑_{w∈W} g(w) − ḡ)^{2} | |W | = s) = γ(s)

and let Pr(|W | = s) = ν(s). Then equation 9 gives us

∑_{s} ν(s)(s − T_{1}A_{1})^{2} ≤ T_{1}A_{1} (1 + \frac{θ^{T_{0}}}{π_{∗}} ) \frac{2(1 + e^{−1})}{(− log θ)}

which implies

∑_{s} ν(s)γ(s)(s − T_{1}A_{1})^{2} ≤ T_{1}A_{1} (1 + \frac{θ^{T_{0}}}{π_{∗}} ) \frac{2(1 + e^{−1})}{(− log θ)}g_{0}^{2}

Now 10 gives

∑_{s} ν(s)γ(s)s^{2} ≤ 1.1(1 − 135α log θ)T_{1}A_{1} (1 + \frac{θ^{T_{0}}}{π_{∗}} ) \frac{2(1 + e^{−1})}{(− log θ)}g_{0}^{2}.

Using the inequality 2(s^{2} + (s − T_{1}A_{1})^{2}) ≥ T_{1}^{2}A_{1}^{2}, we get

∑_{s} ν(s)γ(s) ≤ (11.5 − 900α log θ)\frac{1}{T_{1}A_{1}} (1 + \frac{θ^{T_{0}}}{π_{∗}} ) \frac{1}{− log θ}g_{0}^{2},

which gives us the theorem using proposition [?].

□

Lemma 4 For any real number t (positive, negative or zero), let

K(t) = P′(r, c) ∩ {x : x_{mn} = t} and u(t) = Vol_{N−1}K(t)

and v(t) = Vol_{N−1}(K(t)).

For t_{1} ≤ t_{2} < Min(r_{m}, c_{n}), we have

\frac{v(t_{1})}{v(t_{2})} ≤ (\frac{r_{m} − t_{1}}{r_{m} − t_{2}} )^{n−2} (\frac{c_{n} − t_{1}}{c_{n} − t_{2}} )^{m−2} .

Proof : For a real n+m−4 vector λ = λ_{1n}, λ_{2n} … λ_{m−2,n}, λ_{m1}, λ_{m2}, … λ_{m,n−2} (we assume our λ vector is indexed as above), define K(t|λ) as the set of tables x satisfying

x_{mj} = λ_{mj} for j = 1, 2, … n − 2; x_{m,n−1} = r_{m} − t − ∑_{j=1}^{n−2} λ_{mj};

x_{in} = λ_{in} for i = 1, 2, … m − 2; x_{m−1,n} = c_{n} − t − ∑_{i=1}^{m−2} λ_{in}.

K(t|λ) is the set of tables in K(t) with their last row and column dictated by λ. Let us denote r_{m}−t−∑_{j=1}^{n−2} λ_{mj} by λ_{m,n−1}(t) and c_{n}−t−∑_{i=1}^{m−2} λ_{in} by λ_{n,m−1}(t). Define Λ(t) to be the set of nonnegative vectors λ satisfying :

∑_{j=1}^{n−2} λ_{mj} + t ≤ r_{m}; ∑_{i=1}^{m−2} λ_{in} + t ≤ c_{n};

λ_{mj} ≤ c_{j} for j = 1, 2, … n−2; λ_{m,n−1}(t) ≤ c_{n−1}; λ_{in} ≤ r_{i} for i = 1, 2, … m−2; λ_{m−1,n} ≤ r_{m−1}. Then, we have that K(t|λ) is nonempty iff λ belongs to Λ(t). In general, K(t|λ) is an (m − 2)(n − 2) dimensional set and by the volume of K(t|λ), we mean its (m − 2)(n − 2) dimensional volume. We have

Vol(K(t|λ)) = Vol(P (r(t, λ), c(t, λ))), where,

r(t, λ)_{i} = r_{i}−λ_{in} for i = 1, 2, … m−2; c(t, λ) = c_{j}−λ_{mj} for j = 1, 2, … n−2; r(t, λ)_{m−1} = r_{m−1} − λ_{m−1,n}(t); c(t, λ)_{n−1} = c_{n−1} − λ_{m,n−1}(t).

Consider the linear transformation τ given by

(τ(λ))_{mj} = \frac{r_{m} − t_{2}}{r_{m} − t_{1}} λ_{mj} for j = 1, 2 … n − 2;

(τ(λ))_{in} = \frac{c_{n} − t_{2}}{c_{n} − t_{1}} λ_{in} for i = 1, 2 … m − 2.

It is easy to see that τ is a 1-1 map of Λ(t_{1}) into Λ(t_{2}).

Vol(K(t_{1})) = ∫_{λ∈Λ(t_{1})} Vol(K(t_{1}|λ))dλ = ∫_{α∈Λ(t_{2})} Vol(K(t_{1}|τ^{−1}α))| det(τ^{−1})|dα

= (\frac{r_{m} − t_{1}}{r_{m} − t_{2}} )^{n−2} (\frac{c_{n} − t_{1}}{c_{n} − t_{2}} )^{m−2} ∫_{α} Vol(P (r(t_{1}, τ^{−1}α), c(t_{1}, τ^{−1}α)))dα.

It is easy to check that r(t_{1}, τ^{−1}α) ≤ r(t_{2}, α) and c(t_{1}, τ^{−1}α) ≤ c(t_{2}, α). This implies that the integrand in the last integral is bounded above by Vol(P (r(t_{2}, α), c(t_{2}, α))). This function of course integrates to Vol(K(t_{2})) completing the proof.

□

Lemma 5 For any t_{1} < 0, we have v(t_{1}) (defined in lemma 4) satisfies

v(t_{1}) ≤ \frac{e^{2}}{ρ_{mn}} (\frac{r_{m} − t_{1}}{r_{m}} )^{n−2} (\frac{c_{n} − t_{1}}{c_{n}} )^{m−2} Vol_{N}(QQ(r, c))

where ρ_{mn} = Min(\frac{r_{m}}{n−1}, \frac{c_{m}}{m−1}).

Proof For t_{2} in the range [0 ρ_{mn}], we have

Vol(K(t_{2})) ≥ (\frac{r_{m} − t_{2}}{r_{m} − t_{1}} )^{n−2} (\frac{c_{n} − t_{2}}{c_{n} − t_{1}} )^{m−2} Vol(K(t_{1}))

≥ (\frac{r_{m}}{r_{m} − t_{1}} )^{n−2} (\frac{c_{n}}{c_{n} − t_{1}} )^{m−2} e^{−2}vol(K(t_{1})).

Integrating this over this range of t_{2}, we get the lemma.

Lemma 6 :

∫_{QQ′(r,c)} G(y)dy ≤ e^{2} ∫_{QQ(r,c)} G(y)dy ≤ e^{2}vol(QQ).

Proof : From the last lemma, we have for any t_{1} < 0,

v(t_{1}) ≤ e^{\frac{−(n−2)t_{1}}{rm} + \frac{−(m−2)t_{1}}{cn}} \frac{e^{2}}{ρ_{mn}} Vol(QQ).

Now the lemma follows by integration.

□

Lemma 7 Let T be the union of C(y) over all lcc ’s y. Then we have

∫_{T} G(y)dy ≤ e^{2+ \frac{1.6}{N}} ∫_{QQ(r,c)} G(y)dy.

Proof : We have :

T ⊆ \{y : \frac{1}{N} ∑_{j} ρ_{ij}y_{ij} ≤ r_{i}+\frac{.4}{N} ∑_{j} ρ_{ij}∀i \frac{1}{N} ∑_{i} ρ_{ij}y_{ij} ≤ c_{j}+\frac{.4}{N} ∑_{i} ρ_{ij}∀j y_{ij} ≥ 0\}.

(11)

This implies that T ⊆ (1 + \frac{.4}{N})QQ′. Note also that for any y in QQ′, we have G((1 + \frac{.4}{N})y) ≤ G(y)e^{0.4M ∑_{ij} min(∑_{i} r_{i},∑_{j} c_{j})/(N^{2})} ≤ G(y)e^{1.6/N}. Thus we have the lemma using the last lemma.

□

Proposition 3

π_{∗}^{−1} ≤ e^{5N}(5N)^{N}e^{4}.

Proof The number of states is at most (5N)^{N}. It is easy to check that the ratio of the maximum value of G to the minimum over lcc ’s is at most e^{5N+4} completing the proof.

□

Algorithm to estimate the number of contingency tables

Using the sampling algorithm, we will be able to estimate the number of contingency tables in time polynomial in the data and α(u, v). We will only sketch the method here.

We first estimate by sampling from P (r, c) the following ratio :

\frac{vol(P (r, c))}{vol(∪_{w∈I(u,v)}F (w))}.

Next, we define a sequence of contingency polytopes P_{1}, P_{2}, …, each obtained from the previous one by increasing r_{m} and c_{j} (j = 1 ⋯ n) by ⌊min (\frac{1}{m}c_{n}, \frac{1}{n}r_{m})⌋ until we have c_{j} ≥ ∑_{i=1}^{m−1} r_{i} for j = 1 ⋯ n. Then |I(r, c)| = ∏_{i=1}^{m−1} \binom{r_{i}+n−1}{n−1} and the volume of P (r, c) is read off as the coefficient of β^{N} in the easily calculated polynomial |I(β × r, β × c)| (= ∏_{i=1}^{m−1} (r_{i}^{n−1}/((n − 1)!))).

6 Bound on the spectral gap

We refer the reader to Diaconis and Strook [5] or [12]for a discussion of the eigenvalues of a reversible Markov Chain. The second largest absolute value of such a chain gives us a bound on the “time to converge” to the steady state as described in these references. Here, we use Theorem 2 of [12] to bound the second largest absolute value of our chain. Note that though their theorem does not explicitly state so, the λ_{2} in that Theorem is the second largest absolute value of any eigenvalue of the chain. This section is a technical evaluation of the various quantities needed to plug into the expression for λ_{2} in [12]; we do not redefine these quantities here.

We first calculate the diameter of T (the largest Euclidean distance between two points in T). To do so, we let

I_{1} = \{(ij) : ρ_{ij} = \frac{r_{i}}{n − 1}\} I_{2} = \{(ij) : ρ_{ij} < \frac{r_{i}}{n − 1}\}.

Then, for any y ∈ T, (using the fact that (y − .4\underline{1}) ∈ QQ′,)

∀i, ∑_{\{j:(ij)∈I_{1}\}} y_{ij} ≤ N(n − 1)(1 + \frac{.4}{N}) and

∀j, ∑_{\{i:(ij)∈I_{2}\}} y_{ij} ≤ N(m − 1)(1 + \frac{.4}{N}).

So,

∑_{ij} y_{ij}^{2} = ∑_{I_{1}} y_{ij}^{2} + ∑_{I_{2}} y_{ij}^{2} ≤ N^{2}(1 + \frac{.4}{N})^{2}[(n − 1)^{2}(m − 1) + (m − 1)^{2}(n − 1)].

So the diameter d = d(T) of T satisfies

d(T) ≤ \sqrt{2}N^{1.5}\sqrt{(n + m − 2)}(1 + \frac{.4}{N}).

For every unit (length) vector u ∈ R^{N}, with u_{ij} ≥ 0, let l(u) be the ray {y = λu : λ ≥ 0} from the origin along u. Note that the ray intersects T in a line segment (since y ∈ T iff (0.4)⌊2.5 y⌋ ∈ QQ′ where floor denotes componentwise floor.) . Let R = R(u) be the length of the segment l(u) ∩ QQ′(r, c) and R_{1} = R_{1}(u) be the length of the segment l(u) ∩ T. Then there exists an i_{o} ∈ {1, 2, … m − 1} such that

∑_{j} Ru_{i_{o}j} \frac{ρ_{i_{o}j}}{N} = r_{i_{o}}

or a j_{o} ∈ \{1, 2, … n − 1\} such that ∑_{i} Ru_{ij_{o}} \frac{ρ_{ijo}}{N} = c_{j_{o}}. Assume without loss of generality the first option. Since R_{1}u belongs to T, the vector R_{1}u − \underline{0.4} belongs to QQ′. So, we also have,

∑_{j} R_{1}u_{i_{o}j} \frac{ρ_{i_{o}j}}{N} ≤ r_{i_{o}} + (0.4) ∑_{j} \frac{ρ_{i_{o}j}}{N} .

Using the fact that ∑_{j} ρ_{ij} ≤ r_{i}, we get that

R_{1} ≤ R(1 + \frac{.4}{N}).

The above implies that R_{1} ≤ 2R as required by (7) of [[12]].

Also, Ru belongs to QQ′ implies

R ∑_{j} u_{ij} \frac{ρ_{ij}}{N} ≤ r_{i}∀i R ∑_{i} u_{ij} \frac{ρ_{ij}}{N} ≤ c_{j}∀j

which implies by [2]

(R_{1} − R) ∑_{i,j} u_{ij} \frac{ρ_{ij}}{N} ≤ \frac{R_{1} − R}{R} ∑_{i} r_{i} ≤ \frac{0.4r_{m}}{n − 1} .

Similarly (R_{1}−R) ∑_{i,j} u_{ij} \frac{ρ_{ij}}{N} is also at most 0.4c_{n}/m−1. Thus e^{M(R_{1}−R) ∑_{i,j} u_{ij} \frac{ρ_{ij}}{N}} is at most e^{1.6}. This implies that we may take κ_{1} of (8) of [[12]] to be 25/3.

Also, R_{1} − R ≤ .4R/N. So κ_{2} of [[12]] is given by (since δ, the step size is 0.4 here)

κ_{2} ≤ \frac{50}{3}(\frac{R}{N} + \sqrt{N}) = (50/3)(\sqrt{N} + o(1)).  (12)

Also,

κ_{0} = \frac{1}{.8}N^{3}(n + m − 2)(1 + o(1)).  (13)

We now want to bound (1 + α) of [12]. To this end, we need to prove upper and lower bounds on

\frac{∫_{z∈C(y)} G(z) dz}{G(y)δ^{N}}

where y is an lcc. Since G is a convex function, it is clear that 1 is a lower bound on this ratio. To get an upper bound, we use Heoffding’s bounds [15] on the probability that the sum ∑_{ij} ρ_{ij}[z_{ij} − y_{ij}] deviates from its mean of 0 and (after some complicated integration) arrive at an upper bound of 1.05e^{8/N^{1.5}}.

Plugging all this into the formula for λ_{2}^{−1} in Theorem 2 of [[12]], we get :

Theorem 3

(1 − θ)^{−1} ≤ [1.2e^{24/N^{1.5}}]\frac{N}{.4} [\frac{N^{3}}{.8} (n + m − 2)(1 + \frac{1.5}{N})(28/3)

+25\sqrt{N}\sqrt{n + m − 2}(1+\frac{.8}{3N})+2N^{3/2}\sqrt{n + m − 2}] ≤ (35+o(1))N^{4}(n+m−2),

For N ≥ 10, n+m ≥ 11, a calculation shows that the (35+o(1)) may be replaced by 36.

References

[1] D. Aldous, “Some inequalities for reversible Markov chains,” Journal of the London Mathematical Society (2)25, pp. 564-576, 1982.

[2] A. I. Barvinok, “A polynomial time algorithm for counting integral points in polyhedra when the dimension is fixed”, Proceedings of the Symposium on the Foundations of Computer Science, 1993.

[3] P. Diaconis, B. Efron, “Testing for independence in a two-way table”, Annals of Statistics, 13, pp 845-913, 1985.

[4] P. Diaconis, and A. Gangolli, “Rectangular arrays with fixed margins”, To appear in the Proceedings of the workshop on Markov Chains, Institute of Mathematics and its Applications, University of Minnesota, (1994).

[5] P. Diaconis, D. Stroock, “Geometric bound for eigenvalues of Markov chains,” Annals of Applied Probability 1, pp. 36-61, 1991.

[6] P. Diaconis, and B. Sturmfels, “Algebraic algorithms for sampling from conditional distributions”, Technical Report, Department of Statistics, Stanford University (1993)

[7] M. E. Dyer and A. M. Frieze, On the complexity of computing the volume of a polyhedron, SIAM J. Comput., 17 (1988), 27–37.

[8] A. Schrijver, Theory of linear and integer programming, John Wiley, Chichester, 1986.

[9] R. P. Stanley, Enumerative combinatorics, Vol. I, Wadsworth and Brooks/Cole, Monterey, CA, 1986.

[10] M. E. Dyer, A. Frieze, R. Kannan, A. Kapoor, L. Perkovic, U. Vazirani, “A mildly exponential time algorithm for approximating the number of solutions to a multidimensional knapsack problem,”, Combinatorics, Probability and Computing, (to appear).

[11] J. A. Fill, “Eigenvalue bounds on the convergence to stationarity for nonreversible Markov Chains, with an application to the exclusion process”, The Annals of Applied Probability, Vol. 1, No. 1, pp 62-87, 1991.

[12] A. Frieze, R. Kannan and N. Polson, “Sampling from log-concave distributions”, The Annals of Applied Probability, Vol. 4, No. 3 (1994) pp 812-837.

[13] A. R. Gangolli, Convergence Bounds for Markov Chains and Applications to Sampling, Phd. Thesis, Stanford STAN-CS-91-1361, 1991.

[14] D. Gillman, “A Chernoff bound for random walks on expander graphs,” Proceedings of the 34th Annual IEEE Conference on Foundations of Computer Science, 1993.

[15] Hoeffding, W. Probability Inequalities For Sums of Bounded Random Variables, American Statistical Association Journal, March 1963, pp. 13-30.

[16] N. Kahale, “Large deviation bounds on Markov Chains”, DIMACS Technical Report, 94-39, 1994.

[17] C. R. Mehta, N. R. Patel, “A Network Algorithm for Performing Fisher’s Exact Test in r × c Contingency Tables,” Journal of the American Statistical Association 78, pp. 427-434, 1983.

[18] E. Seneta, “Non-negative matrices and Markov Chains” Springer-Verlag , New York, 1973.

[19] B. Sturmfels, “How to Count Integer Matrices,” manuscript, December 13, 1993.

[20] L. G. Valiant, “The complexity of computing the permanent,” Theoretical Computer Science 8, pp. 189-201, 1979.