Faster Sampling from Log-Concave Densities over Polytopes via Efficient Linear Solvers
Authors: Oren Mangoubi
Nisheeth K. Vishnoi
Venue: ICLR 2024 poster (2024)
Keywords: Logconcave sampling
Dikin walk
Markov chain Monte Carlo
Interior point methods
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)

Abstract: We consider the problem of sampling from a logconcave distribution $\pi(\theta) \propto e^{-f(\theta)}$ constrained to a 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$. The fastest-known algorithm for the setting when $f$ is $O(1)$-Lipschitz or $O(1)$-smooth runs in roughly $O(md \times md^{\omega -1})$ arithmetic operations, where the $md^{\omega -1}$ term arises because each Markov chain step requires computing a matrix inversion and determinant ($\omega \approx 2.37$ is the matrix multiplication constant). We present a nearly-optimal implementation of this Markov chain with per-step complexity that is roughly the number of non-zero entries of $A$ while the number of Markov chain steps remains the same. The key technical ingredients are 1) to show that the matrices that arise in this Dikin walk change slowly, 2) to deploy efficient linear solvers which can leverage this slow change to speed up matrix inversion by using information computed in previous steps, and 3) to speed up the computation of the determinantal term in the Metropolis filter step via a randomized Taylor series-based estimator. This result directly improves the runtime for applications that involve sampling from Gibbs distributions constrained to polytopes that arise in Bayesian statistics and private optimization.

Published as a conference paper at ICLR 2024
FASTER SAMPLING FROM LOG-CONCAVE DENSITIES OVER POLYTOPES VIA EFFICIENT LINEAR SOLVERS
Oren Mangoubi Worcester Polytechnic Institute
Nisheeth K. Vishnoi Yale University
ABSTRACT
We consider the problem of sampling from a log-concave distribution x(f) oc e -((()) constrained to a polytope K == (0 E IRd AO < b),, where A E Rexxd and b Rm The fastest-known algorithm Mangoubi & Vishnoi (2022) for the set- ting when f is O(L)-Lipschitz or -mmoooh runs in roughly O( x mdd=1) arithmetic operations. where the mdd term arises because each Markov chain step requires computing matrix inversion and determinant (here w & 2.37 is the matrix multiplication constant). We present a nearly-optime implementation of this Markov chain with per-step complexity which is roughly the number of non-zero entries of while the number of Markov chain steps remains the same. The key technical ingredients are I) to show that the matrices that arise in this Dikin walk change slowly, 2) to deploy efficient linear solvers that can leverage this slow change to speed up matrix inversion by using information computed in previous steps, and 3) to speed up the computation of the determinantal term in the Metropolis filter step via randomized Taylor series-based estimator. This result directly improves the runtime for applications that involve sampling from Gibbs distributions constrained to polytopes that arise in Bayesian statisties and private optimization.
INTRODUCTION
We consider the problem of sampling from log-concave distribution supported Given a polytope K == {0 E Idd AO < b}. where A E Rmxd and b E IRm and convex function f K - TR, output a sample E K from the distribution x(f) xx e-f(0) This problem arises in many applications, including in Bayesian inference, differentially private optimization, and integration. For instance, in Bayesian Lasso logistic regression, f(0) Σn E(0 x;) where é is
polytope:
the logistic loss and I; are datapoints with |2:1|2 < 1, and K 10 LERE |0ll1 < O(1) is an balbll (see e.g.. Tian al. (2008): Zhang et al. 2017): ; Kim al. (2018)) The closely related optimization problem of minimizing T.-L.ischitz or 3-smooth convex function f constrained to polytope arises in numerous contexts as well. A special case the setting where f is linear on which corresponds to linear programming, since in this case 3 O. In applications where one must preserve privacy when minimizing f constrained to K, sampling from the exponential mechanism of McSheery & Talwar (2007), where (() x e f()) on K, allows one to obtain optimal utility bounds under e-differential privacy (see also Bassily et al. (2014); Munoz et al.) 2221 Kapralov & Talwar (013)). Depending on the application, the polytope could be e.g. the probability simplex, hypereube, or -bll,, where m (d) and nnz(A) the number of constraints defining the polytope is m 2 dl and nnz(A) = ((m) (see e.g. Hsu et al. (0144) Barvinok & Rudelson (2021) for some examples). Here nnz(M) denotes the number of non-zero entries of matrix M.
Od), as well as many applications where
This sampling problem. and its generalization convex K, has been studied multiple prior works.
One line of work Kannan & Narayanan (2012) Narayanan 2016). [ee & Vempala (2018). Cousins Vempala (0018). (Chnn et al. (2017). Jia et al. (2021) focuses on the special case ( 0). In particular, Narayanan (201))))n)0))))))) (2016) gives an algorithm, the Dikin walk Markov chain, which takes Omdx md=-1) x log( ( arithmetic operations to sample within total variation (TV) error d > from the uniform distribution on a polytope K C Rd given by m inequalities from a w-warm start. Here, 2.37 ..... is the matrix-multipliattin constant and distribution y is wo-warm for 2 w.r.t.

---

Published as a conference paper at ICLR 2024
the target distribution π if (x=)/<(=)) S w for every K. Other works consider the more general problem when π may not be the uniform distribution. Multiple works Applegate & Kannan [199]);
Frieze et al. (1994); Frieze & Kannan (1999); Lovász & Vempala (006 2007 Brosse et al. (2017)
have given algorithms that apply in the general setting when K is a convex hody with access given by membership oracle (or related oracles).
Narayanan & Rakhlin (2017) give an algorithm based on the Dikin-walk Markov chain to sample
on K where f is L-Lipschitz or 3-smooth. Building upon the work of Narayannn & Rakhlin (017)) Mangoubi & Vishnoi (2022) presented a Dikin walk which utilizes a thrt-threshold" regularizer, which takes O[ dL " R²) x Xlog(% ) x (md + Ty)) arithmetic operations to sample from an L-log-Lipschitz log-concave distribution π on a polytope K within error& > 0 in the TV distance from a w-warm start. Here Ty the number of arithmetic operations to compute the value of f- Specifically, each iteration in their soft-threshold Dikin walk algorithm requires mdd-1 arithmetic operations to compute the (inverse) Hessian of the log-barrier function for the polytope K and its determinant.
from any log-concave
Our contributions. We show that the per-iteration cost of computing the (inverse) Hessian of the log-barrier with soft-threshold regularizer, and its determinant, at each step of the soft-threshold Dikin walk can be reduced from md= to nnz (A) + de (see Theorem 2.1.. More specifically, our version of the Dikin walk algorithm (Algorithm 1 takes Ó((md + dI. R²) log( 8 )) x (nnz(A) + w
d² + Ty) arithmetic operations to sample within total variation error O()) from an L-log-Lipschitz log-concave distribution on polytope K defined by m inequalities, and ( + +d log( ()) (nnz(A) d TTy arithmetic operations to sample from a 3-log-smooth log-concave distribution on K. Compared to the implementation of the soft-threshold Dikin walk in Mangoubi & Vishnoi 2022), , we obtain an improvement in the runtime for the case of L-log-Lipschig log-concave distributions of at least d 2 in all cases where, c.g., Ty = ((d²). If one also has nnz(A) = O(m), the improvement is d=-1 if m 2 d2 and md= if m S d: see Table for comparison to prior works.
3-log-smooth
The main challenge improving the per-step complexity of the soft-threshold Dikin walk that the current algorithm uses dense matrix multiplication to compute the Hessian of the log-barrier function with soft-threshold regularization, as well as its determinant, which requires mdd arith- metic operations. The fact that the Hessian of the log-barrier function oftentimes changes slowly at each step of interior point-based methods was used by Karmarkar (1844), [aidya (1989). and Nesterov & Nemirovsky (1991), and later by Lee & Sidford (2015), to develop faster linear solvers for the Hessian to improve the per-iteration complexity of interior point methods for linear program- ming. To obtain similar improvements in the computation time of the soft-threshold Dikin walk, we need to show that the regularized barrier function proposed by Mangoubi & Vishnoi (2022) changes slowly as well. The notion of a barrier function changing slowly with respect to the Frobe- nius norm proposed in Laddha et al. (2020) suffices, but we have to show that this notion holds for the soft-threshold regularized log-harrier function. To improve the per-iteration complexity of the soft-threshold Dikin walk, we show how to use the efficient inverse maintenance algorithm of Lee) & Sidford (2015) to maintain a linear system solver for the Hessian of the log-barrier function with soft-threshold regularizer, and then show how to use this solver to efficiently compute random es- timate whose expectation is the Hessian determinant in nnz(A) + d² arithmetic operations. This is accomplished by first computing a randomized estimate for the log-determinant of the Hessian, and then converting this estimate into a randomized estimate for a "smoothed" Metropolis acceptance probability via a a detailed overview of the techniques in Section
series expansion for the acceptance probability. We present
2 MAIN RESULT
Theorem 2.1 (Main result) There exists an algorithm (Algorithm) 7 which, given the following
inputs, 1) 6, R > and either L > 0 or 3 > 0, 2) E Rmxdd be R.m that define a polytope
K = (0 IRd AO < b} such that K is contained in a ball of radius R and has nonempty interior,
3) an oracle for the value of a convex function f K > R. where f is either L-Lipschitz or 3- smooth, and 4) initial point sampled from distribution supported on K which is w-warm with
respect to π x e for some w > O, outputs a point from a distribution i where l|p - π TV S 6. This algorithm takes at most
2

---

Published as a conference paper at ICLR 2024
O(md+ R log( x ( arithmetic operations in the setting where f is L-Lipschitz,
Ö((md d3R²) log( k)) x (nnz(A) + de Ty) arithmetic operations in the setting where f is B-smooth,
where T the number of arithmetic operations to evaluate f.
Theorem 2.1 improves by a factor of (uax) A), arithmetic operations on the previous bound of Mangoubi & Vishnoi 2022) of O((md + dL²R²)= x (md + Ty)) arithmetic operations for e-f where f is L-Lipschitz from a w- w-warm start. Here Ô hides logarithmic factors in 6, L, R2 d. When is instead ß-smooth their bound is O( + d3 dBRR) x (md=-1 + Ty), , and the improvement is the same. When T, S O(d2) (as may be the case when evaluating fequires computing inner products with = d datapoints 12),, the improvement is (at least) d='2. If one also has that A is (mm)-spars,, (unz(A) = (m)), as is the case in applications where each constraint inequality involves O(1) variables), the improvement is d= in the regime [dd2) and md=-3 3 if m S (d?)).
sampling from distribution
Theorem 2.1 also improves on the bound of (Lovász & Vempala 2006); Theorem 1.1) of [m + Ty) arithmetic operations by min(({)² ,ي when e.g. Tf = Mid²) and d
m Odd). In the setting where f is f-smooth, the improvement is min
min(( ( ş2z).
In the Bayesian Lasso logistic regression example considered in Mangoubi & Vishnoi 2022). our algorithm takes O(d*) operations from w-warm start as f is 3-smooth and L-Lipschitz, with == = = d, = L, r 1/va, nnz(AA O(d). This improves by d=-2 on their bound of (() and by d on the bound of (dd) implied by Lovász Vempala (2006).
Theorem 2.1 also improves upon the running time for (s, O)-iifferentially private low-rank approx- imation by a factor of d--2 on the bound (dd log( operations of Mangoubi & Vishnoi 2022). Moreover, it improves by a factor of mmd=-1
their bound of O((md + dn e2) +
mne(A)+o"
(log(nRd((re)) (nddd//rr))))mddd- + T))) arithmetic operations for the problem of finding a minimizer Ô K of convex empirical risk function f(0, z) Σ (, 2,) under (e, 0) d)dfifeerential privacy, when the losses (z( are i--Lipschitz, for any L and dataset & E D See the details in their paper.
Arithmetic operations per iteration
Arithmetic operations if mox(A) O(1)). L,, O(1), R O(4). O(d)
Arithmetic operations if Old), O 3 L Old), 70
Algorithm
Iterations for L-Lipschitz f
lierations for f-smooth f
Prosimal Langevin MC
O(maa=13 Tey
0(d58-912(45)
0(810115-9)
(Brosse et (2017)
Hit-and run
O(yn))
((d"()1))
O(d")
O(P)
(Lovdsz Vompala (2006))
Dikin Walk of
O(ma"=11) Ty
O(d L²RR)) R²))
Off (AR))
Ö(ofthu)
Namyanan Rakhlin (2017)
Soft-Threshool Dikin Wilk
O(ma)-11) Ty
Ô + (L)
O(mn+
Ö(dFttr)
O(dPtwr
(Mangoubi Vishaoi (20221)
paper
O(une(A)++
O(nd+ d1
O(d9)
O(dt)
Table Number
arithmetic operations per-iteratioo) from logeoneave / on polytope K when L-Lipschtiz 3-smooth, from number of operations to compute f f. The Ô notation hides lognaishmic factors of d, 6, R. The fifhh column gives rumtines when A is polytope with on, maz(4) O(d") and R - O(d] that conrains the unit ball (and is thus not well-rouneed and / is Of O())-Lipschiz. The sixth colomn corresponds to sampling a Bayesian Lasso logistic regression posterior distribution Oid' ducapoints,
algorithons which imply
sampling within O())
start. Tj and Tey respectively, the
the 4.-ball.
3 ALGORITHM
For input E Int(K), K),, let S(0) = Diag(A0 - b). and ( == - where
5(0) = (()))
A
and Á ==
For input W E Rddd and t E [o, 1]. let
0
a=n
La
Ta-
I(W.t) := Im+d d + t(W Im+ For a rectangular matrix M, a vector u E Rd, a diagonal matrix D. and of algorithmic parameters P, Solve(v, D, M; P) outputs vector which is (ideally)
3

---

Published as a conference paper at ICLR 2024
solution the system of linear equations M DM Given an initial diagonal matrix Do Initialize Do. M) generates parameters P corresponding to Do. Given a new diagonal matrix D, and set of parameters P. Update( D, M;; P) updates the parameters P and outputs some new set of parameters P corresponding to the new diagonal matrix D. For any C TRd denote the interior of 5 by (S) (0 E S B(0,t) t) for some f 0}.
Algorithm 1: Fast implementation of soft-th Dikin walk
Input: m,, d N, A E R=xd RR,, which define the polytope :== (0 Rd A0 b} Input: Oracle returning the value a convex function Hyperparameters: a > 0: v > 0: E N; E N
+ R. Initial point 0 E Int(K)
Set 0 0o. A
log det (et A) P Initialize Qu Initialize A)
Ta.
for 1, 1.... do Sample point & NNO
Set A S(0)E
Set a a=3Sollveu,, S(0)), A; ( (Computes Gaussian with covariance a -'A' AA" ((0)
if E Int(K) then
Set W S(0) (()
for ......WW do
10
Sample (0, Ia)
11
Sample from the uniform distribution on [0, 1].
12
Set Q = Update(T(W.O) 4) Qo)
Set w = Å (W Av
13
Set Y, = Solve(w, r(w,0) (). Å; Q)+ +
14
(Computes random variable Y, with expectation log det ((2) log det (0))
16
if : < Yi <. 2lgg then
Set X 1 Σ EE
19
Y < 1 then
Set co.
to be the first coefficients of the Taylor expansion of
1:2==7
Set X EwD (ee [[)))))))
22
Y 2 2 log( then Set =1
24
Accept 0 +- with probability & min(max(XX)) min (:-11)) ol= 11= -1:cc -110- -1,,,,,
Set P = Upd))e(S(()). P)
26
27
else
28
Reject
end
30
end
31
Output: 0
In Theorem
we set = 1020 12R5 log 1.02 (w/ss). We set the step size hyperparameters
a
1/(1084) 4) log (1/1) and = 1/au if f is L-Lipschitz, and the number steps to be
109 (2ma + 17' R²2) log(w/)) x log 1.01 (109 (2ma + - R²2) x log(w/s)) When f is 8 smooth (but not necessarily Lipschitz), we instead set = To (md+ log 1.01 (w/s). and
= 1/(10"d) log (1/7) and 1/()) Thus, in either case, we have < 6/(10007), Finally, we set the parameter N, the number of terms in the Taylor expansions, to be N l0
4 OVERVIEW OF PROOF OF MAIN RESULT
We give overview of the proof of Theorem 2.1 Suppose given A Rxdd that define a polytope K == {0 E Rd : A0 < of with non-empty interior and contained in a ball of radius R > 0. and oracle for the value of convex function f K - R, where f is either
xd 5 E R"
4

---

one chooses Published as a conference paper at ICLR 2024
L-Lipschitz or 3-smooth. Our goal is to sample from O e- K within any TV enror å > 0 in a number of arithmetic operations which improves on the best previous runtime for this problem (when, e.g., L, B, R (()), which given for the soft-threshold Dikin walk Markov chain. We do this by showing how to design an algorithm to implement the soft-threshold Dikin walk Markov chain with faster per-iteration complexity.
Dikin walk in the special case f = 0. First, consider the setting where one wishes to sample from the uniform distribution (ʃ 0) on polytope K, where K contained in a ball of radius R and contains a ball of smaller radius > 0.
Kannan & Narayanan 2012) give Markov chain-based algorithm, the Dikin walk, which uses the log-barrier function from interior point methods to take large steps while still remaining inside the polytope K. The log-barrier (p()) -ド log(b, aj 0) for the polytope K defines Riemannian metric with associated norm ||@]| Vu (0)u, where ( V2)(())
is the Hessian of the log-barrier function. At any point A E [nt(K), the unit ball Σ (b, oo-aa00 0ys with respect to this norm, E(O)= {a IR || 0 (0))) S 1), referred to the Dikin ellipsoid, is contained in K. To sample from K, from any point 0 E Int(K), the Dikin walk proposes steps from a Gaussian with covariance (0), where a > 0 is a hyperparameter.
S Og)) this Gaussian concentrates inside F((0) C K, allowing the Markov chain to propose
steps which remain inside K w.h.p.
To ensure the stationary distribution of the Dikin walk is the uniform distribution on K, if = E [n)( the Markov chain accepts each proposed step z with probability determined by the Metropolis acceptance rule min (p(=-1) p(0==]) 1) = (del(H((1) =10-==1feco)-ie -1== -1==000 where p(0 (del(H(=))
=) x (det -(0
denotes the probability (density) that the Markov chain at 0 pro-
poses an update to the point
If one chooses a too large, each step of the Markov chain will be rejected with high probability.
Thus, ideally, one would like to choose as large as possible such that the proposed steps -10-=
cepted with probability (11). To bound the terms Hden(H(s))// (det(H(o)) and
Kannan & Narayanan (2012) use the fact that the Hessian of the log- log-barrier changes slowly w.r.. the local norm |[ (0)- More specifically, to bound the change in the determinantal term, they use the following property of the log-barrier, discovered by Vaidya & Atkinson (1993), which says its log-determinant V V(0)= == log (f) satisfies
(VV (VV)) S Odd
v0 E Int(K).
Since the proposed update (z 0) is Gaussian with covariance (0), if one chooses 4. by standard Gaussian concentration inequalities
S ((1) w.h.p. Thus, log det (H)()) log dct(H(()) (2) V(0)= = 22)1), and (des( (o)) 2211 w.h.p.
To bound the total variation distance between the target (uniform) distribution π and the distribution Vr of the Markov chain after T steps, Kannan & Narayanan (2012) use an isoperimetric inequality of Lovász & Vempala (2003) which is defined in terms of the Hilbert distance on the polytope K. For any distinct points u, v E Int(K) the Hilbert distance
- 10-__1-X||p-0000 1lp-912 ou, r) = p- - x
where p. q are endpoints of the chord through K passing through u, v in the order p, u, U, g.
To apply this isoperimetric inequality to the Dikin walk, which proposes steps of size roughly ((a) Ila- (00), they show that the Hilbert distance satisfies
w.r.t. the local norm
a² (w, v) DE
[a] (u-(2=))) [o, u-
- vlla-llla-000 H(6)
They then use the isoperimetric inequality w.r.t. the Hilbert distance together with standard conduc- tance arguments for bounding the mixing time of Markov chains (see e.g. Lovász & Simonovits (1933), to show that, the Markov chain is initialized with a w-warm start, after steps, the total variation distance to the stationary distribution æ satisfies
||4
F||TV S O(v ma -1)"
Wi N,
5

---

Published as a conference paper at ICLR 2024
provided that each step proposed by the Markov chain is accepted with probability 21(). For (1)), after T = O log( w )) = Omm ()) steps, 1|r - x S
8
a- 1 At each step in their algorithm they compute the proposed update ----==-------- (V--HO00 where & (0, Ia), and the terms /aool//(aa////ooeeeee and e -l==00 ( -100-=11l000 in the acceptance probability. Writing (0) = A 5(0) where S()) = Diag(AR A0 b). these computations can be done in (mdd= ) arithmetic operations using dense matrix multiplication. Their algorithm's runtime is then T x mdd = ((md x md x log( ) arithmetic operations.
Dikin walks for sampling from Lipschitz/smooth log-coneave distributions. One can extend the Dikin walk to sample from more general log-cones log-concave functions on a polytope K, by introducing a term ((()) to the Metropolis acceptance probability. This causes the Markov chain to have stationary distribution (0) x e e-((0), since for every 0, E Int(K),
min p( e (00 100))))) 1) p(0 + 2))(0) = min (100-3) ( 1) p(z - +)T(2)) pp=--0]
However, the term e (0) can become exponentially small with the distance - 0112 (and thus with the step-size o)) even if, e.g., I is T-LLipschitz. Thus, to obtain runtime polynomial in T,, d, one may need to re-scale the covariance (0) to ensure f0)) 2 S2(1).
One approach, taken by MMngoubi & Vishnoi (2022), is to propose steps with (inverse) covariance matrix equal to regularized log-barrier Hessian, (0) = a H0) + + for some > O. Setting 17 at ensures (by standard Gaussian concentration inequalities) that the proposed step satisfies ||= 2|12 S (1) w.h.p. and hence ef(0)- = 52(1) w.h.p.
To bound the total variation distance, in place of I they show lower bound on the Hilbert distance w.r.t. the local norm || [1)))))
o² 2
1 R² 2m Di=1 (a. b.)²
(a (u-e)²² 2 2ma +22 - g2 |/u 011 (l))(w))
1
This leads to a bound of O1ma-1 + 1 122) log(()) = 0 + L22 log( on the number of steps until the Markov chain is within 0(5) of π in the total variation distance. In addition to the operations required in the special case when f = 0, each step requires calling the oracle for f which takes Ty arithmetic operations. Thus, they obtain a runtime of O L R?) (md (ma=-1+Th) +T))
log( 8 arithmetic operations to sample from e
Faster implementation of matrix operations when / = O. When f = 0, one can use that the barrier Hessian H (0) changes slowly each step of the Dikin walk to compute the proposed step, and determinantal terms in the acceptance probability, more efficiently.
The fact that the Hessian of the log- barrier function oftentimes changes slowly at each step of interior poset-based methods was used by Karmarkar (1984). and later by Vaidya (1989), and Nesterov & Nemirovsky (1991) when developing faster línear solvers for (0) to for linear programming. Specifically, they noticed for many interior point methods, which take steps 0,, .2.... [n), the update (0.) = A A at each iteration where D: S(O.). is (nearly) low-rank, since in these interior point methods most of the entries of the diagonal matrix D D.,,D, Dj are very small at each step This allowed them to use the Woodbury matrix formula to compute a low-rank update for the inverse Hessian (0:), reducing the per-step complexity of maintaining a linear solver for H at each iteration i. Specifically, the Woodbury matrix formula says
interior point methods
(M + UCV)-1 M-1 M V -U----C-------U-VV-----
for any M E Rdd U E Bddk V E Rexd C E txk, where l may be much smaller than d.
Lee Sidford (015) obtain an improved algorithm for maintaining linear solvers for sequences of matrices {A D. under the assumption that D. changes slowly at each iteration i with respect weighted Frobenius norm, [|(D++1 D.)D, "IIF < 0(1). This assumption is weaker than requiring most of the entries to be nearly constant at each step. as it only requires a (weighted) sum-of-squares of the change in the entries of D; to be O(1). To obtain their method, they sample diagonal entries of D, with probability proportional to the (change in) the entries' leverage scores, quantity used to measure the importance of the rows of matrix (for matrix M the i'th leverage

---

Published as a conference paper at ICLR 2024
score is defined [M to obtain a spectral sparsifier for the updates A (D,+1 A, allowing them to make low-rank updates to the matrix A D,A at each iteration. Specifically, for any sequence of diagonal matrices D1, D2,....attte satisfying
Using these subsampled diagonal entries they show how
(3) at each i E N, they show that they can maintain a ((nnz(A) A) + +d d-mmt linear system solver for the matrices A D,4, at a cost of ((nz((A) + d²) arithmetic operations per-iteration. (plus an initial cost of ( A) + d²)) to initialize their algorithm).
[1 11) - 0(1
Laddha et al. 2020) apply the fast linear solver of Lee & Sidford (2013) to the Dikin walk in the special case where f O. Specifically, they show the sequence of log- barrierrirr Hessians A D,A where D. = (6,). also satisfies (3) when 0, are the steps of the Dikin walk. They then use the fast linear solver of Lee & Sidford 2015) to evaluate the matrix operations (A (AT)( (())))=A)= -) &gg and compute an estimate for dct(A' (0.) (-2))/2 A) at each iteration. This gives a per-iteration complexity of (n)( + d' for the Dikin walk in the special case = 0.
Our work. To obtain a faster implementation of the (soft-threshold) Dikin walk in the general setting where f is Lipschitz smooth on K, we would ideally like to apply the efficient inverse maintenance algorithm of [ee & Sidford] (2015) to obtain fast linear solver for the log-barrier Hessian with thrthhreehold regularizer (0) a'AA S()) AA + n/a- Specifically, writing
A -----A' S(0) À where À =- (A] and S(0) = ((0)
a21 1a we would like to obtain a fast linear solver for the sequence of matrices ((01). ) ((02).... ...., where 0g, 02,02.-. 02.... are the steps of the soft-threshold Dikin walk, and use it improve the per-iteration complexity of the walk. This poses two challenges:
(()) = a
Ia.
0
1. To obtain an ((nzz(A))) (A) + -time linear solver with the efficient inverse maintenance algorithm, we need show that (00.) a- 'À ((,))² does not change too quickly at each step of the soft-threshold Dikin walk. Specifically, we need to show, at each i w.h.p..
(1($2(0++1) (0++1))) (0+)))))-((+++)) < (1)..
(4)
Once we have a linear solver for (()), this immediately gives us a way to compute the proposed updates of the Markov chain by solving a system of linear equations (see e.g. Section 2.1.1 of Kannan & Narayanan (2012)).
We also need to evaluate the terms det (6)) in the acceptance probability. However, access to a linear solver for (0) does not directly give way to compute its determinant.
Bounding the change in the soft-threshold log-barrier Hessian. From concentration inequalities for multivariate Gaussian distributions, one can show that, with high probability, the Dikin walk
proposes updates = <- 0 which have length ((vd) in the local norm: || 0lla(0) = O(va) log-barrier Hessian matrix 2(0) at
w.h.p. To bound the change in the soft-thar each step of the walk, we would ideally like to show that the Frobenius norm of the derivative of this matrix changes slowly with respect to the local norm || [| The woff-hhreshold-rggulaate
log-barries function & satisfies the self-concordance property [M) 3 T (())v<(()) (0) for any 0 2 E Int(K). and e E Idd, where DM (+) denotes the derivative of valuee function M in the direction h. However, the self-concordance property does not directly imply that (4) holds, as (4) contains a bound in the Frobenius norm. To overcome this issue, we show that (a rescaling of) the regularized log-barrier Hessian, (0) == a ((0) = <2) + a la satisfies the following strengthening of the self-concordance property w.r.t. the Frobenius norm:
(z) - ((0))-1--1()) S
||0 |10==llace:
v0, 2 E Int(()).
(5)
Laddha et al. (2020) show that (5) holds in the special case when is replaced with the Hessian (0) of the barrieri function without regularizer. However, the soft-threshodd log-barrier Hessian, 4. is not the Hessian around this problem, we instead use the fact (from Mangoubi & Vishnoi (2022) that & is the limit of a sequence of matrices (H, ) where each H, is the Hessian of lar-barrier function for an increasing set of (redundant) inequalities defining the polytope K. Combining these two facts
log-barries function for any set of equations defining the polytope K. To get
7

---

Published as a conference paper at ICLR 2024
5) holds for the Hessian & of the log-barries function with soft soft-threshold regularizer (Lemma B.I). Next, we apply (5) to show that (Lemma B.2
allows us to show Equation
$ 1100
(z) ((0))) - (()))))((O)-
I <((a+vv)),
(1-a 11000))))))
w.h.p. As the step size a (V/d), we get that
S 50(). w.h.p
Computing a randomized estimate for the determinantal term. To compute an estimate for the determinantal term in the Metropolis acceptance rule, we use a well-known method for estimating matrix-valued functions using polynomials (see e.g. Barvinok (0117) Rudelson & Zeitouni (2016). A related approach to the one presented here for computing the determinant in the special case when f = was given in Laddha et al. (2020), however contains a number of gaps in the proof and algorithm. We discuss these gaps, and the differences between their approach and the one given here, in Section C
Towards this, one can apply the linear solver for ((0) compute random estimate Y with mean log det( (()) log (z)) Y = (Â (t$(()) (1 ( (m+a) Au + log det ÀT A, Where ~ NO0, Ia) and ~ uniform(0, 1). While this provides an estimate for the log-determinant of 4 it does not give estimate with mean equal to Metropolis acceptance rule for π.
A
We
first replace the Metropolis acceptance
rule min 1 ( 1) =
min (Veura el [ll--llloo -100 - e00))))))11) with a (different) Metropolis acceptance (Veecoron)
rule, 2)), = (1 ( (e)
(-())) : π,y
ell=-01laco -100 -100,188,,, 1). whose
x min
determinantal term is smooth function of det (0)) and det()(=)):
(1 det(((())) ((e)
)ー
sigmoid (1 (ogg)e((()))) - det((((()))))
(6)
where sigmoid(t) := Ite-r. Since a(0, =)p(0 - (0) = (== 0) (2), this acceptance rule preserves the target stationary distribution (0) x e [(0) of the Markov chain.
1
Next, would idcally like to plug i.d. samples Yi.... Yn of mean log det (0) log (2) into a Taylor series for (igmoid(t) to obtain an estimate with mean equal to the l.h.s. of (6) Unfor- tunately, the Taylor series of sigmoid(t) at 0 has finite region of convergence T).
While one can use convexity of the function log det (#))), together with a bound on the gradient of log det (O)(), to show that 2 log (dc((((()))
1 log det( (())) > -I w.h.p., one may have that log (et()())) - 1/ log (0)) > π with probability (2)).
To get around this problem, we use two different series expansions for sigmoid(t) with overlapping regions of convergence:
1. a Taylor series expansion centered at 0 with region of convergence (-),,)).
2. a series expansion "at +oc" with region of convergence (0, +oo), which a polynomial in e sigmoid(t) = Σ oz k= cpe-kt kt
To determine which series expansion to use at each step of the Dikin walk, we show that each random estimate Y concentrates in the ball of radius 3 about its mean E{y] log (et((())))- w.h.p. (Lemma BB3).. As the regions of convergence of the two series expansions have intersection (0, this allows us use the value of the random estimate Y select choice of series expansion that, with high probability, contains E[[] = log det(((z)) - log (e))).
(e(()))))
Bounding the number of arithmetic operations. We have already shown that computing the proposed update 2 (Line 6 of Algorithm I and computing each random estimate Y, for the log- determinant (Line 14) can be done in O ( A) + d)) arithmetic operations at each step of the Markov chain by using the efficient inverse maintenance linear solver.
As both series expansions converge exponentially fast in the number of terms, only log(1/a) terms in the series expansion are needed compute estimate X (Lines 18 and 21 with
8

---

Published as a conference paper at ICLR 2024
mean within error ( of the determinantal term in the Metropolis acceptance probability (e))+ (e))
I O(6). Thus, the number of calls to the linear solver at each Markov chain step is ((log(/8)), taking (((nz((A) + df) log(!/) arithmetic operations.
[[X]
Computing the terms f()) and fz) in the Metropolis acceptance rule requires two calls to the oracle for f. which takes O(Ty) arithmetic operations, at each step of the Markov chain. The other steps in the Algorithm take no more time to compute. In particular, while several "initialization" steps for the linear solver and other computations (Lines 4-6) require O(md==1) arithmetic operations, these computations are not repeated at each step the Markov chain, and thus do not change the overall runtime.
Since the Markov chain run for roughly T O((md + 1,2 122R²) log(m/6) steps, the runtime is roughly O(T ((T(((z((A) + df) log(?/6 + Ty)) )- = O( + L² R²)) (nnz((A)))))))) + dl +Ty log (w/6)))
Bounding the total variation distance. Mangoubi & Vishnoi 2022) show that the soft-threshold Dikin walk outputs a point with distribution frr within TV distance O()) of the target distribution ' after roughly T O(mmd + L²R²) log("/6) steps; the same TV bound holds if one replaces their acceptance rule with the smooth acceptance rule we use.
Our bound the TV distance follows directly from their bound, with some minor adjustments. This is because, while we give a more efficient algorithm to implement the soft-threshold Dikin walk Markov chain, the Markov chain itself is (approximately) the same the Markov chain considered in their work. The main difference (aside from our use of a smoothed Metropolis acceptance rule) is that the (mean of) the randomized estimate we compute for the acceptance probability accurate to within some small error O6)), since we only compute a finite number of terms in the Taylor series.
Thus, can define coupling between the Markov chain 01, 1,,02..... ; Or generated by our algorithm and the Markov chain }/, 1/ YT defined in their paper (with smoothed acceptance rule), such that 0g = Vi every i {1.... T} with probability 2 T'8. This implies that the distribution of the last step of our Markov chain is within total variation distance O(T5) of the distribution of Vr. Thus, || T Tv| S || Pr Prllrv + -|r- ællrv I OTT)) + O6). Pluggining 8/r in place of 6, we get || T| S 6.
CONCLUSION
Our main result improves on the runtime bounds of a line of previous works for the problem of sampling from a log-Lipschitz or log- fogoth log-concave distribution a polytope (Table I). Key to our result is showing that the regularized log-barrier Hessian changes slowly at each step of the soft-threshold Dikin walk, which allows us deploy fast linear solvers from the interior point method literature to reduce the per-iteration complexity of the walk.
The of fast linear solvers allows us achieve per-iteration complexity nnz(A) d2 for the soft-threshold Dikin walk that is nearly-linear in the input complexity of the problem when A is dense matrix; when A is sparse, the dependence of the per-iteration complexity on nnz(A) is nearly-linear as well. On the other hand, it would be interesting to see if the d² term can removed by deploying a different choice of fast linear solver for the Hessian.
Finally, while the fastest-known algorithms for minimizing linear functions on polytope K have runtime (nearly) equal to the matrix multiplication time Cohen et al. (02T); Diang al. (202T), the best runtime bounds for sampling from distribution π x e on a polytope (even when uniform) are larger polynomial the dimension Thus, would be interesting to see the overall runtime (i.e., number of steps times the per-step complexity) of sampling methods can be improved to match the runtime of linear programming methods.
ACKNOWLEEDGMENTS
NV supported in part by NSF CCF-21 award. OM was supported in part by an NSF CCF-2104528 award and a Google Research Scholar award. The authors thank Yin Tat Lee and Sushant Sachdeva for valuable feedback and discussions.

---

Published as a conference paper at ICLR 2024
REFERENCES
David Applegate and Ravi Kannan. Sampling and integration of near log-concave functions. In Proceedings of the twenty-thrrd annual ACM symposium Theory of computing, pp- 156-163, 1991.
Alexander Barvinok. Approximating permanents and hafnians. Discrete Analysis, 2017.
Alexander Barvinok and Mark Rudelson. A quick estimate for the volume of a polyhedron. arXiv preprint arXiv:21 2021.
Raef Bassily, Adam Smith, and Abhradeep Thakurta. Private empirical risk minimization: Efficient algorithms and tight error bounds. In 2014 IEEE 55th Annual Symposium on Foundations of Computer Science, pp. 464 473. IEEE, 2014.
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 Learn- ing Theory, pp. 319-342. PMLR, 2017.
Yuansi Chen, Raaz Dwivedi, Martin J Wainwright, and Bin Yu. Vaidya walk: A sampling algorithm based on the volumetric barrier. In 2017 55th Anmual Allerton Conference Communication, Control, and Computing (Allerton), pp. 1220-1227. 2017.
Michael B Cohen, Yin Tat Lee, and Zhao Song. Solving linear programs in the current matrix multiplication time. Journal of the ACM (JACM), 68(1):1-39, 2021.
Ben Cousins and Santosh Vempala. Gaussian cooling and (n²)) algorithms for volume and gaus- sian volume. SIAM Journal on Computing, 47(3): 1237-1273, 2018.
Alan Frieze and Ravi Kannan. Log-Sobolev inequalities and sampling from log-concave distribu- tions. The Annals of Applied Probability, 91):14-26, 1999.
Alan Frieze, Ravi Kannan, and Nick Polson. Sampling from log-concave distributions, The Annals of Applied Probability, pp. 812-832 12-8877,1 1994.
Justin Hsu, Aaron Roth, Tim Roughgarden, and Jonathan Ullman. Privately solving linear pro- grams. In Automata, Languages, and Programming: 41st International Colloquium, ICALP 2014, Copenhagen, Denmark, July 8-11, 2014, Proceedings, Part I 41, pp. 612-624. Springer, 2014.
He Jia, Aditi Laddha, Yin Tat Lee, and Santosh Vempala. Reducing isotropy and volume to KLS: an o" x2) volume algorithm. In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pp. 961-974, 2021.
Shunhua Jiang, Zhao Song, Omri Weinstein, and Hengjie Zhang. A faster algorithm for solving gen- eral LPs. In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pp. 823-832 2021.
Ravi Kannan and Hariharan Narayaman. Random walks on polytopes and affine interior point method for linear programming. Mathematics of Operations Research, 37(1):1-20, 2012.
Michael Kapralov and Kunal Talwar. On differentially private low rank approximation. In Proceed- ings of the twenty-foure annual ACM-SIAM symposium Discrete algorithmss, pp. 1395-1414, SIAM, 2013.
Narendra Karmarkar. A new polynomial- -imm algorithm for linear programming. In Proceedings of the sixteenth annual ACM symposium Theory of computing, pp. 302-311, 1984.
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 lexi- con for ultrasonography. Uttrasonggaahyy 3761):3, 2018.
Aditi Laddha. Yin Tat Lee, and Santosh Vempala. Strong self-concordance and sampling. In Pro- ceedings of the 52nd annual ACM SIGACT symposium on theory of computing, pp. 1212-1222, 2020.
10

---

Published as a conference paper at ICLR 2024
Yin Tat Lee and Aaron Sidford. Efficient inverse maintenance and faster algorithms for linear pro- gramming. In 2015 IEEE 56th Annual Symposium on Foundations of Computer Science, pp- 230-249, IEEE, 2015.
Yin Tat Lee and Santosh S Vempala. Convergence rate of Riemannian Hamiltonian Monte Carlo and faster polytope volume computation. In Proceedings of the 50th Anmual ACM SIGACT Symposium on Theory of Computing. pp. 1115-1121, 2018.
László Lovász and Mikiós Simonovits. Random walks in a convex body and an improved volume algorithm. Random structures & algorithms,4Ch) 4(4):359-412, 1993.
László Lovász and Santosh Vempala. Hit-and-run is fast and fun. preprint, Microsoft Research, 2003.
László Lovász and Santosh Vempala. Fast algorithms for logconcave functions: Sampling, rounding, integration and optimization. In Foundations of Computer Science, 2006. FOCS'06. 47th Annual IEEE Symposiun on, pp. 57-68. IEEE, 2006.
László Lovász and Santosh Vempala. The geometry of logconcave functions and sampling algo- rithms. Random Structures Algorithms, 30(3):307-358, 2007.
Oren Mangoubi and Nisheeth K Vishnoi. Sampling from log- concave distributions over polytopes via soft-threshol. Dikin walk. arXiv preprint arXhv:2206.09384, 2022.
Frank McSherry and Kunal Talwar. Mechanism design via differential privacy. In 48th Annual IEEE Symposium on Foundations of Computer Science (FOCS'07) 077), pp. 94-103. IEEE, 2007.
Andres Munoz, Umar Syed, Sergei Vassilvtiskii, and Ellen Vitercik. Private optimization without constraint violations. In International Conference on Artificial Intelligence and Statistics, pp. 2557-2565. PMLR, 2021.
Hariharan Narayanan. Randomized interior point methods for sampling and optimization. The Annals of Applied Probability, 26(1):977644 2016.
Hariharan Narayanan and Alexander Rakhlin. Efficient sampling from time-varying log-concave distributions. The Journal of Machine Learning Research, 18(1):4017-4045,2 2017.
Yurii Nesterov and Arkadi Nemirovsky. Acceleration and parallelization of the path-following in- terior point method for a linearly constrained convex quadratic problem. SIAM Journal on Opti- mization, 1(4):548-564 1991.
Mark Rudelson and Ofer Zeitouni. Singular values of gaussian matrices and permanent estimators. Random Structures & Algorithms, 48(1):183 22222,2244 2016.
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. Computational statistics & data analysis, 52(7):3528 -542,, 2008
Pravin M Vaidya. Speeding-up linear programming using fast matrix multiplication. In 30th annual symposium on foundationn of computer science, pp. 332-337. IEEE Computer Society, 1989.
Pravin M Vaidya and David S Atkinson. A technique for bounding the number of iterations in path following algorithms. In Complexity in Numerical Optimization, pp. 462-489. World Scientific, 1993.
J Von Neumann. Some matrix- inttiixiiequalittess and metrization of matric-space. tomsk univ. rev. 1, 286-300 (1937). reprinted in collected works, 1962.
Yuan Zhang, Jonathan N Hofmann, Mark P Purdue, Shili Lin, and Swati Biswas. Logistic bayesian lasso for genetic association analysis of data from complex sampling designs. Journal of Human Genetics, 62(9):8 19-829, 2017.
11

---

Published as a conference paper at ICLR 2024
CONTENTS
Introduction
Main result
3
Algorithm
Overview of proof of main result
Conclusion
Preliminaries
13
Proof of the main result
15
Additional comparisons of approaches to determinant computation
22
12

---

Published as a conference paper at ICLR 2024
A PRELIMINARIES
The following lemma Lee & Sidford (2015) allows us to maintain an efficient linear systems solver for the Hessian of the regularized barrier function t
Lemma A.1 (Efficient inverse maintenance, Theorem 13 Lee & Sidford (2015), Suppose that a sequence of matrices A C(R) A for the inverse maintenance problem, where the se- quence of diagonal matrices C()) C()). ... with ((k) diag(e)
(K) ......())) (k)), satisfies
E (14t))
= O(1) for all E N. Then there is an algorithm that maintains {AT C(k) 1 II
(An)((A) d2) arithmetic operation linear system solver for the sequence of matrices
T rounds in total of ((T(nn)(() d-) d=)) arithmetic operations over all rounds T rounds.
For any symmetric positive definite matrix M E Rddd define ||h]|ax := h hT Mh.
The following lemma is useful in bounding the change in the log-barrier Hessian:
Lemma A.2 (Lemmas 1.2 & 1.5 in Laddha et al. (2200) Let K = {o E R A0 < b}. where A Ruxd and b R such that K bounded with non-empty interior Let H
a
E (o, aa 0yz be the Hessian of the log-barrier function for K. Then for any æ, U E Int(K with [ || - y vll < 1, we have | | (2) H(x)) (x) & (1-1--v)lm)))
The following lemma gives high-probabilltyy lower bound the determinantal term the Metropolis acceptance ratio.
det (=)) in
det()(0))
Lemma A.3 (Lemma 6.9 in Mangoubi Vishnoi (002) (o, Ia). Let = + (l) 2g. Then
Consider any E Int(K), and ~
P det (9(z))
48
98
del (+))
50
100'
and
2
||2 01 - -1 - 0112(00 Σ
2
(8)
50
100"
The following lemma of Mangoubi & Vishnoi 2022) says that the Hessian & of the regularized bar- rier function, while not the Hessian of a log-barrier function for any system of inequalities defining K, is nevertheless the limit of an infinite sequence of matrices H,, where each H, the Hessian of logarithmic barrier function for system of inequalities defining the same polytope K.
Lemma A.4 (Lemma 6.7 in Mangoubi Vishnoi (2022) There exits a sequence of matrix- valued functions {H, 1 where each H, is the Hessian of log- bog-barrier function K, such that for all w In),
lim Hy(w) = (w)) 1
and
lim (Hy(w)) -x ( +
uniformly in w.
The following Lemma (Lemma A.5) gives a mixing time bound for basic implementation of the soft-threshold Dikin walk (Algorithm) below):
13

---

Published as a conference paper at ICLR 2024
Algorithm 2: Basic implementation of soft-threshold Dikin walk
Input: m, d E N, E Rmxd b E Rm, which define the polytope K = {0 E RRd A0 s b)
Input: Vo E (n), 7 N
Hyperparameters: > 0; O:
i . do
Sample ( (0. Ia)
Set H())) Dj= (b - D.
Set (0)) = (HH) + la
Set 2; - -3;
2; E [n)(K) then
Set H(2,) = EØ= (bg 2.
4ja
Set (22,) 1
(O))
+ n La
X min e 112,) (--155)
e 12 -112: 208.-2,113 22,,, 1)
Set ps OnE
(11))) (et(
11
Set Di+1 Z, with probability ps
else
13
Reject 2
end
15
Output: )/7
16 end
Lemma A.5 (Lemma 6.15 of Mangoubi & Vishnoi (2022)) Let Vo. .Y be the Markov chain generated by Algorithm 2 Let ò > O. Suppose that f K - R is either L-Lipschitz (or has B-Lipschitz gradient). Suppose that Vo ~ Vo where vo is w-warm distribution with respect to e with support on K. For any >> 0 denote by denote the distribution of Vt, for S I0 and n S do (or 1) d0 Then for any T 2 109 (2ma 1 R²) x log( 12) have that
|P - क F|l S 6.
Note: Lemma A.5 is given Mangoubi & Vishnoi
(022) for the same Markov chain as in Algo-
rithm 2 but with acceptance rule x min
-((4,) (o
-D: -
F(0s) (det(
in place of the acceptance rule pt used in Algorithm 2 The same proof Lemma A.5 given in
Mangoubi & Vishnoi (2022) holds for the above acceptance rule in Algorithm with minor modifi- cations.
The following trace inequality is needed in the proofs.
Lemma A.6 (von Neumann trace Inequality Von Neumann (1962)) Let M, Z E Rddd be matri-
ces, and denote by σ1 Z ... 2 da the singular values of M and by m 2 ... 2 aa the singular values of Z. Then
Lr(MZ) S Σ om
We need the following fact.
Proposition A.7 Let D E IR ((d+m))(ddm)) be diagonal matrix. Then
det DAA log (AT tr (TÅ
(TĀTDAA DA (1-) DÅ AT
-ATAAT ATA) dr.
Proof:
14

---

Published as a conference paper at ICLR 2024
Let ((T): T - for any O. Then
log (T) tr (J log(((7))) tr (T)
d
के
dr
=tr ttr((TA (1
AT A)).
Thus,
log det(A det(AADA) - log det(A = S/' dr
log (r
(ATAA DÅ ATA)) dr. ff tr (IT DA r) (
Finally, we need the following concentration inequality to show a high probability bound on the error Y - -E] of the estimator Y for log (et log (0).
Lemma A.8 (Hanson-Wrigee Wright concentration inequality) Let M Rdxd Letv N(0, Ia). Then t2
for all t > 0, P(Iu Mul > t) s exp min
B PROOF OF THE MAIN RESULT
To prove Theorem 2.1 we first show that the Frobenius norm of the Hessian of the log-barrier function with soft-threshold regularizer changes slowly with respect to the local norm (Lemmas B.1 and [.2).. We also prove concentration inequality for the random estimate for the log- determinant (Lemma B.3).
The first part of the proof of Theorem 2.1 bounds the running time, by using Lemma B.2 together with Lemma A. in the preliminaries to show that one can deploy the fast linear solver to obtain a bound nnz A) + de on the per-iteration complexity of Algorithm[
We then give proof of Theorem 2.1
The second part of the proof bounds the total variation distance between the output of Algorithm and the target distribution и. Towards this, that the series expansions Lines 18 and cogori w.h.p. and give a good estimate for the Metropolis acceptance probability. To conclude, this fact show that there exists coupling between the Markov chain computed by Algorithm [ and an "exact" implementation of the soft-barrier Dikin walk Markov such that these coupled Markov chains equal w.h.p. This allows us to use the bound on the total variation for the "exact" Markov chain (Lemma A.5 in the preliminaries) obtain a bound on the total variation error of the output Algorithm
the concentration bound in Lemma B.3 to show
Lemma B.1 Let (00) == -m log(b; - aa00) be a log-barrier function for K {0 E Rd Ae < b). where A E pmxd and b E RT. Let 2 O, and define the matrix-valued function (0) = (00) Then for any 0, Int(K) with |10 (10-21lv(o) 1, we have || (z)
(()) (0)ll < ||0 v10)
200)
Proof:
Let 0, E Int(K) be any two points in Int(K) such that ||0 =1] <
From Lemma A.4 we have that there exits a sequence of matrix functions {H,),2,, where each matrix is the Hessian of log-barrier function on K, such that for all w E Int(K),
lim Hw)w)= +(w),
(9)
10000
uniformly in wr and
lim (1) (+(w)) +
(10)
15

---

Published as a conference paper at ICLR 2024
uniformly in w.
Thus we have
lim |10 2||æ,,00) 10000
lim V (0 2) (H) 2) j-x0
(0 2)T
lim H, - 2) 10000
Eq. P
V(O 2) (0) - 2)
110 =|(
(11)
In particular, since |10 2||v(0) < 11 implies that there exists J E N such that
||0 21 1,, <1
Vj 2 J.
(12)
Moreover, since (w) is the Hessian of log-barrier function on K for every j E N, by Lemma A.2 and 12 we have that
||0 2 =|]11,00) H (2) HH,0)H,H (0)lle - (1 100 10)
Vj
(13)
Moreover,
lim |4, (0)(H(()) - (0)0 +0
lim H, (0)) lim H1}())
lim (²) lim Hj 2 (0))l
-=x 8
1+000
TV - (0) (()) 2
(14)
Therefore, plugging
[4] and
into
13 we have that
|| yllvve)
| ( v(0) - (0)1e <
(1 10 = =1(v(0)) )²
Lemma B.2 Let 0 +()) te where & N(0, Ia). Then we have with probability at least y that
[}{=)) Im Imm[F 1000 (og
Proof:
Setting Vi(0) Å ÄT (T((0), Å, have that () (x0(0) H((0) on) -Ta-
16

---

Published as a conference paper at ICLR 2024
Hence, we have by Lemma [.1 that, with probability at least 1 - 7
||((0) (((=)² In = -1000) ImIIF
- | A))) - Ialle
|| -
= | -0== (x)ll
= ||a
-2 - (x)]lr
| (z)
(x)]lr
1 - yll
S (1 | - ))
a =10 =| =1(
(1 $110 2|1o(0))2
S a 28v ddog 1
1000
where the second inequality holds by Lemma B.T
The fourth inequality holds because a
a log-1 (+) The third inequality holds with probability at least - y because
10 =11²(0) (0- = ( 70000)0))000000000000)0000)))))) < &Vdllog(
with probability at least - by the Hanson- Hanson-Wright Wright inequality (Lemma A.8
Lemma B.3 Let W E Bm++dx+ddd be a diagonal matrix, and let t > 0. Define O(1) = AT (m++a I(W Tm+aa) À. Let N(0, Ia), and les
= c e()) - A (W - ImmAv + log del AT
(15)
Suppose that 7m4m 1 W = 21m and that |[ -mm] p S cor some 0. Then
P(YY - E|| 2 c) Se =o
0.
(16)
Proof:
|[0(1)=ATT (W - Am) S 11( (IW + (1 (1 (1-1)))44) |TA((W
(W - All
$ [I(A ( (m) lee T -mm)
-2](AAT (A) (A Im)Alle
=21((ATA)=11) TrA ((Ā (W - Amm)
(I(AAA))~112 (AT (W- -Im)AA' MAATWW-WA (W-Im))
(AĀT (W -Im)AT (W Im))
<2] - V /IAÁTII 112 Ir ((W (WW-m)) -M)AA (W - Im))
$ 2)](ATAA - V VIÁÁTI 1l2 IT (A'A (W 1m))
$ <2)((ATA) A) ||ÁĀ 2/2 I ((W - -mm)
2 - -mlle
2c,
(17)
17

---

Published as a conference paper at ICLR 2024
where the fourth and fifth inequalities hold by the Von Neumann trace inequality (Lemma A.6
Thus, by the Hanson- Wright inequality (Lemmal A.8) we have that for every 2 0,
P(IoTO( (o(tt)-AA A (W Am Elu" ( A (W Im) Avi 2 8)
8²
S 2exp
min
8
[10(t)-1 -AT AA(W - Im AlB
|0(t)- (W - -m))lla>
2exp
(18)
8 |0(t)- A (W - Am)AlF
Plugging
17 into 18). we get that,
P([e T (()) -tt (W lm).Au - ( (W- - Am)Aul] 2 sc) S 2
Vs 0.
Proof: [of Theorem
2.1
Bounding the runtime:
Cost of maintaining and using linear solver: solver:: By Lemma outer "for" loop in Algorithm D
B.2 we have that at every iteration of the
(zz) -
2 = -11$(0)
Im 1m S O(1),
(19)
((0)2
with probability at least
Rewriting the l.h.s. of 19 we also have that with probability at least - 7,
Σ (1S(()))1)5(0) (1 + ()) - (1-
2
($(=) (10)))--1111 -1)² 0(1), (20)
for any t E [0,
To implement Algorithm ((0) A (())-2 2 and A (S(()) (())) t) A (tŠ(()) SS)) (1 (1-t) mm) for
use efficient inverse maintenance maintain linear solvers for
T rounds, where 7 the number Dikin walk steps. Plugging Inequality (99) (for the (()) lin-
ear system solver) and Inequality (0) (for the A (S(())
t) linear system solver) into
Lemma A.I we get that maintaining these linear system solvers for T steps (Lines 2
26 and 12
in Algorithm T.. and applying each linear solver at most A times at each step (Lines 6 and 14] takes
O (nuz (A) + +d) + d²) arithmetic operations.
Since T 2 d=-2, we have that the number of arithmetic operations is O(TNN + d))).
Cost of other steps: The remaining steps require less that ((TN(n))(A) + d2) steps to compute:
Each time they are run, Lines 5 and 13 can be accomplished in ((nzz(A) arithmetic operations by standard matrix-vecte multiplication. Sampling d-dimensional Gaussian vector (Lines and 10) can be done ((d) S O(nz((A) arithmetic operations and sampling from the uniform distribution 11 takes O(1) arithmetic operations. Lines 5 are called most TN times, and thus require (TNNnzz(A) arithmetic operations.
Every time is run, Line 25 requires computing f(0 and f(z), which can be done calls to the oracle for the value of f. It also requires computing || 12-0013 01la and ||0 110-=112 As
||2 0 (2 0)
(0) (T - 0) = (2 - 0)T () -1 ((0) -). 0),
this can be done in (mnz(()) arithmetic operations using standard matrix-vector multi- plication. For the same reason, computing | 2 =1(2()) also takes O(nnz arithmetic operations. Line 25 is called at most times, and thus contributes ( + Ty)) arithmetic operations the runtime.
18

---

Published as a conference paper at ICLR 2024
Each time they run, Lines 18 and
21 require computing Taylor series with at most
( terms, with each term requiring (NN) multiplications. Lines 18 and 21 are called at most T times and thus require O(TNa) arithmetic operations.
Line 2] which computes log (At(AA A), can be accomplished in ((md"=1) arithmetic operations using dense matrix multiplication. Since this step is only run once and O I ((TN(nnzz44 + d))) arithmetic operations, it does not change the overall runtime by more than constant factor.
Runtime for all steps: Thus, the runtime of Algorithm is mnz(A)))) dI.² R x log?.0) 2.01 mdd+12 +logg(((w) log( ())) (T nnz (A) + dB) arithmetic operations in the setting where is L--Lipschitz
O((md+
In the setting where f is B-smooth, the number of arithmetic operations is OTTN(T) + mnz(A) + )² = O((m + ddL2RR) x log2.0) 2.01 +lg(w)))()) log( ())) x (2'y + nnz + df).
Bounding the total variation distance: Plugging Lemma B.2 into Lemma B.3 we get that for every j E (N.....N))
P(V1 -D,l - 100
s 7
(21)
By Proposition A.7 we have that,
E[Y;] Elv" (LÂT S(0)
AT ATA)
+ log det ()
/ tr A -(0))) + (1- 7) (1 (AT -AT))-ATA)] ATA) dr
+ log det (AT(AA
Prop. A.7
log d ((0))))())A)
log det(((z)) log det (0)),,
(22)
where the second equality holds since Algorithm D samples the random vector (O, La) from the standard Gaussian distribution and the random variable ~ unif([0, 1]) from the uniform distri- bution on [0, 1].
Moreover, by Lemma A.3
we have that
B
E
2
log det (+))) log (et( 2
9 10
10
Lemia 4.3
100 7
Next, we show that
2 2Ny
(23)
and that
det (222)
det (00)
E
ST.
(24)
(et(
de)( dct(
dct(((=))
We first show probability at least
23 and
24
when 4
1 Yi $ log(r)- In this have by
21 that with
7
2 210g )
log det (x)()) - log (0) E|Y]
(25)
100"
To show
23
we first note that, by
21 we have that,
< Y, < log -
Vj vIsleN)
21-N7
(26)
19

---

Published as a conference paper at ICLR 2024
Thus, with probability at least NNy,
21
N
X x
Σ
(-1)*
(k)))))
2
j=1
2N-1
NV.
11+
Σ
(-1)*
2k sup Yj)
2
ISISN
k=
t=
2N-1
S1+
Σ
(-1)*
exp(
sup Yj)
2
ISiSN
k=1
1 2N-1
11 + 2
Σ
(-1)* exp(-2k-
sup Yj) ISiSN
k=
1 1 exp(p supisiSN Yj))
1
(27)
where the fourth inequality holds since
1 2N 1 + Σ 1-1 exp(-2kk) s 11+
1 2 k=1
-(-E)) exp(-
A 0,
2
and since the infinite series 1+ 1/2 (0=((-1)) exp(-2kt) = 1/ tanh()t) + Ite of convergence E (0, oo).
has interval
Next, we show that with probability least 1 -No
1 2N
x=
Σ
2kV,)
k=1
Ł=0
j=1
1 2N
2 21+ 1
Σ
(-1)*
Σ
22 2k
inf Vj)
ISJSN
j=1
2N-1
21+ 2 k=1
(-1)* expl 1ap(-2k 2k inf Yj)
ISJSN
=1 exp(e
-4N
inf Yj) + JSPSN
(-1) * exp(-22
inf Yj)
1SSSN
1
1 + exp( infisisN Yg)
exp( (4W inf Yj)
2
I51SN
0,
(28)
where the last inequality holds since infigyev Y3 2 1 with probability at least 1 - No by (26).
Thus, 27 and (88 together imply 23
20

---

Published as a conference paper at ICLR 2024
To show at least
24 %
first note that, since Y1, Y2.... YN are i.d. random variables, with probability
2V-1 N
E
x
= =1
11)
1
II
(-2kB
2
!!
k=1
8=0
==1
2N-1
NV Σ
8
=1+,
Σ (-1)*
2kg dc)((()))))) dog (f)))
2
!!
1=1
1=
j=1
2N
W
=1+
Σ
(-1)*
2k dct( log dct(( ((()))
1 exep
log ( (0)))
log det (+))))
2N
(-1)a
log det (0)))
8=++1
80
(-1)
log 40
2kog
2=
1 + exp(- log (et((())) log det (b)())
21
1-1)
(-2k(logde((0(s) log del (0))) 8=++1
(-1)* exp( 1ap(-22((M)) (log det (D())) log det (())))
2=
where the fourth equality holds since the infinite series + 1 Σ (-1)* eep(11)eppp-----ee = I+c- has
interval of convergence E oo]. and log det( log det ())) [ 100 2l0g( by
25
Thus,
det (+))
E
det (0)
(et(((())
+ (=)
det (=))
(())
B[X]
+ exp( log (c)) log (et((()())
2k log det (()))
k
W+11
1
+
(-1y* - (oog - log (e))))))))
2
2=
() and log dct(((()))) (2)) log dct( E [ 1 logg())) This proves 24
since N 2 10
Thus, we have shown 24 and 233 in the setting when 1 S Yi s 2 log(7).
The proof 24 and 22 when Y 1/ is identical, except that we use the Taylor series expansion
for I+e t t (-), in place of the infinite series 1+ 1 E (-1) ( exp(-2xpe = tanh(}r)+ + I+e-t which has interval of convergence / (0,
tanh (t) + about O, I4e = Σee c,t which has interval of convergence
21

---

Published as a conference paper at ICLR 2024
The proof for (4 and (23 in the case when V1 > log(2) is trivial, as in this case the algorithm
sets = 1 and
1te-t S yy for t 2 log(?))
Thus,
24 and
(33
logether imply that
det (a)))
EJ (ma)(mxx)x)X))))))
det (0))
S 3Ny
(29)
det((b(0)) + (et(()())
det( (2))
det (D()))
Let = 3 min(max(X (), 1)] xmin
(= (-
x - Bb:o - 110--11000, denote the probability
that the proposed step = accepted at any iteration of Algorithm Then for every (1,...;77 .), have by (29) that
dct()(=))
(E-(()) (e)) X min (-)(()) xelle sell--0lae 013:0o -110=21860) 1) 33N7
1
9i 2 de
(30)
(=)) (())
Let Nn, N,... be the Markov chain where at the beginning of each iteration i + of Algorithm X: 0. Moreover, let Jo...... be the Markov chain defined in Algorithm 2 We define a probabilistic coupling between these two Markov chains. First, define a coupling between the random vector ξ & sampled at every iteration + of Algorithm [ and the random vector Ç: sampled at every iteration i of Algorithm 2 such that & = Çi for all 0. Then, by (30). there exists coupling between the Markov chains Ny, X,... and Vo. Vi.... such that with probability at least 1 - 3TN%
X = 3
Vi {0, ,..
(31)
Thus, letting j denote the distribution of the output Mr of Algorithm | and letting or denote the distribution Vr, we have by (31 that
| vrll & 3T < 2
(32)
Moreover, by Lemma A.5 we have that the distribution Mr of Yr satisfies
|PT π T||TV
- 2
1
(33)
Thus, we have that the distribution j of the output of Algorithm Ia satisfies
|| T||TV I ||I pT|lTV + | T||TV S 8.
C ADDITIONAL COMPARISONS OF APPROACHES TO DETERMINANT COMPUTATION
In this section we explain the key differences between our randomized estimator for the determinant term, and the estimator of Laddha et al. 2020), and why these differences necessary.
Approach of Laddha et al. (2020): In Laddha et al. (2020), the determinant ratio is first estimated by plugging independent random estimates for the log-determinant with mean log (e) log (e) (z) into the Taylor expansion for the exponential function e'. This gives an estimate V with mean V.V]
det (8)) (their Lemma 4.3). Roughly speaking, the authors then suggest plugging
dot( ([)))
this estimate (and another estimate V with mean (0))) p(=-9)) pute an estimate for their Metropolis acceptance rule +((0=0)) = g
det ) into a smooth function g(t, s), to com-
det( H(A) dett
det (()) det (o)))
Each time their Markov chain proposes an update f + 2, it is accepted with probability g(V.
However, the randomized estimator g(V, V) is not an unbiased estimator for the Metropolis rule (e)(H(()) det(A(())) even though E[V] det(H(() and [[ det(H(=)) This is because, in (x))
gi (x)) (e))
det(H(()))
22

---

Published as a conference paper at ICLR 2024
general, EV V) # ((E[V)EE[) since 9V, V) is not linear function in V or V. If one utilizes this randomized estimator in the acceeptance rule of the Dikin walk Markov chain, it can cause the Markov chain generate samples from distribution not equal to the (uniform) target distribution. This problem does not seem to be easily fixable, and requires a different approach.
Our approach to estimating the determinantal term: To obtain unbiased estimator for the determinantal term, our algorithm bypasses the use of the Taylor expansion for the exponential function ex used in Laddha et al. (2020), and instead
two different infinite series expan- sions for the sigmoid function (a Taylor series expansion sigmoid(() = Ei=g ojll centered at O with region of convergence (x, R), and series expansion "at +oo" with region of convergence (0, +20) which is polynomial in e e-") to compute an unbiased estimate with mean equal to (different) Metropolis acceptance rule. The Metropolis rule in our algorithm is proportional to sigmoid (2 log det( - 2 log det(((0)))) To generate a randomized estimator for this accep- tance rule, we compute independent samples ,..... Yn
with mean log det(((()) log (z)),, and plug these samples into one of the two series expansions, c.g. Σi= C.Yi small number terms N. To select which series expansion to use, our algorithm chooses a series expansion whose region of convergence contains ball of radius ğ centered at Yı.
......,YY
Y, for some
Since YI,...,Y YN are independent, roughly speaking, the mean of our estimator is proportional to
N Σ
N
N
E
GY
Y Σ
GEM E|G = Σ c (log det( - -log ( (2)))',
1=
(or to a simmilar quantity if the other series expansion is selected). We show that, w.h.p., the choice of series expansion selected by our algorithm converges exponentially fast in N when Y ..... YN are plugged into this expansion. This implies that, if we choose roughly N = log 4. our esti- mator is a (nearly) unbiased estimator with mean within any desired error y > of the correct value sigmoid { log det det(((2)) (+)) log (e)).. We then show that this error y is sufficient to guarantee that our Markov chain samples within total variation error O(8) from the correct target distribution, if one chooses y = (poly(6, 1/d, 1/))).
23


Reviews:
[{"review_id": "lA8JtCVkur", "reviewer": "Reviewer_6Phh", "summary": "This paper studies sampling from a logconcave distribution on a polytope. To this end, it uses a soft-threshold Dikin walk introduced in MV22. The paper signficantly improves upon the per iteration cost by applying the inverse maintenance techniques from LS15 and LLV20. The technical contribution is in showing how to apply these to the soft-threshold logbarrier.", "questions": "-", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "I wasn't able to check the proofs, but the results suggest that the paper overcomes a technical difficulty prior works such as LS15, MV22, and LLV20 hadn't been able to address. Specifically, showing that the soft threshold logbarrier Hessian is slow-changing in a certain norm is a novel contribution of the paper, and it has consequences to improving the runtime of what is clearly an important problem.\n\n-------- \n\nAfter rebuttal: I'm increasing my score and confidence.", "weaknesses": "I think the paper is, currently, not as self-contained as it should be. I believe this is easily fixable by adding relevant background material in the appendix. I also believe perhaps there may not be too much relevance of this paper in ICLR, and the paper would be much better appreciated in the theoretical CS community's conferences such as SODA, COLT, etc. Perhaps it would be helpful to state a more direct connection to ICLR.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "fkXbsqaAa6", "reviewer": "Reviewer_nR3E", "summary": "The paper studies a fundamental problem of sampling from convex sets, the problem of sampling with respect to a log-concave distribution. Previous work [Mangoubi-Vishnoi, NeurIPS'22] has given an algorithm with $O^*(md)$ iteration and each iteration takes $O^*(md^{\\omega-1})$ time. In this paper, the author improves the cost per iteration to $O^*(nnz(A)+d^2)$. This result directly answers the open problem proposed in [Lee-Sidford15, FOCS'15] which asks whether one can achieve such running time for the case $f\\equiv 0$. To achieve this, the authors make use of the inverse maintenance technique in  [Lee-Sidford15, FOCS'15] and show one can compute the estimation of determinant to high accuracy by cleverly constructing an unbiased estimator and making use of the linear system solver as a primitive.", "questions": "-", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "My general evaluation of this paper is very positive.  I did not manage to check the correctness of the proofs. Condition on the new Metropolis update rule and the log-determinant estimator is correct, I think everything goes through. This paper is technically solid and answers an important open problem.", "weaknesses": "I think the primary weakness of this paper is the presentation. I understand that due to the nature of theory papers, it's hard to present everything important within the page limits. The most interesting part of the paper  for me is how to get a good estimation of log-determinant and I think it deserves some space in the main paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qwSYq61QXu", "reviewer": "Reviewer_CJYs", "summary": "The paper presents the current fastest algorithm for sampling from a log-concave distribution over a given polytope. In particular, they implement a more optimized version of the algorithm proposed by Mangoubi and Vishnoi [1], which is to appear in NeurIPS 2023, the arXiv version of which has been around for more than a year at this point. It is important to understand what [1] does, as the current paper heavily builds on it... \n\nThe algorithms considered in these papers are Metropolis-Hastings algorithms, which are a subclass of Markov chain Monte Carlo (MCMC) methods. For each step of the Markov chain, [1] needs time $O(md^{\\omega-1})$ as it involves some matrix inversion and determinant computation steps, which can all be done in $O(md^{\\omega-1})$ time by using algorithms for these problems which work for arbitrary matrices. The main claim in the current paper is that instead of recomputing these matrix operations from scratch in every step of the Markov chain, we can use the information from the previous step to speed up the computations of the current step.\n\nReferences:\n\n1: Sampling from Log-Concave Distributions over Polytopes via a Soft-Threshold Dikin Walk, Mangoubi, Oren and Vishnoi, Nisheeth K, arXiv preprint arXiv:2206.09384. To appear in NeurIPS 2023: https://neurips.cc/virtual/2023/poster/72502", "questions": "No additional questions", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "Significance: The problem of sampling from a log-concave distribution over a polytope is interesting and also has several ML applications.\n\nOriginality: Clearly, the paper is an advance over the state-of-the-art. However, I am not super sure about how original this is in terms of technical contributions...\n\nClarity: Paper is clear but a bit dense and slightly hard to read, especially if you are not already familiar with Mangoubi and Vishnoi [1]. I re-read this paper after reading [1] and it was much clearer the second time. It might have to do with space-constraints, as the arXiv version of [1] obviously has more space so can go over a lot of things in more detail...\n\nQuality: Overall seems to be of good quality. But I have not checked the correctness of the technical details. As it builds quite heavily on [1] and also a bunch of other results which I am not familiar with, I don't have the expertise to ascertain the correctness of the paper. But at least on the surface, it does seem like there are no major technical issues.", "weaknesses": "I don't have too much to point here, except what I already wrote about clarity before. The authors should try to make the paper more accessible to people who might not have read [1]. I know this is a generic remark but honestly I have no clear idea how to improve the paper in terms of clarity... There were a bunch of typos here and there, some of which are:\n\nIn Theorem 2.1, 3) convex function $f:K \\mapsto \\mathbb{R}^d$. Here $f$ should be $K \\mapsto \\mathbb{R}$. I saw this repeated in a bunch of other places. I think it's a typo but since I saw it repeatedly, both in this paper, and also in [1], I tried to dig around to see if it makes sense for the codomain of a convex function to be $\\mathbb{R}^d$ instead of $\\mathbb{R}$ but in that case, $e^{-f}$ would not be real, but it has to be real for it to be a probability density function?\n\nPage 5, line 3: ball \"or\" radius -> ball of radius, again: $f: K \\mapsto \\mathbb{R}^d$\n\nPage 7, two lines below inequality (3): “arithemtic”", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "NjVdlof2ny", "reviewer": "Reviewer_mHXa", "summary": "This paper considers using a soft-threshold Dikin walk to sample from a polytope contained in a box of radius $R$. The sampling distribution is log-concave and specified by a Lipschitz/smooth function $f$. By using a pipeline introduced in [LLV20], they show how to speed up the algorithm introduced in [MV23] so that each iteration, instead of computing the Hessian of log-barrier using fast matrix multiplication, one can resort to the inverse maintenance data structure due to [LS15]. To this end, they achieve a similar per iteration cost as [LLV20], namely $\\mathrm{nnz}(A)+d^2$.", "questions": "See weakness.", "limitations": "", "rating": 3, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "This paper improves the per iteration cost of [MV23] from $md^{\\omega-1}$ to $\\mathrm{nnz}(A)+d^2$. In most regimes where $m\\geq d$ and $A$ is relatively sparse, this is a strict upgrade from prior state-of-the-art.", "weaknesses": "I'm very dubious of the novelty of this paper. What's the major difference between the algorithm in this paper and [LLV20]? While this paper only mentions [LLV20] sparingly, the algorithmic framework is almost identical. Specifically, [LLV20] shows that for sampling uniformly over a polytope with log-barrier, one can use the [LS15] inverse maintenance to compute a new sampling point. The determinant ratio term can be estimated via an unbiased estimator, which can further be estimated using Taylor expansion together with terms that can be quickly computed using the inverse maintenance data structure. The only difference is that this paper also needs to handle the regularization term, but it is neither surprising nor novel the machinery of [LLV20] also works here. \n\nIt is worth noting that [LLV20] does not provide any proof on why using an approximate solver and samples, the algorithm still converges. This paper provides a very simple argument to show it indeed works.\n\nOverall, I think this paper should provide comprehensive comparison with algorithm in [LLV20]. What's the difference here? What's new? Otherwise, it should acknowledge that the algorithm is largely derivative from [LLV20]. In its current writing, the authors acknowledge the fast linear solver part follows from [LLV20]. What about the determinant ratio part? What's the main difference between your approach for determinant and [LLV20]? What's the novelty of your algorithm?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)