On sampling with Markov chains

F. R. K. Chung University of Pennsylvania Philadelphia, PA 19104

R. L. Graham AT&T Bell Laboratories Murray Hill, NJ 07974

S.-T. Yau Harvard University Cambridge, MA 02138

1. Introduction

There are many situations in which one would like to select a random element (nearly) uniformly from some large set. One method for doing this which has received much attention recently is the following. Suppose we can define a Markov chain M on X which has the uniform distribution as its stationary distribution. Then starting at some (arbitrary) initial point x_{0}, and applying M for sufficiently many (say t) steps, the resulting point M^{(t)}(x_{0}) will be (almost) uniformly distributed over X, provided t is large enough. In order for this approach to be effective, however, M should be “rapidly mixing”, i.e., M^{(t)} should be very close to the uniform distribution U (in some appropriate sense) in polynomially many steps t, measured by the size of X.

There are in use a variety of methods for estimating the “mixing rate” of M, i.e., the rate of convergence of M to its stationary distribution, e.g., coupling, path-counting, strong stationary time, eigenvalue estimates and even Stein’s method (cf. [A87], [SJ89], [S93], [DFK89], [FKP94]). However, for many problems of current interest such as volume estimation, approximate enumeration of linear extensions of a partial order and sampling contingency tables, eigenvalue methods have not been effective because of the difficulty in obtaining good (or any!) bounds on the dominant eigenvalue of the associated process M.

In this paper we remedy this problem to some extent by applying new eigenvalue bounds of two of the authors [CY1], [CY2] for random walks on what we call “convex subgraphs” of homogeneous graphs (to be defined in later sections).

A principal feature of these techniques is the ability to obtain eigenvalue estimates in situations in which there is a nonempty boundary obstructing the walk, typically a more difficult situation then the boundaryless case. In particular, we will give the first proof of rapid convergence for the “natural” walk on contingency tables, as well as generalizations to restricted contingency tables, symmetric tables, compositions of an integer, and so-called knapsack solutions. We point out that these methods can also be applied to other problems, such as selecting random points in a given polytope (which is a crucial component of many recent volume approximation algorithms for polytopes). Also, the same ideas can be carried out in the context of weighted graphs and biased random walks (see [CY1] for the basic statements). However, we have deliberately restricted ourselves to some of the simpler applications for ease of exposition.

Before proceeding to our main results, it will be necessary to introduce a certain amount of background material. This we do in the next section.

2. Background

We begin by recalling a variety of concepts from graph theory and differential geometry. Any undefined terms can be found in [SY94] or [BM76].

For a graph G = (V, E) with vertex set V and edge set E, we let d_{v} denote the degree of the vertex u ∈ V, i.e., the number of w ∈ V such that {u, w} ∈ E is an edge (which we also indicate by writing v ∼ w). For the most part, our graphs will be undirected, and without loops or multiple edges. Define the square matrix L = L_{G}, with rows and columns indexed by V by:

(1) L(u, v) = \{ \begin{matrix}d_{u} & if u = v, \\ −1 & if u ∼ v, \\ 0 & otherwise . \end{matrix}

Let T denote the diagonal matrix with (v, v) entry T (v, v) = d_{v}, v ∈ V (and, of course, 0 otherwise).

A key concept for us will be the following operator: The Laplacian L = L_{G} for G is defined by L := T^{−1/2}LT^{−1/2} . (2)

Thus, L can be represented as a |V | by |V | matrix given by

(3) L(u, v) = \{ \begin{matrix}1 & if u = v, \\ −\frac{1}{\sqrt{d_{u}d_{v}}} & if u ∼ v, \\ 0 & otherwise . \end{matrix}

Considering L as an operator acting on {f : V → ℂ} we have:

Lf(v) = \frac{1}{\sqrt{d_{v}}} ∑_{{u∈V \atop u∼v}} (\frac{f(v)}{\sqrt{d_{v}}} − \frac{f(u)}{\sqrt{d_{u}}}) . (4)

In the case that G is d-regular, i.e., d_{v} = d for all v ∈ V, (3) and (4) take the more familiar forms

(3′)  L = I − \frac{1}{d}A

where I is the identity matrix and A = A_{G} is the adjacency matrix for G, and

(4′)  Lf(v) = \frac{1}{d} ∑_{{u \atop u∼v}}(f(v) − f(u)) .

Denote the eigenvalues of L by

0 = λ_{0} ≤ λ_{1} ≤ ⋯ ≤ λ_{n−1}

where we let n denote |V | for now. It is known that λ_{1} ≤ 1 unless G = K_{n}, the complete graph on n vertices, and λ_{n−1} ≤ 2 with equality if and only if G is bipartite.

The most important of these eigenvalues from our perspective will be λ_{1} which we call the Laplacian eigenvalue of G, and denote by λ = λ_{G}. Note that with 1 : V → ℂ denoting the constant 1 function, T^{1/2}1 is an eigenfunction of L with eigenvalue 0. It follows (see [CY1]) that

λ := λ_{1} = inf_{f⊥T 1} ∑_{u∼v} \frac{(f(v) − f(u))^{2}}{∑_{v} d_{v}f(v)^{2}}

(5)

= inf_{f} sup_{c} \frac{∑_{u∼v}(f(v) − f(u))^{2}}{∑_{v} d_{v}(f(v) − c)^{2}} .

For a (nonempty) set S ⊂ V, consider the subgraph of G induced by S. This subgraph, which we also denote by S, has edge set E(S) = {{u, v} ∈ E(G) : u, v ∈ S}. The edge boundary ∂S is defined by

∂S := {{u, v} ∈ E(G) : u ∈ S, u ∈ V \ S} .

The vertex boundary δS is defined by

δS := {v ∈ V \ S : {u, v} ∈ E(G) for some v ∈ V } .

S V \ S ∂S  δS u v

Figure 1: Various boundaries of S

Let E′ := E(S) ∪ ∂S.

We illustrate these definitions in Figure 1.

The next key concept we will need is that of the Neumann eigenvalues for S on G. First, define (parallel to (5))

λ_{S} = inf_{{f \atop ∑_{x∈S} dxf(x)=0}} \frac{∑_{\{x,y\}∈E′}(f(x) − f(y))^{2}}{∑_{x∈S} d_{x}f^{2}(x)}

(6)

= inf_{f} sup_{c} \frac{∑_{\{x,y\}∈E′}(f(x) − f(y))^{2}}{∑_{x∈S} d_{x}(f(x) − c)^{2}} .

In general, define

λ_{S,i} = inf_{f} sup_{f′∈C_{i−1}} \frac{∑_{\{x,y\}∈E′}(f(x) − f(y))^{2}}{∑_{x∈S} d_{x}(f(x) − f′(x))^{2}} (7)

where C_{k} is the subspace spanned by the eigenfunctions ϕ_{i} achieving λ_{S,i}, 0 ≤ i ≤ k. In particular,

λ_{S} = λ_{S,1} = inf_{g⊥T^{1/2}1} \frac{⟨g, Lg⟩_{S}}{⟨g, g⟩_{S}} (8)

where L is the Laplacian of G and ⟨f_{1}, f_{2}⟩_{S} denotes the inner product ∑_{x∈S} f_{1}(x)f_{2}(x). Note that when S = V, then λ_{S} is just the usual Laplacian eigenvalue of G. The λ_{S,i} are eigenvalues of a matrix L_{S} defined by the following steps.

For X ⊂ V, let L_{X} denote the submatrix of L restricted to rows and columns induced by elements of X. Define the matrix N with rows indexed by S ∪ δS and columns indexed by S given by

(9) N(x, y) = \{ \begin{matrix}1 & if x = y, \\ 0 & if x ∈ S, x ≠ y, \\ \frac{1}{d′_{x}} & if x ∈ δS, x ∼ y ∈ S, \\ 0 & otherwise \end{matrix}

where d′_{x} denotes the degree of x to points in S. Then the matrix

L_{S} := T^{−1/2}N^{tr}L_{S}NT^{−1/2} (10)

has the λ_{S,i} as its eigenvalues (together with 0, of course), with λ_{S} = λ_{S,1} being the least eigenvalue exceeding 0. The matrix L_{S} is associated with the following random walk P on S ⊂ G, which we call the Neumann walk on S. A particle at vertex u moves to each of the neighbors v of u with equal probability 1/d_{u} except if v ∉ S (so that v ∈ δS). In this case, the particle then moves to each of the d′_{v} neighbors of v which are in S, with equal probability \frac{1}{d′_{v}}. Thus, P (u, x), the probability of moving from u to x is

P (u, x) = \{ \begin{matrix}1/d_{u} & if u ∼ x \\ 0 & otherwise \end{matrix} \} + \frac{1}{d_{u}} ∑_{{v∈δS \atop v∼u}} \frac{1}{d′_{v}} .

It is easy to see that when G is regular, the stationary distribution of P is uniform for any choice of S connected. The rate of convergence of P to its stationary distribution is controlled by the Neumann Laplacian eigenvalue λ_{S}, which we eventually will bound from below.

The final key concept we introduce here is the heat kernel H_{t} of S. To begin with, let us decompose L_{S} in the usual way as a linear combination of its projections P_{i} onto the corresponding eigenfunctions ϕ_{i} of L_{S}:

L_{S} = ∑_{i} λ_{S,i}P_{i} . (11)

For t ≥ 0, the heat kernel H_{i} is defined to be the operator of {f : S ∪ δS → ℂ} given by

H_{t} : = ∑_{i} e^{−λ_{S,i}}P_{i}
(12)  = e^{−tL_{S}}
= I − tL_{S} + \frac{t^{2}}{2} L_{S}^{2} − \frac{t^{3}}{6} L_{S}^{3} + ⋯

Thus, H_{0} = I.

For f : S ∪ δS → ℂ, set

F = H_{t}f = e^{−tL_{S}}f .

That is,

F (t, x) = (H_{t}f)(x) = ∑_{y∈S∪δS} H_{t}(x, y)f(y) .

Some useful facts concerning H_{t} and F are the following (see [CY1] for proofs):

(i) F (0, x) = f(x)

(ii) For x ∈ S ∪ δS,

∑_{y∈S∪δS} H_{t}(x, y)\sqrt{d_{y}} = \sqrt{d_{x}}

(iii) F satisfies the heat equation \frac{∂F}{∂t} = −L_{S}F

(iv) For any x ∈ δS, L_{S}F (t, x) = 0

(v) For all x, y ∈ S ∪ δ(S), H_{t}(x, y) ≥ 0

The connection of the heat kernel to the Neumann eigenvalue λ_{S} of S is given by the following result:

Theorem [Chung, Yau [CY1]]. For all t > 0,

λ_{S} ≥ \frac{1}{2t} ∑_{x∈S} inf_{y∈S} \frac{H_{t}(x, y)\sqrt{d_{x}}}{\sqrt{d_{y}}} . (13)

This inequality will be an essential ingredient in eventually obtaining lower bounds for λ_{S} (and upper bounds on the mixing rates) for the various Markov chains we consider.

3. A direct application to heat kernel inequality

The main use of the heat kernel inequality (13) in lower bounding λ_{S} will depend on controlling the behavior of H_{t} by connecting it to an associated continuous heat kernel h_{t} for an appropriate Riemannian manifold M = M_{S} containing the points of S ∪ δS. This we do in Section 4.

However, there are situations in which we can get bounds on λ_{S} directly from (13) without going through this process. We describe several of these in this section.

Let us consider the special case for which S = V (so that L_{S} = L) and, further, suppose the graph G has a “covering” vertex x_{0} with the property that x_{0} is adjacent to every g ∈ V \ {x_{0}} (so that d_{x} = n − 1 where n := |V |. We will apply (13) with t → 0. Thus, by (13)

+O(t^{2}) = H_{t} = I − tL + O(t^{2}), t → 0 x ∼ y \frac{t}{\sqrt{d_{x}d_{y}}} · · · · · · x ≁ y′ x (14)  1 − t  y 0 y′ 1 − t

Thus,

+O(t^{2}) x  y 0 y′ 1 − t  x ∼ y \frac{t}{d_{y}} · · · · · · · · ·  \frac{t}{d_{y}} ⋯ x ≁ y′ 1 − t (H_{t}(x, y)\begin{matrix}√ & \overline{d_{x}} \\ √ & \overline{d_{y}} \end{matrix} ) = (15)  1 − t =  \frac{t}{d_{y′}} x_{0}  +O(t^{2}) y y′ 1 − t

Thus, by (13)

λ_{S} = λ ≥ \frac{1}{2t} ∑_{x} inf_{y} \frac{H_{t}(x, y)\sqrt{d_{x}}}{\sqrt{d_{y}}}
(16)  = \frac{1}{2t} ( inf_{y≠x_{0}} \frac{t}{d_{y}} + O(t^{2}))
= \frac{1}{2t} (\frac{t}{δ_{2}} ) + O(t), t → 0

where δ_{2} denotes the second largest degree in G. Thus, we have

λ ≥ \frac{1}{2δ_{2}} (17)

for this situation. For example, for G = P_{3}, the path with 3 vertices, it is true that λ_{0} = 0, λ = λ_{1} = 1 and λ_{2} = 2, while our estimate in (17) gives λ ≥ 1/2.

A similar analysis shows that of G has k > 1 covering vertices then (13) implies

λ ≥ \frac{k}{2(n − 1)} . (18)

Applying this to G = K_{n}, the complete graph on n vertices yields

λ ≥ \frac{n}{2(n − 1)}

while the truth is λ = \frac{n}{n−1} (again off by a factor of 2).

4. The basic set-up

As hinted at in the previous section, the real use of (13) in boundary λ_{S} will proceed along the following lines:

(a) Embed G and S as “convex” sets into some Riemannian manifold M = M_{S} (usually ≅ 𝔼^{N} for our applications).

(b) Relate the heat kernel H_{t} on S ∪ δS to the continuous heat kernel h_{t} on M.

(c) Relate h_{t} to various properties of M (and S), such as the “density” of points of S in M, the diameter of M, the dimension of M, etc.

We next describe a specific situation to which we can apply this procedure.

We begin with an infinite “lattice” graph Γ = (V, E) with V ⊂ M ≅ 𝔼^{N}, with the property that the automorphism group H of G is transitive. Thus, H : V → V so that u ∼ v ⇔ gu ∼ gv, g ∈ H, and for all u, v ∈ V,

v = gu for some g ∈ H .

Suppose Γ has an edge-generating set K ⊂ H, so that every edge of Γ is of the form {v, gv} for some v ∈ V, g ∈ K, and any such pair is an edge (hence the term lattice graph). We will further assume (since we think of Γ as undirected) that g ∈ K ⇒ g^{−1} ∈ K. Finally, we assume that for any x, any lattice point of y ∈ M in the convex hull of {gx : g ∈ K} is also adjacent to x in Γ. Let S be some finite subset of V. Assume S is “convex”, which means that for some submanifold M ⊂ M with nonempty convex boundary, S consists of all lattice points of V which are in M. Let ℓ denote the minimum length of any of the edges {x, gx}, g ∈ K, and assume that for each x ∈ S, the ball B_{x}(ℓ/3) of radius ℓ/3 centered at x is contained in M.

For x ∈ V, let U(x) denote the Voronoi origin for x, i.e., the set of all points of M closer to x than to any other point y ∈ S. Since H is transitive, all U(x) have the same volume, denoted by vol U. Finally, set

diam M := diameter of M
dim M := dimension of M

Under the preceding assumptions and notation, we have the following estimate.

Theorem [Chung, Yau [CY1]]. For convex S,

(19) λ_{S} ≥ c_{0} (\frac{ℓ}{dim M diam M} )^{2} \frac{|S| vol U}{vol M}

for constant c_{0} > 0 depending only on Γ and not on S.

The proof of (19) depends upon a variety of ideas and techniques from differential geometry, linear algebra and combinatorics (as well as (13)), and can be found in [CY1]. The bulk of our results presented here will depend on applying (19) and its generalizations to specific cases of interest.

We remark that the results in [CY1] are similar in spirit to those in [CY2] which also gives lower bounds for λ_{S} for Neumann walks. There, S is required to satisfy more restrictive conditions (e.g., being “strongly convex”) but the bounds on λ_{S} are sharper. In particular, it is shown that under the appropriate hypotheses on S (and Γ),

λ_{S} ≥ \frac{1}{8kD^{2}} (20)

where k = |K| and D := graph diameter of S (see Section 8 for a comparison of (19) and (20) in a specific case).

5. Contingency tables

Given integer vectors \overline{r} = (r_{1}, …, r_{m}), \overline{c} = (c_{1}, …, c_{n}) with r_{i}, c_{j} ≥ 0 and ∑_{i} r_{i} = ∑_{j} c_{j}, we can consider the space of all m × n arrays T with the property that

∑_{j} T (i, j) = r_{i}, 1 ≤ i ≤ m
∑_{i} T (i, j) = c_{j}, 1 ≤ j ≤ n .

Let us denote by T = T (\overline{r}, \overline{c}) the set of all such arrays. The arrays T ∈ T are often called contingency tables (with given row and column sums). These tables arise in a variety of applications, such as goodness of fit tests in statistics, enumeration of permutations by descents, describing tensor product decompositions, counting double cosets, etc., and have a long history. (An excellent survey can be found in [DG].) It seems to be a particularly difficult problem to obtain good estimates of the size of T (\overline{r}, \overline{c}) for large \overline{r} and \overline{c}. In order to attack this problem, a standard (by now — see [A87], [SJ89], [S93], [G91], [DSt]) technique depends on rapidly generating random tables from T = T (\overline{r}, \overline{c}) with nearly equal probability.

To do this, we first consider the following natural random walk P on T . ¿From any given table T ∈ T , select uniformly at random a pair of rows {i, i′} and a pair of columns {j, j′}, and move to the table T′, obtained from T by changing four entries of T as follows:

T′(i, j) = T (i, j) + 1, T′(i′, j′) = T (i, j) + 1
T′(i′, j) = T (i, j′) − 1, T′(i, j′) = T (i, j′) − 1 .

Such a move we call a basic move. The table T′ clearly has the same line sums (i.e., row and column sums) as T. The only problem is that T′ may have negative entries (because we might have T (i, j′) = 0, for example) and so is not in T . To deal with this “boundary” problem, we instead execute the corresponding Neumann walk in T , as described in Section 2.

We first need to place our contingency table problem into the framework of the preceding section. The manifold M will consist of all real mn-triples \overline{x} = (x_{11}, x_{12}, …, x_{mn}) satisfying

∑_{j} x_{ij} = r_{i}, ∑_{i} x_{ij} = c_{j} .

Since ∑_{i} r_{i} = ∑_{j} c_{j} then

dim M = N := (m − 1)(n − 1) .

The graph Γ has as vertices all the integer points in M, i.e., all \overline{x} with all x_{ij} ∈ ℤ. The edge generating set K consists of all the basic moves described above. Thus, |K| = \binom{m}{2}\binom{n}{2}. The set S will just be T = T (\overline{r}, \overline{c}), the set of all T ∈ Γ with all entries nonnegative. Thus,

S = ⋂_{i,j}{T ∈ Γ : x_{ij} ≥ 0} .

Similarly, the manifold M ⊂ M is defined by

M = ⋂_{i,j}\{\overline{x} ∈ M : x_{ij} ≥ −2/3\} .

It is clear that M is an N-dimensional convex polytope and S = T is the set of all lattice points in M, and consequently convex in the sense needed for (19). It is easy to see that T is connected by the basic moves generated by K, and that each edge of Γ has length 2. Our next problem is to deal with the term \frac{|S| vol U}{vol M} in (19). In particular, we would like to show this is close to 1, provided that r_{i} and c_{j} are not too small. To do this we need the following two results.

Claim 1. Suppose L ⊂ 𝔼^{N} is a lattice generated by vectors v_{1}, …, v_{N}. Then the covering

radius of L is at most R := \frac{1}{2} ( ∑_{i=1}^{N} ∥v_{i}∥^{2})^{1/2}.

Proof. The assertion clearly holds for N = 1. Assume it holds for all dimensions less than N. It is enough to prove that any point \overline{x} = (x_{1}, …, x_{N}) in the fundamental domain generated by the v_{i} is at most a distance of R from some v_{i}. Let \overline{x}_{0} be the projection of \overline{x} on either the hyperplane generated by v_{1}, …, v_{N−1}, or a translate of the hyperplane by v_{N}, whichever is closer (these are two bounding hyperplanes of the fundamental domain). Thus, d(\overline{x}_{1}, \overline{x}_{0}) ≤ \frac{1}{2}∥v_{N}∥. By the induction hypothesis,

d(\overline{x}_{0}, v_{j}) ≤ \frac{1}{2} ( ∑_{i=1}^{N−1} ∥v_{i}∥^{2})^{1/2} for some j < N .

Hence

d(\overline{x}, v_{j}) ≤ (\frac{1}{4}∥v_{N}∥^{2} + \frac{1}{4} ∑_{i=1}^{N−1} ∥v_{i}∥^{2})^{1/2} = R

as claimed.

Claim 2. If M is convex and contains a ball B(cRN) of radius cRN, c > 0, then

e^{−1/c} < \frac{|S| vol U}{vol M} < e^{1/c} (21)

where v_{1}, … v_{N} generate Γ, and R = \frac{1}{2} ( ∑_{v=1}^{N} ∥v_{i}∥^{2})^{1/2}.

Sketch of proof: Consider an enlarged copy (1 + δ)M of M expanded about the center of the ball B(cRN), δ > 0. Let L be some bounding hyperplane of M, and let (1 + δ)L be the corresponding expanded copy of L (see Figure 2)

(1 + δ)L L M  xγ B(cRN) cRN  y

Figure 2: A large ball in M

Let x ∈ S ⊂ M and suppose there exists y ∈ U(x), the Voronoi region for x, with y ∉ (1 + δ)M. Thus

d(x, y) ≥ γ > cδRN .

However, by Claim 1, each point in M has distance at most R from lattice point in M. This is a contradiction if we take δ = \frac{1}{cN}.

Thus, for all x ∈ S, U(x) ⊂ (1 + δ)M and so

|S| vol U ≤ vol(1 + δ)M = (1 + δ)^{N} vol M

i.e.,

\frac{|S| vol U}{vol M} ≤ (1 + δ)^{N} = (1 + \frac{1}{cN} )^{N} < e^{1/c} .

A similar argument shows that

\frac{|S| vol U}{vol M} > e^{−1/c}

and Claim 2 is proved.

In order to apply the result in Claim 2, we must find a large ball in M. Let s_{0} denote the smallest line sum average, i.e.,

s_{0} = min (min_{i} \frac{r_{i}}{n} , min_{j} \frac{c_{j}}{m}) .

We begin constructing an element T_{0} ∈ M recursively as follows. Suppose without loss of generality that

s_{0} = \frac{r_{1}}{m} .

Then, in T_{0}, set all elements of the first row equal to s_{0}, and subtract s_{0} from each value c_{j} to form c′_{j} = c_{j} − s_{0}, 1 ≤ j ≤ n. Now, to complete T_{0}, we are reduced to forming an (m − 1) by n table T′_{0} with row sums r_{2}, … r_{m} and column sums c′_{1}, …, c′_{n}. The key point here is that all the line sum averages for T′_{0} are still at least as large as s_{0}. Hence, continuing this process we can eventually construct a table T_{0} (with rational entries) having least entry equal to s_{0}. Consequently there is a ball B(s_{0}) of radius s_{0} centered at T_{0} ∈ M which is contained in M (since to leave M, some entry must become negative). Therefore, if we assume s_{0} > cN^{3/2} then by (19) and (21),

λ_{S} ≥ \frac{c_{0}e^{−1/c}}{N^{2}(diam M)^{2}} (22)

for some absolute constant c_{0} > 0 (since for tables, all the generators have length 2, so that R ≤ N^{1/2}).

In the Appendix, we illustrate how a specific value can be derived here for c_{0} (as well as in several other cases of interest, as well). In particular, for contingency tables, we can take c_{0} = 1/800.

Since

diam M < 2 min {(∑_{i} r_{i}^{2})^{1/2}, (∑_{j} c_{j}^{2})^{1/2}}

then (22) can be written as follows:

For the natural Neumann walk P on the space of tables T (\overline{r}, \overline{c}) where

min \{min_{i} \frac{r_{i}}{n} , min_{j} \frac{c_{j}}{m}\} > c(m − 1)^{3/2}(n − 1)^{3/2}

we have

λ_{S} > [3200e^{1/c}(m − 1)^{2}(n − 1)^{2} min {∑_{i} r_{i}^{2}, ∑_{j} c_{j}^{2}}]^{−1} . (23)

To convert the estimate in (23) to an estimate for the rate of convergence of P to its stationary (uniform) distribution π, we use the following (standard) techniques (e.g., see [S93]).

Define the relative pointwise distance of P^{(t)} to π to be

Δ(t) := max_{x,y} \frac{|P_{y,x}^{(t)} − π(x)|}{π(x)} (24)

where P_{y,x}^{(t)} denotes the probability of being at x after t steps starting at y. It is not hard to show (see [S93]) that

Δ(t) < (1 − λ)^{t} \frac{vol S}{min_{x} deg_{Γ} x}

(25)

≤ e^{−λt} \frac{vol S}{deg Γ} = e^{−λt}|S|

where

vol S : = ∑_{x∈S} deg_{Γ} x,
deg G = deg_{Γ} x for any x

and λ is the eigenvalue of L_{S} which maximizes |1 − λ| for λ ≠ 0.

In order to guarantee that λ is in fact λ_{S}, we can modify P to be “lazy”, i.e., so that the modified walk \overline{P} stays put with probability 1/2, and moves with probability 1/2 times what P did. The eigenvalues for \overline{P} are just 1/2 times those for P, and so, are contained in [0, 1]. Thus, if

t > \frac{2}{λ_{S}} ln \frac{|T |}{ϵ} (26)

then Δ(t) < ϵ. Note that

|T | ≤ min {∏_{i} r_{i}^{n}, ∏_{j} c_{j}^{m}} .

Thus, by (19) and (26), if

(27) t > 6400e^{1/c}m^{2}n^{2} min \{∑_{i} r_{i}^{2}, ∑_{j} c_{j}^{2}\} (ln \frac{1}{ϵ} + min \{n ∑_{i} ln r_{i}, m ∑_{j} ln c_{j}\})

then Δ(t) < ϵ, provided

min \{min_{i} \frac{r_{i}}{n} , min_{j} \frac{c_{j}}{n} \} > c(m − 1)^{3/2}(n − 1)^{3/2} .

As remarked earlier, this shows that the natural Neumann walk on T (\overline{r}, \overline{c}) converges to uniform in time polynomial in the dimensions of the table and the sizes of the line sums (and as the square of the diameter of the space). This strengthens a recent result of Diaconis and Saloff-Coste [DS1] who showed that for fixed dimension the natural walk on T (where any step which might create a negative entry is simply not taken) converges to uniform in time polynomial in the sizes of the line sums and as the square of diameter of the graph T (they use total variation distance instead of relative pointwise distance) but with constants that grow exponentially in mn.

By taking successive relaxations of the line sum constraints (as described in [DG] or [DKM]), it is possible to approximately enumerate T in polynomial time as well.

We also note that Dyer, Kannan and Mount [DKM] have developed a rather different (continuous) random walk which is rapidly mixing on T . They show that the dominant eigenvalue λ for their walk satisfies

λ > \frac{c}{(m − 1)^{4}(n − 1)^{4}(m + n − 2)}

for m, n > 1.

6. Restricted contingency tables

A natural extension of a contingency table is one in which certain entries are restricted, e.g., required to be 0. In this section we indicate how such restricted tables can be dealt with.

Given m, n > 0, let A ⊆ [1, m] × [1, n] be some nonempty index set. By a A-table T we mean an m × n array in which T (i, j) ≡ 0 if (i, j) ∉ A. For given row sum and column sum vectors \overline{r} = (r_{1}, …, r_{m}) and \overline{c} = (c_{1}, …, c_{n}), respectively, with ∑_{i} r_{i} = ∑_{j} c_{j}, we let

T_{A} = T_{A}(\overline{r}, \overline{c}) = \{A-tables T : ∑_{j} T (i, j) = r_{i}, ∑_{i} T (i, j) = c_{j}, 1 ≤ i ≤ m, 1 ≤ j ≤ n\} .

As before, we want to execute a “natural” Neumann walk P on T_{A} and show that it is rapidly mixing. However, several new complications arise from what had to be considered in the preceding section.

To begin with, what will we use for steps in our walk? Let us associate to A a bipartite graph B = B_{A} with vertex sets [1, m] and [1, n], and with (i, j) an edge of B if and only if (i, j) ∈ A. In the case of unrestricted tables, B is the complete bipartite graph on these vertex sets, and the basic steps just occur on the 4-cycles of B. For our more general bipartite graph B, we first normalize it as follows. Let C be a connected component of B. For each cut-edge e = (i, j) of C, the corresponding value T (i, j) is easily seen to be determined by \overline{r} and \overline{c}, say, it is w(e). Then replace r_{i} by r_{i} − w(e) and c_{j} by c_{j} − w(e). Now continue this process recursively until we finally arrive at the 2-connected components C′ of C, with correspondingly reduced row and column sums \overline{r}′ and \overline{c}′. It is not hard to show that there are always feasible assignments satisfying the modified line sum constraints.

Next, we have to describe the basic moves of our Neumann walk on T_{A}. For each 2connected component C of B, let T_{C} denote a fixed spanning tree on the vertex set V (C) = V_{1} ∪ V_{2} of C, where V_{1} ⊂ [1, m], V_{2} ⊂ [1, n]. For each edge i = {i, j} not in T_{C}, with i ∈ V_{1}, j ∈ V_{2}, the addition of e to T_{C} creates some simple even cycle Z(e). We then assign alternating ±1’s on the consecutive edges of Z(e) to generate a possible move ℤ^{±}(e) on the table. It is not hard to show that the set of cycles ⋃_{e} Z(e) form a cycle basis over ℤ for the even cycles on V (C) (in fact, only coefficients of 0 or ±1 are ever needed). It is also not difficult to see that the set of moves K = ⋃_{C} ⋃_{e} Z^{±}(e), C a 2-connected component of B, connects the set T_{A} of all A-tables. Note that K has size O(mn), and so is smaller than what was used in the unrestricted (complete) case.

Of course, if A′ denotes the set of non-cut-edges of B then our tables T ∈ T_{A} only have possibly varying entries T (i, j) where (i, j) ∈ A′. Thus, our A′-tables can now be represented as integer points in 𝔼^{|A′|},. In fact, because of the line sum constraints, the set T_{A′} of all A′-tables actually lies in a subspace M of dimension N = ∑_{C}(m_{C} − 1)(n_{C} − 1) where C ranges over the 2-connected components of B, and m_{C} and n_{C} are the sizes of the vertex sets of C.

The same “minimum average line sum’ technique from the preceding section now applies to each (2-connected) component C of B to produce a “central point” in M_{A′}, the expanded submanifold of M which allows real entries ≥ −2/3 in A′-tables. Of course, S_{A′} = T_{A′} is now the set of all lattice points in A_{A′} (and so, consists of A′-tables with all entries ≥ 0). This shows that there is a ball of radius s_{0} in M_{A′}, where s_{0} is the minimum average line sum occurring in all the components C. Of course, as before, we must restrict s_{0} to be reasonably large. The final result is an estimate for λ_{S} of the form (22) where now the “constant” c_{0} depends on the geometry of the specific generators in K. At worst, c_{0} decreases by a factor of at most O((m^{2} + n^{2})m^{2}n^{2}) from the unrestricted case. We sketch how this comes about in the Appendix.

Another interesting special case we mention here is that of symmetric contingency tables, i.e., with m = n, \overline{r} = \overline{c} and T (i, j) = T (j, i) for all i and j. We can use as basic moves in this case the \binom{n}{2} symmetric transformations

T (i, i) → T (i, i) + 1
T (j, j) → T (j, j) + 1
T (i, j) → T (i, j) − 1

T (j, i) → T (j, i) − 1

and their inverses, for any i ≠ j. Any symmetric table can be transformed to the diagonal table with T (i, i) = r_{i} for all i (and 0 otherwise), so that these moves connect the space S = S(\overline{r}) of symmetric contingency tables for \overline{r}. Easy calculations show that diam M < 2(∑_{i} r_{i})^{1/2} and that M contains a ball of radius \frac{1}{n} min_{i} r_{i}. Thus, by the same arguments that led to (23), we have

λ_{S} > \frac{c_{0}e^{−1/c}}{n^{4} ∑_{i} r_{i}^{2}} (28)

for an absolute constant c_{0}, provided min_{r} r_{i} > c\binom{n}{2}^{3/2}. As usual, the translation of this bound to one for the rate of convergence of the walk to uniform is straightforward.

We remark here that a similar analysis can be applied to the more general problem in which our space of objects on which to walk is now a general graph G with nonnegative integer weights assigned to its edges so that for each vertex v of G, the sum of the weights on the edges incident to v is exactly some preassigned value S(v). Again, the steps of the walk will consist of modifying edge weights by alternating ±1’s on certain simple even cycles (forming a cycle basis over ℤ of all even cycles in G, analogous to what was done for general bipartite graphs). Not surprisingly, the bound on λ_{S} for the corresponding Neumann walk has the same general form as for the bipartite case.

7. Compositions of an integer

An easy application of the preceding ideas concerns compositions of an integer T with a fixed number of parts. These are just ordered partitions (r_{1}, r_{2}, …, r_{n}) with integers r_{i} > 0 so that ∑_{i} r_{i} = T. The basic moves for the walk will be the n(n − 1) transformations of the type: (x_{1}, …, x_{i}, …, x_{j}, …, x_{n}) → (x_{1}, …, x_{i} + 1, …, x_{j} − 1, …, x_{n}). Of course, Γ consists of all points \overline{x} = (x_{1}, …, x_{n}) ∈ 𝔼^{n}, and M = M_{T} = \{\overline{x} ∈ 𝔼^{n} : x_{1} + ⋯ + x_{n} = T \}. We let S = \{\overline{x} ∈ M : x_{i} ≥ −2/3\}. Then diam M ≤ \sqrt{2}(T + \frac{2n}{3}), dim M = n − 1 and M contains a ball of radius T/n. Thus, by (19) if T > cn^{5/2} then

λ_{S} > \frac{c_{0}e^{−1/c}}{n^{2}T^{2}} (29)

for an absolute constant c_{0}. In the Appendix, we show that we can take c_{0} = 1/200. Hence, for T > cn^{5/2}, if t > 400e^{1/c}n^{2}T^{2}(n ln T + ln 1/ϵ) then Δ(t) < ϵ. We point out that compositions can also be treated somewhat more directly by the methods in [CY2]. This is essentially because S in this case is what is called there “strongly convex”, and, in fact, satisfies the stronger condition that if x, y ∈ S, z ∈ ∂S with x ∼ z ∼ y then x ∼ y. As was pointed out in [CY2], we therefore can conclude that

λ_{S} ≥ \frac{1}{8kD^{2}}

where k = |K| and D := graph diameter of S. Thus,

λ_{S} ≥ \frac{1}{8n(n − 1)(T − n)^{2}} > \frac{1}{8n^{2}T^{2}}

for this Neumann walk on compositions (with no restrictions on T). We note that Diaconis/Saloff-Coste [DS1] also treat compositions (by quite different methods) for n and T in restricted ranges. In particular, they conjecture that t = Ω((nT + T^{2})(n ln T + ln 1/ϵ)) steps suffice to guarantee that Δ(t) < ϵ.

8. Knapsacks

In our final example, we will consider the following natural generalization of compositions, which we call “knapsack” solutions (e.g., see [DFKKPV], [DGS]). We are given an integer vector \overline{a} = (a_{1}, …, a_{n}) with a_{i} > 0 and gcd(a_{1}, …, a_{n}) = 1. For an integer T, we consider the set S = S_{T} of all integer vectors \overline{r} = (r_{1}, …, r_{n}), r_{i} > 0, such that

\overline{r} · \overline{a} = ∑_{i} r_{i}a_{i} = T . (30)

As usual, our goal will be to (approximately) select a random element from (and enumerate) the set of “knapsack” vectors in S. We do this by constructing a rapidly mixing Markov chain on S, and estimating the corresponding Neumann eigenvalue λ_{S}. The manifold M is just

\{\overline{x} ∈ 𝔼^{n} : ∑_{i} a_{i}x_{i} = T \}, and M := \{\overline{x} ∈ M : x_{i} > −2/3, 1 ≤ i ≤ n\}. Let e_{i} = (0, …, 1, …, 0)

with 1 in the i^{th} component, and define g_{ij} = a_{j}e_{i} −a_{i}e_{j}. Our generator set will be K = {±g_{ij} : 1 ≤ i ≠ j ≤ n}.

The first problem is to show that S is connected with the generators from K, provided T is sufficiently large. For ease of exposition we are going to make the assumption that gcd(a_{1}, a_{2}) = 1. The general case in which we only assume gcd(a_{1}, a_{2}, …, a_{n}) = 1 is somewhat more technical but offers no real new difficulty. Define A = a_{1} + a_{2} + ⋯ a_{n}, and suppose T > 2A max_{i}{a_{i}}. Then it is easy to see that by repeated application of generators of the form g_{3i}, i ≠ 3, we can reach a vector (r_{1}, r_{2}, …, r_{n}) ∈ S so that r_{3}a_{3} > 2a_{1}a_{2}. Now it is well known (see [EG80]) that any w > 2a_{1}a_{2} can be represented as

w = w_{1}a_{1} + w_{2}a_{2}, w_{1}, w_{2} integers ≥ 0 .

Hence, we can apply g_{13} w_{1} times, and then g_{23} w_{2} times to reach the vector

(r_{1} + w_{1}a_{3}, r_{2} + w_{2}a_{3}, 0, …, r_{n})

which has 0 for its 3rd component. Now, can repeat this argument for each of the other components r_{i}, i ≥ 4, and eventually arrive at a vector (r′_{1}, r′_{2}, 0, 0, …, 0) ∈ S. Finally, by applying g_{12} appropriately we can finally reach (r_{1}^{∗}, r_{2}^{∗}, 0, 0, …, 0) ∈ S where 0 ≤ r_{2}^{∗} < a_{1}. Such a vector is unique, and in reaching it we always remained in S. Thus, S is connected by the generators in K. Also, dim M = n−1 and diam M ≤ T. Since the point (T/A, T/A, …, T/A) ∈ S, then M contains a ball of radius T/A. Hence, by the usual arguments (by now) where R < (2(a_{1}^{2} + ⋯ + a_{n}^{2}))^{1/2} in Claim 2 by using generators g_{i,i+1}, if T > cn^{3/2} (∑_{i} a_{i}^{2})^{1/2} then

λ_{S} > \frac{c_{0}e^{−1/c}}{n^{2}T^{2}} (31)

where c_{0} depends only on the geometry of the generators in K, and not on T. (We estimate c_{0} in the Appendix.) Thus, by (26), if T > cn^{3/2} (∑_{i} a_{i}^{2})^{1/2} then for

t > c_{1}e^{1/c}n^{2}T^{2}(n ln T + ln 1/ϵ) (32)

we have Δ(t) < ϵ (since |S| < T^{n}) where c_{1} depends only on the geometry of the vectors in K.

Appendix

In this section, we provide additional details needed for bounding the constant c_{0} occurring in (19). The arguments will extend (and depend on) those in [CY1]. Briefly, in [CY1] we have

−\frac{2d}{ℓ^{2}}L_{s} = \frac{2d}{|K|} ∑_{g∈K^{∗}} (\frac{μ(x, gx)}{ℓ})^{2} \frac{∂^{2}}{∂g^{2}}
= ∑_{i,j} a_{ij} \frac{∂^{2}}{∂x_{i}∂x_{j}}

where d := dim M, ℓ := min_{g∈K} μ(x, gx), μ denotes (Euclidean) length, and K^{∗} ⊂ K consists of exactly one element from each pair {g, g^{−1}}, g ∈ K.

Suppose C_{1} and C_{2} are constants so that

C_{1}I ≤ (a_{ij}) ≤ C_{2}I

where I is the identity operator on M, and X ≤ Y means that the operator Y − X is positive definite. In particular, we can take for C_{1} and C_{2} the least and greatest eigenvalues, respectively, of (a_{ij}) restricted to M. Now, it follows from the arguments in [CY1] that when M is Euclidean then the constant c_{0} in (19) can be taken to be

c_{0} = \frac{1}{100} min(C_{1}, C_{2}^{−1}) .

Thus, to determine c_{0} in various applications, our job becomes that of bounding the eigenvalues of the corresponding matrix (a_{ij}).

First, we consider m × n contingency tables. With each edge generator g = x_{ij} − x_{i′j} − x_{ij′} − x_{i′j′} we consider \frac{∂^{2}}{∂g^{2}} in terms of the x’s.

Expanding, we have

\frac{∂^{2}}{∂g^{2}} = \frac{∂^{2}}{∂x_{ij^{2}}} + \frac{∂^{2}}{∂x_{i′j}^{2}} + \frac{∂^{2}}{∂x_{ij′}^{2}} + \frac{∂^{2}}{∂x_{i′j′}^{2}} − 2\frac{∂^{2}}{∂x_{ij}∂x_{i′j}}

− 2\frac{∂^{2}}{∂x_{ij}dx_{ij′}} − 2\frac{∂^{2}}{∂x_{i′j′}dx_{i′j}} − 2\frac{∂^{2}}{∂x_{i′j′}dx_{ij′}}

+ 2\frac{∂^{2}}{∂x_{ij}∂x_{i′j′}} + 2\frac{∂^{2}}{∂x_{i′j}dx_{ij′}} .

We can abbreviate this in matrix form as

x_{ij} x_{i′j} x_{ij′} x_{i′j′}
x_{ij} 1 −1 −1 1
x_{i′j} −1 1 1 −1
x_{ij′} −1 1 1 −1
x_{i′j′} 1 −1 −1  1

We need to consider the operator

∑_{g∈K^{∗}} \frac{∂^{2}}{∂g^{2}} .

The corresponding matrix Q has the following coefficient values for its various entries:

\underline{Entry}  \underline{Coefficient}
(x_{ij}, x_{ij})  (m − 1)(n − 1)
(x_{ij}, x_{i′j})  −(m − 1)
(x_{ij}, x_{ij′})  −(m − 1)
(x_{ij}, x_{i′j′})  1

Thus, Q has two distinct eigenvalues: one is mn with multiplicity (m − 1)(n − 1), and the other is 0 with multiplicity m + n − 1. Now, dim M = (m − 1)(n − 1) and the operator corresponding to Q when restricted to M has all eigenvalues equal to mn. So the matrix (a_{ij}) has all eigenvalues equal to 2\frac{mn(m−1)(n−1)}{\binom{m}{2}\binom{n}{2}} = 8, and consequently we can take C_{1} = C_{2} = 8, and c_{0} = 1/800.

Next, in the case of restricted tables, we associate our restrictions with a bipartite graph B (which we assume for now is 2-connected; the general case involves taking the union of such graphs). As usual, let T denote a fixed spanning tree of B, and let C denote the associated (even) cycle basis for B. For each cycle Z in C with edges e_{1}, e_{2}, …, e_{2r} we have the edge generator g = x_{e_{1}} − x_{e_{2}} + ⋯ − x_{e_{2r}}. We consider the matrix (a_{ij}) associated with the operator

\frac{2d}{|K|} ∑_{g∈K^{∗}} (\frac{μ(x, gx)}{ℓ})^{2} \frac{∂^{2}}{∂g^{2}} = ∑_{i,j} a_{ij} \frac{∂^{2}}{∂x_{i}∂x_{j}} .

Clearly, for f : E(B) → ℝ, we can express (a_{ij}) as a quadratic form:

⟨f, (a_{ij})f⟩ = ∑_{Z∈C}(f(e_{1}) − f(e_{2}) + ⋯ − f(e_{2r}))^{2} · r

since \frac{2d∥Z∥^{2}}{|K|ℓ^{2}} = r.

To upper-bound the eigenvalues of (a_{ij}) we have

C_{2} ≤ \frac{2d}{|K|}(m + n)^{2}mn ≤ 2(m + n)^{2}mn

since

\frac{⟨f, (a_{ij})f⟩}{⟨f, f⟩} ≤ ∑_{Z∈C} \frac{(f(e_{1}) − f(e_{2}) + ⋯ − f(e_{2r}))^{2}r}{∑_{i} f^{2}(e_{i})}
≤ \frac{∑_{Z∈C}(f^{2}(e_{1}) + ⋯ + f^{2}(e_{2r}))r^{2}}{∑_{i} f^{2}(e_{i})}
≤ r^{2}|C| ≤ (m + n)^{2}mn .

To establish a bound for C_{1}, more work is needed. We will use the following modified discrete version of Cheeger’s theorem (see [C] ):

For a graph G, suppose f is the eigenfunction associated with some nonzero eigenvalue λ of the Laplacian of G. Then λ satisfies

λ = ∑_{u∼v} \frac{(f(u) − f(v))^{2}}{∑_{v} f^{2}(v)d_{v}} ≥ \frac{h}{2}λ^{2}

where h_{λ} = min_{v} h(v) and

h(v) := \frac{|\{\{u, w\} ∈ E(G) : f(u) ≤ f(v) < f(w)\}|}{min \{ ∑_{{u \atop f(u)≤f(v)}} d_{u}, ∑_{{w \atop f(v)<f(w)}} d_{w}\}} .

In particular, for λ ≠ 0,

h_{λ} ≥ \frac{1}{E(G)} .

Before applying the above result, we will modify A = (a_{ij}). First, choose a root ρ in T. For each tree edge e = {u, v} with d_{T} (ρ, u) < d_{T} (ρ, v) and d_{T} (ρ, u) odd (where d_{T} denotes the graph distance in T), we define

e′ = e − e_{1} + e_{2} − ⋯

where the unique path from u to ρ consists of the edges e_{1}, e_{2}, …. Let X denote the matrix corresponding to this change of coordinates. Thus,

\frac{⟨f, X^{tr}(a_{ij})Xf⟩}{⟨fX^{tr}, Xf⟩} · \frac{⟨fX^{tr}, Xf⟩}{⟨f, f⟩} = \frac{⟨f, A′f⟩}{⟨f, f⟩} (33)

where A′ = X^{tr}(a_{ij})X corresponds to the quadratic form

⟨f, (a′_{ij})f⟩ = ∑_{Z∈C}(f(e′_{1}) − f(e′_{2}) + f(e) − f(e_{1}))^{2}r .

Note that there are just four terms in each of the squared terms in the sum. By (33), the eigenvalues of (a′_{ij}) are products of eigenvalues of (a_{ij}) and eigenvalues of X^{tr}X. So, to lowerbound the nonzero eigenvalues of A′, we apply the above result twice, since the 4-term sum can be interpreted as two 2-term sums. Thus, we have a lower bound in this case of

\frac{1}{(m + n)^{2}} · \frac{1}{(2mn)^{2}} .

This implies that we can take

c_{0} = \frac{1}{400(m + n)^{2}m^{2}n^{2}}

for restricted tables.

For compositions, our generators are of the form

g_{ij} = (0, …, 1, …, −1, …, 0) = x_{i} − x_{j} .

Then

\frac{2d}{|K|} ∑_{\{i,j\}} \frac{∂^{2}}{∂g_{ij}^{2}} = ∑_{ij} a_{ij} \frac{∂}{∂x_{i}} \frac{∂}{∂x_{j}}

where

a_{ij} = \{ \begin{matrix}\frac{2(n−1)}{n} & if i = j \\ −\frac{2}{n} & if i ≠ j \end{matrix} .

Thus, (a_{ij}) has eigenvalues 0 of multiplicity one, and 2 of multiplicity n − 1. This implies that C_{1} and C_{2} are both equal to 2 so that we can take c_{0} = 1/200.

Finally, for the knapsack problem, we have edge generators (…, a_{j}, …, −a_{i}, …). These correspond to g_{ij} = a_{j}x_{i} − a_{i}x_{j}. Therefore

∑_{\{i,j\}} (\frac{a_{i}^{2} + a_{j}^{2}}{2} ) \frac{∂^{2}}{∂g_{ij}^{2}} = ∑_{\{i,j\}} (a_{i}^{2} \frac{∂^{2}}{∂x_{i}^{2}} − 2a_{i}a_{j} \frac{∂}{∂x_{i}} \frac{∂}{∂x_{j}} + a_{j}^{2} \frac{∂}{∂x_{j}^{2}} )
= \frac{|K|}{2d} ∑_{i,j} a_{ij} \frac{∂}{∂x_{i}} \frac{∂}{∂x_{j}}

where the sum is taken over all unordered pairs i, j. The matrix (a_{ij}) corresponds to the quadratic form

\frac{2d}{|K|ℓ^{2}} ∑_{i,j} (\frac{a_{i}^{2} + a_{j}^{2}}{2} ) (a_{j}x_{i} − a_{i}x_{j})^{2} = ⟨\overline{x}, A\overline{x}⟩

where d = n − 1, |K| = n(n − 1), and ℓ = min_{i,j}(a_{i}^{2} + a_{j}^{2})^{1/2} is the minimum edge length. Set

β_{i} = ∑_{j} \frac{1}{2}(a_{i}^{2} + a_{j}^{2})a_{i}^{2}a_{j}^{2}. Then

min_{i,j}(a_{i}^{2} + a_{j}^{2}) · \frac{n}{2} \frac{⟨\overline{x}, A\overline{x}⟩}{⟨\overline{x}, \overline{x}⟩} = \frac{∑_{\{i,j\}} \frac{1}{2}(a_{i}^{2} + a_{j}^{2})(a_{j}x_{i} − a_{i}x_{j})^{2}}{∑_{i} x_{i}^{2}}
= \frac{∑_{\{i,j\}} \frac{1}{2}(a_{i}^{2} + a_{j}^{2})a_{i}^{2}a_{j}^{2}(y_{i} − y_{j})^{2}}{∑_{i} a_{i}^{2}y_{i}^{2}} where y_{i} = \frac{x_{i}}{a_{i}}
≥ \frac{min_{i} a_{i}^{2} ∑_{\{i,j\}} a_{i}^{2}a_{j}^{2}(y_{i} − y_{j})^{2}}{∑_{i} a_{i}^{2}y_{i}^{2}}
≥ min_{i} a_{i}^{2} · min_{i} ∑_{{j \atop j≠i}} a_{j}^{2} \frac{∑_{\{i,j\}} a_{i}^{2}a_{j}^{2}(y_{i} − y_{j})^{2}}{∑_{i} a_{i}^{2}y_{i}^{2} ∑_{{k \atop k≠i}} a_{k}^{2}} .

Now, by the modified Cheeger theorem referred to previously, we have for any eigenvalue λ ≠ 0 of A

min_{i,j}(a_{i}^{2} + a_{j}^{2})\frac{n}{2}λ ≥ \frac{1}{2}h^{2} (∑_{j} a_{j}^{2} − max_{i} a_{i}^{2}) min_{i} a_{i}^{2}

where h is the Cheeger constant for the complete graph K_{n} with edge weights a_{i}^{2}a_{j}^{2}, which is defined by

h = inf_{I⊂V} \frac{∑_{i∈I} a_{i}^{2} ∑_{j∉I} a_{j}^{2}}{∑_{i∈I} a_{i}^{2} ∑_{j≠i} a_{j}^{2}}

taken over all I ⊂ V = V (K_{n}) satisfying

∑_{i∈I} a_{i}^{2} (∑_{k≠i} a_{k}^{2}) ≤ \frac{1}{2} ∑_{i} ∑_{{k \atop k≠i}} a_{i}^{2}a_{k}^{2}

i.e.,

∑_{j∉I} a_{j}^{2} ∑_{{k \atop k≠j}} a_{k}^{2} ≥ ∑_{i∈I} a_{i}^{2} ∑_{{k \atop k≠i}} a_{k}^{2} .

First, suppose ∑_{i∈I} a_{i}^{2} ≤ ∑_{j∉I} a_{j}^{2}. Then

h ≥ \frac{∑_{i∈I} a_{i}^{2} ∑_{j∉I} a_{j}^{2}}{∑_{i∈I} a_{i}^{2} ∑_{{k \atop k≠i}} a_{k}^{2}} ≥ 1/2 .

On the other hand, suppose ∑_{j∉I} a_{j}^{2} < ∑_{i∈I} a_{i}^{2}. Then

h ≥ \frac{∑_{i∈I} a_{i}^{2} ∑_{j∉I} a_{j}^{2}}{∑_{j∉I} a_{j}^{2} ∑_{{k \atop k≠j}} a_{i}^{2}}

≥ \frac{∑_{i∈I} a_{i}^{2} ∑_{j∉I} a_{j}^{2}}{∑_{j∉I} a_{j}^{2} ∑_{k} a_{k}^{2}} ≥ \frac{1}{2} .

Thus, in both cases we get h ≥ 1/2 so that

\frac{⟨\overline{x}, A\overline{x}⟩}{⟨\overline{x}, \overline{x}⟩} ≥ (\frac{2}{n}) · (\frac{1}{2})^{2} \frac{min_{i} a_{i}^{2} min_{i} ∑_{j≠i} a_{j}^{2}}{min_{i,j}(a_{i}^{2} + a_{j}^{2})} = C_{1}

For C_{2}, we compute

\frac{⟨\overline{x}, A\overline{x}⟩}{⟨\overline{x}, \overline{x}⟩} ≤ \frac{2}{n min_{i,j}(a_{i}^{2} + a_{j}^{2})} \frac{∑_{\{i,j\}} \frac{1}{2}(a_{i}^{2} + a_{j}^{2})(a_{j}x_{i} − a_{i}x_{j})^{2}}{∑_{i} x_{i}^{2}}

≤ \frac{2}{n min_{i,j}(a_{i}^{2} + a_{j}^{2})} \frac{∑_{\{i,j\}}(a_{i}^{2} + a_{j}^{2})(a_{j}^{2}x_{i}^{2} + a_{i}^{2}x_{j}^{2})}{∑_{i} x_{i}^{2}}

≤ \frac{2}{n min_{i,j}(a_{i}^{2} + a_{j}^{2})} max_{i} ∑_{{j \atop j≠i}} a_{j}^{2}(a_{i}^{2} + a_{j}^{2}) = C_{2}

and, as usual, we take c_{0} = \frac{1}{100} min(C_{1}, C_{2}^{−1}).

References

[A87] D. Aldous, On the Markov-chain simulation method for uniform combinatorial simulation and simulated annealing, Prob. Eng. Info. Sci. 1 (1987), 33-46.

[BM76] J. A. Bondy and U. S. R. Murty, Graph Theory with Applications, American Elsevier, New York, 1976.

[C] F. R. K. Chung, Laplacians of graphs and Cheeger’s inequalities (to appear).

[CY1] F. R. K. Chung and S.-T. Yau, A Harnack inequality for homogeneous graphs and subgraphs, Comm. Analysis and Geometry, 2 (1994) 628–639.

[CY2] F. R. K. Chung and S.-T. Yau, Heat kernel estimates and eigenvalue inequalities for convex subgraphs (preprint).

[DG] P. Diaconis and A. Gangolli, Rectangular arrays with fixed margins (preprint).

[DGS] P. Diaconis, R. L. Graham and B. Sturmfels, Primitive partition identities (to appear).

[DS1] P. Diaconis and L. Saloff-Coste, Random walk on contingency tables with fixed row and column sums (preprint).

[DSt] P. Diaconis and B. Sturmfels, Algebraic algorithms for sampling from conditional distributions (to appear in Ann. of Statistics).

[DFK89] M. Dyer, A. Freeze and R. Kannan, A random polynomial time algorithm for approximating the volume of convex bodies, in Proc. 21st ACM Symp. on Theory of Computing, 1989, pp. 375–381.

[DFKKPV] M. Dyer, A. Freeze, R. Kannan, A. Kapour, L. Perkovic, U. Vazirani, A mildly exponential time algorithm for approximating the number of solutions to a multidimensional knapsack problem, (to appear in Combinatorics, Probability and Computing).

[DKM] M. Dyer, R. Kannan and J. Mount, Sampling contingency tables (preprint).

[EG80] P. Erdős and R. L. Graham, Old and new problems and results in combinatorial number theory, L’Enseignement Mathematique, Monographe N^{∘}28, Genève, 1980, p. 86.

[FKP94] A. Frieze, R. Kannan and N. Polson, Sampling from log-concave distributions, Annals Appl. Prob. 4 (1994), 812–837.

[G91] A. Gangolli, Convergence bounds for Markov chains and applications to sampling, Ph.D. Dissertation, Dept. of Comp. Sci., Stanford, 1991.

[JS] M. Jerrum and S. Sinclair, Fast uniform generation of regular graphs, Th. Comp. Sci. 73 (1990), 91–100.

[LS90] L. Lovàsz and M. Simonovits, The mixing rate of Markov chains, an isoperimetric inequality and computing the volume, in Proc. 31st IEEE Symp. on Foundations of Computer Science, 1990, pp. 346–354.

[MW] B. McKay and N. Wormald, Uniform generation of random regular graphs of moderate degree (to appear in J. of Algorithms).

[SY94] R. Schoen and S-T. Yau, Differential Geometry, International Press, Cambridge, 1994.

[S93] A. Sinclair, Algorithms for random generation and counting: a Markov chain approach, Birkhauser, Boston, 1993.

[SJ89] A. Sinclair and M. Jerrum, Approximate counting, uniform generation and rapidly mixing Markov chains, Infor. and Comput. 82 (1989), 93–133.