# efficient_sampling_from_timevarying_logconcave_distributions__79beb6c3.pdf

Journal of Machine Learning Research 18 (2017) 1-29 Submitted 6/14; Revised 11/16; Published 10/17

Efficient Sampling from Time-Varying Log-Concave Distributions

Hariharan Narayanan harin@uw.edu Department of Statistics and Department of Mathematics University of Washington

Alexander Rakhlin rakhlin@wharton.upenn.edu Department of Statistics University of Pennsylvania

Editor: Robert Mc Culloch

We propose a computationally efficient random walk on a convex body which rapidly mixes with respect to a fixed log-concave distribution and closely tracks a time-varying log-concave distribution. We develop general theoretical guarantees on the required number of steps; this number can be calculated on the fly according to the distance from and the shape of the next distribution. We then illustrate the technique on several examples. Within the context of exponential families, the proposed method produces samples from a posterior distribution which is updated as data arrive in a streaming fashion. The sampling technique can be used to track time-varying truncated distributions, as well as to obtain samples from a changing mixture model, fitted in a streaming fashion to data. In the setting of linear optimization, the proposed method has oracle complexity with best known dependence on the dimension for certain geometries. In the context of online learning and repeated games, the algorithm is an efficient method for implementing no-regret mixture forecasting strategies. Remarkably, in some of these examples, only one step of the random walk is needed to track the next distribution.1

1. Introduction

Let K be a compact convex subset of Rd with non-empty interior. Let µ0, . . . , µt, . . . be a sequence of probability measures with support on K. Suppose each probability distribution µt has a density

dx = e st(x)

Zt , Zt = Z

x K e st(x)dx (1)

with respect to the Lebesgue measure, where each st(x) is a convex function on K. This paper proposes a Markov Chain Monte Carlo method for sequentially sampling from these distributions. The method comes with strong mixing time guarantees, and is shown to be applicable to a variety of problems. Observe that, by definition, the distributions µt are

1. An extended abstract containing partial results appeared in the proceedings of the NIPS 2010 conference (Narayanan and Rakhlin, 2010).

c 2017 Hariharan Narayanan and Alexander Rakhlin.

License: CC-BY 4.0, see https://creativecommons.org/licenses/by/4.0/. Attribution requirements are provided at http://jmlr.org/papers/v18/14-223.html.

Narayanan and Rakhlin

log-concave, and thus our work falls within the emerging body of literature on sampling from log-concave distributions. The problem of sampling from distributions arises in many areas of statistics, most notably in Bayesian inference (Robert and Casella, 2004). In particular, Sequential Monte Carlo methods (Doucet et al., 2001) aim to sample from time-varying distributions. The need for such methods arises, for instance, in the case of online arrival of data: it is desirable to be able to update the posterior distribution at a low computational cost. If the distributions are changing slowly with time, sequential methods can re-use samples from the previous distribution and perform certain re-weighting to track the next distribution, thus saving computational resources. These ideas are exploited in particle filtering methods (see (Chopin, 2002; Doucet et al., 2001) and references therein). Beyond Bayesian inference, other applications of sampling from distributions include simulated annealing, global optimization, and regret minimization. The main critique of the MCMC methods is, in many situations, the lack of mixing time analysis. In practice, the number of steps of the chain required to obtain an honest sample from a distribution is mostly calculated based on heuristics. There is a growing body of literature that presents exceptions to these heuristic approaches (Diaconis, 2013). Coupling methods, spectral gap methods, as well as the more recent study of positive Ricci curvature, yield geometric decrease of the distance to the desired stationary distribution a property known as geometric ergodicity. The most well-understood cases in this context are those with a finite or countable state space (see (Meyn and Tweedie, 2009; Diaconis, 2009)). In contrast, we are interested in a random walk on a non-discrete set. This paper is focused on a particular circle of problems defined via log-concave distributions. These distributions constitute an important subset of the set of unimodal distributions, a fact that has been recognized within Statistics (see e.g. (Walther, 2009)). We are not the first to study mixing times for such distributions: this line of work started with the breakthrough paper of (Dyer et al., 1991), followed by a series of improvements (Frieze et al., 1994; Lov asz, 1999; Lov asz and Vempala, 2006, 2007). However, the recent advances in (Kannan and Narayanan, 2012) on sampling from convex bodies give an edge to obtaining stronger guarantees. In particular, a variant of the Dikin walk studied in this paper is analyzed in (Narayanan, 2016) for the case of log-linear distributions. Extending the study of this random walk, we show rapid mixing to a log-concave distribution, and further show that we can provably track a changing log-concave distribution with a small number (or even only one step) of a random walk, provided that the distribution changes slowly enough. Such a result seems out of reach with other random walk methods due to the lack of scale-free bounds on conductance. We assume that we can compute a self-concordant barrier (see Section 5 and Appendix 8) for the set K, a requirement that is satisfied in many cases of interest. For instance, the self-concordant barrier can be readily computed in closed form if K is defined via linear and quadratic constraints. While the availability of the barrier is a stronger assumption than, for instance, access to a separation oracle for K, the barrier gives a better handle on the geometry of the space and yields fast mixing of the Markov chain. In Section 5, we illustrate the method within several diverse application domains. As one of the examples, we consider the problem of updating the posterior with respect to a conjugate prior in an exponential family, where the parameter is taking values in a space of a

Efficient Sampling

fixed dimensionality given by the sufficient statistics. The constraints then constitute a prior knowledge about the possible location of the parameter. As another example, we consider sampling from a time-varying truncated distribution, as well as the extension to sampling from mixture models fitted to streaming data. We employ the sampling technique to the classical problem of linear optimization via simulated annealing. The final example concerns the problem of regret minimization where the log-concave distribution arises naturally from the exponential weighting scheme. The paper is organized as follows. In the next section we study the geometry of the set K induced by a self-concordant barrier and prove a key isoperimetric inequality in the corresponding Riemannian metric. The Markov chain for a given log-concave distribution is defined in Section 3. Conditions on the size of a step are introduced in Section 3.1, and a lower bound on the conductance of the chain is proved in Section 3.2. Section 4 contains main results about tracking time-varying distributions given appropriate measures of change between time steps. Section 5 is devoted to applications. Finally, Sections 6 and 7 contain all the remaining proofs.

2. Geometry Induced by the Self-Concordant Barrier

The Markov chain studied in this paper uses as a proposal a Gaussian distribution with a covariance that approximates well the local geometry of the set K at the current point. This local geometry plays a crucial role in the theory of interior point methods for optimization, yet for our purposes a handle on the local geometry yields a good lower bound on conductance of the Markov chain. Further intriguing similarities between optimization and sampling will be pointed out throughout the paper. We refer to (Nemirovskii, 2004) for an introduction to the theory of interior point methods, a subject centered around the notion of a self-concordant barrier. Once we have defined a self-concordant barrier for K, the local geometry is defined through the Hessian of the barrier at the current point. For the reader unfamiliar with the literature on self-concordant barriers, it is useful to think of the simple 2-dimensional example of F(x) = log x1 log x2 defined on the positive quadrant K = [0, )2, with x = (x1, x2). The Hessian of F at x is 1/x2 1 0 0 1/x2 2

, and a unit ball centered at x and reshaped according to this matrix is the

ellipse that touches the axes, as shown below. That is, the ellipsoid corresponds well to the

local geometry of K at the point x. For a function F on the interior int(K) having continuous derivatives of order k, for vectors h1, . . . , hk Rd and x int(K), for k 1, we recursively define

Dk F(x)[h1, . . . , hk] lim ϵ 0 Dk 1(x + ϵhk)[h1, . . . , hk 1] Dk 1(x)[h1, . . . , hk 1]

Narayanan and Rakhlin

where D0F(x) F(x). Let F be a self-concordant barrier of K with a parameter ν (see Appendix 8 for the precise definition and Section 5 for examples). The value of the parameter ν depends on the shape of K and the choice of F. The barrier parameter will enter many of the results in this paper, but it is worth mentioning that there always exists a barrier with ν = O(d). The barrier induces a Riemannian metric whose metric tensor is the Hessian of F (Nesterov and Todd, 2008). In other words, the metric tensor on the tangent space at x assigns to a vector v the length v 2 x D2F(x)[v, v],

and to a pair of vectors v, w, the inner product

v, w x D2F(x)[v, w] .

The unit ball in x around a point x is called the Dikin ellipsoid (Nemirovskii, 2004), and it describes well the local geometry of K at x in the following sense: (i) the unit Dikin ellipsoid at any point is contained in K, and (ii) the Dikin ellipsoid of radius r = 2(1 + 3ν) contains sym(K, x) = {K x} {K x}. The random walk, introduced in the next section, is anisotropic, i.e. the steps change in size and shape from point to point. It is then useful to connect the properties of this random walk directly to the Riemannian distance that is defined in terms of the Hessian of F. More precisely, for x, y K, let ρ(x, y) be the Riemannian distance ρ(x, y) = infΓ R

z dΓ z where the infimum is taken over all rectifiable paths Γ from x to y. Let M be the metric space whose point set is K and metric is ρ, and define ρ(S1, S2) = inf x S1,y S2 ρ(x, y). The first main

ingredient of the analysis is an isoperimetric inequality.

Theorem 1 Let S1 and S2 be measurable subsets of K and µ a probability measure supported on K that possesses a density whose logarithm is concave. Then it holds that

µ((K \ S1) \ S2) 1 2(1 + 3ν)ρ(S1, S2)µ(S1)µ(S2).

The theorem ensures that two subsets well-separated in ρ distance must have a large mass between them. A lower bound on conductance of our Markov chain will follow from this isoperimetric inequality. We remark that convexity of the set K is crucial for the above property. A classical example of a non-convex shape with a bottleneck is a dumbbell. For this body, the above statement clearly fails, and a local random walk on such a body gets trapped in either of the two parts for a long time.

2.1 Connections to Interior Point Methods

Interestingly, the idea of tracking a changing distribution with only one step of a random walk parallels the technique of following a central path in the theory of interior point methods for optimization, as discussed in (Narayanan and Rakhlin, 2010). In the analysis of interior point methods, the local (according to the Dikin ellipsoid introduced in the next section) quadratic convergence of the Damped Newton step counters the slowly changing temperature parameter of the barrier to ensure sufficiency of one optimization step; in our method, the geometric ergodicity of the scale-free random walk (which is based on the

Efficient Sampling

local shape of the Dikin ellipsoid) balances the additive change in the distribution due to the changing temperature. Once again, the scale-free property of the Dikin walk, introduced below, is crucial for drawing this parallel between one-step interior point methods and our one-step random walk. We remark that a different parallel between sampling and interior point methods has been recently outlined in (Abernethy and Hazan, 2016). The authors show that the mean of the log-concave distribution coincides with the central path when the barrier is chosen to be entropic. Since one still needs to sample from the distribution, this does not imply an algorithmic equivalence of the two methods unless one proves a scale-free bound on conductance, as we do in this paper.

3. The Markov Chain

Let B be the Borel σ-field on K. Given an initial probability measure on K, a Markov chain is specified by a collection of one-step transition probabilities

{P(x, B), x K, B B}

such that x 7 P(x, B) is a measurable map for any B B and Px( ) P(x, ) is a probability measure on K for any x K. For x int(K), let Gr x denote the unique Gaussian probability density function on Rd

Gr x(y) exp d x y 2 x r2 + V (x) , V (x) 1

2 ln det D2F(x)

and r is a parameter that is chosen according to a condition specified below. The covariance of this distribution is given by the Hessian of F at point x, and thus the contour lines are scaled Dikin ellipsoids. The Markov chain considered in this paper is based on the Dikin Walk introduced by Kannan and Narayanan (2012). Adapted to sampling from log-concave distributions in this paper, the Markov chain is parametrized by a convex function s and a step size r. Rather than writing out the unwieldy explicit form of the transition kernel Px, we can give it implicitly as the following random walk:

With probability 1/2, set w := x.

With probability 1/2, sample z from Gr x and

If z / K, let w := x.

If z K, let w :=

( z with prob. min 1, Gr z(x) exp(s(x)) Grx(z) exp(s(z))

x otherwise.

The Markov chain is lazy, as it stays at the current point with probability at least 1/2. This ensures uniqueness of the stationary distribution (Lov asz and Simonovits, 1993). Furthermore, a simple calculation shows that the detailed balance conditions are satisfied with

Narayanan and Rakhlin

respect to a stationary distribution µ whose density (with respect to the Lebesgue measure) is proportional to exp( s(x)). Indeed, to see that µ(x)Px(dz) = µ(z)Pz(dx), it suffices to observe that

exp( s(x))Gr x(z) min 1, Gr z(x) exp(s(x)) Grx(z) exp(s(z))

= exp( s(z))Gr z(x) min 1, Gr x(z) exp(s(z)) Grz(x) exp(s(x))

Therefore the Markov chain is reversible and has the desired stationary measure µ. The value of r has a specific meaning: most of the y s sampled from Gr x are within a thin Dikin shell of radius proportional to (E x y 2 x)1/2 = r by measure-concentration arguments. We will therefore refer to r as the effective step size . An important and non-trivial result from the theory of interior point methods is that the unit Dikin ellipsoid is contained in the set K and gives a good approximation to the local geometry of the set (see Figure 1 below). Thanks to this fact, the sampling procedure has in general better mixing properties than the Ball Walk (Lov asz and Simonovits, 1993; Vempala, 2005).

3.1 Step Size Conditions

The analysis of the Markov chain requires the steps r to be not too large to ensure that different enough transition probability functions happen only for far away points. The precise upper bounds on r depend on the convex function s(x) and can be calculated on the fly when we move to the setting of a time-varying function. We give four conditions:

Sufficient Condition 1 (Linear Functions) If s is linear, we may set r = 1/d.

Sufficient Condition 2 (Lipschitz Functions) For a function s that is L-Lipschitz with respect to the Euclidean norm, we may set the step size r = min 1

Sufficient Condition 3 (Smooth Functions) Suppose s has Lipschitz-continuous gradients: there exists σ > 0 such that s(x) s(y) σ x y . We may then set the step size to be min n 1 d, 1 σ

These three conditions can be shown to follow from a more general sufficient step size condition that is based on local information:

Sufficient Condition 4 (General Condition) Fix constants C, C > 0. Given the convex function s(x), the step size r min 1

d, r is a valid choice if there exists a linear function < g, x > such that

r sup n r : z, w K with z w z C r, s(z) s(w) g, z w < C o

The condition says that for two points, with one being inside the O(r)-Dikin ellipsoid around the other point, the function is within a constant of being linear. It follows from the last condition that, for instance, if s(x) = b, x + a(x) is a sum of a linear and a non-linear Lipschitz part, the step size is only affected by the Lipschitz constant of the non-linear part.

Efficient Sampling

It is simple to verify that the step size in Condition 2 satisfies Condition 4. Indeed, for any w such that z w z C r, we have z w C r R (where R is the radius of the largest ball contained in K). Take gz and gw to be any subgradients of s at z and w, respectively. We then have

|s(z) s(w) gw, z w | gz gw, z w 2L z w 2 .

Notice that for Condition 3, the above calculation becomes

gz gw, z w σ z w 2 1 .

In the remainder of this paper, C will denote a universal constant that may change from line to line. The exact value of the final constant in Lemma 4 below can be traced in the proofs; we omit this calculation for the sake of brevity.

3.2 Conductance of the Markov Chain

In order to show rapid mixing of the proposed Markov chain, we prove a lower bound on its conductance

Φ inf µ(S1) 1

S1 Px(K \ S1)dµ(x)

µ(S1) , (2)

where Px is the one-step transition function defined earlier. Once such a lower bound is established, the following general result on the reduction of distance between distributions will imply exponentially fast convergence.

Theorem 2 ((Lov asz and Simonovits, 1993)) Let γ0 be the initial distribution for a lazy reversible ergodic Markov chain whose conductance is Φ and stationary measure is γ. For every bounded f, let f γ q R

K f(x)2dγ(x). For any fixed f, let Ef be the map that

takes x to R

K f(y)d Px(y). Then if R

K f(x)dγ(x) = 0, it holds that

To prove a lower bound on conductance Φ, we first relate the Riemannian metric ρ to the proposed Markov Chain. Intuitively, the following result says that for close-by points, their transition distributions cannot be far apart in the total variation distance d TV .

Lemma 3 If x, y K and ρ(x, y) r C

d for some constant C, then

d TV (Px, Py) 1 1

for some constant C .

Lemma 3 together with the isoperimetric inequality of Theorem 1 give a lower bound on conductance of the Markov Chain.

Narayanan and Rakhlin

Lemma 4 Let µ be a log-concave distribution with support on K whose density with respect to the Lebesgue measure is proportional to exp{ s(x)}, and suppose an appropriate step size condition (Section 3.1) for the Markov chain is satisfied. Then there exists a constant C > 0 such that the conductance of the above Markov chain is bounded below as

We remark that the step size r enters the lower bound on Φ. While we would like the steps to be large, the conditions outlined earlier dictate a limitation on how large r can be. In particular, we always have r 1/d. The step size needs to be even smaller for functions s for which a linear approximation is poor. We also remark that Lemma 4 establishes a scale-free bound on the conductance, that is, a bound that does not depend on the measure of the set S1 in the definition (2). Such a scale-free conductance is needed for the geometric ergodicity of the chain.

4. Tracking the Distributions

Having specified the Markov chain and the step size, we now turn to the problem of tracking a sequence of distributions µ1, . . . , µt, . . .. For each t 1, define a Markov chain with parameters rt and st, and let its transition kernel be denoted by Pt(x, B) for x K and B B. Let Φt denote the conductance of this chain. The chain will be run for τt steps starting from the end of the chain at time t 1. Formally, let the i-th step of the t-th chain be denoted by the random variable Xt,i. Define τ0 = 0 and let σ0,0 be the initial distribution of X0,0. Then Xt,i has distribution

σ0,0Pτ1 1 Pτt 1 t 1 Pi t and we have made the identification Xs,τs = Xs+1,0, gluing the successive chains together. Let the distribution of Xt,i be denoted by σt,i. By the definition of the chain, σt,i is a distribution with bounded density, supported on K.

Figure 1: Steps of the Dikin Walk. The next point is sampled from a Gaussian distribution with a shape (contours depicted with dashed lines) corresponding to Dikin ellipsoids. These ellipsoids approximate well the local geometry.

4.1 Measuring the Change

Let t denote the L2 norm with respect to the measure µt, defined as f t = R

K f2dµt 1/2

for a measurable function f : K R. Further, let K denote the supremum norm

Efficient Sampling

f K = supx K |f(x)| and let

βt+1 = max { dµt/dµt+1 K , dµt+1/dµt K} . (3)

This ratio provides an upper bound on the point-wise change of the density function. A straightforward way to upper bound βt+1 is by writing

K e st+1(x)dx R

K e st(x)dx sup x K e2|st(x) st+1(x)|

and, hence,

log βt+1 2 st(x) st+1(x) K . (4)

Another way to measure the change in successive distributions is with respect to the L2 norm:

αt+1 = dµt/dµt+1 t+1 . (5)

In contrast to the point-wise change, the ratio αt+1 is more difficult to calculate. In this respect, the following result, which follows from the proof of (Lov asz and Vempala, 2006; Kalai and Vempala, 2006), will be useful:

Lemma 5 Let st be a convex function and st+1 = (1 δ) 1 st. Let µt and µt+1 be defined as in (1). Then

αt+1 1 + δ2

In particular, if δ d 1/2 1/3, then αt+1 5.

We remark that the ratio between µt and µt+1 measured in the supremum norm may be exponentially large, while the L2 change is small. As in (Lov asz and Vempala, 2006; Kalai and Vempala, 2006), this fact will be crucial in this paper when we study simulated annealing.

4.2 Tracking the Distributions: Main Results

Denote the error in approximating the stationary distribution at the end of t-th chain by

dµt 1 t (6)

t r2 t Cν2d .

Theorem 6 The errors ξt satisfy the recurrence

ξt (1 t)τt(β3/2 t ξt 1 + p

βt(βt 1)) (7)

for any t 1.

Narayanan and Rakhlin

Proof [Proof of Theorem 6] We iteratively apply Theorem 2 with f = dσt,j

dµt 1 and the stationary distribution γ = µt, and observe that Ef takes σt,j to σt,j+1. Then from Lemma 4, for t 1 and i 1, dσt,i

dµt 1 t dσt,0

dµt 1 t (1 t)i

Using the first part of Lemma 11 (see Section 6) dσt,0

dµt 1 t β3/2 t

dσt,0 dµt 1 1 t 1 + p

concluding the proof. An alternative recurrence, using the second part of Lemma 11, is

ξt (1 t)τt( p

which is better for large βt but worse for βt 1.

We would like to adaptively choose τt to make the right-hand side (7) small. While the value of the error ξt 1 at the previous round is not available for this purpose, let us maintain an upper bound ut 1 on this error. Thus, we may write τt as a function τt(ut 1, st, rt, βt). Suppose at round t = 0 we ensure that ξ0 u0. Then, recursively, we may compute ut as the upper bound in (7):

ut (1 t)τt(β3/2 t ut 1 + p

βt(βt 1)) (8)

Then, given the initial condition, we have ξt ut for all t 0. Let us consider some consequences of Theorem 6. In particular, we are interested in situations when we can track the distributions with only one step of the random walk.

Corollary 7 Let τt = 1 for all t 1 and suppose ξ0 u0 = β0(β0 1)/ 0 with 0 = 1 Cd3ν2 1

2. Assume that βt is non-decreasing and t is non-increasing in t, and suppose

β3/2 t 1 + 2 t 1 t (9)

for all t 1. Then we have

ξt ut = βt(βt 1)

t for all t 0. In particular, (9) is satisfied whenever βt 1 0.4 2 t .

The proof of the above corollary follows from the more general result:

Corollary 8 Fix a sequence ϵ0, . . . , ϵt, . . . of positive target accuracies and assume ξ0 ϵ0. It is then enough to set

t log β3/2 t ϵt 1

ϵt + βt(βt 1)

in order to ensure ξt ϵt for each t 0.

Efficient Sampling

Proof Immediate by writing

ut = (1 t)τt(β3/2 t ϵt 1 + p

βt(βt 1)) ϵt,

solving for τt, and using the approximation log(1/(1 t)) log(1 + t) t .

We now consider the case when one has control on the L2 norm αt of the change between successive distributions. First, observe that closeness of the distributions in the norm t implies closeness in total variation distance as

Z |dσt,i dµt| = Z dσt,i

dµt 1 dµt dσt,i

dµt 1 t . (11)

Proposition 9 Fix a sequence ϵ0, . . . , ϵt, . . . of positive target accuracies and assume d TV (σ0,0, µ0) ϵ0. Suppose we set

Then the total variation distance between σt,τt and µt is bounded as

d TV (σt,τt, µt)

s=0 ϵs (13)

for each t 0.

Proof For any t 1, let us write

σt,τt = µt + γt (14)

with a signed measure γt = σt,τt µt. By way of induction, suppose (13) holds for time t. Consider the operator Et+1 corresponding to the random walk of the t + 1-st chain. The operator acts on a function f by taking f to R

K f(y)d Pt+1(x, y). Then applying Theorem 2 to the function dµt/dµt+1 1, we have

dµt+1 1 t+1 ϵt+1

by the choice of τt+1 and the definition of αt+1. That is, upon the action of Eτt+1 t+1 , µt is mapped to µt+1 within an error of at most ϵt+1 in the L2 sense (and, hence, in the total variation sense). Since the operator Et+1 is non-expanding in the L1 sense, total variation of γt does not increase under the action of Eτt+1 t+1 . In view of the inductive hypothesis for step t, we conclude d TV (σt+1,τt+1, µt+1) Pt s=0 ϵs + ϵt+1, as desired.

Narayanan and Rakhlin

5. Applications

Before diving into the applications of the random walk, let us give several examples of sets K for which the self-concordant barrier F and its Hessian can be easily calculated. In the following examples, assume that K has non-empty interior.

Example 1 Suppose K is given by m linear constraints of the form aj, x bj, j = 1, . . . , m. Then F(x) = Pm j=1 log(bj aj, x ) is a self-concordant barrier with parameter ν = m. The Hessian is easily computable:

aja T j (bj aj, x )2 .

Example 2 Let K = {x Rd : fj(x) 0, j = 1, . . . , m} where each fj is a convex quadratic form. Then F(x) = Pm j=1 log( fj(x)) is a self-concordant barrier with parameter m. As an example, the function log(R x 2) is a self-concordant barrier for the unit Euclidean sphere {x : x 2 1 0}, with parameter ν = 1, and the Hessian is given by

D2F(x) = 2 1 x 2 I + 4 (1 x 2)2 xx T .

Importantly, there always exists a self-concordant barrier with ν = O(d); yet, for some convex sets (such as the sphere) the parameter can even be constant. Self-concordant barriers can be combined: if Fj is νj-self-concordant for Kj, j = 1, . . . , m, then P j Fj is P j νj-self-concordant for the intersection i Ki, given that it has nonempty interior. Thus, closed forms for the Hessian of the barrier, required for defining Gr x in our Markov chain, can be calculated for many sets K of interest. We refer to (Nemirovskii, 2004; Nesterov and Nemirovskii, 1994) for further powerful methods for constructing the barriers.

5.1 Sampling from Posterior in Exponential Families

Suppose data y1, y2, . . . Y are distributed i.i.d. according to a member of an exponential family with natural parameter x:

p(y|x) = exp{ x, T(y) A(x)}h(y)

where A(x) = R h(y) exp { x, T(y) } is a convex function and T : Y 7 Rd is a sufficient statistic. Suppose x K; that is, we have some knowledge about the support of the parameter. We have in mind the situation where data arrive one at a time and we are interested in sampling from the associated posterior distributions. The likelihood function after seeing y1, . . . , yt is

Efficient Sampling

and, together with a conjugate prior πκ1,κ2(x) exp { x, κ1 κ2A(x)} for some (κ1, κ2) Rd+1, we obtain the posterior distribution at time t

pt(x|y) exp

(t + κ2)A(x)

We apply the sampling technique to this scenario by defining

s0(x) = x, κ1 + κ2A(x), st(x) =

+ (t + κ2)A(x) .

It remains to calculate the number of steps required to track the distributions as additional data arrive one-by-one. Let L be the Lipschitz constant of A(x) over K with respect to Euclidean norm, and let us assume L to be finite. Then Condition 4 is satisfied with r = min n 1 (t+κ2)L, 1

d o . Furthermore, we may set

βt = sup x K exp {2| x, T(yt) A(x)|} ,

a quantity that depends on the observed data. Importantly, we do not need to provide an a priori data-independent bound of this type, which might not be finite. Suppose we would like to maintain a constant level ϵ > 0 of accuracy at each step t. Corollary 8 guarantees this accuracy if each chain is run for

τt = O ν2d max{(t + κ2)2L2, d2} + log(1/ϵ) .

One of the features of this bound is a relatively benign dependence on the dimension d, especially if the geometry of the set K allows the parameter ν = O(1), as in the case of a sphere. On the negative side, the number of steps needed after seeing t data points is proportional to t2. Such an adverse dependence, however, is to be expected as the posterior distribution becomes concentrated very quickly. We now demonstrate that stronger results can be achieved under additional assumptions via Condition 3. Suppose that A is smooth: there exists H 0 such that

A(x) A(w) + A(x), w x + (w x)TH(w x)

for any w, x K. This is a natural assumption, as the second derivative of the log normalization function A corresponds to the variance of the random variable with the given parameter; furthermore, A is differentiable. Let λmax be the largest eigenvalue of H. Then the condition yields rt = C

(t+κ2)λmax . To obtain ϵ-accuracy, it suffices to set

τt = O ν2d max{(t + κ2)λmax, d2} + log(1/ϵ) ,

which has only linear dependence on the size of the data seen so far. We remark that each step of the random walk requires evaluation of the log-partition function A(x). If this function is not available in closed form, we may approximate the value A(x) for each query x. In order to do this, we may run an additional sampling procedure with s (x) = x, T(y) . Alternatively, we may appeal to known methods for this problem, such as Hit-and-Run (Vempala, 2005).

Narayanan and Rakhlin

5.2 Examples of Sampling from Drifting Truncated Distributions

In the previous example, we employed the Markov chain to sample a parameter from a logconcave posterior. We now turn to the question of sampling from a log-concave distribution restricted to a convex set. This problem has a long history (see e.g. (Devroye, 1986; Gilks and Wild, 1992)), and it is recognized that sampling from truncated distributions is difficult even for nice forms such as the Normal distribution. One successful approach to this problem is the Gibbs sampling method (Robert, 1995; Damien and Walker, 2001), yet the rate of convergence is not generally available. The MCMC method of this paper yields a provably fast algorithm for such situations. Furthermore, we can track a drifting distribution over K with a small number of steps. For illustration purposes, we study a simple example of a truncated Normal distribution; the same techniques, however, apply more generally. To simplify calculations, suppose the distributions µt are defined to be N(ct, 1

d I) over a convex compact set K Rd and suppose the mean ct is drifting within a Euclidean ball of radius R. With the definition in (1) we have st(x) = 1

2 x ct 2. Define the drift δt = ct ct 1 . In view of (4),

log βt sup x K ct ct 1 2x ct ct 1 CR,Kδt

where CR,K depends on the radius R and the radius of a smallest Euclidean ball enclosing K. In the same manner, the Lipschitz constant of st(x) over K can be upper bounded by LR,K that depends solely on the two radii. We may thus set the step size to be rt = min{ 1

d, 1 LR,K }. If we aim for a fixed target accuracy ϵ for all t, by Corollary 8, it is enough to make

t log β3/2 t + βt(βt 1)

steps. In the case that the drift δt is small enough, only one step is sufficient. To quantify the regime when this happens, observe that βt exp{CR,Kδ} 1 + Cδt, and it is then enough to require

δt = O 2 t = O

min{1/d2, 1/L2 R,K}

in view of (9). It is quite remarkable that the one-step random walk can track the changing distribution up to the accuracy O δt ν2d

, proportional to the size of the drift. Of course, better accuracy can be achieved by performing more steps, as per Corollary 8. Another related application is to modeling with mixtures of log-concave distributions. Such models have been successful in clustering (Mc Lachlan and Peel, 2000; Walther, 2009), with a mixture of normal distributions being a classical example (Fraley and Raftery, 2002). A mixture of parametric log-concave distributions can be written as Pk i=1 αiπi(θi; x); here αi are positive mixing weights summing to one, and πi are a distributions on K parametrized by θi. A classical method for fitting models to data is the EM algorithm. Given that the parameters {θi}k i=1 and the mixing weights {αi}k i=1 have been estimated from data, one may require random samples from this model for integration or other purposes. Given our procedure for sampling from a single log-concave distribution, one may simply pick the mixture according to the weights αi and then sample from the component. The situation

Efficient Sampling

becomes interesting in the case of online arrival of data, when we need to re-compute the EM solution in light of additional data. By the arguments of (Rakhlin and Caponnetto, 2006; Caponnetto and Rakhlin, 2006), the solution to clustering problems (the analysis was performed for square loss) is stable in the following sense: addition of o( n) new data to a sample of size n is unlikely to drastically move the solution (the argument is based on uniqueness of the maximum of an empirical process). This in turn implies that the parameters {θi} are unlikely to change by a large amount, and we may thus use the method of sampling from a drifting distribution described earlier. We also remark that the method can be easily parallelized since the Markov chains for the k components do not interact.

5.3 Simulated Annealing for Convex Optimization

Let f(x) be a proper convex 1-Lipschitz function. The aim of convex optimization is to find x with the property f( x) minx K f(x) ϵ for a given target accuracy ϵ > 0. We consider the special case of linear function f(x) = ℓ, x , known as Linear Optimization. Complexity of an optimization procedure is often measured in terms of oracle calls queries about the unknown function. A query about the function value is known as the zero-th order information, while a query about a subgradient at a point as the first order information. In the case that the oracle answer is given without noise, it is known that the complexity scales as O (poly(d, log(1/ϵ))). The state-of-the-art result here is the method of (Kalai and Vempala, 2006; Lov asz and Vempala, 2006) which attains the d4.5 dependence on the dimension. We now apply our machinery to obtain a O ν2d3.5 log(1/ϵ) method. In particular, this yields an improved d3.5 dependence on the dimension for the case when K has a favorable geometry: there exists a self-concordant barrier with a parameter ν = O(1). We use the annealing scheme of (Kalai and Vempala, 2006). To this end, we set st = 1 d 1/2 t f and observe that the assumption of Lemma 5 is satisfied with δ = d 1/2. Since functions are linear, we may set the step size rt = 1/d for all t. Hence, αt 5 whenever d > 8 (and a different constant can be obtained for smaller d from the proof). By Proposition 9 with a constant accuracy ϵt = ϵ (

d log(d/ϵ)) 1, by making

steps for t = 1, . . . , k, we guarantee

d TV (σk,τk, µk) kϵ(

d log(d/ϵ)) 1 . (17)

According to (Kalai and Vempala, 2006, Lemma 4.1), if X is chosen from a distribution with density proportional to exp{ T 1 ℓ, x }, with ℓ = 1 and some temperature T > 0, then E( ℓ, X ) min x K ℓ, x d T.

Hence, we take the desired temperature to be T = ϵ/d, and the number of chains that permits the annealing schedule to reach this temperature can be calculated as k =

ϵ ). In view of (17), the final output of the procedure is an ϵ-accurate solution to the optimization problem. The complexity of the method is then O(d3.5ν2 log2(d/ϵ)).

Narayanan and Rakhlin

This result can be extended to Lipschitz convex functions beyond linear optimization. However, the step size condition for convex Lipschitz functions requires the steps to be O(1/ϵ) towards the end of the annealing schedule. This in turn implies only a suboptimal O(ν2d/ϵ2) complexity. It is an open question of whether Dikin Walk can handle such annealing schedules in a more graceful manner.

5.4 Sequential Prediction

Another application of the proposed sampling technique is to the problem of sequential prediction with convex cost functions. Within this setting, the learner (or, the Statistician) is tasked with making a series of predictions while observing a sequence of outcomes on which we place no distributional assumptions. The goal of the learner is to incur cost comparable to that of a fixed strategy chosen in hindsight after observing the data. Initially studied by Hannan (Hannan, 1957), Blackwell (Blackwell, 1956), and Cover (Cover, 1965), the problem of achieving low regret for all sequences has received much attention in the last two decades, and we refer the reader to (Cesa-Bianchi and Lugosi, 2006) for a comprehensive treatment. As we show in this section, a strategy that exponentially down-weighs the decisions with large costs is a good regret-minimization strategy, and this exponential form is amenable to the sampling technique of this paper whenever the costs are convex. More specifically, let K Rd be a convex compact set of decisions of the learner. Let ℓ1, . . . , ℓT be a sequence of unknown cost functions ℓt : K R. On round t, the learner chooses a distribution (or, a mixed strategy) µt 1 supported on K and plays a decision Yt µt 1.2 Nature then reveals the next cost function ℓt. For example, in the well-studied problem of sequential probability assignment, the Statistician predicts the probability xt [0, 1] = K of the next outcome {0, 1} and incurs the cost ℓt(xt) = |xt yt| with respect to the actual outcome yt. A randomized strategy Yt then incurs a cost ℓt(Yt). The goal of the learner is to minimize expected regret

Reg T (U) E

with respect to all randomized strategies defined by p U P, for some collection of distributions P. A procedure that guarantees sublinear growth of regret with respect to any distribution p U P and for any sequence of cost functions ℓ1, . . . , ℓT will be called consistent with respect to P. Let Lt(x) = Pt s=1 ℓs(x) denote the cumulative cost functions, and let η > 0 be a parameter called the learning rate. Fix R(x) to be some convex function that defines the prior, let

st(x) = ηLt(x) + R(x), s0(x) = R(x) (18)

and define the probability distributions µt as in (1). It turns out that this choice of µt is indeed a good regret-minimization strategy, as we show next. The method is similar to the Mixture Forecaster used in the prediction context (Yamanishi, 1998; Vovk, 2001; Azoury and Warmuth, 2001; Kakade and Ng, 2005), and for a discrete set of decisions it is known

2. The index t 1 on µt 1 reflects the fact that Yt is chosen without the knowledge of ℓt.

Efficient Sampling

as the celebrated Exponential Weights Algorithm (Vovk, 1990; Littlestone and Warmuth, 1994). Let D(p||q) stand for the Kullback-Leibler (KL) divergence between distributions p and q.

Lemma 10 For each t 1, let Yt be a random variable with distribution µt 1 as defined in (1). The expected regret with respect to U with distribution p U is

Reg T (U) = η 1 (D(p U||µ0) D(p U||µT )) + η 1 T X

t=1 D(µt 1||µt).

Specializing to the case ℓt : K 7 [0, 1] for all t,

Reg T (U) η 1D(p U||µ0) + Tη/8.

If the KL divergence between the comparator distribution p U and the prior µ0 is bounded for all p U P, the second statement of the lemma yields consistency and a O(

T) rate of regret growth. To bound the divergence between a continuous initial µ0 and a point distribution at some x K, the analysis can be carried out in two stages: comparison to a small-covariance Gaussian centered at x , followed by an observation that the loss of the small-covariance Gaussian strategy is not very different from the loss of the deterministic strategy x . The analysis can be found in (Cesa-Bianchi and Lugosi, 2006, p. 326) and gives a near-optimal O( T log T) regret bound. The easy proof of Lemma 10 appeared in (Narayanan and Rakhlin, 2010) and we include it in Section 6 for completeness. Having exhibited a good prediction strategy, a natural question is whether there exists a computationally efficient algorithm for producing a random draw from a distribution close to the desired mixed strategy µt 1. To this end, we use the sampling method proposed in this paper. As a concrete example, consider linear functions ℓ1, . . . , ℓT and let R 0. For simplicity assume boundedness ℓt : K 7 [0, 1]. In this case, we may choose η = O(1/

βt exp {2η ℓt K} 1 + Cη

for large enough T. Further, we set rt = 1/d according to Condition 1, and the requirement (9) is seen to be satisfied for large enough T. With these choices of the parameters, the sequence of distributions µ1, . . . , µt can be tracked with only one step of a random walk per iteration. The quality of this approximation is O ηd3ν2 at each step. Therefore, regret of the proposed random walk method is within O Tηd3ν2 from the ideal procedure of Lemma 10, as can be seen by writing

|Eℓt(Yt) Eℓt(Xt 1,1)| Z

x K |ℓt(x)| |dσt 1,1(x) dµt 1(x)| Cηd3ν2 .

By choosing η = 1 d3/2ν

Reg T (U) Cd3/2νD(p U||µ0)

Narayanan and Rakhlin

A similar results holds for nonzero R, under the assumption that the L2 distance between dµ0(x) exp{ R(x)}dx and the uniform distribution on K is bounded. We now discuss interesting parallels between the proposed randomized method and the known deterministic optimization-based regret minimization methods. First, the statement of Lemma 10 bears striking similarity to upper bounds on regret in terms of Bregman divergences for the Follow the Regularized Leader and Mirror Descent methods (Rakhlin, 2008; Beck and Teboulle, 2003), (Cesa-Bianchi and Lugosi, 2006, Therem 11.1). Yet, the randomized method operates in the (infinite-dimensional) space of distributions while the deterministic methods work directly with the set K. Second, deterministic methods of online convex optimization face the difficulty of projections back to the set K. This issue does not arise when dealing with distributions, but instead translates into the difficulty of sampling. We find these parallels between sampling and optimization intriguing. Third, a single step of the proposed random walk requires sampling from a Gaussian distribution with covariance given by the Hessian of the self-concordant barrier. This step can be implemented efficiently whenever the Hessian can be computed. The computation time exactly matches (Abernethy et al., 2008, Algorithm 2): it is the same as time spent inverting a Hessian matrix, which is O(d3) or less. Finally, as already mentioned, the idea of following a timevarying distribution is inspired by the method of following the central path in the theory of interior point methods (Nesterov and Nemirovskii, 1994; Nemirovskii, 2004). Similarly to the fast convergence of the chain under the lower bound on conductance, one has fast quadratic local convergence of interior point methods. One may therefore make parallels between conductance and local curvature. A further investigation of these connections is needed, especially in view of the recent developments on positive Ricci curvature of Markov chains (Ollivier, 2009).

Lemma 11 For any t and i 0, it holds that dσt,i

dµt 1 t β3/2 t

dσt,i dµt 1 1 t 1 + p

and, alternatively, dσt,i

dµt 1 t β1/2 t

dσt,i dµt 1 1 t 1 + p

Proof Let us use the shorthand dσ = dσt+1,i and β = βt+1. Using (3), we may write dσ dµt+1 1 t+1 p

β dσ dµt+1 1 t

β dσ dµt+1 1 t dσ dµt 1 t + dσ dµt 1 t

By the triangle inequality, dσ dµt+1 1 t dσ dµt 1 t dσ dµt+1 dσ

Efficient Sampling

For any function f : K R, let f+(x) = max(0, f(x)) and f (x) = min(0, f(x)). In view of (3),

dσ dµt+1 dσ

dσ dµt (β 1)1 1 < dµt dµt+1

1 1 dµt dµt+1

(β 1)2 dσ dµt

Therefore, dσ dµt+1 1 t dσ dµt 1 t (β 1) dσ dµt

t (β 1) 1 + dσ dµt 1 t

The first statement follows by rearranging the terms. Alternatively, we can obtain an inequality that is slightly weaker for β 1 0 and stronger for large β by simply writing dσ dµt+1 1

dµt+1 1 2 dµt+1

dµt+1 1 = Z

dµt dµt+1 dµt 1 .

Using β as an upper bound on the one-sided change dµt/dµt+1 K leads to

dµ2 t dµt 1 = β dσ dµt 1

and subadditivity of the square root function concludes the proof.

Proof [Proof of Theorem 1] Given interior points x, y in int(K), suppose p, q are the ends of the chord in K containing x, y and p, x, y, q lie in that order. Denote the cross ratio by

σ(x, y) = |x y||p q|

|p x||q y|,

and for two sets S1 and S2 let

σ(S1, S2) inf x S1,y S2 σ(x, y).

A result due to Lov asz and Vempala (2007) states the following. If S1 and S2 are measurable subsets of K and µ a probability measure supported on K that possesses a density whose logarithm is concave, then

µ((K \ S1) \ S2) σ(S1, S2)µ(S1)µ(S2).

Narayanan and Rakhlin

This is a non-trivial isoperimetric inequality which says that for any partition of the convex set K into S1, S2 and S3, the volume of S3 is large relative to that of S1 and S2 whenever S1 and S2 are separated. Given this isoperimetric result, to prove the theorem it only remains to show that the σ-distance can be lower bounded (up to a multiplicative constant) by the Riemannian metric ρ. The proof of this fact goes through the Hilbert (projective) metric, which is defined by d H(x, y) ln (1 + σ(x, y)) .

Further, for x K and a vector v, let

|v|x sup x αv K α.

The following two relations between the introduced notions hold. The first one (see Nesterov and Nemirovskii (Nesterov and Nemirovskii, 1994, Theorem 2.3.2 (iii))) is

|h|x h x 2(1 + 3ν)|h|x (20)

for all h Rd and x int(K), where ν is the self-concordance parameter of F. The second relation (see Nesterov and Todd (Nesterov and Todd, 2008, Lemma 3.1)) states that

x y x x y 2 x ρ(x, y) ln(1 x y x). (21)

whenever x y x < 1. For any z on the segment xy an easy computation shows that d H(x, z) + d H(z, y) = d H(x, y). Therefore it suffices to prove the result infinitesimally. From (21), limy x ρ(x,y) x y x = 1, and a direct computation shows that

lim y x d H(x, y)

|x y|x = lim y x σ(x, y) |x y|x 1.

Hence, in view of (20), the Hilbert metric and the Riemannian metric satisfy

ρ(x, y) 2(1 + 3ν)d H(x, y).

Using ln(1 + x) x concludes the proof.

Proof [Proof of Lemma 4] The argument roughly follows the standard path, which is explained, for instance, in (Vempala, 2005). Let S1 be a measurable subset of K such that µ(S1) 1

2 and S2 = K \ S1 be its complement. Fix a C > 1 and let

S 1 = S1 {x Px(S2) 1/C} and S 2 = S2 {y Py(S1) 1/C}.

That is, points in the set S 1 are unlikely to transition to the set S2, and S 2 is analogously unlikely to reach S1 in one step. By the reversibility of the chain, which is easily checked, Z

S1 Px(S2)dµ(x) = Z

S2 Py(S1)dµ(y).

Efficient Sampling

For any x S 1 and y S 2,

d TV (Px, Py) = 1 Z

dµ (w), d Py

dµ (w) dµ(w) 1 1

That is, the transition probabilities for a pair in S 1 and S 2 must be dissimilar. But Lemma 3 implies that if ρ(x, y) r C

d, then d TV (Px, Py) 1 1

C . Therefore

ρ(S 1, S 2) r C

We conclude that the sets S 1 and S 2 must be well-separated. Therefore, the isoperimetric result of Theorem 1 implies that

µ((K \ S 1) \ S 2) ρ(S 1, S 2) 2(1 + 3ν) min(µ(S 1), µ(S 2)) r Cν

d min(µ(S 1), µ(S 2)).

First suppose µ(S 1) (1 1

C )µ(S1) and µ(S 2) (1 1

C )µ(S2). Then, Z

S1 Px(S2)dµ(x) = 1

S1 Px(S2)dµ(x) + 1

S2 Px(S1)dµ(x)

2C µ((K \ S 1) \ S 2)

d min(µ(S 1), µ(S 2))

d min(µ(S1), µ(S2)),

proving the result. Otherwise, without loss of generality, suppose µ(S 1) (1 1

C )µ(S1). Then Z

S1 Px(S2)dµ(x) = 1

S1 Px(S2)dµ(x) + 1

S2 Px(S1)dµ(x)

S1\S 1 Px(S2)dµ(x) µ(S1)

concluding the proof.

Proof [Proof of Lemma 5] The proof closely follows that in (Kalai and Vempala, 2006). By definition,

dµt/dµt+1 2 t+1 = Z

2 dµt+1 = Z

dµ2 t dµt+1 = Z

Z2 t Zt+1 exp{ st+1} .

Writing out the normalization terms,

dµt/dµt+1 2 t+1 =

K exp{ st+1} R

K exp{st+1 2st} R

K exp{ st} 2 = Y (1)Y ( 1 + 2(1 δ))

Y (1 δ)Y (1 δ)

Narayanan and Rakhlin

where Y (a) = R

K exp{ ast+1}. As shown in (Kalai and Vempala, 2006, Lemma 3.1), the function ad Y (a) is log-concave in a, and thus

Applying this inequality with a = 1 and b = 1 + 2(1 δ),

dµt/dµt+1 2 t+1 1 + δ2

In particular, if δ d 1/2 1/3 (that is, d > 8), we obtain an upper bound of exp n d d 2

Proof [Proof of Lemma 10] Observe that D(µt 1||µt) can be written as Z

K dµt 1 log qt 1Zt

Zt 1qt = log Zt

K ηℓt(x)dµt 1(x) = log Zt

Zt 1 + ηEℓt(Yt). (22)

Rearranging, canceling the telescoping terms, and using the fact that Z0 = 1

t=1 ℓt(Yt) =

t=1 D(µt 1||µt) log ZT .

Let U be a random variable with a probability distribution p U. Then

t=1 Eℓt(U) = η 1 Z

K ηLT (u)dp U(u) = η 1 Z

K dp U(u) log q T (u)

K dp U(u) log q T (u)/ZT

q0(u) + η 1 T X

t=1 D(µt 1||µt)

= η 1 (D(p U||µ0) D(p U||µT )) + η 1 T X

t=1 D(µt 1||µt).

Now, from Eq. (22), the KL divergence can be also written as

D(µt 1||µt) = log

K e ηℓt(x)qt 1(x)dx R

K qt 1(x)dx + ηEℓt(Yt) = log Ee η(ℓt(Yt) Eℓt(Yt))

By representing the divergence in this form, one can obtain upper bounds via known methods, such as log-Sobolev inequalities (e.g. (Boucheron et al., 2003)). In the simplest case of bounded loss, it is easy to show that D(µt 1||µt) O(η2), and the particular constant 1/8 can be obtained by, for instance, applying Lemma A.1 in (Cesa-Bianchi and Lugosi, 2006). This proves the second part of the lemma.

Efficient Sampling

7. Smooth Variation of the Transition Kernel

In this section, we study the transition x y. For this purpose, it is enough to assume that x is the origin and that the Dikin ellipsoid at x is a unit Euclidean ball. This can be achieved by an affine transformation, leading to no loss of generality since the resulting statement about measures on K is invariant with respect to affine transformations. Hence, in what follows, for the particular x we have < , >x=< , > and x = . Since x is the origin, we have E z 2 x = r2 for z sampled from Gr x. Further, without loss of generality, we may also assume s(x) = 0. Proof [Proof of Lemma 3] In view of the first inequality in Eq. (21),

x y x x y 2 x ρ(x, y) r C

Without loss of generality, assume r C

d 1 8. First, we claim that x y x must be small. For the sake of contradiction, suppose x y x > 1/2 and consider a point y

with x y x = 1/2 and lying on the geodesic path between x and y with respect to the Riemannian metric. Clearly, ρ(x, y ) r C

8, yet by Eq. (21) we have 1

4 ρ(x, y ),

contradicting our assumption. Hence, x y x 1/2, and, therefore, x y x 2r C

d. It remains to show that if x, y K and

then d TV (Px, Py) = 1 1

By definition, we have that

1 d TV (Px, Py) = Ez

min 1, Gr y(z) Grx(z), Gr z(x) exp(s(x)) Grx(z) exp(s(z)), Gr z(y) exp(s(y)) Grx(z) exp(s(z))

where the expectation is taken over a random point z having density Gr x. Thus, it suffices to prove that for some C > 1

P min Gr y(z) Grx(z), Gr z(x) exp(s(x)) Grx(z) exp(s(z)), Gr z(y) exp(s(y)) Grx(z) exp(s(z))

By our assumption, x is the origin and D2F(x) = I, the latter implying that V (x) = 0. Thus, Gr y(z) Grx(z) = exp

d y z 2 y r2 + V (y) + d z 2

Gr z(x) exp(s(x)) Grx(z) exp(s(z)) = exp d z 2 z r2 + V (z) + d z 2

r2 + (s(x) s(z)) ,

and Gr z(y) exp(s(y)) Grx(z) exp(s(z)) = exp d y z 2 z r2 + V (z) + d z 2

r2 + (s(y) s(z)) .

Narayanan and Rakhlin

Thus, it remains to prove that there exists a constant C such that

P h max n d y z 2 y r2V (y), d z 2 z + r2(s(z) s(x)) r2V (z),

d z y 2 z + r2(s(z) s(y)) r2V (z) o < d z 2 + r2C i 1

This fact is shown in technical Lemmas 13 and 14 below.

In proving the technical lemmas, we will use the fact that x y x 2r C

d as shown above, and that x z x (for z sampled from Gr x) is likely to be bounded above by a multiple of r by straightforward concentration arguments.

Lemma 12 There exists a constant C > 0 such that

P [max ( V (y), V (z)) < C] > 0.9

Proof Fix a constant c. First, notice that over a Euclidean ball of radius c/d around the origin, the Hessians D2F(u) are lower-bounded by a factor of (1 c/d)2 from the Hessian at the origin (the identity) by (24). Hence, the determinant function can decrease from 1 by at most a constant factor. Thus V (u) < C for some constant C for any u with x u x c/d. Now recall that y is deterministically within the 1/d ball, while z is in the ball of radius c/d with high probability.

Lemma 13 Under step size Condition 4, for any

P h max n s(z) s(x), s(z) s(y) o < C i > 0.32.

Proof Since with large enough probability x y x < C r and x z x < C r, we also have z y x < 2C r. Then, by (24), the norms at z and x are within a multiplicative constant, and thus the pairs (z, x) and (z, y) are subject to the step size choice specified in the condition. That is, there exists a g such that

s(z) s(x) = s(z) s(x) g, z x + g, z x C + g, z x

and similarly

s(z) s(y) = s(z) s(y) g, z y + g, z y C + g, z y

Then, assuming (without loss of generality) x = 0,

P [max { g, z x , g, z y } < 0] = P [ g, z min {0, g, y }] .

Observe that g, z is a Gaussian random variable whose standard deviation is larger than g y . Therefore,

P [ g, z min {0, g, y }] erfc 1/

where erfc(x) 2 π R x e t2dt is the usual complementary error function.

The following probabilistic upper bound completes the proof.

Efficient Sampling

Lemma 14 There exists a constant C > 0 such that

P max n y z 2 y, z 2 z, z y 2 z o z 2 < Cr2

Proof [Proof of Lemma 14] Since y < Cr

d, y y and y z are less than Cr

d. So it suffices to show that

P max n z 2 y z 2, z 2 z z 2, y, z y , y, z z o < Cr2

We proceed to do so by proving probabilistic upper bounds on each of the terms

(a) z 2 y z 2 , (b) z 2 z z 2 , (c) y, z y , and (d) y, z z

separately, and finally applying the union bound. We first prove an upper bound on (a) and (b). Note that r 1

d and thus r3 r2

d . It suffices to observe that by (24)

whenever z < 1/2. Similarly, for y < 1/2,

There exists a constant C such that the quantity z 3 is bounded by Cr3 with probability at least 0.99. We now turn to bounding (c) and (d). Let [0, u] denote the line segment between the origin and u. By the mean-value theorem,

y, z y = y, z + ( y, z y y, z ) y, z + sup y [0,y] D3F(y )[y, y, z]

y, z z = y, z + ( y, z z y, z ) y, z + sup z [0,z] D3F(z )[y, z, z]

Observe that

d with probability at least 0.99 by a measure-concentration argument. Indeed, most of the vectors z are almost perpendicular to the given vector y. Now, using (23),

sup y [0,y] D3F(y )[y, y, z] sup y [0,y] 2 y 2 y z y Cr2

sup z [0,z] D3F(z )[y, z, z] sup z [0,z] 2 y z z 2 z Cr3

Narayanan and Rakhlin

with probability at least 0.99. Therefore, there exists a constant C > 0 such that

P y, z y < Cr2

and the same statement holds for y, z z. We also have that

d + sup z [0,z] 2 y z z 2 z Cr2

P y, z z < Cr2

8. Self-concordant barriers

Let K be a convex subset of Rd that is not contained in any (d 1)-dimensional affine subspace and int(K) denote its interior. Following Nesterov and Nemirovskii, we call a real-valued function F : int(K) R, a regular self-concordant barrier if it satisfies the conditions stated below. For convenience, if x int(K), we define F(x) = .

1. (Convex, Smooth) F is a convex thrice continuously differentiable function on int(K).

2. (Barrier) For every sequence of points {xi} int(K) converging to a point x int(K), limi f(xi) = .

3. (Differential Inequalities) For all h Rd and all x int(K), the following inequalities hold.

(a) D2F(x)[h, h] is 2-Lipschitz continuous with respect to the local norm, which is equivalent to D3F(x)[h, h, h] 2(D2F(x)[h, h]) 3 2 .

(b) F(x) is ν-Lipschitz continuous with respect to the local norm defined by F,

|DF(x)[h]|2 νD2F(x)[h, h].

We call the smallest positive integer ν for which this holds, the self-concordance parameter of the barrier.

The following results can be found, for instance, in (Nesterov and Nemirovskii, 1994; Nemirovskii, 2004; Nemirovski and Todd, 2008). First,

|D3F(x)[h1, . . . , hk]| 2 h1 x h2 x h3 x . (23)

Second, if δ = h x < 1, then

(1 δ)2D2F(x) D2F(x + h) (1 δ) 2D2F(x) . (24)

Efficient Sampling

J. Abernethy, E. Hazan, and A. Rakhlin. Competing in the dark: An efficient algorithm for bandit linear optimization. In Proceedings of The Twenty First Annual Conference on Learning Theory, 2008.

J. D. Abernethy and E. Hazan. Faster convex optimization: Simulated annealing with an efficient universal barrier. In Proceedings of the 33nd International Conference on Machine Learning, ICML 2016, New York City, NY, USA, June 19-24, 2016, pages 2520 2528, 2016.

K. S. Azoury and M. K. Warmuth. Relative loss bounds for on-line density estimation with the exponential family of distributions. Machine Learning, 43(3):211 246, June 2001.

A. Beck and M. Teboulle. Mirror descent and nonlinear projected subgradient methods for convex optimization. Oper. Res. Lett., 31(3):167 175, 2003.

D. Blackwell. An analog of the minimax theorem for vector payoffs. Pac. J. Math., 6:1 8, 1956.

S. Boucheron, G. Lugosi, and P. Massart. Concentration inequalities using the entropy method. Annals of Probability, 31:1583 1614, 2003.

A. Caponnetto and A. Rakhlin. Stability properties of empirical risk minimization over Donsker classes. Journal of Machine Learning Research, 6:2565 2583, 2006.

N. Cesa-Bianchi and G. Lugosi. Prediction, Learning, and Games. Cambridge University Press, 2006.

N. Chopin. A sequential particle filter method for static models. Biometrika, 89(3):539 552, 2002.

T. Cover. Behaviour of sequential predictors of binary sequences. In Proc. 4th Prague Conf. Inform. Theory, Statistical Decision Functions, Random Processes, 1965.

P. Damien and S. Walker. Sampling truncated normal, beta, and gamma densities. Journal of Computational and Graphical Statistics, 10(2), 2001.

L. Devroye. Non-uniform random variate generation (1986). Springer Verlag, 1986.

P. Diaconis. The markov chain monte carlo revolution. Bulletin of the American Mathematical Society, 46(2):179 205, 2009.

P. Diaconis. Some things we ve learned (about Markov chain Monte Carlo). Bernoulli, 19 (4):1294 1305, 2013.

A. Doucet, N. De Freitas, N. Gordon, et al. Sequential Monte Carlo methods in practice. Springer New York, 2001.

M. Dyer, A. Frieze, and R. Kannan. A random polynomial-time algorithm for approximating the volume of convex bodies. Journal of the ACM (JACM), 38(1):1 17, 1991.

Narayanan and Rakhlin

C. Fraley and A. Raftery. Model-based clustering, discriminant analysis, and density estimation. Journal of the American Statistical Association, 97(458):611 631, 2002.

A. Frieze, R. Kannan, and N. Polson. Sampling from log-concave distributions. The Annals of Applied Probability, pages 812 837, 1994.

W. Gilks and P. Wild. Adaptive rejection sampling for gibbs sampling. Applied Statistics, pages 337 348, 1992.

J. Hannan. Approximation to Bayes risk in repeated play. Contributions to the Theory of Games, 3:97 139, 1957.

S. Kakade and A. Ng. Online bounds for Bayesian algorithms. In Proceedings of Neural Information Processing Systems (NIPS 17), 2005.

A.T. Kalai and S. Vempala. Simulated annealing for convex optimization. Mathematics of Operations Research, 31(2):253 266, 2006.

R. Kannan and H. Narayanan. Random walks on polytopes and an affine interior point method for linear programming. Mathematics of Operations Research, 37(1):1 20, 2012.

N. Littlestone and M. K. Warmuth. The weighted majority algorithm. Information and Computation, 108(2):212 261, 1994.

L. Lov asz. Hit-and-run mixes fast. Mathematical Programming, 86(3):443 461, 1999. ISSN 0025-5610.

L. Lov asz and M. Simonovits. Random walks in a convex body and an improved volume algorithm. Random Structures and Algorithms, 4(4):359 412, 1993.

L. Lov asz and S. Vempala. Simulated annealing in convex bodies and an o (n4) volume algorithm. J. Comput. Syst. Sci., 72(2):392 417, 2006.

L. Lov asz and S. Vempala. The geometry of logconcave functions and sampling algorithms. Random Struct. Algorithms, 30(3):307 358, 2007.

GJ Mc Lachlan and D Peel. Finite mixture models. 2000.

S. Meyn and R. L. Tweedie. Markov chains and stochastic stability. Cambridge University Press, 2009.

H. Narayanan. Randomized interior point methods for sampling and optimization. The Annals of Applied Probability, 26(1):597 641, 2016.

H. Narayanan and A. Rakhlin. Random walk approach to regret minimization. In Advances in Neural Information Processing Systems, 2010.

A. S. Nemirovski and M. J. Todd. Interior-point methods for optimization. Acta Numerica, pages 191 234, 2008.

A. S. Nemirovskii. Interior point polynomial time methods in convex programming, 2004.

Efficient Sampling

Y. E. Nesterov and A. S. Nemirovskii. Interior Point Polynomial Algorithms in Convex Programming. SIAM, Philadelphia, 1994.

Y.E. Nesterov and M. J. Todd. On the Riemannian geometry defined by self-concordant barriers and interior-point methods. Foundations of Computational Mathematics, 2(4): 333 361, 2008.

Y. Ollivier. Ricci curvature of markov chains on metric spaces. Journal of Functional Analysis, 256(3):810 864, 2009.

A. Rakhlin. Lecture notes on online learning, 2008.

http://stat.wharton.upenn.edu/~rakhlin/papers/online learning.pdf.

A. Rakhlin and A. Caponnetto. Stability of K-means clustering. In Advances in Neural Information Processing Systems 19, pages 1121 1128. MIT Press, 2006.

C. P. Robert. Simulation of truncated normal variables. Statistics and computing, 5(2): 121 125, 1995.

C. P. Robert and G. Casella. Monte Carlo statistical methods, volume 319. 2004.

S. Vempala. Geometric random walks: A survey. In Combinatorial and computational geometry. Math. Sci. Res. Inst. Publ, 52:577 616, 2005.

V. Vovk. Aggregating strategies. In Proceedings of the Third Annual Workshop on Computational Learning Theory, pages 372 383. Morgan Kaufmann, 1990.

V. Vovk. Competitive on-line statistics. International Statistical Review, 69:213 248, 2001.

G. Walther. Inference and modeling with log-concave distributions. Statistical Science, pages 319 327, 2009.

K. Yamanishi. Minimax relative loss analysis for sequential prediction algorithms using parametric hypotheses. In COLT 98, pages 32 43, New York, NY, USA, 1998. ACM.