Sampling from Structured Log-Concave Distributions via a Soft-Threshold Dikin Walk
Authors: Oren Mangoubi
Nisheeth K Vishnoi
Venue: NeurIPS 2023 poster (2023)
Keywords: Logconcave sampling
Dikin walk
Markov chain Monte Carlo
interior point methods

Abstract: Given a Lipschitz or smooth convex function $f:K \to \mathbb{R}^d$ for a bounded polytope $K:=${ $\theta \in \mathbb{R}^d: A\theta \leq b$}, where $A\in  \mathbb{R}^{m\times d}$ and $b \in \mathbb{R}^m$, we consider the problem of sampling from the log-concave distribution $\pi(\theta) \propto e^{-f(\theta)}$ constrained to $K$. Interest in this problem derives from its applications to  Bayesian inference and differential privacy. We present  a generalization of the Dikin walk  to this setting that requires at most $O((md + d L^2 R^2) \times md^{\omega-1} \log(\frac{w}{\delta}))$ arithmetic operations to sample from $\pi$ within error $\delta>0$ in the total variation distance from a $w$-warm start. Here $L$ is the Lipschitz constant of $f$, $K$ is contained in a ball of radius $R$ and contains a ball of smaller radius $r$, and $\omega \approx 2.37$ is the matrix-multiplication constant. This improves on the running time of prior works for a range of structured settings important for the aforementioned inference and privacy applications. Technically, we depart from previous Dikin walks by adding a soft-threshold regularizer derived from the Lipschitz or smoothness properties of $f$  to a barrier function for $K$ that allows our version of the Dikin walk to propose updates  that have a high Metropolis acceptance ratio for $f$, while at the same time remaining inside the polytope $K$.

Sampling from Structured Log-Concave Distributions via a Soft-Threshold Dikin Walk
Oren Mangoubi Worcester Polytechnic Institute
Nisheeth K. Vishnoi Yale University
Abstract
Given a Lipschitz or smooth convex function f Rd for a bounded polytope K {0 E R2 AO b}, where A Ruxd and 6 E g" we consider the problem of sampling from the log-concave distribution (0) & e- constrained to K. Interest in this problem derives from its applications to Bayesian inference and differential privacy. We present generalization of the Dikin walk to this setting that requires at most ((Im + RL² x md log( / )) arithmetic operations to sample from i within error 5 > 0 in the total variation distance from w- w-warm start. Here I. is the Lipschitz constant of f, K is contained in a ball of radius R and contains ball of smaller radius and ss 2.37 is the matrix-multiplicatinn constant. This improves on the running time of prior works for a range of structured setsings important for the aforementioned inference and privacy applications. Technically, we depart from previous Dikin walks by adding a soft-thresholl regularizer derived from the Lipschitz or smoothness prop- erties of f to a barrier function for K that allows our version of the Dikin walk to propose updates that have high Metropolis acceptance ratio for f, while at the same time remaining inside the polytope K.
Introduction
We consider the problem of sampling from a log-concave distribution supported on a poly- tope: Given polytope K = {0 R AA < b}, where A Rmxd and E Rm and convex function f K + R, output sample 0 E K from the distribution T00) e- f(0), Our in- terest in this problem derives from its applications to Bayesian inference and differentially private optimization. In Bayesian inference, the ability to sample from (00) -(()) allows one to compute Bayesian confidence intervals and other statistics for the Bayesian posterior distribution of many machine lcarning models (scc c.g. [7
20 In differentially private optimization, sampling from the "exponential mechanism" 36 allows one to get optimal utility bounds for the problem of minimizing under e-differential privacy
The instances of the volytope-constraie are more general than the two well-studied special cases- the uniform density case (f 0) and the unconstrained case (K Rd). However, they still have more structure than the case of general log-coneave function supported on an arbitrary convex body. For instance, in Bayesian Lasso logistic regression, f(0) > ((0;;;), where & is the logistic loss and are datapoints with ||xح||2 < 1, and K= (0 E Rd ||0|| O(1)); see 44 the logistic function both 0(1)-smooth and Lipshhi,, I is both Oin)-hischhi and ((n)shmooth, and K is defined by O(d) inequalities and contained in ball of radius O(1).
sampling problem that arise in these applications
43
26
47
Since
To obtain an s-differentially private mechanism for the Lasso logistic regression problem, = using the exponential mechanism, the goal is to sample from exp(- n. Σ ((0;x;)) where l is the logistic loss and is contained in ball of radius R. Thus, the log-density is both
37th Conference on Neural Information Processing Systems (NeurIPS 2023).

---

8-smooth and L-Lipschitz for 3 = L = ne/n= = (dd) n n d and E < 1, since R = O1). Another example is a result of [28] that reduces the problem of e-differentially private low- rank approximation of a symmetric p x p matrix to a constrained sampling problem where f linear (and thus 0-smooth) of dimension d d==p p² and K is the Gelfand-Tsetlee polytope (a generalization of the probability simplex). Here K has d incqualities and diameter O(vd).
Importantly, when sampling from the exponential mechanism in privacy applications, sam- pling with total variation (TV) bounds-thi the case that has received the most attention- insufficient to guarantee s-differential privacy, the strongest notion of differential privacy; see 14. Instead, one requires bounds in the stronger infinity-distance metric d (V, SUPaEK [log((()/=(0)))- A recent work 35 showed how to convert samples within O(8)-TV distance from continuous log-Lipscti Lig-hiisc densities n, into samples with D(e)-infnity-distanee bounds, but requires the TV distance 5 to be very small-roghlyy 5 = O(Ee (isse-d-AA, RR), rais- ing the question of designing Markov chains whose runtime bounds also have a low-order dependence on log l/.. Thus, for the aforementioned applications to Bayesian inference and privacy which give rise to structured instances of sampling from log-concave densities over polytopes, is desirable to design sampling algorithms that have a low-order polynomial dependence not only on the parameters d, T., R. B3 but also on log 1/8.
Main related works. A line of work has developed algorithms for sampling in the general sctting when K is an arbitrary convex body given by membership oracle T T 15 [3 use the "hit-and-run" framework to give an algorithm to sample from a log- distribution x e e-f on a convex body K which also contains a ball of radius r with TV-error 8 > O in (dd))))) logg(wdd/((p) (wdR((r5)) log" (w/o)) + [k)) arithmetic operations from w-warm start (their Theorem 1.1) and Old (n) (6r) + fx)) arithmetic operations from a cold start (Corollary 1.2). Here Ty is the time required to evaluate f and Tx is the time for a membership oracle query. Here, a distribution v is w-warm for w 2 1 w.r.t. the stationary distribution ı if supzek "(*)/T(s) S w.
33
3el.
38 give "Dikin-walk -based algorithm to sample from any log-concave a e- on K where f is T--Lisschizz or 8-smooth. The Dikin walk Markov chain was in- troduced in 24 in the special case where f O (see also (7)) Their runtime is O((d - T3, ST2RR R2) log("/s)(T) + md)-1)) arithmetic operations, where w. = 2.37 is the matrix-multiplication constant. From a cold start, their runtime is (( + d LRRRR (dlog((/))+ + M log()/)))(T + +md-1), where M := log(maX@eK (()))) Their bounds when f is 8-smooth are the same, but with each I term replaced with B. We discuss additional related work in Appendix A
Our contributions. We prosent Markov chain sampling algorithm that generates sam- ples from an L-log-Lipschitz or 8-log-smooth log-concave distribution π oc e on an R- bounded polytope K given by m incqualities, with an error bounded in the TV distance (Al- gorithm I and Theorem 2.1 Our algorithm requires O (mmd+dL1 R22) log( ( arithmetic operations to sample with TV error (65) from π x e-f when f is L-Liipscitz, and O(md + d2R2) log( (TT) mdd-1) arithmetic operations in the setting where f is 8-smooth, where Ty is the number of arithmetic operations to compute the value of f and md=-1 is the number of arithmetic operations to compute the Hessian of the log-barrier of K.
In comparison to 33 we improve dependence of the running time on the parameters R/r and log !/s. In comparison to BS,, , we improve the dependence on d while retaining the same dependence on 1/8. Our result directly implies faster runtimes with improved dependence on the dimension d for structured inference problems such as Bayesian Lasso logistic regression where c.g. R/r ((Vd) (Corollary B.1 Moreover, plugging our algorithm into the TV-to- infinity distance bound converter of 155 we obtain an algorithm to sample from a logconcave density constrained to polytope with error bounded by infinity distance that improves upon prior work of 2 see Corollary C.1 and the subsequent discussion. Corollary C.1 along with the exponential mechanism 36 allows to obtain faster runtime bounds with improved
2

---

dependence on d for applications to differentially private empirical risk minimization 2 and matrix approximation 28 see Corollary D.I and the subsequent discussion.
Technically, our algorithm is a Markov chain inspired by the Dikin walk [2 whose steps are determined by barrier function that generalizes the log-barrier function by adding "soft-threshold" f2-norm regularizer. The regularized barrier allows our Markov chain to take larger steps, while still retaining high acceptance probability on Lipschitz or smooth log-densities f allowing our Markov chain to sample from these distributions with faster runtime. key technical stop in obtaining our results is show that our self-concordant barrier function is the limit of an infinite sequence of log- reg-barrier functions for K "padded" with additional redundant inequalities. This allows us to leverage well-known properties of the log-barrier to bound the acceptance probability and mixing time of our Markov chain.
While la regularization is optimal for classes of functions f which are Lipschitz smooth in the lz-norm, for other classes of functions (e.g., Lipschitz in the (-norm), &2 regularization may not be optimal. Moreover, it remains open obtain runtime bounds for the Dikin walk that do not require f to be Lipschitz or smooth, and /or depend polynomially on log R. This leads to the question of whether one can design other tractable self-concordant barriers to obtain further runtime improvements for sampling log-concave distributions on a polytope; we discuss this in Appendix F
2 Results
Our main result (Theorem 2.1 and Algorithm 1 is Markov chain algorithm that generates TV-error bounded samples from L-log- Lipschits or 3-log-smooth log-concave distributions on a polytope. As explained later in this section, Theorem 2.1 often results in the fastest known algorithm for some of the applications to Bayesian inference and differentially private optimization mentioned in the introduction.
Notation. the following, I is L-Lipschitz or 3-smooth for some L, 3 > O. Tr denotes the number of arithmetic operalions to evaluate J, Tk is the number of arithmetic operations to compute the Hessian of the log-barrier function, Tx the operations for membership oracle query, and Tk the operations to compute a projection oracle for K. When K a polytope K (0 E Rd AO < b} given by A E Rmxdd and b E RM one has Tk 0(md==1), Tk O(md), and Tx Ommd"1)). For every j E [1....,m), we denote the j'th row of A by ay and the j'th entry of b by by. When comparing runtimes, we often assume for simplicity that Tj = O(f)) unless otherwise stated, which is the case, e.g., in logistic regression with n = (dd) datapoints. For any two distributions μl, u on Rd, we denote their total variation distance by || p|/TV == supsCR" (p(S) v(8)). For @ E Idd t> O, denote the ball of radius at f by B(0, t) == {2 E Rd ||2 -12--122- S t) where || 1|2 is the Euclidean norm. Denole the interior of any C Rd by [nt(() == (0 ( B(0, C S for some 1 0}. > 0 is step-size hyperparameter shared by our algorithm and the original Dikin walk of 24 37 y > 0 is the hyperparameter in 33 and 17 > a hyperparameter for our algorithm's regularizer.
Theorem 2.1 (Sampling with TV bounds via a soft-threshold Dikin Walk)
There erists algorithm (Algorithm D which, given 8, R 0 and either L > 0 or 3 0, A E Igaxd, b E Rm that define polytope K = (0 E Rd A0 < b} such that K is contained in a ball of radius R and has nonemply interior, an oracle for the value of a conver function f: K + Rdd, where I is eilher L-Lápschitz or B-smooth, and an inélial point sampled from distribution supported on K which is w- wa-arrm with respect to π. x e- for some w > O, outputs point from distribution μ where || T||TV 8. Moreover, this algorithm takes 0(mnd + dL RL2 (og(//)) x (Ty + +K) arithmetic operations in the setting where f is L--Lipschitz, or O((md= + d3R²) log(w/s)) (T, + +T) arithmetic operations when f 3-smooth, where Ty the number of operations to evaluate f and Tr- 0(md--1)
Theorem 2.1 improves on the previous bound of 33 Theorem 1.1) of log (wdt/(()) (w) arithmetic operations, for sampling from x e
3

---

with TV error (()) from a w- waraarm start, by a factor of min min(d2-<(//)),, ( d(-2(())))) log (w/6) when f is L-Lipschitz on a polytope K defined by m = O(d) inequalities (see Table 1). Thus, Tk = O(md) and if, e.g., we also have LR = O(v)) and B/r = vd, the im- When f is 3-smooth, the improvement is
provement is d arithmetic operations P min(d-- ((/)), o' (/(r))) log" (*%o) arithmetic operations.
When a warm start is not provided, Algorithm takes at most O((md+dL3 R²) )²)xx (ogoo( M + Log( (oo)()))) + Tx) arithmetic operations when f is L-Lipschitz (or O((md + dBR²) x (d R/r) + M + (T) + Tkk arithmetic operations when f is f-smooth), since an dlog("/-)+ -warm start can be obtained by sampling uniformly from the ball B(a,r) S K, where M = log(mnxxx e(())-1(0)) < LR. In comparison, the work of [3] for the hit-and- run algorithm gives bound of ((d))) log" to sample with TV error (6) without a warm start. Thus, when we are not given a warm start and क xx e- is constrained to a polytope K defined by m m= O(d) incqualities, Theorem 2.1 improves on the bounds of [3 for the hit-and-run algorithm by a factor of min(d (R/+)² " (xL=)) log (1/6) arithmetic operations in the setting when f is L- Lipschitz, and min(d (B./))2 d² ((((()))) log' when f is 3-smooth. On the other hand, we note that the results of [8 apply more generally when и is a log-concave distribution on convex body, while our bounds for the soft oftttthreshol Dikin walk (Theorem 2.1 apply to the setting where π is a log-Lipschitz (or log-smooth) log-concave distribution on a polytope. In the example of Bayesian Lasso logistic regression, our algorithm takes O(d++w [og@+wlog(///))) arithmetic operations without a warm start since f is both 3-smooth and I--iisscitz, with
+ Tk) arithmetic operations
== = n = m d, and R- 1 and r 3/va (Corollary B.1 in Appendix B this (ogg((/s) arithmelic operations for the hil-and-run algorithm of [3] This is because, while it only takes fx= (d) arithmetic operations to establish membership in the 4 ball, computing the function f in logistic regression requires T = nd = d² arithmetic operations.
improves by a factor of dZ log (4/5) arithmetic operations on the bound of d8 log (4/6)
Arithmetic
Biraaiinns
Ierations if
Algorishm
Iocrations L-Lipschitz J
Deraaionns f-amooth /
openation per Aleration
O(d), I,r - 0(1.
Oldi, Or Old), R= OI
R O(d)
Proximal Langevia MC
((()6
O(vabr-1 1) Tes
6(86-1)
0(855
Hit- ani-ru
O(d((1)))
O(md)
oid")
6(d)
[32]
Dikin Walk of [37]
old xPaxe)
O' 19
ō(d))
Soft Thresheld Dikdn
O(vid d12
Olend daRe)
O'nd=-1)-T,
O(P))
O(F))
Walk (this paper)
Table
Number itemations (and arithmetic operations per-iteration) of different algorishms which imply bounds for sampling within TV emor O(f) from logooncave c- on polytope K when I is L-Lipcchiiz f-smooth, from w--wrm start. Ty and Tey are, respoctively, the mumber of operations to compute / V. The Ó notation hides logarillmic faciors of d. 6.r, R, The titth column gives runtimes when K is a polytope with = O(d) that contains the unit ball (and thus not well-rounddd), and f O(1) Lippshii, The sixth column corresponds to sampling Bayesian Lasso logistic regression posterior distribution with O(d) datapoints, where the unit ball.
We also compare to the bound of O((d' + dd R²) log("/)))TT + Tx) arithmetic operations to sample from an TV error O()) from a w-warm start of 38 dependence on d (and L) than the bounds for the hit-and-run algorithm, the dependence on log 1/s is much smaller- by a factor of log' 1/8- which can allow for faster runtimes when the goal is to sample with infinity-distacce error. We explain this in detail below. Theorem 2.1 improves on the bound of [38] bound by factor of d² arithmetic operations if m 0 (d),
density on a polytope defined by m incqualitics with We note that while this bound has larger
AAy convex body K contained in a ball of radius R > 0 and containing a ball of smaller radius > 0 satisfies the lower bound f/-- 2 1. Moreover, for any convex body K, there always exists a linear transformation T' for which TK is contained in a ball of radius At and contains a ball of radius f such that R/+ S (ddd) (T is referred to as a "well-rounded" convex body).
Whee setting their scalar step size hyperparameter a²/² to (min(1/a 1/r.))) [S get bound of Φ 2 a f'(2/( on the conductance of their Dikin walk (their Lemma 4). Here is the self- concordance parameter of the barrier function; for the log-barrier λ m, although there are other barriers for which k = Od). Plugging their conductance bound into Corollary 1.5 of []] implics a bound of o (olog((/5))) (I)( +i) = O((d' + dL L22 H2) log( (T) + Tk)) arithmetic operations from w-warm start for their walk to sample with O(8) TV error from n. See Appendix A for details.
4

---

and dd if we also have that LR = O(vd). When f is instead f-smooth their bound is O((d + +d log (T, + Te); and our improvement on this bound is d2 arithmetic operations if m = O(d), and d8 if we also have that 31² Odd).
S give faster bounds for the Dikin walk and its variants than the bounds in B8 these only apply in the special case when is the uniform distribution on K. The proof of Theorem 2.1 appears in Appendix - We give
We also note that while many works, e.g. 24 37
an overview of the main ideas in the proof of Theorem 2.1 in Section 3 In Appendix F we
give an axiomatic approach to arrive at our barrier function and discuss possible extensions.
Infinity-distance sampling. In applications of sampling to differentially private opti- 18 23,, bounds in the total variation (TV) distance are insufficient to guarantee "pure" e-differentia privacy, and one instead requires bounds in the infinity- distance do (v, T) = s [log((()/x()))] see e.g. [ 35 give an algorithm that con- verts samples from TV bounds to those bounded in infinity-distance. Namely, given >> 0 and a sample from a distribution μ within TV distance d S (((((((logg()))))))) /er)- 4 of π,, this post-processing algorithm outputs sample within infinity-distance O(E) from т. Plugging the TV bounds from our Theorem 2.1 into their Theorem 2.2 gives a faster algorithm to sample from a log-concave and log-Lipschitz (or log-smooth) distribution strained to a polytope K, with O(e) error in dx (Corollary C.I in Appendix C In partic- ular, Corollary C.1 gives bound of ((md + d R$) x [LR + d log( (Rd [T, + +T)) arithmetic operations to sample within OEe) error in dx from a log-concave distribution x e constrained to a polytope K from a cold start when f is L-Lipschitz (when f is also 3-smooth, the bound is O(mmd + d X [LR + d (R++ + Tx)). This improves on the bound of O( +m +m2 R2) x [T.R + d log( Rd- Rd- Rd-LR)/ee))(T,+ +Tk) arith- metic operations in Theorem 2.1 of [35] by a factor of dº. Moreover, it further improves on the bound of ((mdd + =md3/ L R2 (1/=2) x polylog(Yee,4) 1/r R, (T,) + Tx) operations implied by 2 Corollary C.1 improves on this bound by a factor of dl "/2² when, e.g, each function evaluation takes Ty (dd) operations and m = O(d) may be the case in privacy applications.
mization B0 23 2
Differentially private optimization. A randomized mechanism h D R is s- dilferentially privale (e-DP) if for any datasets I,, x' E D which differ by a single datapoint, and any S C R, we have P(h(x)) E S) S PP((xx)) E S); see [ a-disferential privacy is the strongest notion of differential privacy, holds several advantages over weaker notions of differential privacy (see Remark 2.2 and has been widely studied in the literature [4 25 2
In the application of the exponential mechanism to c-DP low-rank approximation of a p pxp symmetric matrix M [28] (sce also 19 one wishes to sample within infinity distance O(s) from a log-linear distribution xx e on the Gelfand-Tsetlin polytope K C IRd (which generalizes the probability simplex), where d d= d²,², and where K has m = d inequalities with diameter = O(vd). In this application, the log-linear density f is (trivially) 0-smooth and d'or-Lipschitz where 0n := |[M]]2 is the spectral norm of M. Thus, when applicd to the mechanism of [ our algorithm takes d4.b+ a log(1/=) arithmetic operations. This improves by factor of d on the bound of O(d arithmetic operations implied by
35 and improves by a factor of oi/e on the bound of ofd")//)) implied by 2
Consider the problem of finding an (approximate) minimum i of an empirical risk function f:Kx K x D + R under the constraint that i is e-differentially private, where f(0, x) == EL (0, x.). We assume that the (x(xx)) are L-Lipschitz for all x E D", i N, for some given i > O. In this setting, _ show the minimum ERM utility bound under the constraint that 0 is pure e-DP, Eg[f(0,z)] x)] mineek f(0,x) x) = ((di.R/s) is achieved if one samples from the exponential mechanism æ εx exp(- ( with infinityddistance error at most O(s). Plugging Corollary C.I into the exponential mechanism, we obtain faster algorithm for a pure e-DP mechanism which achieves the minitnum expected risk (Corollary D.1. Specifically, the runtime bound implied by Corollary D.1 is O(md + dn²e?) (en + (od((n///rrr))TTTytt arithmetic operations each é, is Í-Lisschitz (or ((md+ Rs) (en + dlog(nndd/re)))( + Tk)) if f is also B-Lipschitz). This improves upon the bound of
5

---

+ +=2n'd (((2) of max(d" "/le²m) m), nd (V/e), when the & are L-Lipschitz
dd) x polylog(nRtffe))((T + Tk) arithmetic operations in 2 by factor
polytope K and f can be evaluated in T = ((nd) operations. And it improves by a factor of (at least) md on the bound of O((m² d + + d dnn (n x md= operations obtained in [5].
For instance, when applying the exponential mechanism to Lasso logistic regression, if e.g. n = d and ε < 1, our algorithm requires ((d+t)) arithmetic operations, improving by d3 on the bound of dlow operations implied by 35 and by roughly dd-w the bound of O(d operations implicd by 2 In another example, when training a support vector machine with hinge loss and Lasso constraints under e-DP,
algorithm requires O(d--4) arithmetic operations, improving by a factor of de on the bound of do+w operations implied by 35 and by d=-w on the bound of O(d2) implied by 2 See Appendix D for details.
Remark 2.2 (Weaker notions of differential privacy) E-DP holds several practical
advantages over weaker notions of differential privacy (DP), including (E.8)-DP- notion of differential privacy where the privacy of the mechanism is allowed to foil with probability O()). E.g., when group privacy- privacy of subsels of k individuals- must be preserved, any pure e-DP mechanism is also ke-DP. In contrast, (e,6)-DP only implies (c, ke(k-1))8)- DP for subsets of k individuals- the failure probability grows exponentially with k.
Overview of proof of Theorem 2.1
Suppose we are given any polytope K = (0 E TRd AO < b} defined by m inequalities, and a convex f K + Rd which is L-Lipschitz (or B-smooth) and given by an oracle which returns f(0) at any 0 K. Our goal is to sample from π x e- on K within any TV error 5 >> 0, in number of arithmetic operations and oracle calls that has dependence on the dimension that a lower-order polynomial than currently available bounds for sampling from log- Lipschitt (or log-smooth) log-concave distributions, and is logarithmic in 1/8.
Extending the Dikin walk to sample from log-concave distributions on K. As a first step, we begin by attempting to generate samples from π via the (Gaussian) Dikin walk, by extending the standard analysis given in e.g. 24 37 for the special case when i is uniform on K to the more general case where π is a Lig-hipsciitz log-concave on K.
In the special case where π is the uniform distribution interior of K, the Dikin walk proposes updates and (0) = (P) is the Hessian of the barrierrier function -(0)= = Σ}= log(b, 0 [)) for K = (0 E Rd AO < 6}, and a > 0 is a scalar hyperparameter. To ensure that the stationary distribution of the Dikin walk is the uniform distribution on K, if proposed update falls in the interior of K, it is accepted with probability
K, from any point A in the 0 + VaH--10)) where & (N, Ia)
min ((voe(((0) ell=- -1l=0002m00 10 10-=11
determined by the metropolis rule; otherwise, it is rejected. The use of the log-barrier is to ensure that the steps proposed by the Dikin walk remain inside the polytope K w.h.p.
The hyperparaeeter is chosen as large as possible while still ensuring the proposed steps remain in K and are accepted w.h.p. On the one hand, since the covariance matrix aH - of the proposed updates is proportional to a, larger choices of a allow the walk to propose larger update steps. On the other hand, if a is too large, the proposed steps may fall outside the polytope and be rejected w.h.p To see how to choose a, note that for any 0 E Int(K), the Dikin ellipsoid Dg {w (w 0)T (0) 0) S 1} is contained in K. Thus, standard Gaussian concentration inequalities which guarantee [|Ella O(Vi) w.h.p. imply 0 + VaHl+t(0) is in K w.h.p. if a < O(/a). Moreover, using properties of log-barrier functions, one can show the term (e((H()))) in the acceptance ratio is also 2(1) for O1/4),, as is done in 24 B7. To see why, Lemma 4.3 of [b implies that, if (0) is the Hessian of a log- barrier function for K, its log- degedrmeemnnn V(0) log(det(H (0))) satisfies
(VV(O)) (H() < O(d)
V0 Int(K).
(1)

---

Thus, i a I 1/d, the proposed update= : = -++ / aH) has variance (11) in the direction VV(0), and (by Gaussian concentration), (0 - 2) vV(V I O(1) w.h.p. This implies V(2) - (0) = log (=)( - log Q()), and hence del((()))/dd((((0)) = 22)).
=
In [38] the Dikin walk is applied to the more general problem of sampling from a L-log- Lipschitz (or 3-log-smooth) log-concave т x e on K (the problem of interest in this paper). To guarantee the walk has the correct stationary distribution π, the Metropolis acceptance probability of the proposed updates 2 0 + v
xe
where Y is a hyperparameter,
gains an additional factor e ---00) To ensure this acceplance probability remains I(1), they modily the scalar step size 7 such that w.h.p. the walk lakes steps where ʃ changes by O(1). To see how to choose 7, note that since f is L-Lipschitz, (=(=))((0) = 2111) if the Euclidean distance || - 0|12 is (1)2). This can be shown to occur w.h.p. if = 0(//LR))) since the fact that the Dikin ellipsoid is in K C (0, R) implies that the eigenvalues of (() are all S R and hence the variance "ye H (6)u of the proposed step is I ((/((11)) in any given direction E Rd (where v is a unit vector). Thus, it is sufficient for them to choose y = min( \/(LR)) to ensure the proposed step remains in K and is accepted w.h.p.
On the one hand, to ensure the Markov chain proposes steps that change f by an amount at most ((1) for any I--Lipschitz f. it necessary and sufficient to ensure that from any point 0 E Int K), the Markov chain makes updatos which fall w.h.p. inside Euclidean ball B(0, 1/z) of radius 1/z centered at 9. This because the Lipschitz condition on f: llf - -f())]2 S L|00 - =l12 for all 0, 2 E K holds w.r.t. the Euclidean norm || 1l22- On the other hand, to ensure that the Markov chain remains inside the polytope K, it is sufficient to propose steps that lie inside the Dikin ellipsoid D, {w: (w - 0) H (0) -0) 0) S 1} centered at. 0. Roughly, the scalar step size y is chosen such that this ellipsoid is contained inside the Euclidean ball B(0 1/2,, this guarantees that w.h.p. the steps proposed by the walk will both remain in K and will also not change the value of f by more than (11).
However, at many points 0 the Dikin ellipsoid Do is such that the ratio of the largest to smallest eigenvalues of (0) may be very large (this ratio can grow arbitrarily large as approaches face of the polytope). Thus, roughly speaking, modifying the covariance matrix of the Dikin walk by a scalar constant (yH-((0) can cause the Dikin walk to propose steps whose variance in some directions is much smaller than is required for either of the two goals: staying inside K and slaying inside the ball B(0, 1/z) defined by the Lipschitz condition on f. This suggests that modifying the log-barrier function for K by scalar multiple may not be the most efficient way of extending the Dikin walk to the problem of sampling from a general L-log- Liog-Lipschitz (or 3-log-smooth) log- concave distribution on K, and that one may be able to obtain faster runtimes by making other modifications to the barrier function.
A soft-threshold regularized Dikin walk. Before we introduce our soft-threshold Dikin walk, we first note that even in the special case where is the uniform distribution on K, the analysis in BS does not recover the bounds given in 24 Br| for this special case, as [8 use a different runtime analysis geared to lime-varying distributions studied in that paper. Namely, BN imply a bound of O(m²d log(</8)) steps to sample from a uniform distribution on K, while 24 337 show a bound of O log( (=) For this reason, first extend the analysis of the Gaussian Dikin walk given in 24 37 for the special case of uniform т, to the more general problem of sampling from an E-log-Lipschin
5-log-smooth log- concave density. The analysis in 24 37 uses the cross-ratio distance metric. More specifically, if for any distinct u, E Int(K) (K) we let p. g be the endpoints of the chord in K passing through and u such that the four points lic in the order p, u, v. q. the cross-ratio distance is
ou, v) = |u u112 x 1pp gll2
(2)
pp x gll2*
1 One can show that for any u, v E Int(K), 0² (u, v) 2 (4/mm>>))ll - v - H()) (see e.g.
24
37 ). Thus, as the usual Dikin walk takes steps that have roughly identity covariance matrix
Ta with respect to the local norm ||w|,=>=0(0== == (uTn-IH(0),, for
min(1/a, (L(())))),
the bound we would obtain on the number of steps until the Dikin walk is within TV error
5 from is O(A-1 -log(=))) 0(mmd+ + mL2 R²) log(w/6) steps from a w-warm start.
7

---

To obtain even faster bounds, we would ideally like to allow the Dikin walk to take larger steps by choosing a larger Y closer to the value 1/a that is sufficient to ensure an Q(1) acceptance probability in the special case when π is uniform. Unfortunately, if e.g. LR 2 d, reducing y from 1/a to (V/mn))) may be necessary to ensure the variance of the Dikin walk steps is (1/(at2)) in every direction, and hence that the acceptance probability is 2(1).
To get around this problem, we introduce a new variant of the Dikin walk Markov chain for (or 8-log-smooth) log-concave distributions on polytope K, which generalizes the Dikin walk introduced in [2 for sampling from æ in the special case when is the uniform distribution on K. The main difference between our Dikin walk and the usual Dikin walk of [4 (and of [8]] is that our Dikin walk regularizes the regularization term -a proportional to the identity matrix, where 1 is a hyperparameter and a the same hyperparameter appearing in the original Dikin walk of [24] Since the log-barrier Hessian H0) and regularization term 17 - have different scalar hyperparameters a, 7,, we can set and 1 independently from each other: roughly, a= is chosen to be the largest value such that the Dikin ellipsoid defined by the matrix a -1 (0) remains inside K, while n is independently chosen to be the largest value such that, with high probability, the steps proposed by our Markov chain remain inside the ball B(0,1/2 defined by the Lipschitz condition on f. Roughly, the addition of the soft-threshole regularization term to the log- barrier Hessian allows us to reduce the variance of the proposed steps of the Dikin walk only in those directions where a choice of = 1/a would cause the variance to be greater than 1/(ar") while leaving the variance in other directions unchanged.
sampling from any
Hessian a= (0) of the log- barrier for K by adding a
More specifically, the steps proposed by our soft-threshold Dikin walk are Gaussian with + mean 0 and covariance matrix & --1(0)= == (a)
The addition of the soft- to ensure that the largest eigenvalues of the
threshold regularization term 7) 7-' La allows covariance matrix & (0) are S 0(/(a1²)) without reducing (by more than a constant factor) the eigenvalues which were already S O(/(a13)). This allows our Dikin walk to take larger steps, while still ensuring these steps are accepted w.h.p. by the Metropolis accept/reject rule for f. Taking larger steps allows our Dikin walk to converge more quickly to the target distribution (0) e --((0) on K.
The (inverse) covariance matrix ((0) of the steps proposed by our Dikin walk is the Hessian of the function (0) a (0)) + barrier for K. The modified function (0) can be seen to also be a self-concordant barrier for K. In the special case where æ is the uniform distribution, L = 0 and 7 = 0, and our Dikin walk recovers the original walk of 24. Thus, our walk generalizes the original Dikin walk to the problem of sampling from general smoothh log-concave distribution on a polytope.
11 where p(0) -ド log(b, - az 0) is the log-
(or 8-log-
Bounding the number of Markov chain steps. Setting n = 1/(a13) ensures the vari- ance u. affff11000 in any given unit-vector direction v of the proposed update z - 0 of our Markov chain is at most 0(/(a12)), and hence the term (=()))(()) in the Metropolis accep- tance rule is (((1) with high probability (Lemma E.5p. Moreover, we also show that, if we choose a = 1/d, the other terms in the Metropolis acceptance rule are also 2((1) (Lemmas E.8 E.9 While the proofs of these lemmas follow roughly the same outline as in the spe- cial case of the original Dikin walk where л is uniform (e.g., 24 37) ) our bound on the determinantal term (et(())/da0(0)) must deal with additional challenges, which we discuss in the next subsection.
To bound the number of steps required by our Markov chain, we first bound the cross-ratio distance o au, v) between any u, a E [n)( by the local norm I|u - vll (Lemma E.2
1
(a (u u))² (a a b.)²
1
Iu 2
1
0² v)
+
2ma 1 + 2n 122 1|u u U122())
(3)
2m
2
R²
Using (3] Logether with the isoperimetric inequality for the cross-ratio distance (Theorem 2.2 of 32 ). we show that, if the acceptance probability of our Markov chain is 2(1) at each step, then the number of steps for our Markov chain to obtain a sample within a TV distance of 5 from is O(A-1 -1log(4/6) = O((md+ + dL R²) log(w/6) from an w- warm start,
8

---

where A = (2//22a--42y11 In particular, if m = O(d) and LR > d, this improves on the bound we would get for the basic Dikin walk by factor of
Bounding the determinantal term in the acceptance probability. For our mixing time bound to hold, we still need to show that dot(0())/aat(000)) is ((11) w.h.p. We would idcally like to follow the general approach previous works 24 37 use to show the term de(((())aaa((((()) in the basic Dikin walk is ((1) w.h.p., which relies on the property of log- barrier functions in Inequality ① Unfortunately, as ((0) is not the Hessian of a log-barrier function for any system of inequalities defining K, we cannot directly apply [ to (()).
To get around this problem, we show that, while (()) is not the Hessian of log-barring function, it is in fact the limit of a sequence of matrices H,( E N, where each matrix (0) in this sequence is the Hessian of a (different) log-barrier function for K. Specifically, for every j E N. we consider the matrices A = [AAT, II Ia......4. Ial where A is concatenated with (m,-1)/a copies of the identity matrix Ia, and m, = m + Lam ] And we consider the vectors b = (b" .j1 T "...,.""""..................." )TT where is concatenated with (m--11/a copies of the vector j1, where
(1, (1................ I)T R2 is the all-ones vector. Then, for large enough j,
K (0 E Rdd A < }}, and the Hessian of the corresponding log-barrier functions is
H2( = Em d((a)) (a)))T00))2 H(0) + lan-
ese, (o) 0-)))
(4)
Using this fact 4 we show that, for every 0 E Kn)(K), every E Do, and every sequence { C & Da such that 2j 2, we have (Lemma E.7
lim (et()) (=))) 8-000 det(
(5)
det(
Moreover, since each H, is the Hessian of log- barrier function for K, we have that [1] does hold for H, and hence (from the work of 24 37) that (=)(Hi,)))/aau (0)) = 52(1) w.h.p. for all N, if we set 2j = 0 H 1/2 (0), h~ (((0 Ia), and choose = 1/a. Thus, (5) implies that (el(((=))/ee(()()) 2(1) w.h.p. as well (Lemma E.9 and hence the acceptance probability of the soft-threshold Dikin walk is ((1)) at each step.
Bounding the number of arithmetic operations. Since the acceptance probability is Q(1), from the above discussion, the number of steps for our Markov chain to obtain sample within TV distance > 0 from π is 0((md + R²))) log((())) from a w-warm start.
Each time our Markov chain proposes step 0 + (00)) it must first sample a Gaussian vector & ~ NO, Ia) which takes (d) arithmetic operations. It must then compute the log- log-barrier Hessian H (0), and invert the matrix (()) = a (() + +17 -
Since (()) = (O)) where C(0) is a d m matrix with columns cy(0) = as/(a] 0-b;) for all j E [m], we can compute (f)) in md=-1 arithmetic operations. And since (()) is a dxd matrix, computing & (0) can be accomplished in d" arithmetic operations by computing the singular value decomposition (SVD) of Q. Next, we compute the acceptance probability
min (- /det(<(()) 11=0011 100= 1): The determinants can be computed in ((d)) 100 (det(((0)))
arithmetic operations via the SVD. Evaluating [(), takes two calls to the oracle for J.
Thus, from w-warm start, the soft-threshold Dikin walk takes at most 0((md + dT.2 x (g(()))) Markov chain steps to obtain a sample from æ with total variation error ô > 0, where each step takes Ommd=11) arithmetic operations, and one function evaluation.
Remark 3.1 (Optimality of eg-regularization) For the class of functions considered in our paper, we conjecture that an lz-regularizer that does not depend on f is optimal. This is because we consider the class of functions which are L-Lipschilz or B-smooth with respect to the t2-norm, and our bound on L. or 3 does not depend on 0. Moreover, we only have access the function through an oraele which returns the value of f at any given point but does not tell us how ff changes at nearby points.
9

---

Algorithm 1: Soft-threshold Dikin walk
Input: m, E N, A E Rm=x be R", which define K := (0 E IR" : b).
Input: Oracle returning the value of convex f K -> RR. Initial point 8o E [nt(()).
Hyperparameters: a > 0, 17 0, and T N.
Set 8 4 0o
3 i I.... do
Sample point & ~v ((,
Set (() + 2-=1 [u] 2sa
b-1
Set (0) +- - ((0) +1 La
Set 0 +0) (00)-112
E Int(K) then
Set (z) 25=1 (a b-b
10
Set ((z) + a (z)) +
Accept 2 with probability 1 min
-1()) / (a))
x -|2 (l=-01ãc00 -110- =locco, 1)
11
5(6) (o))
12
else
Reject
end
end Output 8
In Theorem 2.1 we set the step size hyperparameters a 2/(10/ and 77 1/(10 if f is L-Lipschitz, and the number of steps to be T = 10 (2ma R2) log(w/6). When is 8-smooth (but not necessarily Lipschitz), we instead set = 1/(10')) and 1 1/(10ta3).
In many applications, bound on L or can be calculated analytically, allowing one to set a, as above. This includes, e.g., applications to training Bayesian or differentially private logistic regression models (or other generalized linear models such support vector machines). When a bound on I or B is not known, one can in practice set a, by hand such that the average acceptance probability 2 2(1).
Conclusions, limitations, and future work
Our result improves on the runtime bounds of line of previous work, for the problem of sampling from several classes of log-Lipschitz or log-smooth log-concave distributions on polytope (see Table 1). To the best of our knowledge, this is the first result to introduce regularized barrier functions that simultancously take into account the geometry of both the constraint polytope and the Lipschitz or smoothness property of target logconcave func- tion. These barrier functions may be of independent interest for sampling or optimization.
On the other hand, we note that 27 give an implementation of the Dikin walk in the special case where is constant, where the (average) cost of computing the Hessian matrix of the log-barrier of the polytope K (0 TRd AO < b} at each step of the walk is improved to roughly (d (nz((A) arithmetic operations, where nnz(A) (A) denotes the number of non-zero entries of A. Whether this improvement. in the per-step computation time be achieved for the problem of computing the regularized barrier functions used in our algorithm, in the more general setting where L-Lipschitz or 3-smooth, is an interesting open problem.
Moreover, we note that our bounds are polynomial in L or 3, yet there are algorithms for sampling from log-concave distributions x e which do not assume f is L-Lipschitz 3- smooth. Thus, another interesting open problem whether one can obtain runtime bounds for version of the Dikin walk which do not require f to be Lipschitz smooth.
Our results have applications to Bayesian inference and differentially private optimization. Bayesian inference can lead to algorithms with better generalization properties and quan- tification of uncortainty, and differential privacy guarantecs are important to protecting the privacy of individuals in medical and other sensitive datasets. Thus, we believe our results will have positive societal impacts, and do not anticipate any negative impacts to society.
10

---

Acknowledgments. NV was supported in part by an NSF CCF-2112665 award. OM was supported in part by an NSF CCF-2104528 award and Google Research Scholar award.
References
[1] David Applegate and Ravi Kannan. Sampling and integration of near log- concave functions. In Proceedings of the twenty-third onnual ACM symposium on Theory of computing, pages 156 163, 1991.
[2] Raef Bassily, Adam Smith, and Abhradeep Thakurta. Private empirical risk mini- mization: Efficient algorithms and light error bounds. In 2014 IEEE 55th Annual Symposium on Foundations of Computer Science, pages 464 473. IEEE, 2014.
3] Nicolas Brosse, Alain Durmus, Éric Moulines, and Marcelo Pereyra. Sampling from a log-concave distribution with compact support with proximal Langevin Monte Carlo. In Conference on Learning Theory, pages 319 342. PMLR, 2017.
[4] Sébastien Bubeck and Ronen Eldan. The entropic barrier: simple and optimal univer- sal self-concordant barrier. In Conference on Learning Theory, pages 279 279. PMLR, 2015.
[5] Sebastien Bubeck, Ronen Eldan, and Joseph Lchec. Finite-time analysis of projected Langevin Monte Carlo. In Advances in Neural Information Processing Systems, pages 1251- Citescor, 2015.
[6] Bradley P Carlin and Siddhartha Chib. Bayesian model choice via markov chain monte carlo methods. Journal of the Royal Statistical Society: Series B (Methodological), 57(3):473-44) 1995.
[7] Apostolos Chalkis, Vissarion Fisikopoulos, Marios Papachristou, and Elias Tsigaridas. Truncated log-concave sampling for convex bodies with reflective Hamiltonian Monte Carlo. ACM Transactions on Mathematical Software, 2023.
[8] Yuansi Chen, Raaz Dwivedi, Martin J Wainwright, and Bin Yu. Vaidya walk: A sampling algorithm based on the volumetric barrier. In 2017 55th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 1220 1227. IEEE, 2017.
[9] Arnak S Dalalyan and Lionel Riou-Durand. On sampling from a log-concave density using kinetic Langevin diffusions. Bernoulli, 26(3):1956-1988, 2020.
[10] John Duchi. Derivations for linear algebra and optimization. Berkeley, California, 3:2325 5870, 2007.
[11] Alain Durmus, Szymon Majewski, and Blazej Miasojedow. Analysis of Langevin Monte Carlo via convex optimization. J. Mach. Learn. Res., 20:73-1, 2019.
[12] Alain Durmus and Eric Moulines. Nonasymptotic convergence analysis for the unad- justed Langevin algorithm. Annals of Applied Probability, 27 1751111111 2017.
[13] Raaz Dwivedi, Yuansi Chen, Martin J Wainwright, and Bin Yu. Log-concave sampling: cattropolishhastnngs algorithms are fast! In Conference on Learning Theory, pages 793 797. PMLR, 2018
[14] Cynthia Dwork and Aaron Roth. The algorithmic foundations of differential privacy. Foundations and Trends in Theoretical Computer Science, 9(3-4)2211440) 2014.
[15] Alan Frieze and Ravi Kannan. Sobolev inequalities and sampling from loncave distributions. The Annals of Applied Probability, 9(1):14 26, 1999.
[16] Alan Frieze, Ravi Kannan, and Nick Polson. Sampling from log-concave distributions. The Annals of Applied Probability, pages 812-837 1994.
11

---

[17] Dani Gamerman and Hedibert F Lopes. Markov chain Monte Carlo: lation for Bayesian inference. CRC press, 2006.
stochastic simu-
[18] Arun Ganesh and Kunal Talwar. Faster differentially private samplers via Rényi di- vergence analysis of discretized Langevin MCMC. Advances in Newral Information Processing Systems, 33, 2020.
[19] Rong Ge, Holden Lee, Jianfeng Lu, and Andrej Risteski. Elficient sampling from the Bingham distribution. In Vitaly Feldman, Katrina Ligett, and Sivan Sabato, editors, Algorithmic Learning Theory, 16-19 March 2021, Virtual Conference, Worldwide, vol-
ume 132 of Proceedings of Machine Learning Research, pages 673 685. PMLR, 2021.
[20] Andrew Gelman and Donald B Rubin. Markov chain Monte Carlo methods in bio- statistics. Statistical methods in medical research, 5(4):339-35) 3535 1996.
[21] Sivakanth Gopi, Yin Tat Lec, and Daogao Liu. Private convex optimization via ex- ponential mechanism. In Po-Ling Loh and Maxim Raginsky, editors, Proceedings of Thirly Fifth Conference on Learning Theory, volume 178 of Proceedings of Machine Learning Research, pages 1948 1989. PMLR, 02 05 Jul 2022.
[22] David Lee Hanson and Farroll Tim Wright. A bound on tail probabilities for quadratic forms in independent random variables. The Annals of Mathemalioal Statistics, 42(3))1079-1100 1083, 1971.
[23] Moritz Hardt and Kunal Talwar. On the geometry of differential privacy. In Proceedings of the forty-second ACM symposium on Theory of computing, pages 705-714, 2010.
[24] Ravi Kannan and Hariharan Narayanan. Random walks on polytopes and an affine interior point method for linear programming. Mathematics of Operations Research, 37(1):1-20, 2012.
[25] Michael Kapralov and Kunal Talwar. On differentially private low rank approxima- tion. In Proceedings of the twendy--ourth annual ACM-SIAM symposium on Discrete algorilhms, pages 1395 1414. SIAM, 2013.
[26] Sun Mi Kim, Yongdai Kim, Kuhwan Jeong, Heeyeong Jeong, and Jiyoung Kim. Logistic lasso regression for the diagnosis of breast cancer using clinical demographic data and the bi-rads lexicon for ultrasonography. Ultrasonography. 37(1))36, 2018.
[27] Aditi Laddha, Yin Tat Lee, and Santosh Vempala. Strong self-concordance and sam- pling. In Proceedings of the 52nd annual ACM SIGACT symposium theory of com- puting, pages 1212-1222, 2020.
[28] Jonathan Leake, Colin S McSwiggen, and Nishceth K Vishnoi. A polynomial-time Saanddrr-tyykon--222 Zuber
algorithm and applications for matrix sampling from densities. In ACM symposium on theory of computing STOC, 2021.
[29] Yin Tat Lee and Aaron Sidford. Solving linear programs with (qrt(rank) linear system solves. arXiv preprint ar arXiv:1910.0810 2019.
[30] Yin Tat Lee and Santosh S Vempala. Convergence rate of Riemannian Hamiltonian Monte Carlo and faster polytope volume computation. In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing, pages 1115 1121, 2018.
[31] László Lovász and Miklós Simonovits. Random walks in a convex body and an improved volume algorithm. Random structures & algorithms, 4(4):350-412, 1993.
[32] László Lovász and Santosh Vempala. Hit-and-run is fast and fun. preprint, Microsoft Research, 2003.
[33] Lásuló Lovász and Santosh Vempala. Fast algorithms for logconcave functions: Sam- pling, rounding, integration and optimization. Tn Foundations of Compuler Science, 2006. FOCS 4 7th Annual IEFE Symposium on, pages 68. IEFF, 2006.
12

---

[34] László Lovász and Santosh Vempala. The geometry of logconcave functions and sam- pling algorithms. Random Structures & Algorithms, 30(3):307-355 2007.
[35] Oren Mangoubi and Nisheeth K. Vishnoi. Sampling from log-concave distributions with infinity-dislance guarantees. In Advances in Newral Information Processing Systems, 2022.
[36] Frank McSherry and Kunal Talwar. Mechanism design via differential privacy. In 48th Annual IEEE Symposium Foundations of Computer Science (FOCS'07), pages
94 103. IEEE, 2007.
[37] Hariharan Narayanan. Randomized interior point methods for sampling and optimiza- tion. The Annals of Applied Probability, 26(1):597 641, 2016.
[38] Hariharan Narayanan and Alexander Rakhlin, Efficient sampling from time-varying log-concave distributions. The Journal of Machine Learning Research, 18()):401744445, 2017.
[39] Yurii Nesterov and Arkadi Nemirovskii. Interior-point polynomial algorithms in conver programming. SLAM, 1994.
[40] Mark Rudelson and Roman Vershynin. Hanson- inequality and sub-gaussian concentration. Electronic Communications in Probabilily, 8(82):1 9, 2013.
[41] Sushant Sachdeva and Nisheeth K Vishnoi. The mixing time of the Dikin walk in polytopo-a simple proof. Operations Research Letters, 44(5):630-634, 2016.
[42) Adil Salim and Peter Richtárík. Primal dual interpretation of the proximal stochas- tic gradient langevin algorithm. Advances in Neural Information Processing Systems, 33:3786 3796, 2020.
[43] Mervyn J Silvapulle and Pranab Kumar Sen. Constrained statistical inference: Order, imequality, and shape constraints. John Wiley & Sons, 2011.
[44] Guo-Liang Tian, Man-Lai Tang, Hong-Bin Fang, and Ming Tan. Efficient methods for estimating constrained parameters with applications to regularized (lasso) logistic regression. Compulational slatistics & dala analysis, 52(7):3528 3542, 2008.
[45] Pravin M Vaidya and David S Atkinson. A technique for bounding the number of iterations in path following algorithms. In Complezity in Numerical Optimization, pages 462-489. World Scientific, 1993.
[46] Andre Wibisono. Sampling as optimization in the space of measures: The Langevin dynamics as composite optimization problem. In Conference on Learning Theory, pages 2093-300300 3027. PMLR, 2018
[47] Yuan Zhang, Jonathan N Hofmann, Mark Purduc, Shili Lin, and Swati Biswas. Logistic baycsian lasso for genetic association analysis of data from complex sampling designs. Journal of Human Genelics, 62(9))819 829, 2017.
13

---

Contents
Introduction
2 Results
3
Overview of proof of Theorem 2.1
Conclusions, limitations, and future work
10
Additional discussion with related work
B Application to Bayesian Lasso logistic regression
16
Infinity-distance sampling
D Differentially private optimization
17
Proof of Theorem 2.1
E.
Bounding the number of arithmellic operations
18
E.2
Bounding the step size
19
................................
E.3 Bounding the acceptance probability
Bounding the conductance
29
Extension to general barrier functions?
33
Proof Lemmall
34
2
Lower bounds for sell-concordance parameler
35
14

---

A Additional discussion with related work
The problem of sampling from the uniform distribution when K is polytope given by a set of inequalities, and the distance to π is measured in total-variation (TV) distance has been Also widely studied is the problem of sampling from an unconstrained log-concave distribution. One direction of study includes works that require the log-density π to be hipg-Lipp or 3-log-smooth on all RK for some L, > O, and give bounds on the distance to in terms of TV distance
widely studied; see, e-g-; 24
[3 Wasserstein distance
12 9
E2
and Kullback-Leible (KL) divergence 4 IT.
5 3 provide versions of the Langevin dynamics for sampling from log- concave distribution л x e- on K where ʃ is I-Lipschitz and f-smooth and K is given by a projection oracle- when, e-g.- K is contained in a ball of radius O(1), it contains a ball of radius (21), and I, = / = O(1), 5 give a bound of roughly O(d /822) gradient and projection oracle calls sample from т with TV error d O, while [ give bound of roughly 0(4%/8) gradient and projection oracle calls. [2]] provide an algorithm for sampling from distributions x e- where f is both L-Lipschitz and p-strongly convex (that is, f(0) is the sum of a convex function f(0) and the quadratic function (0/2)]101² for some / > 0) on a convex body K in roughly (((²2/x) (4/)())) membership and function oracle calls (see also 77 who require f to be μ-strongly convex for μ > and B-smooth). However, the algorithms of [2I. 7 do not apply to the more general setting where f is convex and L-T.ipschitz (or B-smooth), but not necessarily p-strongly convex for any / > 0.
Further,
Detailed comparison to 38 Comparod to in our version of the Dikin walk, 2 is sampled from Gaussian with covariance matrix
38
(a- (1(2))+1111 (==)
where a- and 1) are hyper- pnpeer-prameterrs chosen to be a ~ d² and 17 dL2 if is L-Lipschitz (or n7 ~ d3 if f is L-smooth). The "het-thrrehold" regularization term n La prevents the Markov chain from taking steps where the value of f decreases by more than O(1) w.h.p., ensuring that the term The soft-threshold regularizer chosen to be a multiple of the identity matrix la since the Lipschitz condition on f' is rotationally invariant- it bounds the derivative of f by the same amount L in each direction (the same true for the sccond derivative of f if f is f-smooth). This in turn allows our choice of scaling a- for the "Dikin ellipsoid" term Hz2)- which is not in general rotationally invariant, and determined only by the geometry of the polytope rather than the geometry of the function ʃ to be independent of L (or 3). This is in rontrast to the Dikin walk in BSJ where the scaling parameter for H(z) must depend on (or B) to ensure an ?(1) acceptance probability, which allows our Markov chain to propose steps with a larger variance than the Dikin walk in 38 in directions which are not the largest eigenvector of H (s).
in the acceptance probability is I(11..
The (inverse) covariance matrix a (-) of our soft-threshold Dikin walk updates is the Hossian of the function v0()) = p(0) + 17 ||0113 This barrier function can be seen to be a Hessian of a self-concordant barrier. On the other hand, it is not the Hessian of a logarithmic-barrier function for any polytope defined by any set of inequalities. This prevents us from directly applying the analysis of the Dikin walk in the special case where == 0 B7 which relies on properties of log- log-barrier functions- to our soft-threshold Dikin walk on Lipschitz or smooth f. To get around this problem, we show that, while eb() is not log- barrier function of any polytope K, it is the limit of a sequence of log- barrier functions E.... where (0) + (0) uniformly in 0 as + co. See Section 3 for
detailed overview of the proof.
An open problem is to obtain runtime bounds for the Dikin walk that do not require f to be Lipschitz or smooth, and/or that depend polynomially on log R. This leads to the question of whether one can design other tractable self-concordant barriers to obtain further runtime improvements for sampling log-concave distributions on a polytope. We discuss possible extensions in Appendix
15

---

Explanation of bounds in 38 38 consider a version of the Dikin walk for sampling from a log-concave distribution x x e on a polytope K contained in a ball of radius > (and also more generally on any convex body with self-concordant barrier function). Their version of the Dikin walk has proposed updales = a² H 2 (0)S where & ~ N(, Ia), and acceptance probability given by min c (-1 [dec[H(=)) :) Here H is the hessian of /-
barrier function for the convex body K, which may be the log-barrier or another barrier function. For a step size of a (min (2 HR)) (Condition 4 in Section 3.1 of their paper; note that is considered a constant in [38] and therefore does appear explicitly in their bounds 3 they obtain bound on the conductance & of their version of the Dikin walk (Lemma 4 of their paper):
Q2
(6)
RVd
where x the coloocanner parameter of the barrier function; for the log- log-barrier A = m, although there are other barrier functions for which к Od)) Plugging Inequality 6 into Corollary 1.5 of 3T implics that their Dikin walk can obtain a sample from a with ( TV error in roughly & log( 3) = O((d di L R2) log( g steps from a w-warm start.
The number of arithmetic operations per step is md=-1 and the number of arithmetic oper- ations to obtain a sample with TV error O6)) from π is O(mm++# + mdd+wL3 log (og())) )))
If one also has that K contains ball of radius > 0, from cold start initial point which is ed B 4 LR by sampling from the uniform distribution this ball, and hence their bound O(mm + mdd-w R² (d + LR +
find
(
B Application to Bayesian Lasso logistic regression
The following corollary follows directly from Theorem 2.1
Corollary B.1 There erists an algorilim which, given 8, > O and dalasel (aili= C Rd S.L. || 2| < 1 V i E [n], and lelting [(0) Ei= E(0Ta), x:) where l(x) ITe-= is the logistic loss, K {0 E Rd [10]]1 S c} the Lasso constraint, and π e has support K, outputs point from a distribution μ s.t. ||e (|p-m ~V I 6, in O(dd + nnde) + nd?) log(()) arithmetic operations.
C Infnitt-distance sampling
In applications of sampling to differentially private optimization 30 23 2
18
2 bounds
in the total variation (TV) distance are insufficient to guarantee "pure" differential privacy, and one instead requires bounds in the infinity-distance dx (1, == supge k [log( (()())((()))) ; see e.g. [4. 35 give an algorithm that converts samples from TV bounds to those bounded in infmnit--disance. Namely, given & 0 and sample from a distribution p within TV distance ö O (E [Ridllog((/)++ -d - of π, this post-processing algorithm outputs
a sample from a distribution a with infinity distance doe (v, T) < from Plugging the TV bounds from our Theorem [2.1 into their Theorem 2.2 gives a faster algorithm sample from a log-concave and log-Lipschitz distribution constrained to a polytope K, with ((e) error in d In particular, Corollary C.1 improves the bound of O(m2dd + m R²) x [LR + d (d+=Ld/(r=)))(TT) identical to that of how their Theorem 2.2 implies their Theorem 2.1, and we refer the reader to their paper.
Tk) arithmetic operations in their Theorem 2.1. The proof
Corollary C.1 (Sampling with infinity-distane guarantees) There exists an algo- rithm which, given E, L,,, R > 0, E Rmmxd b E Rmm (and possibly B > 0), that de- fine a polytope K = {0 E Rdd A0
< b} contained in a ball of radius R, a point
There is a typo page 7 of 38): "R" should be the radius of the smallest ball containing K, not the radius of the largest ball contained in K.
16

---

a E Rd such that K contains a ball B(a,r) of smaller radius T, and an oracle for the value of a conver function f K - R2, where I is I.TTipschitz (or is both L-- Lipschitz and 3-smooth), and defining π to be the distribution x & e-1, outputs a point from distribution v such that dx (V, < E. Moreover, this algorithm takes O ((md + RR) x LR + dlog ( (Ra+LRd) )]) (T, + TK) arithmetic operations, where Tj is the number of arithmetic operations to evaluate f and Tк O(md=-1). If is also B- Lipschitz, the number of operations is O ((md d x [R + d dllog (dd+LRd
)l) (Ty+TK).
re
Corollary C.1 further improves the dependence on the dimension d over 2 35 Specifically, 2 implies bound of ((mdd md L (t/=)) polylog(//s, 1/r, R, L, (T) + Tk) arith- metic operations, where Ty is the number of operations to compute and Tx the operations to compute a projection oracle for K. Corollary C.1 improves on this bound by factor of roughly d'h */(²mm when, e.g, each function evaluation takes T, = =Od²) operations. For example, when m Od),, as may be the case in privacy applications, the improvement is d' "/2. Moreover, it improves by d on the bound of (((m² + min RI, R x [LR + d (oog((((t)R))/rr)))
+ Tx) arithmetic operations obtained in 35
D Differentially private optimization
A randomized mechanism h : D + R said to be e-differentially private (s-DP) for any datasets x E D which differ by a single datapoint, and any S S RR we have that P(h(x) E S) S (h(('z) E S); see 1.. In the application of the exponential mechanism to c-DP low-rank approximation of a p x p symmetric matrix M 28 (sce also [9] one wishes to sample within infinity distance O(E) from log- linear distribution x e on the Gelfand-Tsetlin polytope K C Rdd (which generalizes the probability simplex), where = p², and where K has m = d inequalities with diameter R = O(vØ). In this application, the log-linear density is (trivially) 0-smooth and dio--iips,,,,,, where 01 := ||M]|2 the spectral norm of M. Thus, when applied to the mechanism of 28), our algorithm takes d45++ 01 log(1/e) arithmetic operations. This improves by a factor of d on the runtime bound of (dd a)) arithmetic operations implied by [5 B8, , and improves by a factor of ofte on the bound of O (d"")t/+) arithmetic operations for the bound in [2].
Consider the problem of finding an (approximate) minimum è of an empirical risk function f: K x D + R under the constraint that Ô is s-differentially private, where f00 := Σx (0 .r)). We assume that the (()xx) are î-T.iischitz for all x E DD, i E N, for some given L > O. In this setting, [2] show the minimum ERM utility bound under the constraint that i is pure eDDP, [=] - minsek (00,x) = ((dLR/))) is achieved if one samples 0 from the exponential mechanism & 1ap(-2i.Rd1 with infinity-distance error at most O(e). Plugging Corollary C.1 into the exponential mechanism, we obtain a faster algorithm for a pure DP mechanism which achieves the minimum expected risk (Corollary D.I.
Corollary D.1 (Differentially private empirical risk minimization) There erists
algorithm which, given L, r, R > 0, A E Rmxdd b E RT (and possibly B > 0) that define a polytope K == {0 E IR A0 < b} contained in ball of radius R and containing ball of smaller radius r, and an empirical risk function f(0,, x) = : Ei= (, xr), where each 4 : K + R is L-Lipschitz (and possibly also 3 - smooth), outputs random point Ô E K which is pure e-differentially privale and satisfies Eg x:) mingek f(0,a:) S O(dLn/e). Moreover, this algorithm takes (((dd + d (en + dlog (Rdd/(r)))( + Tk)) arithmetic Tx)) is f is also B- L-pischhil where T3 is the number of operations to evaluate f and Tx = 1(mdw-1).
operations if each E, is L-Lipschitz (or O((md - dn & RE) x (en +
Corollary D.1 improves on 2 35.. In particular, it improves upon the bound of (1E2))
d x (olylog(nRd/rre)))) + Tx) arithmetic operations in [2] by a factor of roughly max(d'0 -/(2m), nd (1/=) arithmetic operations, in the setting where the l are L-pisshiizz
17

---

on a polytope K and f can be evaluated in Tÿ = ((nd) operations (Tk is the number of operations compute a projection oracle for K). And it improves by a factor of (at least) md on the bound of ( + m (en + (vRd(v)))) operations obtained in R5. The proof is identical to that of how their Theorem 2.2 implics their Corollary 2.4, and we refer the reader to their paper. For instance, when applying the exponential mechanism to the Lasso logistic regression problem, each loss l is both 3-smooth and i-Lipschitz for 3 = f. 1 and R O(1). Thus, if n = d and < 1, for this problem our algorithm requires O(d+t)) arithmetic operations, which improves by factor of d on the bound of d arithmetic operations implied by Corollary 2.4 in 35 and by roughly d=- on the bound O(d")) arithmetic operations implied by [2] In another example, when training a support vector machine model with hinge loss and Lasso has that & is L-Lipschitz but not smooth, for i= I, and R O(1). Thus, if n d and < our algorithm requires O(d++*) arithmetic operations, which improves by a factor d² on the bound of dotw operations implied by Corollary 2.4 in 35 and by roughly d=~ the bound of O(d'2) implicd by 2
mdd-1) arithmetic
constraints under e-differential privacy,
E Proof of Theorem 2.1
E.1 Bounding the number of arithmetic operations
In the following, we assume the hyperparameters a, satisfy a S log( for some choice 8 > 0, and either v S T (in the setting where f is I-Lipschitz) or S 1o (in the setting where f is 8-smooth). To obtain our main result (Theorem 2.1 we set 8 100 To obtain improved runtime bounds in the setting where A is sparse, will set
8 100 aL" TT 8 when f is L-Lipschitz, and 8 100(mdd33R"))og 1o when ʃ is
8-smooth.
Lemma E.1 Each iteration of Algorithm can be implemented in O(md==1) arithmetic operations plus O(1) calls to the oracle for the value of f.
Proof: We go through each step of Algorithm and add up the number of arithmetic operations and oracle calls for each step:
1. Line 4 samples a d-dimensional Gaussian random vector & No, Ia), which can be performed in (dd) arithmetic operations.
2. At each iteration, Algorithm 1 computes the matrix Hww) at = 0 (line 5 and
w = (line 9):
aja H(w) = Σ (aJ w b₂)
Computing H(w) at any w e gd can be accomplished in md=-1 operations as follows:
Define C(w) be the d m matrix where each column cy(w) aj w b
for all
j [m]. Then
Hww) = C(w)
Since ((w) is a m matrix the product C(w) can be computed in mdd-1 arithmetic operations. Thus, Lines 5 and g of Algorithm [ can each be computed in md arithmetic operations.
3. Since Lines and 10 compute sum of two x matrices, Lines and 10 can each be performed in d arithmetic operations.
4. Line 7 computes the proposed update
= 0 (() igg.
18

---

Computing $(6) can be performed by taking the singular value decomposition of the matrix { and then inverting and taking the square root of its eigenval- ues. This can be accomplished in d arithmetic operations since Φ is d x d matrix. Once ((0)- i is computed, the computation 0 1b(00)+g &gg can be performed in d² arithmetic operations. Thus Line can be computed in ((d') s O(md--1) arithmetic operations.
5. Line @ requires determining whether 2 E K. This can be accomplished in ( arithmetic operations, by checking whether the inequality Az S is satisfied.
6. Line 11 requires computing the determinant det det( and f(w) for 0 and w 2. Computing det(@ (w() can be accomplished by computing the singular value decomposition of dee (W()),, and then laking the product of the resulting singular values to compule the determinant. Since ((w) is a x matrix, computing the singular value decomposition can be done in d arithmetic operations. Computing f(w) for any w 5 I can be accomplished in one call to the oracle for the value of f. Thus, Line 11 can be computed in (dd)) I O(md==1) arithmetic operations, plus two calls to the oracle for the value of f.
Therefore, adding up the number of arithmetic operations and oracle calls from all the different steps of Algorithm [ we get that each iteration of Algorithm [ can be computed in 0(md=11) arithmetic operations plus ((1) calls to the oracle for the value of f.
E.2 Bounding the step size
Definition E.1 (Cross-ratio distance) Let u, E Int(K). If + v, let p. be the end- points of the chord in K which passes through u and such that the four points lie in the order p, u, v.g. Define
|/ 2/12 2 x
1p
4112
ou, v) =
1p
2
1
4112
if + u, and o(u, if u.
For convenience, we define the cross-ratio distance between any two subsets Si, S2 C K as
oSi, S2) min(o(u, Si E S2).
And for any S C Rdd and any density = d + R we define the induced measure:
v*(5) / u(2)d2.
ZEB
Definition E.2 For any posiline-defintte matrix M E IRdd x IRd, and any u E Rdd, we define
||u|| = Vu Mu.
Lemma E.2 For any u, E Int(K), we have
(u, 2
11 v v((a(w)-
/ 2mm -1 + 1---1110 122
19

---

Proof:
Let p. q be the endpoints of the chord in K which passes through u and v such that the four points lie in the order p. u, v, q. Then
u 2112 x p 4/12
2
o² u) =
1p
21/2 X 00 2122
1 2 2 ||u - v 2
12 |p 4 12
u 2
max
|Ip
||u 2
(a. (u 2
|| 2/ 1
max
max (E[m]
(0, u b,)²
|pp 4
1 max
(a. (u 0))
U 21
2
+
2 [c(m]
[a, u - 6,)2
1p q
m (a? (u - u))
||u m1
+ 2
2m i=1
(a b;)2
R²
m
(a:a? T. )2
(u - U)T
2 x 11 TTa Ta (u u)
x a
2ma
(a. 6:)²
2ma 1 + 21) RR (u v) [T)((u(( v)
uu v u1(2(u)-
2ma + 2n
Thus, we have
1
ou, 2
| 1 v)(o(u).
v 2ma 1 + 1-nn-11
R
Lemma E.3 For any u, E (K) such that ||u (ell S 1 2= we have that
(1-a" - ((()) = (uu) 1 (1 + a // - 2 ulll(w))) 2 (z))
Proof:
(a (u - 2))
||u - 1 v1(((u) = a- Σ
+ 7 x//
(a u 6,)²
i=1
(a (u 2))²
2 a
max i6 iE(m]
(a. u b,)2
Thus,
[(. u ba) (a. v - b:)l S a /2/l|/uu v1( la u 6il
Vi E [m].
Therefore, for all Rdd we have
a T
T
(1 a 1/2 lillll-olll - -olla())) )aa
+ 17 'I Id w
w
(a v b.)2
1=1
< w T
a1 1
a;ap
+
1 Id
i=1 (0)10 62)2
S w T
(1 + a /21 w|
12cx
a:a
+1 -1 Ta
k(w))
a
b:)²
Thus,
(1 -aa 1/2 ||u v = (w)) = (1 a 2/2l|u 2/ )() ) 2
20

---

E.3 Bounding the acceptance probability
Lemma E.4 Let 0 In)( Then the acceplance ratão of satisfies
n (=) del (=))
x exp (11 112-01 - -10 110-=llee =Ile(=)) x1
3
(=-NN(00--100))
π T det((()))) (())
10
3
Proof:
By 33 of Lemma E.9 we have
P P2-N(0,--1000
[|2 01 - ||2 01 2
2
2
98
(7)
50
100
By Lemmas
R.5
F.9
and
E.8 we have that
P ((0
TT()) det (+())
x exp (11= 01 - 110
2112(s))
x1
det)
10
((=tt)
1 vdct(D(2))
48
P P2NN(0.0- ((01
dee
50
n
n exp (11= 01
100 - 2 =
0.96
n ((z=ek k}
({770) (1)0)
1 1
v /el(((()))
48
c
P. N000(0---10))
2
/det (())
50
u
{ 111=
01 - 110 =1(2(0))
2
u UfzeRy
- 50
( Vdet(((2)
(T)(E)
1 }1) PINN(0.0 1(0)
48
2 N00,0 ((01
(T))
2 2
(del(((1))
50
2
N0000 2(01)
||z - 01b:(0)
- 110
2
PP 1(0) (12 K}9
50
Lemmas
E.8
Eq.
2
2
10
100
100
100
1
3
Lemma E.5 Let 0 E int(K). Then
((2)
99
99
P P=~N(0,0-1(0))
2
(()
100
100"
Proof: Since ((0,4(0)-1) we have that
z 0 (0) bg = 0 (a- (0) + 17 1a) kg
for some E~ N(O,
Since a (1 + 17 La 2 n la, and H0)) and Is are both positive definite, we have that nla - (a H (0) + -
(8)
Thus,
|| 0 11 H(()
+ + ||(() Ta) E
= =VST (a- + 1a)
Eq.
VET
(9)
21

---

Now, since ç Ia), by the Hanson-wright concentration inequality for the X- distribution
40
we have that
P(IE t) S بيف
Vt
2d.
(10)
Thus, Equations and
10
together imply that
P( 011 > vTV40d) se
2dd
(11)
100"
Now, in the setting where f is L-Lipschitz, we have
π (2)
-((())))((0) 2 -L11--00,,
π (0)
Therefore,
P. P=-N(10-110))
[(2) 99
P2-N((,0-1(0) (--ill=0012
>
99
(0)) 100
100
100 log( L
P=~N(0.0-1(0)
[|2 - 0112 5
99
99
100
where the last incquality holds by Incquality Moreover, in the setting where differentiable and B-Lipschitz, we have that, since is multivariate Gaussian random variable,
11 since 1j S 104d12
P((=- -1==0) V) 0) 2
If (z ()) O, we have that
f(2) f(0) S (z 0) TV f(0) + B|1 011
< 0 01
Therefore,
(<2) 2
99
P N(=
((A)
100
[T(2)
99
{(2 V s $0))
2 (P=NN(00-1(0)
(0))
100
-012-10011
99
0) (v((0) 50) x PP((= VV(==)))V/1(0- <
2 P==N((.--100)
100
Vlog(2)
(z - 0) (v) 50) -P((=-0) ()(()) 5
||2 0|12 5
v3
99
1
100
X
2'
where the last Inequality holds by Inequality 11 since - 10 and since the Gaussian
distribution is symmetric about the 1-dimensional hyperplane (2 0) f
Lemma E.6 For any 0, E In)(K) such that ||0 - 2||0(0) - 72 we have that
1 1 2 =122000
[[N 4 (-..4~1(())) (()) N(2, t (2)))l I 3dall0 =11ă,00 +
22

---

The proof of Lemma E.6 is an adaptation of the proof Lemma 3 in AI to our Markov chain's "soft" barrier function, and follows roughly along the lines of that proof.
Proof: Since ||0 - (0) S u by Lemma E.3 we have that
(1 - a /211 (1) (0) (1 + a01100 21(010)) (z)
(12)
From Inequality
12 we have that
(1 (1-a a //2 0 =||0(0)) La = (0) (1+ a'2110 =1(0(0)) 2d.
Thus, denoting by (M) the ith-largest eigenvalue of any matrix M E R x R we have
by Inequality (12 that
(1-a 1310 - = S
=| )
Vi [d. (13)
Now the KI--divergence between any to multivariate Gaussian distributions with any means M1, M2 Rd and any covariance matrices Σ1, Σn E Idd satisfies (sce c.g. Section 9 of 10 or
Fact in AT)
DKL. (N(u,, E), N (u2, E))
(dt((En) deet(E2)
+ (m (2) T n2)) -
2 Tr( - d + log
(14)
Therefore we have that
Pinsker's Inequality
[0, + (0)) -
20 (N(0, $ (0)), N(z, $ (=)))
Eq.
14
Tr(O(0))- - d log
det( det())1(00)
+ + - =T)())()00 - =)
dot(
-[2))
dct((-1(())
Tr( (2)000)))) - d + log
(let((-1(2))) +(0 2) T (() - 2)
- 2
1
1
1
1
2 + log
2 110 =
det( (2)00()) 2
Σ
\x 3 (2000 (2)) + log
A: (=))))))))))=+( (2)) ))
lle 21%
2
\1( 3 (=) (3) (x))
2
+ 2 00 = =112(00
Eq.
max
1+ - 2
+ 110 =1(300
L8 (1- slloto)))
1
+ 1
max
(1 1)2 +
(1 t)2 - 2 2)
0 =1
[E| =1looo0 a -//2110==1
2
IM
1
161
max
612
+ 2 10 =1|2(00
[E| a0t-0-/ = =loo00 a//2110- =11oro])
d
convexity S
1
6all0- - =21(2(0)
+ 10100 0- =1(ō,0)
= 3da||
= 210000
_1100 = =12
where the first inequality is Pinsker's inequality, the second inequality holds because
log(t) < } 1 for all t>0, the fourth inequality holds because (1 + t)2 + (I-15z - 2 < 612
for all / E |- 1. }) and the fifth inequality holds since 1 is convex for l E R.
23

---

For every j E N, let m, = m + lan j2]d. Consider the matrices A = [AT, Id] where A is concatenated with mh=- copies of the identity matrix L And consider the vectors = (bT, "=1 ........'' ), where b is concatenated with - copies of the vector ji, where 1 = (l...., 1)T E R is the all-ones vector. Then K = {8 E IR A - b)) and the Hessian of the corresponding log-barrics functions is
m c((1))
Hy(w)= = ((a)) b,)2
Lemma E.7 For all w E int(K) we have that
lim Hz(w) ab(w), 10000
(15)
uniformly in w, and that
lim 5
(
(16)
uniformly in w. Moreover, for any 0 E inl(K), any 2 E 1 Do and any sequence {zi15=1 Da such that lim 2j uniformly in 2, have that lim Hz(2;) - ( 500
S
(17)
uniformly in and that
det( 5000 (0))H;(0) det (e)))
det (et(H,(z)))
lim
(18)
uniformly in
Proof:
Hyrww
0;0;
Lann 102]] Σ
(a. w
6,)2
(e] w - j)}
d e.c,T
= ((w) + Lanj 1/21
Σ
(19)
(e] w - j}
1=1
Now, since w E K S B(0, Rt,, we have that
(R j)² S (e. w - ;)² $ (-R- ;)3.
(20)
Thus,
[or) lam-11;3] 1/22 lim
Lanj aan--13²
lim
(1)
(21)
J j-000
(e. w - 5)²
3 000
R - j)B
and
lom -122 lim
Lan 112
lim j
a3)
(22)
1-200 (e: w - و)}
(-R R - j)
Thus,
19
and 22 together imply that
Cie lim 11,-
lim H2(w) H(w) + 1+000
Σ
fo
(e] w - j)z
1=1
H(w) + Σ ege
az7
-H(w) + am -1
adD(w).
(23)
24

---

where Inequalities 21 and 22 imply that the convergence to the limit is uniform in w.
This proves Equation Moreover, since {=x [= S 2 Do, and Do C K, we also have that
|a - bil 2 2 | 0 b,/.
(24)
Therefore,
Lan aan-11 ]
lim H,(2)) = (2;) + Σ 500
lim
10000
(e, 2j j)²
i=1
0.0
d
aan lan- 12a
ese, x lim
0 23
b,)2
1000
(c, 23 5)²
e x an-
[g 2j
b,)2
Eq.
24
m
a:a;
d Σ
city
x on
(a. 6.)2
H() + azj IId
a=(22)
where the convergence of the limit in the fourth equality holds uniformly in 2 by 24 and the fact that lim 2j = 2. Thus, we have that
lim Hz(22) (2)
(25)
uniformly in 2. This proves Equation
17
Moreover, since the determinant is polynomial in the entries of the matrix, Inequality implics that
23
lim dot(H, (w))
(26)
10000
uniformly in
We also have, plugging Inoquality
20 into Equation 19
Lan 1an-113 (R
|12
H(w) + E x
Σ
3 Hy( = H(w) +
c x (-R R j}
j]
1=1
i=1
and, hence, that
H(w) Lan (am-1j²) (R j)²
122
Lar lam- Ia
T Thus, Incquality implics that
Ia U T H2(W)u S T H(w) +
Yu Rd.
(-R j)²
(27)
H(w) + lon 1j2]
Ia = H,(w)) 3 H(w)+ + [an Ia YjEN.
(28)
(R j)²
(-R j)²
Thus, Inequality 28 implies that
H(w)= + lam 1am1 la < (H,(w)) 3 [(())=
aan 122
-1
Is
VjeN. (29)
(-R [R-jjy j²=
(R - j)2
Thus, since H(w) is posilive definite, Inequality imply that
29
logether with inequalities and
22
lim (Hr(w))-1 (H(w) am -aa 0
a (')(w)
Vw E Int(K)
25

---

uniformly in w. This proves Equation 16
Moreover, Inequality
28 implies that
(et z det H(ww) + aU7 }
(R -R2 La
umax H(w)+ + Lan j-]
|
d
(R j)² Ia
[vnj |an
Amaax (H(w)) + Amax
(R- - j)2 Ia)
max (H(w)) 3anj 1)d
ci Vj 3R, for some 61 which does not depend on j.
(30)
Inequality
27 also implics that
Lam- ] (R - j)} In
del (H)(W)) 2 det H(w)
Lan- (R j}² Ld
d
Al
Amin
H(w)=
|an |an /2]
max
Amin (Hin) Amiu
(R - j)² Ia)
3 am
Vj
2 3R,
(31)
for some 02 0 which does not depend either j
Thus, Inequalities that
26
30
and
31
together imply that for any 0, E int(K) we have
lim min 10000
det (z)())
1
Eq.
30
3
min
lim del (H,(H)(z)))
lim (1,((1;(0)) 1)
Eq.
26
25
det(a
min
det(a
=min (et()
1),
de)))
where the limit converges uniformly in This proves Equation 18
Lemma E.8 Let & N(O, Ia) and let 0 E int(K). Then with probability at least 1 100 we have that
1 Do
0 + o² H 1/2
and
|Kll2 < 10vd.
Proof:
Let 2 = 0 + at H - Then since & ~ N(0, Ia), by the Hanson-Wright
concentration inequality
22 40 have that
P(| f) < R_d
Wt 2dd
And hence, since a a-Vaaavv V/1|2 0 []] = |1 (0) (0)Ell2 ||E12,, that
- 0 (() > a /2+) s
2-4
Vt > 2d.
26

---

Hence,
1
P(I 01 (11//
S
100
Ing( 100d
Thus, since S
we have that
1 P E Do) = IP (11z 01li10) 2/
1
1
21
8.
100
Lemma E.9 Consider any 0 E int(K), and _ N(0, Ia).
Let 2 0 (8(0))
Then
IP det (+(2))
48
1 8
2 50
2
(32)
(e)
100
and
2
98
||2 (1 -1 01
I
(33)
50
100"
Proof:
Let 2j = - a H (0)E for all jE E N. Since H2 = H()) for all j E N, we have
that 2j 0 0 11 (0)E E De whenever 0 a (00- 1 Dg. Let E be the event that
Ikll & 10vd and that {=,] S Do. Thus, by Lemma E.8 we have that
99
P P(E) >
(34)
100
Moreover, by Equation 16 of Lemma E.7 we have that lim (0) a0) which implies that
lim 2j lim a² Hj (A)E
1000
10000
0 $ (()E
(35)
uniformly g (and hence uniformly in = 0 x (0)E) whenever the event E occurs.
Therefore, by Equation (18] of Lemma E.7 we have that
det det (zz)) det (x)()
lim
(36)
1 det(H,
det(
uniformly in z whenever the event E occurs.
Since, for each j N, Hy is the Hessian of log- barrberrier function for K, by last line in the proof of Proposition 6 in 411, we have that, for all E (0, }) and all y E (0. I],
log leg (0, 7 H,) -,(0)-1() 10) log det (H,()) -m 71/200((2) () Vd
PENNN(O,Ia)
2 -
Setting L == 8 we have
PENN((..La) (logdel(H, at H;(0)) log del (e,((((())))
1
>
1
1
8.
100
100
since 100g
Thus,
de + at -())
49
1
PEWN((,T
2
2
8.
(37)
det(H,
50
100
27

---

Inequalities 34 and
337
together imply that
({
det [ce((H)(z)) (H,
49
n 2
8.
(38)
50
100
Moreover, Equation 36
implies that there exists some number
E N such that
det( det(HxxzzN)
49 nES
det(((z))
48
nE
(39)
det det(HNx(0)
50
det(
Hence,
Eq.
38
{ det(HNN(0)
det( (eN(HN((N))
49 2
1
P (
n E
100
50
Eq.
39
det(((()))
48
FeoN(o.Laa)
n
(40)
det(((())
50
This proves Incquality 132
Moreover, since, for each j E N, H, is the Hessian of log-barrier function for K, by Proposition 7 in 4T, for all (0, g) all a 0 such that vod I 20 log - and all j E N, we have that
PP (1=) 012) || - 00 2 21a) 51-6 VjEN.
(41)
Thus, Equations and 134 imply that
P {1 |z 01
98
01 H2)
|| - 012,, S 50
n 2
VjEN,
(42)
100
since a S 10'd
By Equation
17 of Lemmna
E.7 Equation
335 implies that
lim Hy(2)) = (2), 5000
(43)
uniformly in 2, whenever the event E occurs. Thus, Equation
43 implies that
lim ||= - (s))) - ||> - 0022,,00 = lim (zj - T (H)H((2)){ - 0) - (z) - 0) ())))(z - 0)
Eq. 43
lim (z 0)T ()) - 0) - (z 0)T (0) 0)
a
01 all - 01
(44)
uniformly in z (and hence in & 2b (0)(z 0)) whenever event E occurs. Thus, Equation 44h implies that there exists a number M E N such that
0ll
1| 0l1
nEW c fall -01 all= 0lla S
2
NE
50
50
(45)
Thus,
98 Eq.
42
IF (11Ea -01 1|EM 0 -012 (0)))(0)) S a (E)
100
50
Eq.
45
aall= - 0112(=) -allz -012,00 - 50
nE
This proves Inequality
33
28

---

E.4
Bounding the conductance
Lemma E.10 For any 0, z int(K) have that
Po P2 IV -
whenever
110 - =0((00)) S
30'
Proof: First, we note that,
P P. I P NO, (0))llTV
|N + (0)) - N(2, (x)) -1: ((2,5-112)) (x))(lv
(46)
By Lemma E.6 the middle term on the RHS of 46 satisfies
[| < (0)) N(z, (=)) V
3do/|0 = =lūx(o) 1 = ())
(47)
Plugging in our choice of hyperparamoter a= 105d Incquality
1
147 simplifies to
53
[0, 1b (0)) N(2, Φ (2))lav $ V
|| 00))
(48)
100
Thus, we can show that | NN 4 (0))llVV S 1/ for all E int(K), we will have that [Pa P2 S 29 whenever |0 =llæ(0) S 1.4 as desired.
To bound the other two terms the RHS of 16 observe that
[Po N(0, (
(49)
where
((0,)) min 722) v /et (z) T(0) det
exp(1I= 011b(0) -110- 100=2112:m =1(2(e)) {(
the acceptance ratio.
By Lemma E.5 we have that IP (==N(0,0-100))
[(x)) ((क)
2
99 100
2
99 Therefore, T0
PI~N(0,0 100))((0,2) 2
10
(2))
99
(dct(((2))
48
22~N(0.0-100
<
P
(0)
100
det(((())
50
P
(ell= 11=-0112.0 110=21 2 < 0.96)
-P=-N(00-000) K)
2
1
2
100
100
100
100
9
10'
where the term with the exponent is bounded by Incquality other two terms are bounded by Lemmas E.8 and E.9
33
of Lemma E.9 and the
Therefore,
Po NO0, + (0))ltv = -E.-N000-
2
(50)
10
10
10
5
Plugging Inequalities
48 and
50
into
46).
we obtain that P Palltv
whenever
||A 2 ())) S 2 as desired.
We recall the following isoperimetric inequality for log-concave distribution on a convex body, which uses the cross-ratio distance:
29

---

Lemma E.11 (Isoperimetric inequality for cross-ratio distance (Theorem 2.2 of 32 Let π Rdd + be a log-coneave densiby, with support on a convex body K. Then for any partition of gº into meosurable sets S S2, S3, the induced measure n¹ satisfies π (Sa) 2 (S1. (S1)
For any 0 E Int( K), define the random variable Zg to be the step taken by the Markov chain in Algorithm 1 from the point 8. That is, set = 0+ ( & where & ~ N Ia). If 2 K,
set Zo with probability min d(t) exp [|| 01ū(()) 10 - = =12a()) det)())
Else,
For any 0, S C Rd, define the onc-stop transition probability of the soft toft-threshold Dikin walk Markov chain to be
Pa(S) = P(Z E S).
The next proposition shows that the soft-threshold Dikin walk Markov chain is reversible and æ is stationary distribution of this Markov chain:
Proposition E.12 (Reversibility and stationary distribution) For any Sr, S2 S K have that
(())Po(S))d0 = -f (S))ds
JES
JEES2
Proof: Let po(2) = (det(((o))
(2p1) e e- ((00)=0)))000-)) for any 0, E Int(K) be the density of
the N(0, ()) distribution.
1 (0)))(S2)d0- = Joesh T() ff pa(z)x
DES
=E
<det((()((((3)
min
exp [|z -1 01 110 - -=1(((s)) dæd0
T(0) (e)
(=))p=(0)
T(0)po min
1 dzd6
(0))po(2) p0))
BES,
JES,
T(0)p= min
(T))p=(0)' d0ds
22S2 JES.
(O))oo
JES
Define the conductance of the Markov chain to be
Φ (9)5 T*()) / Po(K))))7(()d)
inf
Lemma E.13 The conductance o satisfies
10 12ma-1+17 I
Proof:
The proof follows the general format for conductance proofs for geometric Markov chains. Let S, C K and let S2 KYS,, Without loss of generality, we may assume that T(i)) < (since otherwise we could just swap the names "S," and "S²"). Let
1
8
0 S Po(S2)
70
1
8 E S. P.(Si) < 70
(51)
30

---

and let
S/ = (K\S)))S
By Proposition 12 we have that
((0)Po(S2)dB
JoES2 (()))
(52)
Jess
Moreover, by Lemma E.10 for any u, Int(K) we have that
29
1
Pa Pollrv S
whenever ||u- vl) $
(53)
30
2
Thus, on the one hand, by Lemma
E.2
Inequality
53 implies that for any u, 0 E Int(K),
29
1
[P, - PolITV S
whenever
a ou, v) 2/2ma 1 R2
(54)
30'
On the other hand, Equations
51
imply that, for any Si, E S% we have that
68
Pr, - Po TV 2 Pa - Pa(Si)
(55)
70
Thus, Inequalities [54] and 55 together imply that
σ (S S2)
(56)
2 2m2x + 112
Moreover by Lemma E.11 we have that
n" (S(g)) 2 (z)
(57)
First, we assume that both n² (St) 2 (S,) and x*(a)) 1 (S²).. In this have
/ (o(S))=(B)de
Eq. 52
(o(S)) + 2 2 Po(S))=(O)de
2 S,
140 " (S')
Eq. 57
1
O S/) (S))
140
Eq. 56
1
1
(8))=((()
280
22na + 1Rk R
1
1
min (S)), T*(s))
560
2ma + R2
2024 /2ma + 11-11 min (T" (S1), **(2)).
(58)
Now suppose that instead oither T'(s) < 1r1) or T' (5)) < (S2). If (S') < (Si) then we have
J (o()))=(0)d0
Eq.
62
(o((2) + J. (o(S))=(0)d0
2 Si
1
Po(S2) (0)d0
2
S.VS
5,
Eq.
61
3
34
π (Si)
70
1
2
min(=*(s))
T (S2)).
(59)
10
31

---

Similarly, if T" (S2) < T (S2)
we have
/ S. (o(S))=(O)dB
Eq.
52
2 2.s.
/
Po(S))T(O)dB + Js.
Po(S))=(0)dB
(o(S))T(O)dB
Eq.
51
34
(S2)
70
1
2
min(s=((1)), n (S2))
(60)
10
Therefore, Inequalities 58
59 and
60
together imply that
min( n² (S2)) / Po(52)=(0)dl) 10 2mma-1 + - R2
(61)
for every partition S S2 = K. Hence, Inoquality 66 implies that
inf Sc (8) / (o()))))T(0)d
1
60 2ma +1 R
Definition E.3
We say that a distribution v is w-warm for some w 2 with respect to the u(z)
stationary distribution ı if supseK T(3) S w.
Lemma E.14 (Corollary 1.5 of [11) Suppose that ll0 is the initial distribution of lazy reversible Markov chain with conductance & > and stationary distribution. n, and let μr be the distribution of this Markov chain after i 2 0 steps. Suppose that Ho is w-warm with respect to n for some w 2 1. Then for all 2 we have
62
l| - x S V 1
2
is either L-Lipschitz (or has B-Lipschitz gradient). Suppose that do po where po is a w-warm distribution with respect to π xc e- with support on K. Let μ denote the distribution of the output of Algorithm ] with teperpprrmmeerss Io and ToaL (or nj S r if it is run for T iterations. Then T 2 109 +n IR²) x log( =)) we have that
Lemma E.15 Let 8 O. Suppose that f:
|| - TV S 6.
Proof:
By Lemma E.13 we have that the conductance & of the Markov chain in Algorithm
satisfies
2 104 2ma -1 + 1
and hence that
= 10 (2ma + + 17- 182)
log( ) 20- log( w
x)
Thus, by Lemma
13 we have that
132
|e T/ITV S V/
1
2²
2Φ
log( ))
w
V WE
log( )
6.
Proof: [of Theorem 2.1
32

---

Total variation bound. Recall that we have set the step size parameters = 106d and
either 77 T0-ar3 (if f is L-Lipschitz) or n 10423 (if is 8-smooth). Thus, after rumning
Algorithm I for T = 10 (2ma-1 1 + 17 - x log( () iterations, by Lemma E.15
we have
that the distribution j of the output of Algorithm satisfies
π TV S.
(62)
Bounding the number of operations. Moreover, by Lemma E.1 we have that each iteration of Algorithm 1 can be implemented (mmd=11) arithmetic operations plus O(1)
calls to the oracle for the value of f. Thus, the number of iterations T taken by Algorithm is (((mmddd²RR)) x l))( iterations in the setting where f is L-Lisschitz and O((md
d3R²) x log( )) iterations in the setting where f is 3-smooth, where each iteration takes one function evaluation and md= arithmetic operations.
F Extension to general barrier functions?
From any point 0, our algorithm proposes step with Gaussian distribution (0, (V where g is the following barrier function
g(0) = (0) +
(63)
where (00) is a barrier function for the polytope K, and the parameter f- 22 in the setting where f is guaranteed to be L-Lipschitz and i = (2)) in the setting where f is guaranteed to be f-smooth. Thus, most of the probability mass of the Gaussian distribution concentrates in the Dikin ellipsoid Dd = 0 + {w : w (V)) -' S 1} for the barrier function g. For simplicity, in our main result, we assume that (o()) == -, log(b, - a (which has self-conccrdaace parameter m), however, we can in principle choose any barrier function to for the polytope K, such as the entropic barrier 4 the Lee-Sidfer Sic-Sidford Barrier 29 which have self-concordance parameter roughly v =dd.
To arrive at our barrier function from a more axiomatic approach, we first consider the definition of v-self concordant barrier function:
Definition F.1 (j-self-concordant barrier function for K) We soy that is a v-self- concordont barrier function g Int(K) +R IR and g satisfies the following conditions:
1. Convex and differentiable barrier: is convex and third order differentiable, and g(z) - +00 as + ŌK.
2. Self-concordance: Vg(z))h, S ((V"g(=)[,,h))) for all h Ed (this en- sures that the Hessian of the barrier function does not change too much each time the Dikin walk takes step from x to 2, that is, [V2((=) g²g(z) g2))
3. g is v-self concordant: h (gVa)) S vuh V2g(z)h for every x E [nt(K), h Rd
The fact that our barrier function 63 satisfies parts (1) and (2) of Definition
F.1 follows from the fact that & satisfies Definition F.1 and that (² 0) 0. We discuss the self-
concordance parameter v for which our barrier function satisfies 63 below.
To ensure that the steps 2 ~~ N(O, ( 2 ((0)==) proposed by Dikin walk Markov chain arising from the barrier function g has an 2(11) acceptance ratio e we require that the function g satisfies the following property. This property says that at least of the volume of the Dikin ellipsoid Do contained the sublevel set {= e J(() < f(0) + 2} where the value of ʃ does not increase by more than 2.
Property F.1 Dikin ellipsoid mostly contained in sublevel set: At every 0 E Int(K), the Dikin Ellipsoid, salisfies
Vol D Donf :())) = f())
23) 2 1
Vol(DD).
33

---

When designing a barrier function 9. there is a trade-off between choosing such that the self-concordance parameter V is small, while at the same time ensuring that Property F. holds. Roughly, to make the parameter y as small as possible, we would like the Dikin ellipsoid to be large relative the Hilbert-distance metric for the convex body K (This is because, by Proposition 2.3.2(1ii) of 39. any v-self concordant function 9 satisfies for any h E I& where |h/ := supfa 0: ath E {= E K)). On the other hand, if we make the Dikin ellipsoid too large (with respect to the sublevel set {= E K : f(z) < (xx + 2]) then Property F.1 will not be satisfied, and the steps proposed by the Dikin walk will have a very low acceptance probability. Setting + when f is I-Tipschitz or r when f is B smooth ensures that our barrier function in (63 satisfies Property
(h (z)))) - < |/l S (1+ V29(z)h)
the hyperparameter
In the special case where f is the log-density of the uniform distribution =f=0) or when f any linear function, we have that f is 8-smooth for 8 = O. Thus, our barrier function [63) we use to encode the geometry of f on K is the same as the barrier function & for the polytope K. This is because, since the level sets of linear functions f are half-planes, any ellipsoid centered at & has at least half of its volume in the sublevel set {z E K f(=) S f()) S f= E K f(z) < f(0) + 2), satisfying Property F.1
The following lemma shows that our barrier function in 663 is u-self concordant with v = Olv' + +n)),, where s' is the self-concordance parameters of 4p (which is v = -m if we choose to be log-barrier function). Thus, our barrier function in 633 is v = O(v + L R²) self- concordnnt in the setting where f is L-Lipschitz, and v = 0 BR in the setting where f is 8-smooth:
Lemma F.2 Suppose that Φ(x) is a v -self concordant barrier function for a convex body K C Rd where B(0,r) S KG B(0, R) for some R > r > 0. Let g(x) = o(x) + qx' x for some a > 0. Then is u-self concordand for v = 4v + 4.
The polynomial dependence of v on LR or 312 is a necessary feature of any objective function satisfying Definition [.1 and Property F.I In appendix F.2. we give explicit examples of classes of objective functions f and polytopes K for which the minimum value of depends polynomially on LR (and classes of smooth functions f where v depends polynomially on VBR.
The proof of Lemma is given in Appendix
An open problem is whether one can design versions of the Dikin walk which sample from log-concave distributions with a runtime that depends only on the dimension d, and is independent. of L., R or B and which are invariant to affine transformations. The difficulty in achieving this lies in the fact that (e.g., in the setting where I is L-Lipschitzz, on the one hand, the level sets of f where most of the probability mass of x e concentrates may have diameter roughly RL times smaller than the diameter R of K. Thus, to have an Ω(1) acceptance probability the Dikin walk may need to take steps that are roughly RL times smaller than the diameter of K. On the other hand, the isoperimetric inequality currently used to bound the mixing time of the Dikin walk uses a metric- the cross-ratio distance for K -w-hhch, roughly speaking, defines distances between steps by how quickly these steps approach the boundary of K. Thus, measured in the cross-ratio distance, the steps proposed by the Dikin walk are of size roughly proportional to πτ, and mixing time bounds obtained with this isoperimetric inequality thus depend polynomially on RL. To obtain mixing time bounds independent of the Rt L, B one may need to show a new isoperimetric inequality which is based on a different metric that encodes the geometry of all of the level sets of f-rather than just the geometry of its support K.
F.1 Proof of Lemma F.2
Proof: [Proof of Lemma F.2
For any x E int(K), we have
Vg(z) = VΦ(z) + ax
Thus, for any h E Rd.
A [V = h Vo(a)) + ah x S h Vo(x)) + ||b||||=|||| S SA Vo)()) + R|AA
(64)
34

---

Thus,
Eq.
(h T Vg(x)))
(h (VVo)))+aRl)
(65)
4(h 44hh"
4(h Vo(==))+ 44 Lah
=4(h" V
(ala)h
44/h vlo 4alt (alaj)
(4v' + 4a/t") (b (alaa)) (ala)h)
(4v' + 4ul (v +alahh
(4v' + 4allt) (
where the third inequality holds since & is v' self-concoddant. Thus, plugging in 4o' + get that
4oR2 to 65
Vg(z)) & (u(hT(V'g())))
F.2 Lower bounds for self-conncordance parameter
Lower bounds on for worst-case L-Lipschitz f. Consider the L-Lipschitz function
f(0) = L| constrained to the convex body K 2v 2val [ which is contained in ball of radius R. In this case, any barrier function has Dikin Ellipsoid Da which is contained in the ball B(0, $). Thus, we have that 2 M(LL). (This is because, by Proposition 2.3.2(iii) of 89 any v-self concordant function y satisfies (h (²'g(0)))= < [hla S (1 + +3v)(h (²g(0)h) -3, for any I&d where Ihle- := sup(a > 0: 0 + ah E {= E Kj. Thus, at 0 = 0 and choosing h = (1,-.", ,1) we have Ih a. = R and (h T ('g())))= < ) Thus, // 2 Ω (A (BTT2g(0)6)
R
which satisfics Property F.1
Ikle
LR).
[hlla
Since thore exists a convex body for which any self-concordant barrier function satisfy- ing Definition F.1 has self-concordance parameter at least 2 d, for any L, 0 there exists a function f and convex body K C B(0, R) such that the self-concordance param- eter of every barrier function satisfying both Definition F.1 and Properly E.I is at least v 2 [(max(d, LR)).
Lower bounds for worst-case 8-smooth f. Consider the f-smooth function f(0) 38 constrained to the convex body K Rt [-,, At 0 0, any ellipsoid Do
2/a
satisfying Property F.1 is contained in the ball
8
((0 1). Thus, we have that > R vB
VB
((V3R) (This is because at 0 = 0 we have
Whle h
= R and (hTT(g(())))) < 0 and
thus, by Proposition 2.i.2(iii) of 39 we have Ω
hls
2 RR). Thus, for
(AT
any R > there exists a function f and convex body K C B(0, R) such that the self- concordance parameter of every barrier function satisfying both Definition F.land and Property F.1 is least 2 max(d, 2(RR)).
35


Reviews:
[{"review_id": "T51hDlR0jt", "reviewer": "Reviewer_5ZqR", "summary": "The paper studies the question of efficient sampling from a log-concave distribution from a constrained polytope. The main idea is to use the well known Dikin's random walk with a soft threshold. The threshold is important to ensure a far more efficient convergence to the desired distribution. This soft threshold ensures that the acceptance ratio is high while also ensuring that the walk remains inside the polytope. ", "questions": "N/A", "limitations": "N/A", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The paper gives a definite improvement over the best-known algorithm in the run time for both Lipschitz and smooth loss function. This has direct implications in making DP algorithm for convex optimization more efficient. Just on the strength of the result of the algorithm, I favor acceptance. I have read the proof and did not find any issue with any of it. ", "weaknesses": "The only weakness of the paper is that the writing can be improved. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable comments, and for taking the time to review our paper. We are glad that you appreciate the strength of our results and we thank you for supporting our paper."}, {"review_id": "Gw8jvJcLt9", "reviewer": "Reviewer_wurv", "summary": "This paper studies the problem of sampling from distribution of the form $\\pi(\\theta) \\propto e^{-f(\\theta)}$, restricted to a polytope $K$. Here, $f$ is either Lipschitz or smooth convex. To this end, the authors propose to use _Dikin walk_ Markov chain (Kannan and Narayanan, 2012), which was originally proposed as a sampler for the uniform distribution on polytopes. Given that $K = \\{\\theta : A\\theta \\leq B\\}$ with $A=(a_1,\\ldots, a_n), B=(b_1,\\ldots,b_n)$, the Dikin walk proposes the update $z = \\theta + \\sqrt{\\alpha H^{-1}(\\theta)}\\xi$ where $\\xi\\sim N(0,I_d)$, $H$ is the Hessian matrix of the log-barrier function $\\varphi(\\theta) = -\\sum_{j=1}^m \\log(b_j-a_j^{\\perp}\\theta)$, and $\\alpha$ is a hyperparameter. If the proposed update is in the interior of $K$, it is accepted with a certain probability, otherwise it is rejected. To ensure that the acceptance probability is $\\Omega(1)$ while allowing the Dikin walk to make sufficiently large steps, the authors propose to add the scaled Hessian $\\alpha H(\\theta)$ by a regularizing term $\\eta^{-1}I_d$, in order to \"round up\" the set of proposals; then, we can adjust $\\alpha$ and $\\eta$ so that the set of proposals can be fitted in $K$ with high probability. The authors have done a runtime analysis and show that their method can be run in $\\tilde{O}(md^{\\omega+1})$ in order to sample with small TV error, where $m$ is the dimension of the polytope and $\\omega$ is the matrix multiplication constant. In particular, when $m = O(d)$ (e.g. $A$ is full rank), it can be run in $\\tilde{O}(d^{\\omega+2})$.", "questions": "- There are some notations that should be introduced in the Notation section. For example, $a_j$ and $b_j$. $\\Phi(\\theta)$ and $\\Phi(z)$ are used at the beginning of Section 3 but I could not find their definitions anywhere before this section.\n- Maybe I have missed it, but there should be a mention at the beginning on a lower bound of $R/r$ so that it is easier to make comparisons of the numbers in Table 1.", "limitations": "The authors have mentioned a limitation of the proposed method that it can only be applied to Lipschitz or $\\beta$-smooth $f$. Also, finding a better regularizing term is also a possible research direction.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "- The authors did an excellent job on literature review, with complete runtime analysis of previous methods. The time for warm-start and the evaluation of $f$ are also taken into consideration.\n- The authors provide detailed explanation on Dikin walk and its limitations.\n- The algorithm seems simple enough to implement with a few lines of code.\n- The authors provide some motivations for sampling on polytopes in the area of differential privacy.\n- I appreciate the overview of the proof of the sampling guarantee in Section 3.\n\n", "weaknesses": "- The authors mention that the hyperparameter $\\alpha$ alone is insufficient for the Dikin walk to make large steps due to the differences in the geometry of the proposals at each point $\\theta$. So I think it is better to add a regularizer that depends on $\\theta$ that allows us to remove the hyperparameter $\\eta$. Have the authors considered such approach? (In other words, can we somehow make $\\eta$ depend on $\\theta$?)\n- The comparison of the runtimes are nice, but I think some experiments on simulated data would make it more convincing that the proposed method is better than the previous ones.\n- The authors might want to also discuss the following paper by Chalkis, Fisikopoulos, Papachristou and Tsigaridas: \"Truncated Log-concave Sampling for Convex Bodies with Reflective Hamiltonian Monte Carlo\".\n- In the introduction, the authors need to be more explicit about the matrix multiplication constant $\\omega$. Outside readers might not know that $\\omega \\leq 2.373$, making it harder to compare the numbers in Table 1.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable comments and suggestions. We are glad that you appreciate our proof overview, and thank you for supporting our paper. We answer your specific questions below.\n\n**\"regularizer that depends on $\\theta$...\"**\n\nThank you for the suggestion.  While there are settings where a regularizer that depends on the position $\\theta$ may lead to a faster runtime, this would likely require additional access to the function $f,$ or different assumptions on the structure of $f$, beyond what we assume in our paper.\n\nFor the class of functions considered in our paper, we conjecture that an $\\ell_2$ regularizer which does not depend on $\\theta$ is optimal. This is because we consider the class of functions $f$ which are $L$-Lipschitz or $\\beta$-smooth with respect to the $\\ell_2$ norm, and our bound on $L$ or $\\beta$ does not depend on $\\theta$. Moreover, we only have access to the function $f$ through an oracle which returns the value of $f$ at any given point $\\theta$, but does not tell us how $f$ changes at nearby points.\n\nWe will add a remark about this in the final version.\n\n\n**“The authors might want to also discuss the following paper…”**\n\nThank you for pointing us to this reference, we will discuss it in the related work section.  In particular, we note that the bounds in [Chalkis, Fisikopoulos, Papachristou, Tsigaridas, ACM Transactions on Mathematical Software, 2023] assume that $f$ is $M$-strongly convex for some $M>0$ (we assume $f$ is convex, but not necessarily $M$-strongly convex for $M>0$), and thus are not directly comparable to our bounds.\n\n**“matrix multiplication constant $\\omega$.”**  We will add the value of $\\omega$ to the introduction.\n\n**“notations…”** Thank you for pointing this out. We will add these definitions to the notation section.\n\n\n**“lower bound of $R/r$...”**\n\nAny convex body satisfies the lower bound $R/r \\geq 1$ ($R/r =1$ for a ball, and one can construct a polytope, which approximates a ball, for which $R/r$ is arbitrarily close to $1$).\n\nMoreover, for any convex body $K$, there always exists a linear transformation $T$ for which $TK$ is contained in a ball of radius $R$ and contains a ball of radius $r$ such that $R/r$ satisfies the upper bound $R/r \\leq O(\\sqrt{d})$ ($TK$ is referred to as a “well-rounded” convex body).\n\nWe will add this information to the table caption."}, {"review_id": "6iNQW30dWP", "reviewer": "Reviewer_ACUh", "summary": "The paper proposes a sampling method for log-concave distribution with bounded polytope. The methods can be used in privacy preservation. The theoretical results show the guarantee of the sampling methods with some accepted error magnitude. The comparisons for other methods are also detailedly proposed, showing the improvement of the proposed one.", "questions": "1. can you provide some insights about how your algorithm can speed up the sampling. I just think that is because compared with the Dikin walk, the designed regularizer supports fast computing, but I am not sure.\n2. Also, I can not figure out how much improvement for the dimension $d$ in the real problem. Some small examples can be provides to both show the constant in the magnitude and $d$'s improvement in some cases.\n3. I am not sure about the adjustment of the line height is allowed.\n4. Maybe we also need to consider the time cost for computation of the Lip/smooth constant for $f$", "limitations": "Yes", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 1, "contribution": 2, "strengths": "Pros:\n1. The work focuses on an important problem in the privacy field--sampling and also provides some ways for us to adopt their methods.\n2. The theoretical results show the effectiveness of the proposed methods in magnitude.\n3. Some interesting discussion and delicate application in the appendix is provided for better practical usage.", "weaknesses": "Cons:\n1. Work in the main context is not friendly to the readers. Maybe authors can use more paras to provide different information. The current version is not so friendly for me to read at least. Some more important Lemma/Corollary can also be organized, I believe. I think the design of the func for the barrier should also be clearly stated at first, but not in the Algorithms.\n2. Some technical contributions can also be provided to better differentiate yours and the Dikin walk with simply applying a new regularizer.\n3. Some numerical studies can be provided to better show the effectiveness of the methods.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable comments and suggestions. We are glad you appreciate our theoretical results and the applications to differential privacy, and thank you for supporting our paper. We are sorry for any difficulty understanding our presentation. We answer your specific questions below.\n\n**“technical contributions can also be provided…”**\n\nOur work makes the following novel technical contributions (discussed in Section 3).\n\n**(1)** We introduce self-concordant barrier functions which simultaneously take into account the geometry of both the constraint polytope $K$ and the Lipschitz or smoothness property of the target log-density $f$.\n\n**(2)** The main technical challenge in bounding the mixing time of our Dikin walk is to prove that the determinantal term in the Metropolis acceptance probability $\\frac{\\textrm{det}\\Phi(z)}{\\textrm{det}\\Phi (\\theta)}$ is $\\Omega(1)$ with high probability (w.h.p.), where $\\Phi$ is the Hessian of our regularized log-barrier function. In previous works on the Dikin walk, which use the log-barrier without regularizer, the determinantal term can be bounded using the following inequality of [Vaidya, Atkinson, 1993] which holds for the Hessian $H$ of the log-barrier:\n$$(\\nabla V(\\theta))^\\top[H(\\theta)]^{-1}\\nabla V(\\theta)\\leq O(d)\\qquad\\forall\\theta\\in\\textrm{int}(K),$$\n\nwhere $V(\\theta):=\\log\\textrm{det}H(\\theta)$. Unfortunately, this inequality does not hold for every self-concordant barrier function. We prove that it *does* however hold for our regularized barrier functions. To see why, we first show that the Hessian of our regularized barrier function, $\\Phi(\\theta)=\\alpha^{-1}H(\\theta)+\\eta^{-1}I_d$, can be viewed as the limit of an infinite sequence {$H_j(\\theta)$}$_{j=1}^\\infty$ of matrices, where each $H_j$ is the Hessian of a log-barrier obtained by representing $K$ by an increasing set of (redundant) inequalities. Roughly, this allows us to show that the above inequality, which holds for any log-barrier, must also hold for our regularized barrier function (if we replace $H(\\theta)$ with $\\Phi(\\theta)$ in the above inequality and $V(\\theta)$ definition).\n\n**“…insights about how your algorithm can speed up the sampling…”**\n\nThe regularizer in our algorithm speeds up the runtime by allowing the Dikin walk to take larger steps, while still ensuring these steps are accepted w.h.p. by the Metropolis accept/reject rule for $f$. Taking larger steps allows our Dikin walk to converge more quickly to the target distribution $\\pi\\propto e^{-f}$.\n\nTo see why, note that from any point $\\theta$, the original Dikin walk proposes updates $z=\\theta+y$ where $y$ is normally distributed with covariance matrix $\\alpha H(\\theta)^{-1}$, where $H$ is the Hessian of, e.g., the log-barrier function for $K$ and $\\alpha>0$ is a hyperparameter. If one applies the Dikin walk to sample from a non-constant distribution $\\pi\\propto e^{-f}$, one needs to ensure the stationary distribution it converges to is equal to $\\pi$. This can be done by accepting each proposed step with probability proportional to the Metropolis rule $e^{f(z)-f(\\theta)}$ and rejecting it otherwise. However, the acceptance probability may be very low (e.g., if half the eigenvalues of $\\alpha H(\\theta)^{-1}$ are $>\\frac{c}{dL^2}$ for some $c>\\Omega(1)$, the acceptance probability may be exponentially small in $c$). \n\nOne approach (used in [Narayanan, Rakhlin, JMLR 2017]) to ensuring the acceptance probability is high is to choose a smaller $\\alpha$ such that all the eigenvalues of $\\alpha H(\\theta)^{-1}$ are $\\leq\\frac{1}{dL^2}$, which ensures the acceptance probability $e^{f(z)-f(\\theta)}=\\Omega(1)$ w.h.p. Unfortunately, this approach can lead the Dikin walk to propose steps with covariance matrix that has many of its eigenvalues unnecessarily small. This is because some eigenvalues of $H(\\theta)^{-1}$ may be much larger than other eigenvalues (e.g., if $K$ is much wider in some directions than in others).\n\nTo overcome this, our Dikin walk proposes steps with covariance $(\\alpha^{-1}H(\\theta)+\\eta^{-1}I_d)^{-1}$, where we set $\\eta=\\frac{1}{dL^2}$ (and set $\\alpha$ to the same value $\\frac{1}{d}$ used in prior works which apply in the special case where $f$ is constant). This ensures the largest eigenvalues of the covariance matrix are no larger than $\\frac{1}{dL^2}$, *without* reducing (by more than a constant factor) the eigenvalues which were already $\\leq\\frac{1}{dL^2}$.\n\nWe will add this discussion to Section 3.\n\n**\"…Some small examples…\"**\n\nThank you for the suggestion. We will add one or two concrete examples, in addition to the examples given in Section 2 and in the last two columns of Table 1, to better illustrate the runtime improvement.\n\n**“…computation of the Lip/smooth constant”**\n\nOftentimes, a bound on the Lipschitz or smoothness constant can be calculated analytically. This includes, e.g., applications to training Bayesian or differentially private logistic regression models (or other generalized linear models such as support vector machines (SVM)). In these applications, $f(\\theta)=\\sum_{i=1}^n\\ell(\\theta^\\top x_i)$ where $\\ell:\\mathbb{R}\\rightarrow\\mathbb{R}$ is a convex loss and {$x_1,…,x_n$} $\\subset\\mathbb{R}^d$ is a dataset. The loss $\\ell$ may be $O(1)$-Lipschitz (e.g., if $\\ell$ is the logistic loss $\\ell(s)=\\log(1+e^{-s})$, or the loss $\\ell(s)=\\max(0,s)$ used to train SVMs) or may be $O(1)$-smooth (in e.g. logistic regression).\n\nWhen the Lipschitz or smoothness constant is not known, one can in practice set our algorithm's hyperparameters by hand such that the average acceptance probability is, e.g., $>\\frac{1}{2}$. One can then run the Markov chain until an easily-computed heuristic convergence metric (e.g., the autocorrelation time) is lower than some desired value (see, e.g., [Durmus, Moulines, “High-dimensional Bayesian inference…” Bernoulli 2019], who use a similar approach to choose hyperparameters of a different Markov chain).\n\nWe will add a remark in the final version."}, {"review_id": "X9fefcYIIx", "reviewer": "Reviewer_D286", "summary": "This paper studies the problem of sampling from a log-concave distribution, parametrized by $\\exp(f)$ where $f$ is $L$-Lipschitz, over a $d$-dimensional polytope specified by $m$ halfspaces. Assuming the polytope is contained in a ball of radius $R$, the authors develop an algorithm that generates a sample from a distribution at most $\\delta$ TV distance away from the target distribution, in time $O((md+dL^2R^2)\\cdot md^{\\omega-1}\\log(w/\\delta))$ where $\\omega\\approx 2.37$ and the initial point is $w$-warm start. Their approach is an extension from Dikin walk by introducing a soft-threshold regularization term. This result has a wide range of applications in differential private optimization. \n\nIt is worth noting that Dikin walk was previously adapted for sampling from a uniform distribution over a polytope, and the iteration count and arithmetic operations per iteration are nearly tight [LLV20]. For sampling from a log-concave distribution, hit-and-run has its runtime depends on the dimension $d$, but under certain regimes, it is slower than the proposed algorithm in this paper. \n\n[LLV20]: Laddha, Lee and Vempala. Strong self-concordance and sampling. STOC'20.", "questions": "See weakness. I'm willing to raise score if authors can address my concerns.", "limitations": "Yes.", "rating": 7, "confidence": 5, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "Under many regimes, this paper provides the state-of-the-art result for sampling from a log-concave distribution over a polytope. Moreover, their algorithm is conceptually easy to understand and should be easy to implement, as it is a variant of the standard Dikin walk for sampling from a uniform distribution. The soft-threshold regularization term is not completely novel, but prior work gets a highly sub-optimal dependence on the self-concordance parameter $\\nu$, at least for log-barrier. \n\nConsequently, this algorithm also leads to many applications, such as differential private empirical risk minimization and Bayesian Lasso logit regression.", "weaknesses": "There are several weaknesses of this paper to highlight.\n\n1. Restriction to log-barrier. Note that the $md$ iteration bound comes from the log-barrier is $m$-self-concordance. It is well-known that if $m\\gg d$, better barrier functions with better parameters exist, such as the hybrid barrier and Lewis weights barrier. However, if one carefully examines the techniques employed in this paper, they are tuned very specified towards log-barrier. In particular, the authors prove that in the limit, sampling from the regularized Gaussian distribution over a polytope whose constraint contains infinitely many copies of identities will converge to sample from the uniform distribution over another polytope. This enables them to utilize prior rich works on Dikin walk for uniform sampling over polytopes to conclude the proof. Unfortunately, this approach does not scale with more complicated barrier functions. Take the volumetric function as an example, the Hessian itself reweights all constraints based on their leverage scores, and duplicating a constraint infinitely many times will make the leverage score goes to 0, invalidating the construction completely. I think it will be much more interesting if the proposed framework can generalize to other barrier functions.\n\n2. Tightness of arithmetic operations. The algorithm takes $O(md^{\\omega-1})$ to form the Hessian matrix. This does not seem to be tight. I think the results in this paper will be much stronger if, at least for log-barrier, the arithmetic operations per iteration are tight (say nearly linear in $md$). \n\n3. Dependence on $L$ and $R$. The Dikin walk is inherently a second-order method, so one would expect a better dependence on $L$ and $R$, hopefully polylog dependence. Notably, the Dikin walk for sampling from uniform distribution over a polytope does not depend on $R$, as the walk exploits the local geometry of the Dikin ellipsoid. Hit-and-run can achieve such a result only if the body is in isotropic position [CE22]. The quadratic dependence on $R$ weakens the result significantly in my opinion. Authors should discuss whether the dependence on $R$ can be improved.\n\nDespite the above concerns, I think this paper is technically solid and has many interesting applications.\n\n[CE22]: Yuansi Chen and Ronen Eldan. Hit-and-run via localization schemes. 2022.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable comments and suggestions.  We are glad that you appreciate our state-of-the-art results, and its many applications to training differentially private and Bayesian ML models, and we thank you for supporting our paper.  We answer your specific questions below.\n\n**“…generalize to other barrier functions”**\n\n\nIt is an interesting question whether one can extend our results to more general barrier functions. In Appendix $F$ we show that, if any $\\nu’$-self-concordant barrier function for the polytope $K$ is used in place of the log-barrier function in our algorithm, then our regularized barrier function is self-concordant with parameter $\\nu = \\nu’ + L^2 R^2$ (in, e.g., the setting where $f$ is $L$-Lipschitz). In the special case where $f$ is constant, [Laddha, Lee, Vempala, STOC 2020] show that their Dikin walk Markov chain has mixing time $O(\\bar{\\nu} d)$, but they require the barrier function to satisfy a stronger condition,  “strong self-concordance with symmetric self-concordance parameter” $\\bar{\\nu}$ (in particular, they show that the Lewis weights barrier satisfies this stronger condition with parameter $\\bar{\\nu}=d$). Showing that our regularized barrier function is *strongly* self-concordant with symmetric self-concordance parameter $\\bar{\\nu} = d + L^2 R^2$ when, e.g., the Lewis weights barrier is used, is thus a natural direction for future work.\n\n\n**\"The algorithm takes $O(m d^{\\omega-1})$ to form the Hessian matrix. This does not seem to be tight.”**\n\nInvestigating whether one can reduce the cost of computing the log-barrier Hessian matrix at each step of our algorithm is an interesting open problem. In particular, we note that [Laddha, Lee, Vempala, STOC 2020] show, in the special case where $f$ is constant, that the (average) cost of computing the Hessian matrix of the log-barrier of the polytope $K:= $ {$\\theta \\in \\mathbb{R}^d : A \\theta \\leq b$} at each step of their Dikin walk can be improved to roughly $O(d^2 + \\textrm{nnz}(A))$ arithmetic operations, where $\\textrm{nnz}(A)$ denotes the number of non-zero entries of $A$. Whether their result can be extended to the problem of computing the regularized barrier functions used in our algorithms, in the more general setting where $f$ is $L$-Lipschitz or $\\beta$-smooth, is an interesting direction for future work.  We will discuss this in the Conclusions, limitations, and future work section.\n\n**“Dependence on $L$ and $R$...”**\n\nIt may be possible to eliminate the polynomial dependence on $L$ and $R$, but it is outside the scope of the current paper.  This is a challenging problem, and we discuss it in the Conclusions, limitations, and future work section, and in more detail in Appendix $F$.\n\nOne challenge in obtaining bounds which are independent of $L$ and $R$ is that the isoperimetric inequality used (in our paper, and in many prior works on the Dikin walk) to bound the mixing time of Dikin walk Markov chains relies on a metric—the cross-ratio distance for the polytope $K$— which, roughly speaking, defines the distances between Markov chain steps by how quickly these steps approach the boundary of the polytope $K$. Measured in the cross-ratio distance, the steps proposed by our Dikin walk (or, more generally, by versions of the Dikin walk which take steps that are small enough such that the term $e^{f(z)-f(\\theta)}$ in the Metropolis acceptance probability is $\\Omega(1)$ when $f$ is $L$-Lipschitz on a polytope contained in a ball of radius $R$) may be as small as roughly $O(\\frac{1}{RL})$ with respect to the cross-ratio distance metric, and the mixing time bounds one would obtain via the aforementioned isoperimetric inequality would depend polynomialy on $RL$."}]

Author rebuttal:

Meta-review: This work studies the algorithmic problem of sampling from a log-concave distribution over a bounded polytope. The main contribution is a probably fast method using the Dikin random walk with a soft threshold regularizer. The reviewers agreed that this is a technically solid contribution that should be accepted to NeurIPS.
Decision: Accept (poster)