A New Approach to Polynomial-Time Generation of Random Points in Convex Bodies^{∗}

Russ Bubley Martin Dyer School of Computer Studies University of Leeds, Leeds LS2 9JT, England

Mark Jerrum Department of Computer Science University of Edinburgh, Edinburgh EH9 3JZ, Scotland

March 28, 1996

Abstract

In this paper we describe a new method for proving the polynomial-time convergence of an algorithm for sampling (almost) uniformly at random from a convex body in high dimension. Previous approaches have been based on estimating conductance via isoperimetric inequalities. We show that a conceptually simpler coupling argument can be used to give a similar result.

1 Introduction

Dyer, Frieze and Kannan [7] gave the first fully polynomial randomized approximation scheme for approximating the volume of a convex body in ℝ^{n}, for large n. The algorithm was based on a reduction to sampling uniformly at random in such a body. The necessary sampling was achieved, to close enough approximation, by simulating a polynomial number of steps of a random walk in the body.

The approach of [7] to proving “rapid mixing” of the random walk was to estimate the conductance [22]. (This methodology was previously used in [13] to prove rapid convergence of an edge-matching random walk proposed in [4].) Conductance estimation requires some form of isoperimetric inequality. In [7] the appropriate inequality was geometric in origin.

∗Research supported by EXPRIT Working Group RAND

The result of [7] has been dramatically improved in a sequence of papers. For example, see [2, 6, 18, 8, 19, 16, 15] and [14] for a recent survey. These papers have introduced ideas which are important in their own right, for example new isoperimetric inequalities. We exploit some of these ideas here, specifically the idea of using a surrogate log-concave distribution, and the use of a “Metropolised” rotationally symmetric random walk. These techniques were introduced in [2, 19]. However, while these advances have been made, all approaches to the problem have essentially similarities to [7].

Here we explore a different approach, coupling, which previously seemed difficult to apply to this problem. We show that, by viewing the problem from the correct perspective, it is in fact rather easy to show polynomial convergence to the uniform distribution of a certain random walk. To achieve this, we first transform the uniform sampling problem to that of generating a random point from the log-concave distribution of [19]. Then we apply our coupling arguments to a suitably chosen random walk for the transformed problem.

The most important application of the sampling algorithm is to the above-mentioned volume computation problem. We do not consider this here, but refer the interested reader to the papers already cited. The best time bound currently known is given in [15].

Our analysis contains many trade-offs between constants. We do not attempt to make these choices “optimally”, and no doubt our estimates can be improved upon. Rather, for ease of exposition, we employ a single constant c which must be chosen large enough for various inequalities to hold.

Similarly, we do not address issues of precision in computation. While there is no essential difficulty here, to do so introduces a level of tedious detail which we prefer to avoid. We generally assume that computations can be carried out to the required precision, and measure the complexity of the algorithm simply by the number of steps of the random walk.

The plan of the paper is as follows. In section 2 we give the necessary backgound, and some initial estimates. In section 3 we analyse the random walk using a coupling argument. In section 4 we describe some improvements, and finally in section 5 we discuss the relationship of our work to the literature cited above.

2 Background and preliminaries

2.1 Convex bodies and gauge functions

Let K ⊆ ℝ^{n} be a convex body with 0 ∈ int K and let ∂K be its boundary. We will assume that K is well-rounded, i.e. for some “small” R, B ⊆ K ⊆ RB where B is the unit ball in ℝ^{n}. A polynomial-time computation will yield R = O(n^{3/2}) [10] or, for almostuniform generation purposes R = O(n) [19]. By a recent result of Kannan, Lovász and Simonovits [16], we may even achieve R = O(n^{1/2}). However, the algorithm producing R = O(n^{1/2}) requires generating almost uniform points in K, so here we assume only R = O(n).

We define the gauge function of K by f_{K}(x) = inf_{t≥0}{x ∈ tK}, and an associated density function^{1}

F_{K}(x) = e^{−f_{K}(x)} (x ∈ ℝ^{n}).

It is well known that f_{K} is convex and hence F_{K} is log-concave. Note that K ⊆ L implies F_{K} ≤ F_{L} pointwise. Also f_{B}(x) = ∥x∥, so if K is well rounded, we have

\frac{∥x∥}{R} ≤ f_{K}(x) ≤ ∥x∥.  (1)

Note that, since f_{K} is convex, it possesses at least one subgradient ∇(x) at each point x. Then, for all y ∈ ℝ^{n} and any subgradient ∇ at x,

f_{K}(y) ≥ f_{K}(x) + ∇.(y − x).  (2)

It is also easy to show that, for any v ∈ ∂B, the directional derivative δ_{v}(x), of f_{K} at x in direction v, satisfies

|δ_{v}(x)| ≤ 1.  (3)

2.2 Sampling equivalence

We will show that sampling uniformly from K is equivalent to sampling with density F_{K}. Specifically, we show that there is a measure-preserving bijection between ℝ^{n} and int K.

Let Γ_{n} the distribution with probability density e^{−t}t^{n−1}/(n − 1)!. We will assume this to be defined for non-integer as well as integer n by using the obvious integral for (n − 1)! in such cases. Let γ_{n}(t) = ∫_{0}^{t} e^{−u}u^{n−1} du/(n − 1)! be its distribution function. Let P_{n} be the distribution on [0, 1] with probability density nt^{n−1} and distribution function t^{n}.

Let us write w(x) = ∥x∥ and b(x) = x/∥x∥. A point x ∈ ℝ^{n} \ {0} can be represented uniquely by the polar coordinates (w, b), where b ∈ ∂B and w ∈ (0, ∞). Let ξ be the uniform probability measure^{2} on ∂B. Then the volume element in ℝ^{n} is proportional to w^{n−1}dwdξ.

Let h(b) = 1/f_{K}(b) be the distance from the origin to ∂K in the direction b. Note that f_{K}(x) = w/h. Then the uniform probability measure μ_{1} on K has element

dμ_{1} = C_{1}w^{n−1}dwdξ (b ∈ ∂B, w ∈ (0, h]),

where clearly C_{1} = 1/ Vol(K). Similarly, the measure μ_{2} determined by F_{K} has element

dμ_{2} = C_{2}e^{−w/h}w^{n−1}dwdξ (b ∈ ∂B, w ∈ (0, ∞)).

Integrating with respect to w, the marginal measures ξ_{1}, ξ_{2} on ∂B induced by μ_{1} and μ_{2} have elements (C_{1}/n)h^{n}dξ and C_{2}(n − 1)!h^{n}dξ. Comparing, we see that ξ_{1} = ξ_{2} = ξ′ (say), and hence in particular C_{2} = 1/n! Vol(K) (c.f. [19]).

1We will use the term density when the integral is not necessarily 1. If we wish to insist on integral 1, we will write probability density.

2Measure will always mean probability measure.

It follows further that, given a point from μ_{2} we can generate a point from μ_{1}. (The reverse is also true, but we will not elaborate the details.) Note, from the expressions for dμ_{1} and dμ_{2}, that the measures ν_{1}, ν_{2} of f(x) = w/h are P_{n} and Γ_{n} respectively. If z is a random point from ν_{2}, then γ_{n}(z)^{1/n} is a point from ν_{1}.

Thus, if we have a point x_{2} from μ_{2}, we can construct a point x_{1} from μ_{1} as follows. If x_{2} = 0 then x_{1} = 0. Otherwise, let b = x_{2}/∥x_{2}∥ ∈ ∂B be drawn from ξ′. Also, in obvious notation, w_{1} = hγ_{n}(w_{2}/h)^{1/n}. Then x_{1} = w_{1}b. Let us write x_{1} = χ(x_{2}). Then clearly χ is a bijection between ℝ^{n} and int K.

In practice, we can only sample approximately from F_{K}. We will take variation distance as the measure of distributional approximation. Let μ_{1}, μ_{2} be two measures defined on the Borel sets in ℝ^{n}. (All measures we employ satisfy this assumption, and henceforward we presume the necessary measurability qualification on sets.) Then their variation distance is defined by

d_{TV}(μ_{1}, μ_{2}) = max_{A}{μ_{1}(A) − μ_{2}(A)},

where A ⊆ ℝ^{n} is a Borel set. See, for example, Appendix A.1 of [3] for more information.

Variation distance is a metric on measures, but we will extend the notation to random variables. If X, Y are random variables, we will write d_{TV}(X, Y) for the variation distance between the measures P(X ∈ ·), P(Y ∈ ·) associated with X and Y,

Thus, let us denote the approximating measure to μ_{2} by μ̂_{2}, and let μ̂_{1} be the resulting approximation to μ_{1}. Then, since χ is a bijection, we have from the definition of total variation distance that

d_{TV}(μ_{1}, μ̂_{1}) = d_{TV}(μ_{2}, μ̂_{2}).

Thus the transformation preserves the quality of the approximation.

Note that, in view of this transformation, we can regard a random walk on ℝ^{n} converging to the F_{K} distribution as a random walk on K converging to the uniform distribution. Thus our development below could be described equivalently in terms of K rather than F_{K}. For example, we define below a distance metric on ℝ^{n} to gauge the convergence of the random walk for F_{K}. The bijection χ then induces a metric on K. We could equally well work with this as a measure for convergence of the random walk in K.

In practice, it is not necessary to determine w_{1} = hγ_{n}(w_{2}/h)^{1/n}. We can merely generate an independent uniform [0, 1] variable U and let w_{1} = U^{1/n}. The additional generator ν̂_{1}, approximating ν_{1}, introduces further error. We have to check that the effect of this on the overall approximation is small. Letting μ̃_{1} be the resulting approximation to μ_{1}, this follows from

Lemma 1

d_{TV}(μ_{1}, μ̃_{1}) ≤ d_{TV}(μ_{2}, μ̂_{2}) + d_{TV}(ν_{1}, ν̂_{1}).

Proof: Let ξ̂′ be the measure resulting from projecting μ̂_{2} onto ∂B. Then it is clear from the definition of variation distance that

d_{TV}(ξ′, ξ̂′) ≤ d_{TV}(μ_{2}, μ̂_{2}).

Let η = d_{TV}(ν_{1}, ν̂_{1}). Now, for any A ⊆ ℝ^{n} and given b ∈ ∂B, let A_{b} = A ∩ (0, h)b. We have μ_{1}(A) = ∫_{∂B} ν_{1}(A_{b})dξ′ and

μ̂_{1}(A) = ∫_{∂B} ν̂_{1}(A_{b})dξ̂′ ≤ ∫_{∂B}(ν_{1}(A_{b}) + η)dξ̂′.

Hence

μ̂_{1}(A) − μ_{1}(A) ≤ ∫_{∂B} ν_{1}(A_{b})(dξ̂′ − dξ′) + η ≤ ∫_{B′}(dξ̂′ − dξ′) + η,

where B′ = {b ∈ ∂B : dξ̂′ ≥ dξ′}. But

∫_{B′}(dξ̂′ − dξ′) = ξ̂′(B′) − ξ′(B′) ≤ d_{TV}(ξ′, ξ̂′),

from which the Lemma follows. □

Thus, if we can generate approximately from F_{K}, we can sample almost uniformly from K.

2.3 Modifying F_{K}

We need to be able to sample from the log-concave distribution defined by F_{K}. It will be easier to show that we can do this for a slightly modified function, which differs from F_{K} only at “large” distances from the origin. We will first show that effect on the sampling distribution will be negligible.

Let ρ = cnR, for some constant c ≥ 2. For convenience, we will assume without loss that R ≥ c^{4}. Then, if X has our target distribution, we have using (1),

P(∥X∥ ≥ ρ) = P(∥X∥/R ≥ cn) ≤ P(f_{K}(X) ≥ cn)

But we know from section 2.2 that f_{K}(X) has the Γ_{n} density. Therefore, using Stirling’s approximation and Lemma 2 of section 2.7 below, we have

P(∥X∥ ≥ ρ) ≤ 2e^{−cn}(cn)^{n−1}/(n − 1)! < (c/e^{c−1})^{n}.  (4)

Thus, by suitable choice of c, ∥X∥ < ρ with very high probability.

In view of this, we may replace f_{K}(x) outside ρB by a larger function without significantly altering the distribution determined by f_{K} for any polynomial-time sampling scheme. We will choose

f′_{K}(x) = max{f_{K}(x), 2(∥x∥ − ρ)}.

Note that this is convex. If ∥x∥ ≤ ρ, then f′_{K}(x) = f_{K}(x), and if ∥x∥ ≥ 2ρ, we have f′_{K}(x) = 2(∥x∥ − ρ) using (1) . For our purposes below, the advantage of this function is that it behaves “essentially” like ∥x∥ at large distances from the origin. Note that f′_{K} has (to within the required approximations) the same properties as f_{K}, except that (3) must be modified to |δ_{v}(x)| ≤ 2.

Henceforward we will simply write f for f′_{K} and F = e^{−f}, assuming only a condition of the form

|δ_{v}(x)| ≤ λ,  (5)

for some given constant λ ≥ 1.

2.4 Metropolis random walks

The following is sufficient for our purpose, but can be generalised. (See, for example, [14, 19].) Consider a lazy (discrete time, continuous space) Markov random walk X_{t} on ℝ^{n}. Denote the (nonnegative) transition density by p(x, x′) where, for each x ∈ ℝ^{n}, we suppose p(x, x′) > 0 in some neighbourhood of x and ∫_{ℝ^{n}} p(x, x′) dx′ ≤ 1. Then we assume P(X_{t+1} ∈ x′ + dx′ | X_{t} = x) = p(x, x′)dx′ for x ≠ x′, and P(X_{t+1} = x | X_{t} = x) = 1 − ∫_{ℝ^{n}} p(x, x′) dx′. (Here dx′ is an infinitesimal neighbourhood of x′, and we abuse notation by employing the same symbol for its volume.)

The random walk is time-reversible if there exists a probability density π(x) on ℝ^{n} such that the conditions

π(x)p(x, x′) = π(x′)p(x′, x) (for all x, x′ ∈ ℝ^{n} with x ≠ x′),

are satisfied. Then it follows that π(x) is the asymptotic density of X_{t}. These conditions are called detailed balance.

The random walks in which we are interested arise in the following way. We have a function q(x, x′) : ℝ^{n} × ℝ^{n} → ℝ, symmetric in x and x′. We further assume q(x, ·) is a probability density for each x, positive in some neighbourhood of x. Suppose F : ℝ^{n} → ℝ is a log-concave function, integrable on ℝ^{n}, such that the derivative of f = − ln F is bounded in absolute value by λ in all directions at every point. Let

M(x) = e^{−λ∥x∥}, and A(x, x′) = \sqrt{M(x′ − x)F (x′)/F (x)}.

Clearly 0 ≤ A(x, x′) ≤ 1 everywhere, since |f(x′)−f(x)| ≤ λ∥x′−x∥ is a consequence of the derivative condition. Thus A(x, x′) is a probability. Also A(x, x′)/A(x′, x) = F (x′)/F (x), since M(x′ − x) = M(x − x′). Now define a lazy random walk as above by p(x, x′) = A(x, x′)q(x, x′). Clearly this satisfies

F (x)p(x, x′) = F (x′)p(x′, x) (x, x′ ∈ ℝ^{n}),

and hence the asymptotic density of the random walk is π(x) = F (x)/ ∫ F (x). This is a modification of the usual Metropolis random walk (see [14]), but is convenient for our proofs. This random walk is easily implemented as follows. If X_{t} = x, we choose a trial point x′ with density q(x, ·). We accept this point with probability A(x, x′). If we accept the point, X_{t+1} = x′, otherwise X_{t+1} = x. We will let A(x, x′) be the 0-1 indicator variable of the acceptance event, so P(A(x, x′) = 1) = A(x, x′). We will denote the complementary quantities (1 − A), (1 − A) by Ā, Ā respectively.

2.5 The random walk

Given a log-concave function of the type discussed at the end of section 2.4, we define a random walk as in section 2.4. We choose q(x, x′) to be Gaussian. Let ϕ be the unit Normal density on ℝ, and G_{n} be the standard Normal density on ℝ^{n}, i.e.

ϕ(x) = e^{−x^{2}/2}/\sqrt{2π} and G_{n}(x) = ∏_{i=1}^{n} ϕ(x_{i}).

Let U be a random variable with density G_{n}. We take q(x, x′) so that x′ = x + σU, for some constant σ, to be determined. Then the trial step is to x′, and it is accepted with probability A(x, x′).

In fact, for convenience, we will modify this random walk slightly. Let r = σ\sqrt{n}. Then, for large n, r is the approximate step size for X_{t}. Now, with c as in (4), let us redefine the acceptance multiplier M(x) to be

M(x) = exp (−λ max {cr, ∥x∥}) .

Now P(∥x′ − x∥ > cr) = P(∥U∥^{2} > c^{2}n). But it is well known that \frac{1}{2}∥U∥^{2} has the Γ_{n/2} density. Hence, using Lemma 2 and Stirling’s approximation,

P(∥x′ − x∥ > cr) ≤ (c^{2}/e^{c^{2}−1})^{n/2} ≤ (c/e^{c−1})^{n}.

Thus, with a suitable choice of c, we may assume that with very high probability there is no step of size greater than cr during any polynomial number of steps of the random walk. Thus M(x) = exp(−λcr), where henceforth we will assume without loss that λr ≤ c^{−4}. In section 4.2 below, we will describe a further modification to M(x) which asymptotically improves the number of steps required, but for ease of exposition here we use the expression above.

2.6 Coupling

The proof technique we use below is coupling, and has been employed in similar contexts by Jerrum [12] and Matthews [20]. An overview of the use of couplings for Markov chains can be found, for example, in [1] and [5]. Details of the extension of the coupling method to more general processes can be found in [9], and [23].

To analyse the convergence rate of our random walk X_{t} to its stationary distribution, we consider a second random walk, Y_{t}. The “Coupling Lemma” (which is very easy to prove) then states that

d_{TV}(X_{t}, Y_{t}) ≤ P(X_{t} ≠ Y_{t}).

This Lemma is usually attributed to Aldous and holds for any (not necessarily Markovian) processes X_{t}, Y_{t}. In order to apply it to convergence of Markov processes, we will require that marginally, X_{t} and Y_{t} are both faithful copies of the process and that Y_{0} is chosen randomly from the stationary distribution (and hence Y_{t} has this property throughout).

The Coupling Lemma will give us an upper bound on the variation distance of the random walk from its stationary distribution. To get good bounds, however, we have to construct the joint process (X_{t}, Y_{t}) so that we quickly get X_{t} = Y_{t}. In order to achieve this here, we use two different couplings depending on ∥X_{t} − Y_{t}∥. The first encourages the two random walks to come closer to each other when they are distant, and the second encourages the random walks to meet if they are sufficiently close.

We wish to have d_{TV}(X_{t}, Y_{t}) ≤ ε for arbitrary ε. However, it is sufficient to have τ, of size polynomial in n, such that P(X_{τ} ≠ Y_{τ}) ≤ \frac{1}{2}. Then we may regard t as divided into blocks of length τ, each constituting a coupling “trial” with “success” probability at least \frac{1}{2}. The blocks are independent. and hence if k = ⌈lg(1/ε)⌉,

d_{TV}(X_{kτ}, Y_{kτ}) ≤ P(X_{kτ} ≠ Y_{kτ}) ≤ (\frac{1}{2})^{k} ≤ ε.

Thus we can guarantee variation distance ε by increasing the execution time only by a factor O(log(1/ε)). Therefore, we aim below simply to attain a probability \frac{1}{2} of coupling.

2.7 Technical Results

We collect here some simple lemmas which are used in the analysis of section 3. All but the most tenacious reader may wish to skip this material, at least on first reading. We first give the proof of a result used in section 2.3.

Lemma 2 Let Z have the Γ_{k+1} density, and let μ ≥ 2k, then

e^{−μ}μ^{k}/k! ≤ P(Z ≥ μ) ≤ 2e^{−μ}μ^{k}/k!

Proof: The left hand inequality is easily proved by proved by parts integration. Also, since μ ≥ 2k,

P(Z ≥ μ) = \frac{1}{k!} ∫_{μ}^{∞} e^{−t}t^{k}dt
≤ \frac{1}{k!}e^{−μ/2}μ^{k} ∫_{μ}^{∞} e^{−t/2}dt
= 2e^{−μ}μ^{k}/k!

□

We now state a very simple bound on ln n, which we use repeatedly without comment. It may be proved by elementary calculus.

Lemma 3 If α > 0, then ln n/n^{α} ≤ 1/(eα).

We next prove two simple approximation results.

Lemma 4 Let a ≥ \frac{1}{4}, and b > 0. Then \sqrt{a + b} ≤ \sqrt{a} + b.

Proof: This follows from a + b ≤ a + 2\sqrt{a} b ≤ (\sqrt{a} + b)^{2}, on taking square roots. We require only 2\sqrt{a} ≥ 1, i.e. a ≥ \frac{1}{4}. □

Corollary 1 Let a ≥ \frac{1}{4}, and b > 0. Then ^{4} \sqrt{a + b} ≤ ^{4} \sqrt{a} + b.

Proof: Applying Lemma 4 twice

4 \sqrt{a + b} ≤ \sqrt{\sqrt{a} + b} ≤ ^{4} \sqrt{a} + b,

where the first application is valid for a ≥ \frac{1}{4}, and the second for a ≥ \frac{1}{16}. □

Lemma 5 If |z| ≤ \frac{1}{2}, then \sqrt{1 − z} ≤ 1 − \frac{1}{2}z − \frac{1}{10}z^{2}.

Proof: Since 1 − \frac{1}{2}z − \frac{1}{10}z^{2} is positive for |z| ≤ 1, we may obtain an equivalent inequality by squaring both sides. Simplifying this gives z^{2} + 10z + 5 ≥ 0. This is satisfied for all |z| ≤ \frac{1}{2}. □

We now prove a simple bound on the tail of the Normal distribution.

Lemma 6 Let v ∼ N(0, σ^{2}) and ℓ = cσ ln n, then P(|v| > ℓ) < n^{−c^{2} ln n/2}.

Proof:

P(|v| > ℓ) = \frac{2}{σ} ∫_{ℓ}^{∞} ϕ(t/σ) dt ≤ \frac{2}{σℓ} ∫_{ℓ}^{∞} tϕ(t/σ) dt
= \frac{2σ}{ℓ} ϕ(ℓ/σ) < n^{−c^{2} ln n/2}.

□ Next we prove two easy results to deal with mild conditioning in the analysis.

Lemma 7 Let Z be a real valued random variable, symmetrically distributed about the origin. For t ≥ 0, let E_{t} be the event |Z| ≤ t. Then if f(z) is an

(i) odd function, E (f(Z)|E_{t}) = 0;

(ii) even convex function, E (f(Z)|E_{t}) ≤ E (f(Z)).

Proof: The first statement follows by symmetry, the second by noting that f(z_{1}) ≤ f(z_{2}) for |z_{1}| ≤ |z_{2}|. □

Lemma 8 Let Z ∼ N(0, σ^{2}), and for t ≥ 3σ let E_{t} be the event |Z| ≤ t, then E (Z^{2}|E_{t}) ≥ 0.97σ^{2}.

Proof: It clearly suffices to prove the Lemma for σ = 1. Then E (Z^{2}|E_{t}) = 1−tϕ(t)/(Φ(t)− \frac{1}{2}). This expression is increasing for all t > 0, and putting t = 3 gives (1−0.0133/0.4987) > 0.97. □

Finally we give a bound on a simple function.

Lemma 9 For z < \frac{1}{10}, ze^{−z} is increasing and bounded above by z − \frac{19}{20}z^{2}.

Proof: The first claim may be proved by elementary calculus. If z ≤ 0, the second follows from e^{|z|} ≥ 1 + |z|. If z > 0, note that the series for e^{−z} is alternating with deceasing terms for z < 1, and hence

e^{−z} < 1 − z + \frac{1}{2}z^{2} = 1 − z(1 − \frac{1}{2}z) ≤ 1 − \frac{19}{20}z,

from which the result follows. □

3 Analysis of the random walk

We assume n ≥ 2. We present the coupling analysis in three sections. In 3.1 we prove boundedness, in 3.2 we show approximate coupling, and in 3.3 we show that exact coupling will eventually occur.

3.1 Boundedness of the walk

It will be necessary for us to be sure that with high probability the random walks do not stray too far from the origin. For Y_{t} (started in the stationary distribution) this can be established easily by a calculation like that leading to (4). However for X_{t} (started at an arbitrary point) it is not so straightforward. To bound this we consider the induced random walk Ξ_{t} = ∥X_{t}∥. We assume only that Ξ_{0} < 2ρ, so in fact the analysis applies to Y_{t} also.

Let us call any set of consecutive values of t for which Ξ_{t} > 2ρ an excursion. We wish to show that any excursion does not include values of Ξ_{t} > 3ρ. Suppose X_{t} = x with Ξ_{t} = ξ ≥ 2ρ, and let s = Ξ_{t+1} − Ξ_{t}. Then we know that f(x) = 2(ξ − ρ), and hence we will assume λ ≥ 2 for the remainder of this section. (The argument may be modified, with a weaker conclusion, for any positive λ.) Also, if x′ = x + σU is the trial point as in section 2.5, then letting ξ′ = ∥x′∥ and s′ = ξ′ − ξ, we have F (x′)/F (x) = e^{−2s′}.

Now write u = x/ξ, and let v = σU.u and w = \sqrt{σ^{2}∥U∥^{2} − v^{2}}. Then

ξ′ = \sqrt{(ξ + v)^{2} + w^{2}} = ξ\sqrt{(1 + v/ξ)^{2} + (w/ξ)^{2}}.

From Lemma 6 we may assume that |v| ≤ ℓ = cσ ln n ≤ \frac{1}{21}. Denote the event |v| ≤ ℓ by E_{1}, and the event w < 2cr/3 by E_{2}. Let E = E_{1} ∩ E_{2}. Then with high probability, E occurs at every step for polynomially many steps, provided c is large enough. Note that E implies σ∥U∥ ≤ cr\sqrt{4/9 + (ln n)^{2}/n} < cr. Thus, if E occurs, the acceptance indicator A has

P(A = 1) = A = \sqrt{e^{−λcr}e^{−2s′}} = e^{−(c′+s′)},

where c′ = λcr/2 ≤ \frac{1}{16}.

Assuming the same event, v/ξ ≥ −ℓ/2cnR = −r ln n/2Rn^{3/2} ≥ −1/4000. Hence we may apply Lemma 4 to give

ξ′ ≤ ξ(1 + v/ξ + (w/ξ)^{2}) = ξ + v + w^{2}/ξ ≤ ξ + v + ϵ,

for some constant ϵ, where 0 ≤ ϵ < 4(cr)^{2}/18cnR ≤ σ^{2}/144. Also it follows that ξ′ ≥ ξ +v, and hence

v ≤ s′ ≤ v + ϵ.

Since ϵ < ℓ/1000, we have |s| ≤ |s′| ≤ ℓ+ϵ < 1.001ℓ ≤ \frac{1}{20}. Now, using Lemmas 7, 8, and 9,

E (s|E) = E (s′e^{−(c′+s′)}|E)

≤ e^{−c′}E ((v + ϵ)e^{−(v+ϵ)}|E_{1})
≤ e^{−c′}E ((v + ϵ) − 0.95(v + ϵ)^{2}|E_{1})
= e^{−c′}E (ϵ(1 − 0.95ϵ) + (1 − 1.9ϵ)v − 0.95v^{2}|E_{1})
≤ e^{−c′}(ϵ − 0.92σ^{2})
≤ −0.85σ^{2},

using ϵ ≤ σ^{2}/144 and c′ ≤ 1/16. Hence |Ξ_{t+1} − Ξ_{t}| < 1.001 cσ ln n and E (Ξ_{t+1} − Ξ_{t}) ≤ −0.85σ^{2}. Thus, suppose the excursion starts at t = a + 1, so Ξ_{a} ≤ 2ρ. Let k_{0} = ⌊c^{4}(ln n)^{4}/2.1σ^{2}⌋. Then, since the process Ξ_{t} is Markovian, Hoeffding’s inequality for Martingales [11] (see also [21]) gives,

P(∃k, 1 ≤ k ≤ k_{0} : Ξ_{a+k} − Ξ_{a} + 0.85σ^{2}k > 0.4c^{4}(ln n)^{4}) ≤ n^{−c^{2} ln n/6}.  (6)

However, (6) implies

P(Ξ_{a+k_{0}} − Ξ_{a} > 0) ≤ n^{−c^{2} ln n/6}.

Thus, with high probability, the excursion must have already ended at step (a + k_{0}). For k < k_{0}, (6) implies

P(∃k : Ξ_{a+k} > 2ρ + 0.4c^{4}(ln n)^{4}) ≤ n^{−c^{2} ln n/6}.

But 0.4c^{4}(ln n)^{4} < 2c^{4}n ≤ 2nR ≤ ρ. Thus each excursion lies entirely in the ball 3ρB, except with probability n^{−c^{2} ln n/3}. Since the algorithm can only perform a polynomial number of excursions, we can choose c so that 3ρB includes the whole walk with high probability.

3.2 Bringing the random walks close

We have two random walks X_{t}, Y_{t} as in section 2.5, with Y_{0} chosen randomly from the equilibrium density π. We will couple them in the following way. Suppose σU is the trial step. Let u_{1}, u_{2}, …, u_{n} be an orthonormal basis for ℝ^{n} such that u_{1} = (Y_{t} −X_{t})/∥Y_{t} −X_{t}∥. Then U = ∑_{i=1}^{n} U_{i}u_{i}, where U_{i} has the density G_{1}. Let us write U′ = ∑_{i=2}^{n} U_{i}u_{i}, so U′ has the density G_{n−1}. Then the trial points are generated by

X′_{t} = X_{t} + U_{1}u_{1} + U′, Y′_{t} = Y_{t} − U_{1}u_{1} + U′.

We accept the trial points independently with probabilities A(X_{t}, X′_{t}), A(Y_{t}, Y′_{t}) respectively. If we accept X′_{t}, then X_{t+1} = X′_{t}, otherwise X_{t+1} = X_{t} and similarly for Y′_{t}. Thus the two random walks are coupled so that, if the two steps are accepted, their steps will be reflections of each other in the hyperplane of points equidistant from X_{t}, and Y_{t}. This has similarities to the coupling given in [17].

We will show that, with this coupling, X_{t}, Y_{t} have a tendency to converge. Specifically, we show that, for some metric d on ℝ^{n}, E (d(X_{t}, Y_{t})) decreases. The metric we choose is the square-root of Euclidean distance, d(X_{t}, Y_{t}) = \sqrt{∥X_{t} − Y_{t}∥}. The rationale for this metric is to downweight larger values of the “neutral” Euclidean distance. We will write d_{t} = d(X_{t}, Y_{t}), and let D denote d_{t}^{2} = ∥Y_{t} − X_{t}∥.

In notation similar to that in 3.1, let us write v = U_{1}, and w = ∥U′∥. Also, abbreviating the indicator variables for acceptance to A_{X}, A_{Y}, let

I = A_{X} + A_{Y} = 1 + A_{X}A_{Y} − Ā_{X} Ā_{Y}
J = |A_{X} − A_{Y} | = |Ā_{X} − Ā_{Y} |

where I ∈ {0, 1, 2} and J ∈ {0, 1}. Then it follows that

d_{t+1} = ^{4} \sqrt{(d_{t}^{2} − vI)^{2} + w^{2}J}.

We will bound E (d_{t+1}|d_{t}) as a function of d_{t}.

As in section 3.1, denote the event |v| ≤ cσ ln n by E_{1}, the event w < 2cr/3 by E_{2}, and E_{1} ∩ E_{2} by E. Then, with high probability, E occurs at every step for sufficiently many steps, and no step size is larger than cr. Now E implies

D − vI ≥ D − 2cσ ln n ≥ \frac{1}{2}D,

provided we do not have the event d_{t} < 2\sqrt{cσ ln n}. Let us call this event C. We will say the walks are “close” when C occurs. Now if C does not occur, using Corollary 1 and Lemma 5,

d_{t+1}/d_{t} = ^{4} \sqrt{(1 − vI/D)^{2} + w^{2}J/D^{2}}.
≤ \sqrt{1 − vI/D} + w^{2}J/D^{2}
≤ 1 − \frac{1}{2}vI/D − v^{2}I^{2}/10D^{2} + w^{2}J/D^{2}  (7)

We have to bound the terms on the right side of (7) conditioned by E.

Note that if E occurs then, for A_{X}, A_{Y},

e^{−λcr} ≤ A ≤ 1, 0 ≤ Ā ≤ 1 − e^{−cλr} ≤ λcr.

Thus, using Lemmas 7 and 8,

E (w^{2}J|E) ≤ E (w^{2}(Ā_{X} + Ā_{Y})|E)
≤ 2λcrE (w^{2}|E_{2}) ≤ 2λcrE (w^{2}) < 2λcr^{3},  (8)

and, providing c ln n ≥ 3,

E (v^{2}I^{2}|E) = E (v^{2}(A_{X} + A_{Y} + 2A_{X}A_{Y})|E)
≥ 2e^{−cλr}(1 + e^{−cλr})E (v^{2}|E_{1}) > 3r^{2}/n,  (9)

since cλr < \frac{1}{8} and E (v^{2}|E_{1}) > 0.97σ^{2}.

It remains to consider E (vI). Now, since E (v|E) = E (v|E_{1}), using Lemma 7 we have

E (vI|E) = E (v + vA_{X}A_{Y} − vĀ_{X} Ā_{Y} |E) = E (vA_{X}A_{Y} |E) − E (vĀ_{X} Ā_{Y} |E) ,

We will bound these two contributions separately. We do this by arguing conditionally on both E_{1} and the value of w (given that this satisfies E_{2}). Since the bounds we obtain are independent of w, we may infer the same bound for conditioning on E. In this respect, for either A_{X} or A_{Y}, the notation A(v) will mean “A considered as a function of v for a given fixed w”.

Thus let p(v) be the conditional density of v given E, and first consider the second term. Then

E (vĀ_{X} Ā_{Y} |E_{1}, w) = ∫_{−∞}^{∞} vĀ_{X}(v)Ā_{Y} (v)p(v)dv
= ∫_{0}^{∞} v(Ā_{X}(v)Ā_{Y} (v) − Ā_{X}(−v)Ā_{Y} (−v))p(v)dv
= 2 ∫_{0}^{∞} \frac{Ā_{X}(v)Ā_{Y} (v) − Ā_{X}(−v)Ā_{Y} (−v)}{2v} v^{2}p(v)dv
= 2 ∫_{0}^{∞} [\frac{d(Ā_{X}(v)Ā_{Y} (v))}{dv}]_{θ} v^{2}p(v)dv,

where −v ≤ θ ≤ v. Note that if Ā_{X}(v)Ā_{Y} (v) is not differentiable everywhere we may consider a suitably close approximant. Now

\frac{d(Ā_{X}(v)Ā_{Y} (v))}{dv} = − (Ā_{X}(v)A_{Y} (v)\frac{d(ln A_{Y} (v))}{dv} + Ā_{Y} (v)A_{X}(v)\frac{d(ln A_{X}(v))}{dv}) ,

and hence, using the bound on the directional derivative of f,

|\frac{d(Ā_{X}(v)Ā_{Y} (v))}{dv}| ≤ 2λ^{2}cr.

Thus, using Lemma 7,

E (vĀ_{X} Ā_{Y} |E_{1}, w) ≤ 4λ^{2}cr ∫_{0}^{∞} v^{2}p(v)dv
= 2λ^{2}crE (v^{2}|E_{1})
≤ 2λ^{2}crσ^{2} = 2λ^{2}cr^{3}/n.  (10)

Now let g(v) = \sqrt{e^{−λcr}F (X_{t} + U′ + vu_{1})}. Then

g(D − v) = \sqrt{e^{−λcr}F (X_{t} + U′ + Du_{1} − vu_{1})} = \sqrt{e^{−λcr}F (Y_{t} + U′ − vu_{1})}.

Note that g(v) is a log-concave function (of one argument), and

A_{X}(v) = \frac{g(v)}{\sqrt{F (X_{t})}}, A_{Y} (v) = \frac{g(D − v)}{\sqrt{F (Y_{t})}}.

Therefore, writing C = (F (X_{t})F (Y_{t}))^{−1/2}, we have

E (vA_{X}A_{Y} |E_{1}, w) = ∫_{−∞}^{∞} vA_{X}(v)A_{Y} (v)p(v)dv
= C ∫_{−∞}^{∞} vg(v)g(D − v)p(v)dv
= C ∫_{0}^{∞}(g(v)g(D − v) − g(−v)g(D + v))vp(v)dv,
≥ 0,  (11)

since the integrand is non-negative by the log-concavity of g(v). To show this, let t(v) = − ln g(v), so t(v) is convex. We have now to show that t(v) − t(−v) ≤ t(D + v) − t(D − v). But this is clear since

t(v) − t(−v) = ∫_{−v}^{v} t′(v)dv, t(D + v) − t(D − v) = ∫_{D−v}^{D+v} t′(v)dv,

for some nondecreasing function t′(v).

We note in passing that it is only to establish (11) that we require log-concavity of F. All our other estimates simply require suitable “smoothness” of F. On the other hand, for something like (11) to be true, it appears we need a property close to log-concavity. It is also the proof of (11) which necessitates our departure from the “usual” Metropolis process, since we have to factor out the different contributions from the points X_{t}, Y_{t}.

From (10) and (11) we now have

E (vI|E) ≥ −2λ^{2}cr^{3}/n.  (12)

Therefore, putting (8), (9) and (12) into (7), we have

E (d_{t+1}/d_{t}|d_{t}) ≤ 1 + \frac{λ^{2}cr^{3}}{nD} − \frac{3r^{2}}{10nD^{2}} + \frac{2λcr^{3}}{D^{2}}
= 1 + \frac{r^{2}}{10nD^{2}} (3 − 10λ^{2}crD − 5λcrn)
≤ 1 − \frac{r^{2}}{10nD^{2}} ,

provided 10λ^{2}crD + 5λcrn ≤ 2. This will be true if 20.1λ^{2}crρ ≤ 1, using R ≥ c^{4} and the fact that, with high probablity, D < 4ρ from sections 2.3 and 3.1. Thus let us take

r = \frac{1}{20.1cλ^{2}ρ} = \frac{1}{20.1c^{2}λ^{2}nR}.  (13)

Thus if R = O(n), r = Ω(1/n^{2}). Now, using D < 4ρ,

E (d_{t+1}) ≤ E (d_{t}) (1 − \frac{1}{65000nc^{2}λ^{4}ρ^{4}} ) = E (d_{t}) (1 − \frac{1}{65000n^{5}c^{6}λ^{4}R^{4}} ) .  (14)

The reduction factor is (1 − Ω(1/n^{9})), assuming R = O(n). Therefore, since d_{0} = O(n) and σ = Ω(n^{−5/2}), in O(n^{9} log n) steps we must either have encountered C or E (d_{t}) ≤ \frac{1}{10} \sqrt{cσ ln n}. In the latter case, by Markov’s inequality, P(C̄) ≤ \frac{1}{20}. Hence, in a sequence of O(n^{9} log n) steps, there is probability at least \frac{19}{20} that C holds at some step.

3.3 Making the random walks meet

When C occurs, we switch to a different coupling. Note that X_{t}, Y_{t} are still relatively far apart compared to the expected movement in any direction in a single step. Therefore we choose instead to correlate a sequence of steps in order to have sufficient probability that X_{t} = Y_{t} at the end of the sequence.

Thus suppose C occurs at step τ, but C̄ held at (τ −1). Let s_{X,t} = X′_{t}−X_{t}, s_{Y,t} = Y′_{t}−Y_{t} denote the trial steps in the two walks. Let k = ⌈(4c ln n)^{2}⌉, and

S_{X} = ∑_{t=τ}^{τ+k−1} s_{X,t},  Q_{X} = X_{τ} + S_{X},  (15)

and similarly S_{Y}, Q_{Y}. Note that X_{τ+k} = Q_{X} and Y_{τ+k} = Q_{Y} if all trial steps are accepted for t = τ, τ + 1, …, τ + k − 1.

Now S_{X}, S_{Y} have density G_{n}(x/σ\sqrt{k}). We will couple S_{X}, S_{Y}. Then s_{X,t}, s_{Y,t} are generated “independently” with density ϕ(x/σ), subject to the sum condition (15).

As in section 3.2, let u_{1}, u_{2}, …, u_{n} be a basis such that u_{1} = (Y_{τ} − X_{τ})/D, where D = ∥Y_{τ} −X_{τ}∥. Suppose S_{X} = (v, v_{2}, …, v_{n}) in this coordinate system, where v, v_{2}, …, v_{n} are independent with probability density ψ(x) proportional to ϕ(x/σ\sqrt{k}). We now take S_{Y} = (−v′, v_{2}, …, v_{n}), where v′ will be defined in terms of v. Note that we have Q_{X} = Q_{Y} if and only if v′ = D − v. Thus, for given v, define v′ by

v′ = \{ \begin{matrix}D − v, & with probability min\{1, ψ(D − v)/ψ(v)\}, \\ v, & otherwise. \end{matrix}

Intuitively, we are reflecting the tail of the density ψ(v) about the point v = \frac{1}{2}D. Then v′ has probability density function

ψ(D − v′) min{1, ψ(v′)/ψ(D − v′)} + ψ(v′)(1 − min{1, ψ(D − v′)/ψ(v′)}) = ψ(v′).

Hence, since ψ(v′) = ψ(−v′), −v′ has density ψ, as required. Also

P(Q_{X} = Q_{Y}) = P(v′ = D − v) = ∫_{−∞}^{∞} min{ψ(D − v), ψ(v)}dv
= 2 ∫_{−∞}^{−D/2} ψ(v)dv
= 2Φ (\frac{−D}{2σ\sqrt{k}} )
≥ 2Φ(−\frac{1}{2}) > 0.616.

We must now consider the probability that all steps accept. To maximise this, we choose to make the random variables A_{X}, A_{Y} dependent at each step. Their values are generated as follows. Let w be a random variable, independent for each step, with distribution uniform on [0, 1]. Then A_{X} = 1 if w ≥ A_{X}, otherwise A_{X} = 0. The value of A_{Y} is chosen similarly, using the same w. The advantage of this is that now

P(A_{X} = A_{Y} = 1) = min{A_{X}, A_{Y} } ≥ e^{−cλr}.

Now if G is the event that all trial steps are accepted in both walks for t = τ, τ + 1, …, τ + k − 1, we have

P(G) ≥ (e^{−cλr})^{k} ≥ e^{−1/18} > 0.945,

using the value of r given in (13). Hence we have

P(X_{τ+k} = Y_{τ+k}) = P(Q_{X} = Q_{Y})P(G) > 0.616 × 0.945 > 0.582.

The number of steps, O((log n)^{2}), is asymptotically negligible compared with the number required in section 3.2 to achieve P(C) > \frac{19}{20}. Hence, after O(n^{5}R^{4} log n) steps, the probability that the walks have failed to couple is at most \frac{1}{20} + \frac{19}{20}(1 − 0.582) < 0.45.

The other possible sources of failure to couple, for example that some step is too large, can be made very small by appropriate choice of c. We may arrange that their combined probability is less than 0.05 for all n. Hence, the probability of coupling will be at least \frac{1}{2}, as required.

This completes the analysis of the algorithm. Putting R = O(n), our “mixing time” is O^{∗}(n^{9}), where O^{∗}(·) is the notation which hides factors of log n. This is larger than is known to be possible [15]. Note that the mixing time is O^{∗}(nρ^{2}/r^{2}) = O^{∗}((ρ/σ)^{2}), By comparison with one-dimensional random walk, this is of the optimal form. The problem is that the step size, r, is too small. We attempt to offset this difficulty in the next section.

4 Improvements

In this section we describe two refinements which reduce the required number of steps of the random walk from O^{∗}(n^{9}) to O^{∗}(n^{6.5}).

4.1 A faster simulation of the random walk

An idea similar to that used to couple the walks in section 3.3 can be employed to simulate the random walk with fewer actual steps. A difficulty with our walk is that r = Θ(1/ρ) is smaller than values known to be achievable [19]. However, the fact that our trial steps are so small means that almost all will be accept, and hence we may simulate many steps simultaneously in a larger “step”.

Suppose at step t, we attempt to perform k steps, where k is a power of 2. We will call this a k-step. Denote the current k-step by S. Let S be the vector sum of the k trial steps in S, as in section 3.3. We generate the components of S as independent N(0, kσ^{2}). If a is the probability that all these trial steps are accepted, then we have

a = \sqrt{e^{−kλcr}F (X_{t} + S)/F (X_{t})},

since all other terms cancel. Thus we accept S with probability a. If we accept, then we update t to t + k and consider the next k-step.

If however, with probability (1 − a) ≤ kλcr, we do not accept S, we know that at least one constituent step was rejected. We will say S rejects. We subdivide S into two \frac{1}{2}k-steps S_{1}, S_{2} (in that order). We generate the corresponding vector sums S_{1}, S_{2} subject to the condition S_{1} + S_{2} = S. Suppose s_{1} is the ith component of S_{1}, and s the ith component of S. Then s_{1} has conditional density

\frac{(e^{−s_{1}^{2}/kσ^{2}}/\sqrt{πkσ^{2}})(e^{−(s−s_{1})^{2}/kσ^{2}}/\sqrt{πkσ^{2}})}{e^{−s^{2}/2kσ^{2}}/\sqrt{2πkσ^{2}}} = \frac{1}{\sqrt{\frac{1}{2}πkσ^{2}}} e^{−2(s_{1}−s/2)^{2}/kσ^{2}}.

Thus s_{1} ∼ N(\frac{1}{2}s, \frac{1}{4}kσ^{2}). Hence we can generate S_{1}, S_{2}. Now let

a_{1} = \sqrt{e^{−kλcr/2}F (X_{t} + S_{1})/F (X_{t})}

be the unconditioned acceptance probability for S_{1}. We must condition on the event that S rejects. Let q be the conditional probability that S_{1} does not reject. Then clearly q = (a_{1} − a)/(1 − a). We use q to decide (by random generation) if S_{1} rejects. If not, we can set X_{t+k/2} = X_{t} + S_{1} and consider S_{2} (with vector sum S_{2}), conditional that it rejects. The base case of the recursion, k = 1, is clearly a (single) step.

If S_{1} rejects, then we (recursively) halve it conditional on this event. Eventually, using this recursion, we compute X_{t+k/2}. Then we consider S_{2}, with vector sum S_{2}, now as an unconditioned \frac{1}{2}k-step. Hence we ultimately compute X_{t+k}.

To analyse this procedure, let t_{u}(k) denote the expected total number of k′-step s which are accepted within S, for all k′ ≤ k. This effectively determines the execution time for S. Let t_{c}(k) be this number conditional that S rejects. For notational convenience, let β = λcr. Then, from the above,

t_{u}(k) ≤ a + (1 − a)t_{c}(k) ≤ 1 + kβt_{c}(k)
t_{c}(k) ≤ q(1 + t_{c}(\frac{1}{2}k)) + (1 − q)(t_{c}(\frac{1}{2}k) + t_{u}(\frac{1}{2}k)) ≤ 1 + t_{c}(\frac{1}{2}k) + t_{u}(\frac{1}{2}k).

Hence we have

t_{c}(k) ≤ 2 + (1 + kβ)t_{c}(\frac{1}{2}k) ≤ 2 + e^{kβ}t_{c}(\frac{1}{2}k).  (16)

Now (16) has solution t_{c}(k) ≤ 2e^{2kβ} lg k, as is easily shown by induction. (Note t_{c}(1) = 0.) Hence

t_{u}(k) ≤ 1 + 2e^{2kβ}kβ lg k.

Let k = 2^{⌊lg(1/β)⌋}. Note that k ≤ 1/β and k = Ω(1/r). Then

t_{u}(k) ≤ 1 + 2e^{2} lg(1/β) = O(log(1/r)).

Thus we trade a factor 1/r in the number of steps for log(1/r). This reduces the number from O^{∗}(n^{9}) to O^{∗}(n^{7}).

4.2 An even faster simulation

We use the factor M = e^{−λcr} in the calculation of A to ensure that the result is at most 1. However this factor is smaller than necessary. From the subgradient inequality (2), we have

f(X′_{t}) ≥ f(X_{t}) + ∇.(X′_{t} − X_{t}),

where ∇ is a subgradient at X_{t}. Now ∇.(X′_{t} − X_{t}) ∼ N(0, σ^{2}∥∇∥^{2}) and ∥∇∥ ≤ λ. Thus, using Lemma 6,

P(∇.(X′_{t} − X_{t}) > cλr ln n/\sqrt{n}) < n^{−c^{2} ln n/2}.

Hence F (X′_{t})/F (X_{t}) ≤ e^{−cλr ln n/\sqrt{n}} = M′, say, with high probability. Let us write A′ = \sqrt{M′F (X′_{t})/F (X_{t})}, A_{0} = \sqrt{M/M′}, and note that A_{0} > \frac{3}{4} is a constant. Then the acceptance probability A = A_{0}A′. We write A = A_{0}A′ for the corresponding indicator variables. Hence we may simulate a step of the original algorithm by first choosing A_{0} and generating A′ only if A_{0} = 1.

Now consider the algorithm which uses only A′. To each step of the modified algorithm there will be a geometrically distributed number of steps of the original algorithm with expectation 1/A_{0}. However, observe that it does not matter if we stop at a random time in the original algorithm provided that it is independent of the random walk and it exceeds the bound on the coupling time (at least with high probability).

Thus we need not generate the sum of geometrics if the number of steps of the modified algorithm exceeds the coupling time. Note that the saving in the number steps from generating this sum (as a negative binomial) would be only a small constant factor, and we will not consider it here.

Now consider a k-step of the modified algorithm. We use the same notation as in section 4.1. We have, with high probability, ∥S∥ ≤ cr\sqrt{k}, and hence f(X_{t} + S) − f(X_{t}) ≤ λcr\sqrt{k}. Thus

a ≥ exp(−\frac{1}{2}λcr(k ln n/\sqrt{n} + \sqrt{k})) = exp(−\frac{1}{2}β(k ln n/\sqrt{n} + \sqrt{k})),

and hence (1 − a) ≤ β(k ln n/\sqrt{n} + \sqrt{k}). With this modification (16) becomes

t_{c}(k) ≤ 2 + e^{β(k ln n/\sqrt{n}+\sqrt{k})/2}t_{c}(\frac{1}{2}k),

from which it follows that

t_{u}(k) ≤ 1 + e^{2β(k ln n/\sqrt{n}+2\sqrt{k})}β(k ln n/\sqrt{n} + \sqrt{k}) lg k.

Hence we may choose k = 2^{⌊lg(\sqrt{n}/β ln n)⌋}, and simulate Ω(\sqrt{n}/r) steps of the modified algorithm in O(log(\sqrt{n}/r)) k-steps. This reduces the required number of “steps” to O^{∗}(n^{6.5}).

5 Conclusions

We have given a coupling proof for the polynomial-time convergence of a random walk for generating an approximately uniformly distributed point in a convex body. The argument is completely elementary, and does not require any auxiliary proof of an isoperimetric inequality, as with the usual approach. On the other hand, the running time of the algorithm is inferior to the best conductance methods. We obtain mixing time O^{∗}(n^{6.5}), whereas Kannan, Lovász and Simonovits [15] obtain O^{∗}(n^{5}) (and an amortized time for multiple points which is even smaller). However, the present conductance techniques are the product of years of refinement. We have drawn on some of these developments, but our proof is a first attempt using a different approach. It would be surprising if it cannot be improved further. The present difficulty is the small step size r. We have described modifications to mitigate this problem, but genuine progress would result from an improvement in the arguments to allow r to be increased.

References

[1] D. Aldous. Random walks on finite groups and rapidly mixing Markov chains. In Séminaire de Probabilitès XVII 1981/1982, volume 986 of Springer-Verlag Lecture Notes in Mathematics, pages 243–297, 1981.

[2] D. Applegate and R. Kannan. Sampling and integration of near log-concave functions. In Proceedings of the Twenty Third Annual ACM Symposium on Theory of Computing, pages 156–163, New Orleans, Louisiana, 1991.

[3] A. D. Barbour, L. Holst, and S. Janson. Poisson Approximation. Oxford Science, 1992.

[4] A. Z. Broder. How hard is it to marry at random ? In Proceedings of the Eighteenth Annual ACM Symposium on Theory of Computing, pages 50–58, Berkeley, California, 1986.

[5] P. Diaconis. Group Representations in Probability and Statistics, volume 11 of Lecture Notes-Monograph Series, chapter 4e, pages 84–87. Institute of Mathematical Statistics, Hayward, California, 1988.

[6] M. E. Dyer and A. M. Frieze. Computing the volume of convex bodies: A case where randomness provably helps. In Symposia in Applied Mathematics, volume 44, pages 123–169. American Mathematical Society, 1991.

[7] M. E. Dyer, A. M. Frieze, and R. Kannan. A random polynomial-time algorithm for approximating the volume of convex bodies. Journal of the ACM, 38:1–17, 1991.

[8] A. M. Frieze, R. Kannan, and N. Polson. Sampling from log-concave distributions. Annals of Applied Probability, to appear.

[9] D. Griffeath. Coupling methods for Markov processes. In Gian-Carlo Rota, editor, Studies in Probability and Ergodic Theory, volume 2 of Advances in Mathematics Supplementary Studies, pages 1–43. Academic Press, 1978.

[10] M. Grötschel, L. Lovász, and A. Schrijver. Geometric Algorithms and Combinatorial Optimization. Springer-Verlag, New York, 1988.

[11] W. Hoeffding. Probability inequalities for sums of bounded random variables. Journal of the American Statistical Association, 58:13–30, 1963.

[12] M. R. Jerrum. A very simple algorithm for estimating the number of k-colourings of a low-degree graph. Random Structures and Algorithms, 7:157–165, 1995.

[13] M. R. Jerrum and A. J. Sinclair. Approximating the permanent. SIAM Journal on Computing, 18:1149–1178, 1989.

[14] R. Kannan. Markov chains and polynomial time algorithms. In 35th Annual IEEE Symposium on Foundations of Computer Science, pages 656–671, Santa Fe, New Mexico, 1994. IEEE.

[15] R. Kannan, L. Lovász, and M. Simonovits. Random walks and an O^{∗}(n^{5}) volume algorithm for convex bodies. Preprint, January 1996.

[16] R. Kannan, L. Lovász, and M. Simonovits. Isoperimetric inequalities for convex sets. Discrete and Computational Geometry, to appear.

[17] T. Lindvall and L. C. G. Rogers. Coupling of multidimensional diffusions by reflection. Annals of Probablilty, 14:860–872, 1986.

[18] L. Lovász and M. Simonovits. Mixing rate of markov chains, an isoperimetric inequality and computing the volume. In Proceedings of the 31st IEEE Symposium on Foundations of Computer Science, pages 346–355, St. Louis, Missouri, 1990.

[19] L. Lovász and M. Simonovits. Random walks in a convex body and an improved volume algorithm. Random Structures and Algorithms, 4:359–412, 1993.

[20] P. Matthews. Generating a random linear extension of a partial order. Annals of Probability, 19:1367–1392, 1991.

[21] C. McDiarmid. On the method of bounded differences. In J. Siemons, editor, Surveys in Combinatorics, LMS Lecture Note Series 141, pages 148–188. Cambridge University Press, 1989.

[22] A. J. Sinclair and M. R. Jerrum. Approximate counting, uniform generation and rapidly mixing Markov chains. Information and Computation, 82:93–133, 1989.

[23] H. Thorisson. Coupling methods in probability theory. Scandinavian Journal of Statistics, 22:159–182, 1995.