Optimal Underdamped Langevin MCMC Method
Authors: Zhengmian Hu
Feihu Huang
Heng Huang
Venue: NeurIPS 2021 Poster (2021)
Keywords: sampling
underdamped Langevin diffusion
Bayesian methods

Abstract: In the paper, we study the underdamped Langevin diffusion (ULD) with strongly-convex potential consisting of finite summation of $N$ smooth components, and propose an efficient discretization method, which requires $O(N+d^\frac{1}{3}N^\frac{2}{3}/\varepsilon^\frac{2}{3})$ gradient evaluations to achieve $\varepsilon$-error (in $\sqrt{\mathbb{E}{\lVert{\cdot}\rVert_2^2}}$ distance) for approximating $d$-dimensional ULD. Moreover, we prove a lower bound of gradient complexity as $\Omega(N+d^\frac{1}{3}N^\frac{2}{3}/\varepsilon^\frac{2}{3})$, which indicates that our method is optimal in dependence of $N$, $\varepsilon$, and $d$. In particular, we apply our method to sample the strongly-log-concave distribution and obtain gradient complexity better than all existing gradient based sampling algorithms. Experimental results on both synthetic and real-world data show that our new method consistently outperforms the existing ULD approaches.

Optimal Underdamped Langevin MCMC Method
Zhengmian Hu, Feihu Huang, Heng Huang Department of Electrical and Computer Engineering University of Pittsburgh, Pittsburgh, PA 15213, USA
(huzhengmian, huanh wateeeeehhhhh11gmmma com, heng www.g.huangpiit..cm edu
Abstract
In the paper, we study the underdamped Langevin diffusion (ULD) with strongly- convex potential consisting of finite summation of N smooth components, and propose an efficient discretization method, which requires O + d / gradient evaluations to achieve e-error (in /BII-1 distance) for approximating d-dimensional ULD. Moreover, we prove lower bound of gradient complexity 9(N + d' N /Ea), which indicates that our method is optimal in dependence of N, e. and d. particular, we apply our method to sample the strongly-log-oo logeccaavee distribution and obtain gradient complexity better than all existing gradient based sampling algorithms. Experimental results on both synthetic and real- world data show that our method consistently outperforms the existing ULD approaches.
Introduction
Sampling is important research problem in statistics learning with many applications such as Bayesian inference [Р, multi-arm bandit optimization [2]. and reinforcement learning [3. One of the fundamental problems in these applications is sample from high-dimensional strongly- log-concave distribution. Recently, several Markov chain Monte Carlo (MCMC) based methods proposed to solve this problem based underdamped Langevin diffusion (ULD). This con- tinuous diffusion process converges to the target distribution exponentially fast. Thus, the methods approximating a ULD process could be used sample from the target distribution within certain accuracy.
Multiple discretization methods have been proposed for approximating ULD. Among them. the Euler-Maruyana. discretization [4] is the simplest one but generates the largest error. Recently the left point method (LPM) [ was introduced to fix the gradient term in ULD to be the gradient at k-th iteration, and then integrate the new linear stochastic differential equation (SDE) with a small time-interval. Subsequently, Shen and Lee (6] proposed randomized midpoint method (RMM) with smaller error. There are also discretization schemes based on splitting (7) or Runge- Kutta complexity lower bound for simulating d-dimensional ULD. Under the assumption that the full gradient oracle fx()) is evaluated at most n times, they show lower bound for worst-case error by perturbation analysis, which matches the discretization eтror upper bound of RMM in the dependence of d and n.
method g 9 More recently, Cao al. [1] derived information-
Although the ULD-MCMC methods with full gradient oracle are largely understood, many real-worke applications involve summation form of potential function and large-scale data, which leads to the need of stochastic gradient methods. The vanilla stochastic gradient methods have been used to replace full gradient [5). Albeit the computational cost for each iteration reduced, the variance of
Althhough this method mostly just denoted as ULD-MCMC we adopt the name LPM to distinguish from other discretization methods. The name comes from the fact that gradient is evaluated at the left point of the time interval.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

Summary of gradient complexity of sampling methods, which is defined as number of gradient evaluation of (x) needed to sample from m-strongly-log -stoongy-log-cconcvee concave distributions up to evd/m m accuracy in 2- 2-Wasserstein distance where S 1 is the target accuracy, d is the dimension. ULA, LPM, RMM, and ALUM are full gradient methods, therefore, gradient complexities for them sample size N times the iteration complexities. Only dependence on d, E, and N are shown below. The dependence of batch size is made clear in Table The dependence of condition number k is discussed in Section 7
Table 1:
Algorithms
Gradient complexities
Unadjusted Langevin Algorithm (ULA) 15
2(Ne-2)
LPM 17
O(Ne-1)
RMM 6
Ö(Ne $)
ALUM (Ours)
Ö 9)
Stochastic Gradient LPM (SG-LPM) 5
Ö( 3)
SVRG-LPM
O(N + N
CV-ULD IB
O(N +6-
SVRG-ALUM (Ours)
Ö(N + N
SAGA-ALUM (Ours)
Ö(N NiE
stochastic gradient is much larger than the discretization error and therefore degenerates the overall
performance. Previous works [2] used stochastic variance reduced gradient (SVRG) [3] and SAGA [14] instead, but the gradient complexities of these methods are still worse than the full gradient RMM in terms of dependence on accuracy Thus, there exists a natural question:
What is the optimal ULDD-MCC method with sum-decompossable potential?
In this paper, we focus on optimal dependence of dimension d, components number N and accuracy e in gradient complexity for estimating a ULD process. We answer this question by two parts. We first provide novel ULD-MCMC method and derive the corresponding complexity upper bound in Sections 4 and After that, we analyze the worse case error and show that the lower bound matches the upper bound in Section 6 The major contributions of our paper can be summarized as follows.
1. We propose new full gradient LLD-MCMC method, called AcceLerated ULD method (ALUM), whose discretization error has the same order dependence on dimension d, step size RMM. Although RMM already has optimal asymptotic complexity in full gradient setting, ALUM still of practical interest. Compared with RMM, which uses two gradient evaluations each iteration, ALUM gradient less frequently and only requires one gradient each iteration to achieve constant speedup.
2. We further propose VR-ALUM methods, including SVRG-ALUM and SAGA-ALUM, which utilize the unbiased variance reduction techniques in ALUM under sum-decomposable setting. We show that these methods achieve better gradient complexity than all existing gradient based MCMC approaches. These gradient complexities for sampling from a strongly- crronglylog-cconcve concave distribution are compared in Table
3. We derive information- based lower bound worst-case error for estimating ULD process with only gradient oracle and weighted Brownian motion oracle. We show that order to achieve approximation accuracy, Q(N + d' N single component gradient evaluations needed. This lower bound matches the upper bound for VR-ALUM in terms of dependence of dimension d, sample size N, and accuracy e. Therefore, our VR-ALUM methods are indeed optimal for estimating ULD process under sum-decomposable setting.
Pn Zou et al. [1]. the authors call their method SVR-HMC. However, their method is not based on Hamiltonian Monte Carlo (HMC), based on ULD. Their method just applying SVRG to replace full gradicnt in I.M.
Further explanation and comparison is shown in Appendix A.2

---

2
Related work
Though we mainly use the gradient oracle in this paper, there is also a large body of sampling algorithms which leverage zeroth order potential value oracle. For example, Metropolis-Hastingg accept-reject step could be used to ensure an MCMC converges to a stationary distribution equal to the target distribution
[9 20 21 22 23 and the linear convergence can be obtained. For example, ULA, can achieve : eπror in total variation distance within Odd (1/))) steps under certain initialization [26] Note that purely gradient based sampling algorithms typically only converge sub-linearly because the stationary distribution is different from the target distribution.
MALA (24] 25), as Metropolis-adiie
Other diffusion processes could also be used for constructing MCMC. ULA, a discretization of Langevin diffusion (LD) [27 24), is the first commonly used gradient-based MCMC. We will provide an overview for some LD-based and ULD-based methods in Appendix A. Later, higher order diffusion process is also discussed. Although the path is even smoother than ULD, currently it is not clear whether such higher order smoothness of path could be leveraged to accelerate the convergence without extra assumptions potential function. Mou et al. [28] studied third-order diffusion, but their acceleration requires special structure or high order smoothness of potential. The more general diffusion process, which converges the target distribution, studied in Ma al. [29]
Preliminary
The problem of sampling from a strongly-log-coneave distribution involves a probability density function p" (z:) defined on real vector space IR A corresponding potential function f(x) - -log(p" (z)) can be defined such that f(x) is strongly convex. One way solve this sampling problem is constructing a Markov chain that converges to a stationary distribution that is the same as or similar to the target distribution.
Appendix [ we introduce several such Markov chains as discretization of certain continuous stochastic processes, and we roughly analyze discretization error of these methods with pointing out the bottleneck.
We use Ō(f) Off) log (f) to omit logarithm factor. I|-I12 means the Euclidean norm. We define a norm of random vector as ||-I|1, = / E1-11 Next, we list our assumptions on potential fax
Assumption 1 (Sum-decomposablee, f(x)= - M fi(z), where integer N is the sample size.
Assumption 2 (Smoothness). Each function fi is twice differentiable on Rd and there exists constant L >> 0, such that f((xx & t for any Id where I is the identity matrix. It can be easily verified that f(x) L-smooth.
Assumption 3 (Strong Convexity). There exists a constant m such that: f(x) - f(y) N (Vf x y) + me
x - yll.. We define the condition number == L/m. 2
We finally define the 2- Wasserstein distance between distributions. For any pair of probability measures μ and y the same parameter space, transference plan & between and v is joint distribution such that the marginal distributions on two sets of coordinates are μ and v. respectively. We denote I(u, v) the set of all transference plans, and define the 2- Wasserstein distance between and y as follows:
(2, inf cer(u.u) flle: -0llay(( v)-
4 AcceLerated ULD-MCMC (ALUM) methods
In this section, propose class of AcceLerated ULD-MCMC (ALUM) methods based approxi- mation of the continuous ULD process:
dxx: Vidt, dV. f - Vidt + 2
The gradient term in the above SDE could be highly non-linear such that closed form solution not available. Thus, propose estimation of SDE solution at time point h which only requires
3

---

gradient evaluation at one single point:
x (o) = Xo + 1/1 (h)Vo har (h-ihhaahVif (0))
an cx.
V = vo(h) - hop(h Vf + ev.(o.A)
(2)
h
X Xo + 1 (ah)Vo + + (z.(0.a)]. ah
where the random variable a is uniform random variable in [0, 11. ah(+) and noise terms [x.[o.N]
Co,(0,A]+ (0,hh defined as follows !
(h) = - (1) 7 "1 (1--
1
cx. 22 / 1h1h - s)dB,, cv. = 227 1/' voln wo(h-8))1 - s)dB..
4.1 Full gradient ALUM
In this subsection, we propose full gradient ALUM method, which is shown in Algorithm Specifically, the update rule is obtained by setting x # (0)) to X V(() where the ULD starts from X, 2
(0) Vo okk (0) and h is small step size.
Next, we compare our method with both RMM and Nesterov's acceler- ated gradient method (NAG) which has been shown optimal for gradient based optimization. Compared to RMM: By comparing the formula in (19) and (2). it is easy to see that ALUM and RMM differ in one single term
Algorithm 1: Full gradient ALUM Method
Input: Initial point (a vn (o)), parameter 7, iteration number K, and step size h > 0. for to K - do
Randomly sample uniformly from [0, 1];
Generate c. 10, Co.(0.N), (z, according to
Appendix A.5
x/ = x + th (axh) (0)
+ (0.,0..o)).s
Our motivation of dropping this term to reduce the computations. By
Calculate full gradient V (ax
deleting this term from the algorithm, only one gradient evaluation a ran- domized midpoint is needed at each iteration. This compares favorably to the RMM, which requires two gradi- ent evaluations at each iteration.
+++
+ th (h)uk -
= vo(h)e 2 (f
k
end for
Output: zx
However, could this method still estimate ULD accurately with only half gradient evaluations?
We firmly answer this question with rigorous error analysis in Theorem 2 Roughly speaking. if we only consider the error's dependence on step size h, the answer is yes. In this section, we only give an intuition how that possible.
In Appendix I. we show that RMM has very low bias and high variance. Therefore, the bottleneck the variance, and increasing the bias slightly would not degenerate the overall performance much if the bias is not larger than the eтror introduced by variance. The dropped term - 2/2 (ah) (Xx) has norm Ohh2), and the coefficient of gradient in Vn, Ohh). Thus, the bias introduced in single step is O(h3) when we use step size h. After accumulating for T/h O(h iterations, the bias is O(h²), which still has better dependence on h than the square root of variance (hh/22)
With increased bias, the complexities of ALUM and some other variants have the same depen- dence on as RMM, the highest order dependence on k could deteriorate as discussed Section 7
Compared to Nesterov's accelerated gradient (NAG) method: NAG is gradient based optimiza- tion method 33 Based on [31]. NAG could be formulated as the following momentum method:
[2,[0o.a]] is obtained by simply substituting &s with ah the definition ex.(o.]]- *e can drop (z.[o.A).A and A), in (2) to derive other variants of ALUM. The bias of these variants increases to ( which is still no larger than the square root of variance. The maximum order dependence on h is the same, thus the final iteration complexity has the same dependence
RMM.

---

#+++ at + Cg 4 V
V (c) ai + Cg (4)
1k+1 CA
We can see that the main differences between Algorithm and (4) are coefficients and additional Gaussian noise terms. The coefficients in NAG are set in deterministic way. However, the coefficients in ALUM could be random values from fixed distributions. The additional noise terms in ALUM come from the Brownian motion term in ULD and are necessary for sampling methods.
Despite the differences, both NAG and ALUM use one gradient at each iteration, and both of them take big jump along current momentum direction to calculate the gradient instead of directly computing the gradient at current iterate.
4.2 Variance-reduced ALUM (VR-ALUM)
In the subsection, we propose variance-reduced stochastic ALUM (VR-ALUM) methods based on two common unbiased o cariance-reduced techniques: stochastic variance reduced gradient (SVRG)
B3] and SAGA [1]. VR-ALUM comes from simply replacing all full gradient in Algorithm with
gradient estimation V.. We show Algorithm 2 in the Appendix A.6 due to the limit of space.
Next, briefly introduce these variance reduction techniques. SVRG utilizes the follow- ing gradient estimation (B& is the batch of k-th iteration and b is the batch size): v
b Eien, (v (x) fz() EE Vf((z). The full gradient Σ V((()) and point updated after every evaluations of VVRGG We call the hyperparameter
epoch length.
SAGA estimates the gradient in the following way, where Φ is set as zk (c)V if and only if E By, otherwise d/+1 = $ A -/- EieB, (VS(aa V)i(E)) + Σ V((())) SAGA does not re-compute but stores the latest gradient information ((((tt for each fi. Therefore, SAGA does not introduce extra gradient evaluation except for initialization, but has much larger storage requirements.
Both SVRG and SAGA are unbiased, which means En, V. Vi(a (c)V expectation over random batch at k-th iteration. Moreover, both SVRG and SAGA reduce the mean- squared error of gradient estimation and satisfy bounded MSE property proposed in Appendix B.3
where EB, means
Theoretical analysis
We provide non-asympiotic upper bounds on sampling error and discretization error for our methods, including full gradient ALUM and VR-ALUMS. The proof shown in Appendix B Throughout this section, we assume W2 (po. p²) ) = O(1) /d/m when deriving asymptotic results. This means the initialization is not too far away from the target distribution, and can be achieved under multiple setting as discussed in Appendix A.7
5.1 Convergence analysis of full gradient ALUM
Recall that ALUM can be used for solving two different but related problems: strongly-log-conced sampling and approximating the ULD. We show the upper bound for sampling error Theorem and the upper bound for approximation error in Theorem separately.
Theorem Suppose Assumptions 1ioB LoBhool hold. Given an initial distribution po po(x)), we initialize ALUM
with random x. 10) based on probability po and random v/ from standard Gaussian distribution.
Assume L = 1 and let 7 = 2, pk be the distribution of ((0) and p" be the target distribution. Assume use step size h < 22 After running the ALUM for k iterations, have the following upper
bound of sampling error in 2- 2-Wasserstein distance:
mh & (2 p²) + 12
hdd
W2(Pk, p²) < 2(1 4
(5)
'The reason we choose unbiased variance reduction instead of biased one is that bias accumulates quicker than variance, therefore generates higher overall error.
5

---

Table 2: Iteration complexities for full gradient ALUM on both sampling and approximation problem.
Problem
Accuracy
Step size h
Iteration complexity
Sampling
evd/m in W2
h == (minne m)) in L2 h O(min( m (m3edd d m))
(maax(x/e3, x22) (max(TT<B 3d TR))
Approximating
Theorem 2.
With same assumptions and setup in Theorem consider Xth which comes from continuous ULD starting from the same position as an))) with probability We have the following upper bound of discrelization error:
h W2(Po.p") 31 m
hdd
| (0) Xah S 38/
(6)
m
Remark 1. The assumption L 1 actually does not limit the availability of the algorithm. For any which satisfies L
function ((x) with L + 1, we can define another function - [ ALUM with step size S 22L and initial momentum from the standard Gaussian distribution can be used to sample x' from the distribution (a') x exp(-f' (x')). The sample x from distribution p(z) & exp(-f(z)) could be obtained by transform æ x'. This is essentially the same as
m
directly incorporating L into ULD process in Appendix A.3
Thus, we can derive the iteration complexity for both problems. For sampling problem, we define the iteration complexity as the number of iterations needed to achieve W2(PK. .p) S sv d/m with certain step size h. For approximating the ULD at given time point T, we define the iteration complexity as the number of iterations needed achieve || Xx l1 S with step size h = T/K. The results are shown in Table Detailed derivations can be found at Appendix A.8
5.2 Convergence analysis of variance-reduced ALUM (VR-ALUM)
We show the upper bound for sampling error in Theorem and the upper bound for discretization error in Theorem separately.
Theorem 3. With same assumptions in Theorem 7 we use SVRG-ALUM with epoch length
[N/b] or SAGA-ALUM. We introduce an extra assumption h I 2304 b'mN-2. We have the following upper bound of sampling error in Wasserstein distance.
mh W&2 p²) y
h
hdd
W2(PE. p²) s 2
92
Vc
W2(Po + +12
D)
(7)
m b
b
m
The constant is defined as 1 for SVRG-ALUM and 2 for SAGA-ALUM.
Theorem 4. With the same assumptions in Theorem
we consider X which comes from continuous ULD starting from the same position as F with probability 1. We have the following upper bound of discretization error.
||
( 1
Xanlle, S <38 92VC
N
h W2(Po, +0) (31 57 m
N
hid
b
m
Remark 2. Theorems 7 F are specializations of more general results in Appendix gradient, SVRG and SAGA are unified under the framework of bounded MSE property that defined in Appendix B.3] A unified approach not only simplifies the proof. but also indicates that our analysis could easily generalize to other gradient estimations that satisfy bounded MSE property.
B. where full
Similar to the full gradient case, we derive the iteration complexity K in Appendix A.9 Moreover, we define gradient complexity as the number of single component gradient evaluation ()((z) needed to achieve certain accuracy. We show the results in Table 3 and add the derivations in Appendix A.9 We finally simplify the result by only considering the dependence of d, N, b, and E.
Corollary 1. When b < O(N5), the gradient complexity of SAGA-ALUM and SVRG-ALUM for sampling problem is Ö(N + N g- 8) and their gradient complexity for ULD approximation problem O(N + dt N +))

---

Table 3: Gradient complexity for SAGA-ALUM and SVRG-ALUM.
Problem
Accuracy
Gradient complexity
Sampling
6 e/d/m W2
O + (bs + (1 + bn?
Approximating
in Lg
O(N + T(nb Trade
6
Information-based complexity
We first declare our setup for the problem class and accessible information. Then we show information- based lower bound for approximation error and oracle complexity. The proof is shown in Appendix C
6.1 Setup
Class of possible potential functions: We denote the class of possible potential functions as U N,,,N,m.b L) = ((......,f. V2(x)) & (L/N)I,m ml & V D fi(() and
|| [X]]l2 S
vd/m3 where m < L and (X) x exp(- En fi(X)..
Besides Assumptions to used in the previous sections, we introduce a new assumption that the target distributions have mean in a ball of constant radius around the origin. This is necessary because the lower bound could be arbitrarily large if the mean is far away from the initialization point.
Solution mapping: We denote the probability space for Brownian motion in the ULD process as (M, Σ, P). Then, the true solution of ULD process starting from the origin at time can be denoted the solution mapping XT: (wi, U) (MI uf) + XT (w, U) Rdd
Gradient oracle: For a given set of potential functions = . fN), the single component (x) E Rd.
gradient oracle Ty (1) E [N] TRd
Brownian oracle: We assume the Brownian motion at given time t > 0 could be evaluated with oracle (w) Rd for any event w E Ml. We further assume that the weighted Brownian motion is also admissible. (w) = Jo" e (B.(w).
Deterministic algorithm: deterministic algorithm starts from empty information Iu O. At i-th step, one oracle and corresponding parameters are picked by certain procedure. If the gradient oracle picked, the algorithm will generate an index and a point æ. If the weighted gradient oracle is picked, the parameters are order 0 for weighted Brownian motion and time L. The picked oracle and parameters are represented by Φ) (T:). The evaluation result represented by Y (c)(T,) w,w, U' where E MI and U U. The picked oracle, parameters, and result will be stored as new information, hence 4++1 = (4, (1,), ((,(L,,,,.)) U)).
We consider deterministic algorithms that stop at n-th step. The final estimation is generated with mapping (T,) E Rdd. We expect (I,) to be as close as possible to true solution (, U).
We denote deterministic algorithm mapping A from MI x l to Rd and A(w; U) (I,) for some & and Y. The family of all such algorithms denoted by Aet
We consider another probability space (M, Σ, P) as the source of ran- domness. A randomized algorithm A with steps is mapping from M x M x l to Rd such that A(, E Apes for any W E M. The family of all such randomized algorithms is denoted by An. Worst error of algorithms: We care about the following worst-case error for any possible algorithms:
Randomized algorithm:
Exu
inf sup Ewer (x( U) A(w, A(w,5U)) U)2 AEA Ueu
We always assume T >>> 0 to avoid trivial cases.
6.2 Lower bounds
We first provide lower bounds worst-case estimation error.
7

---

Theorem 5. When n < N which means that gradient evaluation number is less than components number, we have An 2 dO,, where C₁ is positive and independent of N, and n.
The above theorem is based on the fact that no algorithm can accurately estimate the global minimum point of the sum of quadratic potentials E Sc(æ)) with information from only N - components.
Theorem 6. When the gradient evaluation number is a multiple of N, we have Exn 2 dC2 N² n² where C2 is positive and independent of d, N, and
The method number of gradient evaluations cannot distinguish a class of perturbed quadratic functions.
to prove the above theorem is to show that algorithm that uses limited
Next, we use the above lower bounds on error derive a lower bound on gradient complexity.
Corollary For small enough target accuracy E such that 2² < <CC,, in order to achieve C.A,, u1 S E. the minimum number of single component gradient oracle evaluations is SN(N + -
This lower bound matches the upper bound in Corollary D of pariance-reduced ALUM for estimating a ULD is optimal in the dependence of d, components number N, and approximation accuracy E.
which indicates that gradient complexity
7 Optimality
In this section, we discuss in what sense our ALUM is optimal (or not), and point out possible improvements for future work.
Optimal for approximating problem: ALUM can be used for two tasks: (I) estimating a ULD process, (2) sampling from a strongly-log concave distribution. ALUM is only optimal for the first task, and is not necessarily optimal for the second one. For example, we see in Section 5 that there exists logarithm factor in the gradient complexity for sampling problem. We believe some proper adaptive step size method could cancel that extra term.
Dependency on K: ALUM is only optimal on the dependence of d. N, e, but not optimal in k dependence. Actually both our method and our analysis may not be optimized for the dependence of k. Currently, there is no information- k, therefore no matter how good the dependence on K is, it not sufficient to say is optimal. For sampling problem, it is believed that Ommin(k. (K. d')) is the "natural barrier" [6] for iterations needed for achieving W2 < e s/dd d/m m with
lower bound with clear dependence on condition number
1/2. Currently, the best dependence on is achieved by RMM with 7. For ALUM, the maximum dependence is O(k²). This dependence is worse for two reasons. First, the analysis is not optimized for k. Based on the analysis of bias and variance in Appendix A. 10 we conjecture that the dependence could be improved to (x33/2) Second, the method is not optimized for k. We save one extra gradient compared to RMM, at the price of slightly increased bias. This degenerates the dependence K.
Despite the max dependence order being larger than 1, in the segh-precision regime, when e is small enough, the term k</$ is dominant, therefore, both full-gradient ALUM and full-gradient RMM achieve (xx) dependence.
We note the complexity for VR-ALUM increases the dependence of K 4/3 in high-precision regime. We conjecture that this is an artifact of error-based analysis compared to momentum-based analysis, and could be improved by some tighter analysis.
Assumptions: Finally, we point out that our method and analysis are based on Assumptions Ito and the gradient oracle. It is natural to obtain better algorithm by introducing new assumptions or new oracle.
Many higher order integrators g leverage higher order smoothness assumption to reduce discretiza- tion error. Assumptions on structure of potential (28) have also been shown to make acceleration possible together with high order diffusion process. Apart from the gradient oracle, many Metropolis- adjusted algorithms 24 26 leverage zeroth-order oracle to achieve línear convergence. Second order
oracle [17] has also been incorporated into estimating ULD.
8

---

Experiments
10
10/
Trijetsc Eive
Trijetury Eiror
10-1
10
10-4
PG
TPM
10-9
F6
LPM RMM
SG SVRG SAGA
RMM
SG SVRG SAGA
ALUM
ALUM
10
10-5
Gradiet Evaluation Number
Gradient Evalaztion Number
(a) Gaussian model
(b) Logistic regression on australian dataset
Figure 1: Trajectory error for different algorithms" full gradient, SG means stochastic gradient.
The vertical lines are the error bar. FG means
In this section, we compare different algorithms for estimating a ULD process on both Gaussian model and Bayesian logistic regression model. For the Gaussian model, the potential is defined as:
1 dd - x- x),
[1 = 2N
(10)
where d; and > is generated randomly to satisfy = 5.N 100,m and 10. For logistic regression model, the potential is:
N
f(()) =
172
x
+
log(1 + exp(( x)),
(11)
2N
where m is set such that x = 10 and Vi. a, are data points in australian dataset from LIBSVM
32
We calculate the trajectory error which is defined as R Σ VT - 24,12 12 U 2. where 24 and 8i are generated by ALUM, VR-ALUMS or other algorithms. zyg and v/. come from a reference path that is very close true solution. We specify how generate this reference path in Appendix [..]] The reason for us to average the error along the path instead of just reporting error at the final iterate that we need more data points to reduce the variance.
Figure shows the error ALUM, LPM, and RMM with full gradient, stochastic gradient, SVRG and SAGA B The detailed setup can be found in Appendix D.2 We summarize the messages in Figure as follows:
SAGA-ALUM achieves the best efficiency in the sense that with same number evaluations of single component gradient fzx) SAGA-ALUM has smaller discretization error than any other algorithms. Results on more dataset are shown in Appendix [D.
Variance reduced algorithms constantly outperform full gradient algorithm of the same type large margin.
For full gradient method, the discretization of ALUM and RMM has similar asymptotic dependence on gradient evaluation number, and they are hetter than LPM. This phenomenon will be shown clearer Figure 2
ALUM achieves constant acceleration compared to RMM by saving one gradient evaluation per iteration.
"Most of the bars small be visible. The SG-ALUM and SG- SG-LPM highly overlap. *Currently there is no theoretical result for RMM with stochastic gradient, SVRG or SAGA. We also didn't provide theoretical result for ALUM with stochastic gradient. However, this doesn't 't prevent us from evaluating them experimentally.

---

We also show the relationship between discretization error and step size in Figure 2 Only full gradient method are shown here due to the limit of space, and more results for VR-ALUMS are shown Appendix [D.4
10-1
10-3
Trajectory
LPM RMM ALUM O(A
10
Our analysis in Theorem 2 gives upper bound of discretiza- tion error of ALUM as ( Figure 2 shows that our analysis is tight. Moreover, the discretization error of ALUM is almost the same as RMM. Due to the fact that ALUM requires only half gradient evaluations per step than RMM, ALUM achieves better efficiency.
10-T
((A)
10
10
Step size
Figure 2: Discretization error of full gra- dient methods on australian dataset
Next, discuss the effect of batch size. According to our theory in Corollary [l our method is not sen- sitive to this hyperparameter intro- duced by SAGA SVRG when batch size is relatively small, as the gradient complexity remains the
Taajotor Erve
10
Bajotun Eace
same for b O(N2/3). We verify that in Figure 3 where the sampling efficiency is almost same for small batch sizes, and only deteriorate for very large batch size. We give fur- ther discussion on why our method not sensitive step size in Ap-
b-330
10m2
Gradient
Gradient
(a) SVRG-ALUM
(b) SAGA-ALUM
Figure 3: Discretization eπror for SVRG-ALUM and SAGA-
ALUM with different batch sizes on australian dataset.
pendix D.5
Finally, we apply ALUM and VR-ALUMS to sample from target distribution in Appendix D.6
9 Conclusion
In this paper, we propose class of MCMC methods for finite sum form of strongly-conver potential. Our methods are proven to be optimal in the sense that the discretization error has the best possible asymptotic dependence on dimension, number potential summands, and number of gradient evaluations. Experiments on both synthetic and real data verify the superior performance of our algorithm. We also discuss possible improvements for future work.
Acknowledgments and Disclosure of Funding
This work partially supported by NSF IIS 1845666, 1852606, 1838627, 1837956, 1956002, OIA 2040588.
References
[1] Andrew Gelman, John B Carlin, Hal Stem, David B Dunson, Aki Vehtari, and Donald B Rubin. Bayesian data analysis. CRC press, 2013.
[2] Daniel Russo, Benjamin Roy, Abbas Kazerouni, Osband, and Zheng Wen. A tutorial on thompson sampling. arXiv preprint rrXiv:1707.02035, 2017.
Deepak Ramachandran and Eyal Amir. Bayesian inverse reinforcement learning. In IJCAI, volume pages 2586-2591, 2007.
[4] Peter E Kloeden and Eckhard Platen. Numerical solution of stochastic differential equations, volume 23, Springer Science Business Media, 2013.
[5] Xiang Cheng, Niladri S. Chatterji, Peter L. Bartlett, and Michael I. Jordan, Underdamped langevin meme: A non-asymptotic analysis. In Proceedings of the 31st Conference On Learning Theory, volume 75, pages 300-323, 2018.
The ALUM and RMM highly overlap.
10

---

[6] Ruoqi Shen and Yin Tat Lee. The randomized midpoint method for log-coneave sampling. In Advances in Neural Information Processing Systems, volume 32, 2019.
[7] Benedict Leimkuhler, Charles Matthews, and Gabriel Stoltz. The computation of averages from equilibrium and nonequilibrium langevin molecular dynamics. IMA Journal of Numerical Analysis, 36(1):13 13-17-3333 79, 2016.
GN Milstein and Michael V Tretyakov. Quasi-symplectic methods for langevin-type equations. IMA journal of numerical analysis, 23(4):593- 626, 2003.
[9] James Foster, Terry Lyons, and Harald Oberhauser. The shifted ode method for underdamped langevin meme. arXiv preprint 10Xiv:2
2021.
[10] Yu Cao, Jianfeng Lu, and Lihan Wang. Complexity of randomized algorithms for underdamped langevin dynamics. arXiv preprint arXiv:2003.0 09900, 2020.
[11] Difan Zou, Pan Xu, and Quanquan Gu. Stochastic variance-reduced hamilton monte carlo methods. In Proceedings the 35th International Conference Machine Learning, volume 80, pages 6023-6033 2018
[12] Zhize Li, Tianyi Zhang, Shuyu Cheng, Jun Zhu, and Jian Li. Stochastic gradient hamiltonian monte carlo with variance reduction for bayesian inference. Machine Learning, 108(8):1 1011))170111127) 701-11277 2019.
[13] Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In Advances in neural information processing systems, pages 315-323, 2013.
[14] Aaron Defazio, Francis Bach, and Simon Lacoste-Julien. Saga: A fast incremental gradient method with support for non-strongly convex composite objectives. In Advances in neural information processing systems, pages 1646-1654, 2014.
[15] Alain Durmus and Eric Moulines. High-dimensional bayesian inference via the unadjusted langevin algorithm. Bernoulli, 25(4A):2854 28822 2019.
[16] Amak Dalalyan. Further and stronger analogy between sampling and optimization: Langevin monte carlo and gradient descent. In Conference on Learning Theory, pages 678-666 778-689, 2017.
[17] Arnak Dalalyan and Lionel aoounDd On sampling from a log con-cone density using kinetic langevin diffusions. Bernoulli, 26(3):1956-13 1988, 2020.
[18] Niladri Chatterji, Nicolas Flammarion, Yian Ma, Peter Bartlett, and Michael Jordan. On the theory of variance reduction for stochastic gradient monte carlo. In Proceedings of the 35th International Conference Machine Learning, pages 764-773, 2018.
[19] Claude JP Bélisle, Edwin Romeijn, and Robert I. Smith. Hit-and-run algorithms for generating multivariate distributions. Mathematics of Operations Research, 18(2):555-666, 1993.
[20] László Lovász and Miklós Simonovits. Random walks in convex body and an improved volume algorithm. Random structures algorithms, 4(4):359-412, 1993.
[21] Gareth Roberts and Richard L Tweedie. Geometric convergence and central limit theorems for multidi- mensional hastings and metropolis algorithms. Biometrika, 83(1):95-111 110, 1996.
[22] Kerrie L Mengersen and Richard L Tweedie. Rates of convergence of the hastings and metropolis algorithms. Annals of Statistics, 24(1):101-12) 1211 1996.
[23] László Lovász and Santosh Vempala. The geometry of logconcave functions and sampling algorithms. Random Structures & Algorithms,, 30(3)3007-358, 2007.
[24] Gareth O Roberts and Richard Tweedie. Exponential convergence of langevin distributions and their discrete approximations. Bernoulli, 2(4):341-363, 1996.
[25] Gareth O Roberts and Jeffrey S Rosenthal. Optimal scaling of discrete approximations to langevin diffusions. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 60(1):255-26 268., 1998.
[26] Raaz Dwivedi, Yuansi Chen, Martin Wainwright, and Bin Yu. Log-concave sampling: Metropolis- hastings algorithms are fast! In Conference Learning Theory, pages 793 733777, 2018.
[27] G Parisi. Correlation functions and computer simulations. Nuclear Physics B. 180(3)378-338 1981.
11

---

[28] Wenlong Mou, Yi-An Ma. Martin J Wainwright, Peter L Bartlett, and Michael I Jordan. Highh-rder langevin diffusion yields accelerated meme algorithm. Journal Machine Learning Rexearch, 22(42):1-44 2021.
[29] Yi-An Ma, Tianqi Chen, and Emily B. Fox. A complete recipe for stochastic gradient meme. In Proceedings of the 28th International Conference Neural Information Processing Systems Volume page 2917-2925, 2015.
[30] Yurii Nesterov. Lectures on convex optimization, volume 137. Springer, 2018
[31] Ilya Sutskever, James Martens, George Dahl, and Geoffrey Hinton. On the importance of initialization and momentum in deep learning. International conference on machine learning, pages 1139-1147, 2013.
[32] Chih-Chung Chang and Chih-Jen Lin. Libsvm; library for support on intelligent systems and technology (TIST), 2(3) 27, 2011.
machines. ACM transactions
[33] Simon Duane, Anthony D Kennedy, Brian Pendleton, and Duncan Roweth. Hybrid monte carlo. Physics letters B, 195(2):11-222, 1987.
[34] Yi-An Ma, Niladri Chatterji, Xiang Cheng, Nicolas Flammarion, Peter Bartlett, and Michael T Jordan. Is there an analog nesterov acceleration for meme? arXiv preprint urXiv:I902.00990 2019.
[35] Arnak S Dalalyan and Avetik Karagulyan. User-friendly guarantees for the langevin monte carlo with inaccurate gradient. Stochastic Processes and their Applications, 129(12):5278-5311 2019.
12


Reviews:
[{"review_id": "cuJyuB2phA", "reviewer": "Reviewer_VSLc", "summary": "This paper proposes a randomized algorithm for approximating a trajectory of underdamped Langevin dynamics, similar to the RMM method cited in the paper. It further uses variance reduction method to improve the error over sum-decomposable problems. An information theoretic lower bound is established, which the proposed method matches.", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "This paper is interesting and technically sound. Even tough the lower bounds are only in terms of approximating a specific trajectory---instead of for sampling a posterior distribution---they are still useful for numerical methods with Monte Carlo strategies. I quickly skimmed through the proofs, which seem to make sense. \n\nMore specific comments below:\n\nIn remark 1, the authors stated that the function can be rescaled for arbitrary L. It is worthwhile to point out how does the step size h scale with L as well.\n\nThe notation of equations (2) and (3) causes some confusion: I first thought that the authors are analyzing an algorithm with a shriking step size. It turns out that the step size is a constant. I would suggest chaning t to something else (perhaps \\eta). \n\nIn variance reduction, m has been double used as both the epoch length and the strong convexity. For brevity, the authors can also choose a suitable (if not optimal) epoch length and define it in the algorithm.\n\nIt seems that the improvement of N dependence provided by the current method is due to the variance reduction method. It seems to also work on the RMM method. Is that the take away from the experiments? Also it is mentioned that the error for ALUM, LPM, and RMM with full gradient, stochastic gradient, SVRG, and SAGA are compared. I can't seem to find the curve for LPM with stochastic gradient.\n\nThe overall structure of the paper is clear. However, the literary quality of the paper can be improved. The authors should also be aware of grammatical errors. For example, in the abstract, \"consists of finite summation\" should be changed to \"consisting of finite summation\".", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "aanQMzbFwoG", "reviewer": "Reviewer_aKpK", "summary": "The paper studies the underdamped Langevin diffusion for strongly-convex potential consists of finite sum of smooth components. The authors propose ALUM, which achieves optimal asymptotic complexity in the full gradient setting, and only requires one gradient at each iteration compared to RMM. For the stochastic gradient setting, VR-ALUM methods are proposed which improve over the previous methods in gradient complexity. The authors also accompany with a lower bound showing that the proposed methods are optimal in related parameters for the approximation problem. \n", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "Overall, I believe this work is a nice complement to existing underdamped Langevin MCMC literature. The presentation is good with clear upper and lower bounds and fair experiments. My comments are summarized as below:\n\n1. The paper studies two problems: sampling in $W_2$ distance and approximation in Euclidean distance. It is known that in the full gradient setting, the algorithm matches RMM, which achieves optimal asymptotic complexity for sampling problem, and in the stochastic setting, the proposed algorithm achieves optimal complexity for the approximation problem. Since sampling problem is considered an important target in MCMC, it would be great if the authors can further comment on the optimality of the proposed algorithm for sampling problem in the stochastic setting beyond the paragraph around line 252: is the gap conjectured to be only logarithmic? What are the lower bounds that one can achieve so far? \n\n2. Given that there are two different problems (sampling vs. approximation) and two different settings (full gradient vs stochastic gradient), the authors are encouraged to be explicit on the optimality w.r.t. which problem in the introduction and be careful with overclaim in e.g. line 57. Furthermore, the rate $\\Theta(N+d^{1/3}N^{2/3}/\\epsilon^{2/3})$ for approximation never appears in the introduction, which to me is one of the core contributions of this paper.\n\n3. The theoretical analysis relies on the assumption that the initialization is not too far away from the target distribution. How does this assumption compare to that in other algorithms? Is this the exact same assumption used in other algorithms? Appendix A.6 also appears to be too brief to make the initialization condition clear.\n\n\n** Update on the reply:\n\nThanks for the detailed reply. I'm satisfied with the responses to the comments and remain my score.\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have adequately addressed the limitations and potential negative societal impact of their work.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "Nw_pdloSL5_", "reviewer": "Reviewer_gwq5", "summary": "The paper considers methods for approximating underdamped Langevin diffusion processes with sum-decomposable strongly convex potential. The main contribution in the paper is the proposal of the Accelerated ULD-MCMC algorithm and its variance-reduced variants. This is followed by a detailed analysis of convergence for each of the algorithms, and the derivation of an information-based lower bound for gradient complexity in the task of estimating ULD that matches the upper bound (in dimensionality, component number, and target accuracy). The authors conclude with a discussion of the optimality (or lack thereof) of ALUM with respect to all the variables. Comprehensive experiments are presented to compare the performance of various algorithms on estimation ULD processes, showing that ALUM and its variants achieve better performance than previous algorithms.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "** Update: The additional discussion and context provided in the rebuttal addressed my concerns and I have adjusted the score accordingly. **\n\nThe paper is interesting and well-written overall with good contribution to the algorithms for ULD processes. The authors offer a clear context about ULD and existing algorithms such as LPM and RMM (mostly in appendix, some of which I believe might better fit in the intro if space allows). I find it interesting and a bit surprising that dropping the $-\\psi_2(at)\\nabla f(X_0)$ term in RMM does not hinder the accuracy by much. The authors also offer a very comprehensive and detailed discussion on the experiments on both synthesized and real-world data.\n\nHowever, my main concern is that I’m not fully convinced that the newly proposed algorithm and its performance improvement are significant enough. The difference between ALUM and RMM might seem subtle: I did not completely get the motivation behind ALUM or its intuition (e.g., why the extra gradient term can be drop - such intuition could potentially be of value in other problems); the factor-two reduction in gradient evaluation could also use more justification regarding why it is significant (e.g., is there some applications that can critically benefit from the factor-two reduction?). For the experiment results, while I agree that the proposed algorithms outperform the existing ones (e.g., Fig 1), I would again suggest that the authors argue its significance by, e.g., relating the improvement to some application.\n\nComparison of ALUM with RMM and NAG is also interesting, but I would hope to see deeper insights, e.g., whether there may be underlying intuition for why ALUM might be desirable in its current form. I would expect this theme to re-appear in later discussion if the authors believe there may be more interesting insights.\n\nI did not go over all the proofs. Overall, the theoretical analysis is sound and I believe all the claims are correct. As a minor point I believe some of the proofs in the appendices could use more polishing with explanations around the equations (for easier reading).\n\nTo sum up, I believe the paper can be stronger if the authors could demonstrate better the significance of the improvement and/or offer some motivation or intuition behind the ALUM and RMM distinction.\n\n### Other comments:\n\nShould probably mention $I$ is identity matrix with dimension inferred from context (since $I_n$ is also used in some discussion for other things).\n\nGiven the tiny error bars (e.g., in Fig 1), I would expect the curves to be smoother/straighter. Is there certain artifact that contributes to the curviness or maybe I missed something in my thinking?\n\nApart from the main assumptions in Section 3 (Assumptions 1-3), there seem to be a number of sporadic assumptions in the theorems/lemmas, e.g. bounded 2-Wasserstein distance of initialization (Section 5, L156, and I did not fully get how Appendix A.6 addresses this), and assumption on step size $h$ in Theorem 3 (L185). I would hope that a quick comment can be added to address how strong/realistic these assumptions are.\n\nThe paper could benefit from proof-reading for typos and minor English style/grammar errors.\n\nSome typos spotted:\n\nL74: studies -> studied\n\nL132: comes -> come, is -> are\n\nL184: delete “and”\n\nL521: Eqn (20), don’t think $\\gamma=\\bar{\\gamma}$ is needed\n\nL609: Eqn (35) seems incomplete (I think a summation is missing)", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "N/A", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "18"}, "extra_text": {}}, {"review_id": "DXlQDGX7qV", "reviewer": "Reviewer_NwCb", "summary": "This paper studies sampling from strongly convex potential in the form of finite summation of N smooth components. The authors show that using ALUM to simulate ULD needs $O(N+d^{1/3}N^2/3/\\epsilon^{2/3})$ gradient evaluations.", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "This paper proposes a new algorithm ALUM and combines it with SVRG and SAGA. The algorithms and theorems are clearly written. The claims are supported by theoretical analysis and experimental results.\n\nHowever, this paper does not discuss some of the important works studying the same problem, e.g. [1,2].  The significance of this work is unclear given previous works [1,2]. [1] shows it is possible to achieve better $\\kappa$ and $N$ dependence than ALUM using much simpler algorithms based on similar ideas. The only advantage of ALUM over [1] is that ALUM can improve eps dependence from $1/\\epsilon$ to $1/\\epsilon^{2/3}$, but [2] shows that with a simple MH step, it’s possible to achieve $\\log 1/\\epsilon$ and much better $\\kappa$ and $N$ dependence at the same time. For this reason, it’s unclear why improving $1/\\epsilon$ to $1/\\epsilon^{2/3}$ is important. \n\nMoreover, the algorithm ALUM seems very tailored to sampling from potential in the form of finite summation. It’s based on the idea of RMM, a previous sampler, but its performance on sampling general log-concave distributions is worse than that of RMM. As a result, the contribution of ALUM to understanding broader sampling problems is limited as well. \n\n[1] Chatterji, Niladri, et al. \"On the theory of variance reduction for stochastic gradient Monte Carlo.\" International Conference on Machine Learning. PMLR, 2018.\n\n[2] Lee, Yin Tat, et al. \"Structured Logconcave Sampling with a Restricted Gaussian Oracle.\" arXiv preprint arXiv:2010.03106 (2020).\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Some limitations are discussed in sec 7.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "9OSFZsjh9IR", "reviewer": "Reviewer_5bz1", "summary": "In the context of sum-decomposable, smooth, and strongly convex potential functions the authors propose algorithms to improve the upper bound on the approximation error for estimating a Underdamped Langevin Diffusion (ULD) process and on the sampling error of the strongly-log-concave sampling problem. The continuous time ULD process has the strongly-log-concave distribution as its invariant distribution so approximating this process and sampling from this distribution are closely related. In the full gradient setting, the authors present a novel algorithm, ALUM which has optimal iteration complexity in $d$, $N$ and $\\epsilon$ for the ULD approximation problem. Here, $N$ is the number of terms in the sum-decomposable function, $d$ is the input dimension and $\\epsilon$ is the accuracy parameter. In comparison to the RMM algorithm (Lee et al), the dependence on $d$, $N$ and $\\epsilon$ in the iteration complexity are matched, but the number of gradient evaluations are reduced by a factor of 2.\nNext the stochastic gradient oracle setting is considered: here the authors propose algorithms (SVRG-ALUM and SAGA-ALUM) and analyze the corresponding sampling/discretization error. Finally an improvement in the lower bound for the gradient complexity of the ULD approximation problem to $\\Omega(N + d^{1/3}N^{2/3}/\\epsilon^{2/3})$ is presented. This shows that the gradient complexity of the presented Variance reduction modifications of ALUM match the lower bound in $d$, $N$ and $n$ in the sum-decomposable setting with a stochastic gradient oracle.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "Overall the authors do a good job in conveying the content of the paper. The body of the paper abstracts the lower level details of the paper fairly effectively, and the relationship of the presented algorithms with the rest of the literature is comprehensive (case in point: section 4.1).\n\nHowever, I should point out that I find the main contribution of the paper in the gradient oracle setting questionable. In each iteration, the algorithm saves a single gradient computation, but this is at the cost of worse dependence on $\\kappa$ which can often be $10^2$ or larger in practice. Thus, overall the results seem to suggest that the gradient complexity may even be worse when the error parameter $\\epsilon$ is moderately large.\n\n- Space permitting, I would recommend moving the first paragraph from appendix A.1 into the preliminaries section of the paper. I think it sets up a nice general introduction to the setting and the relation between the sampling and ULD approximation problems.\n\n- It should be made more clear that the gradient complexity upper bounds proved in the paper are not optimal in all parameters. The dependence on the condition number $\\kappa$ can be improved (as with the RMM algorithm of Lee et al).\n\n- In Section 7, the authors mention that the dependence on the condition number $O(\\kappa^2)$ can be reduced to $O(\\kappa^{3/2})$. It would help to be more explicit about the approach here and what are the challenges in improving the guarantee.\n\n- The intuition in Section 4.1, lines 118-123 can be made more clear. Namely, it would be more clear to write the bias accumulated over $T/h$ steps as $O(Th^2)$ (instead of $O(h^2)$). Further, it is mentioned that the variance term is $O(h^{3})$, so I fail to see how the bias is of smaller order than the induced standard deviation. It would also help to point out how the calculation for the variance is carried out.\n\n- According to the first two plots in the experiments section, the ALUM algorithm appears to have the same curve as the RMM algorithm but shifted by what appears to be a constant factor for all ranges of the error parameter. Seeing that the dependence on $\\kappa$ ($\\gg 1$ here) is presumably worse for ALUM compared to RMM in the proved results, does this plot seem to indicate that the analysis of with respect to $\\kappa$ can be improved for ALUM? If the precision level $\\epsilon$ is increased further, would the plots eventually diverge?\n\n- An important point to address in the paper is the case when Assumption 2 is relaxed so that only the total potential function $f$ is L-Smooth and the uniform smoothness assumption on the $f_i$'s is dropped. Is there any intuition about whether ALUM / SVRG-ALUM continue to achieve a similar iteration complexity under the gradient and stochastic gradient oracles?\n\n- This is an important point to address at least in the special case when $f$ takes the form of an empirical risk under L2 loss: the smoothness of each $f_i$ corresponds to $\\| x_i \\|_2^2$ where $x_i$ is the $i^{\\text{th}}$ input point. The current smoothness assumption implies that all the $x_i$'s lie in an L2 ball of radius $\\sqrt{L}$. Without the uniform smoothness assumption, but assuming that $f$ is $L$ smooth relaxes the condition on the points lying strictly in an L2 ball and imposes an eigenvalue constraint on the covariance matrix which is of more practical significance. Do the guarantees carry over here?\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes, I believe the authors have.", "review_text": "", "extra_scores": {"ethical_concerns": "None to the best of my knowledge.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "8yFPzR46dFp", "reviewer": "Reviewer_BMbv", "summary": "The paper proposes a new way of discretizing the underdamped Langevin dynamics, by modifying the Randomized midpoint method, proposed by Shen and Lee 2019, which reduces the number of gradient queries by a constant factor. In the case when the log potential of the target distribution is a sum of many terms, the paper produces a variance-reduced variant of their discretization method, which they show achieves the optimal query complexity for approximating the underdamped Langevin dynamics.", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "My main criticism for the paper is that the problems investigated are not well motivated. The setting is overly specific to be interesting to a practitioner, so the main contribution of the paper is a technical one. However, although the analysis is impressive, there does not seem to be genuinely new proof ideas introduced, so it seems unlikely that the paper will generate much interest among theoreticians. The lower bound that the paper proves is also not the most interesting one: it proves a lower bound for the point-wise approximation error; however, in sampling it is much more natural to consider distributional approximation error.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "86AOIG2WH2", "reviewer": "Reviewer_S5uq", "summary": "The paper makes new contributions to the literature on sampling from a strongly log-concave smooth measure. This is done by approximating the Underdamped Langevin Diffusion (ULD), and so, a fortiori, the paper also contributes to this approximation problem.\n\nThe main result is a new algorithm to approximate the ULD, which is well suited to a stochastic setting, where the gradient of the potential is replaced by certain estimators. \nIn the full gradient (i.e. without using an estimator for the gradient) setting the new algorithm matches the asymptotic performance of previous works while cutting the number of Gradient queries by about half.\nIn the stochastic setting, if the potential can be represented as a sum of $N$ smooth functions, then the authors also show how to modify their algorithm to accommodate estimators to the gradient and reduce the overall query and iteration complexity.\n\nThe results of the stochastic setting are also complemented by a lower bound which essentially show that the dependence on the parameters $N$, and $d$, the dimension are essentially optimal.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "======POST REBUTTAL======\nThe authors addressed my concerns and I have raised my score accordingly.\n\nThe new algorithm, which the authors dub ALUM (with stochastic variant VR-ALUM), is based on the RMM method for discretizing the ULD, introduced by Shen and Lee in 2019.\nThe RMM method makes two gradient queries per iteration, and the gradient cleverly identified a way to remove one query so that the overall bias does not increase too much.\nHowever, this does come at the cost of worsening the dependence on the condition number, which, I believe, was the main motivation in introducing the RMM algorithm.\n\nThe stochastic gradient analog introduced by the authors is a nice addition to the literature and, to my knowledge, outperforms previous algorithms. It is based on thoroughly analyzing variance reduction techniques to stochastic gradients coupled with the discretization scheme of the ULD. The supplied lower bounds further strengthen this.\n\nWhile the authors supply a clear motivation and an explanation of their analysis, I tried to go over the proofs. I found them extremely technical, even when compared to previous works. This is the reason why my current score is not higher, and I am hoping that the authors can comment on that and explain the main technical difficulties and their own contributions in this sense so that I can raise my score.\n\nOther than that, the paper has many typos and English grammar mistakes. I would ask the authors to go over the entire paper and fix such issues carefully. Such a thing could greatly improve the readability of the paper. (I want to emphasize that my score is not based on that point).\n\nOther comments (I do not mention here typos):\nTable 1: Actually, the dependency on $d$ is not shown below ($d$ is not even properly defined by this point)\nLine 36: I agree that the question is natural, but perhaps there is some room for discussion about the notion of optimality here.\nOptimality can be taken with respect to the number of iterations, the total number of queries, and total operations.\nAlso, which parameters are taken into account (this discussion appears later in the paper, but some small mention could be appropriate here.\nLine 63: I could not understand which divergences are discussed. Total variation vs. something else?\nLine 95: I'm not sure what the $h$ in the subscript $X_h$ stands for.\nLine 126: I'm not sure what the comparison to Nesterov's momentum algorithm is supposed to demonstrate. Yes, the algorithms are similar because both of them add a momentum term to accelerate the vanilla gradient method. This similarity was already observed before (for example [CCBG2018]).\nLine 211: I would suggest reconsidering the use of $U$ for a potential function and $\\mathcal{F}$ for a function class. An uncareful reader can think that $\\mathcal{F}$ stands for a sigma-algebra and $U$ is an open set.\nLine 239: Can the authors say something about the requirement that n is a multiple of N. At first glance, it looks very artificial.\nLine 727: I could not find a Lemma 7 in the paper of Dalalyan and Riou-Durand, in any of its versions.\nLine 371: On the same topic, there are only two authors to this paper, no et al.\n\n[CCBG2018] Cheng, X., Chatterji, N. S., Bartlett, P. L., & Jordan, M. I. (2018, July). Underdamped Langevin MCMC: A non-asymptotic analysis. In Conference on learning theory (pp. 300-323). PMLR.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors adequately addressed the limitations and potential negative societal impact of their work.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "8"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes a randomized algorithm for approximating a trajectory of underdamped Langevin dynamics. It also uses variance reduction to improve the error over decomposable problems. It also establishes an information-theoretic lower bound. There have been extensive discussions between the authors and the reviewers, and I have gone through them carefully. I agree with the majority of the reviewers that the paper has made interesting technical contributions to the field. I am happy to recommend acceptance. 


Decision: Accept (Poster)