Langevin Quasi-Monte Carlo
Authors: Sifan Liu
Venue: NeurIPS 2023 poster (2023)
Keywords: Completely uniformly distributed; log-concave sampling; low-discrepancy; MCMC;
TL;DR: This paper proposes a method for applying quasi-random numbers to Langevin Monte Carlo algorithms.

Abstract: Langevin Monte Carlo (LMC) and its stochastic gradient versions are powerful algorithms for sampling from complex high-dimensional distributions. To sample from a distribution with density $\pi(\theta)\propto \exp(-U(\theta)) $, LMC iteratively generates the next sample by taking a step in the gradient direction $\nabla U$ with added Gaussian perturbations. Expectations w.r.t. the target distribution $\pi$ are estimated by averaging over LMC samples. In ordinary Monte Carlo, it is well known that the estimation error can be substantially reduced by replacing independent random samples by quasi-random samples like low-discrepancy sequences. In this work, we show that the estimation error of LMC can also be reduced by using quasi-random samples. Specifically, we propose to use completely uniformly distributed (CUD) sequences with certain low-discrepancy property to generate the Gaussian perturbations. Under smoothness and convexity conditions, we prove that LMC with a low-discrepancy CUD sequence achieves smaller error than standard LMC. The theoretical analysis is supported by compelling numerical experiments, which demonstrate the effectiveness of our approach.

Langevin Quasi-Monte Carlo
Sifan Liu Department of Statistics Stanford University Stanford, CA 94305 sfliu@stanford. edu
Abstract
Langevin Monte Carlo (LMC) and its stochastic gradient versions are powerful algorithms for sampling from complex high-dimensional distributions. To sample from a distribution with density π (f) x exp(-U(0)), LMC iteratively generates the next sample by taking step in the gradient direction VU with added Gaus- sian perturbations. Expectations w.r.t. the target distribution ir are estimated by averaging over LMC samples. In ordinary Monte Carlo, it is well known that the estimation error can be substantially reduced by replacing independent random samples by quasi-random samples like low-discrepancy sequences. In this work, we show that the estimation error of LMC can also be reduced by using quasi- random samples. Specifically, we propose to completely uniformly distributed (CUD) sequences with certain dis-eiscyyyancy property to generate the Gaussian perturbations. Under smoothness and convexity conditions, we prove that LMC with a low-discrepancy CUD sequence achieves smaller error than standard LMC. The theoretical analysis is supported by compelling numerical experiments, which demonstrate the effectiveness of our approach.
Introduction
Sampling from probability distributions crucial task in both statistics and machine learning. However, when the target distribution does not permit exact sampling, researchers often rely on Markov chain Monte Carlo (MCMC) methods. These techniques simulate a Markov chain that converges to the target distribution as its stationary distribution. Recently, MCMC samplers based on discretizing the continuous-timg Langevin diffusion have become popular, due its ease of implementation and ability to handle stochastic gradients (Welling and Teh 2011).
The primary focus of this work is the quality of samples generated by Langevin Monte Carlo (LMC) algorithms in terms of estimating the expectation Ex-z [f)) for some integrand f by sample averages. In the context of Bayesian inference, the target distribution is typically the posterior distribution, and computing the posterior expectation, posterior variance, or confidence intervals are of great interest. the context of post-selection inference, the target distribution и is the probability distribution conditioned on the selection event, and computing the selection -adjusted p-value is the main task. LMC has been widely used in this problem as well (Markovic and Taylor 2016 Shi et al. 2022). In all these situations, the accuracy of the sample average estimator is critical and affects the downstream data analysis.
In traditional Monte Carlo sampling, it is well known that using quasi-Monte Carlo (QMC) samples, instead of independent and identically distributed (i.i.d) random samples, lead significant error reduction. So it is natural to ask whether we can apply QMC techniques to improve Langevin Monte Carlo sampling as well. In this work, we introduce the Langevin quasi-Monte Carlo (LQMC) algorithm, which replaces the i.i.d. random inputs in the LMC algorithm with quasi-random numbers.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

251 points from MT1 19937
251 points from LCG
Figure 1: Scatter plots of 251 points generated from Mersenne Twister 19937 (left) and 251 points generated from linear congruential generator (LCG) of period 251, Points from entire period of pseudo-randon random number generator (right) fill the unit square more evenly than the same number of points from a PRNG with a larger period (left).
These quasi-randon numbers are carefully designed to sample from the target distribution more evenly and more balanced, leading to improved estimation accuracy.
Not all quasi-Monte Carlo point sets suitable for simulating Markov chains. Suppose the Markov chain is driven by sequence of uniform random vectors in the unit cube. A sufficient condition for the sequence is known as completely uniformly distributed (CUD). In our implementation of the driving sequence, we use an entire period of pseudo-random number generator (PRNG). While modern computer simulations often use PRNGs with a large period, such as Mersenne Twister with a period of 2 199377 - 1, our approach runs through the entire period of a PRNG with a relatively small period in the LMC algorithm. The advantage of using an entire period of PRNG is that the points are more evenly distributed, which is more desirable for numerical integration. We illustrate the balancing property of entire PRNG Figure
The main contributions of this paper are threefold. First, we propose a novel technique of using quasi-random numbers in Langevin-type algorithms, which can be applied wide range of such algorithms by substituting i.i.d. random numbers with a sequence of quasi-random numbers. The quasi-random numbers are constructed similarly as usual PRNGs, therefore no extra computational complexity is required. Second, we evaluate the performance of the proposed LQMC algorithm in a variety of numerical experiments, demonstrating that it can significantly reduce the mean squared error (MSE) of traditional LMC by a factor ranging from 2 to 500, depending on the problem. Finally, we provide theoretical analysis showing that LQMC can reduce the Monte Carlo part of the error from O(n-1/2 1/2) On 1+6) for any 8 in situations where the Markov chain strongly contracting and the integrand function ff is sufficiently regular. This eтror reduction is consistent with the usual improvement achieved by using quasi-Monte Carlo place of plain Monte Carlo.
The rest of the paper is organized as follows. Section 2 we provide some background on LMC and QMC, followed by review of related work. Section 3 describes the LQMC algorithm and its implementation details. In Section 4 we present theoretical guarantees for the proposed method. Finally, in Section we provide empirical results evaluate the performance of LQMC and compare with the standard LMC algorithm.
2 Backgrounds
This section provides some background on Langevin Monte Carlo and quasi-Monte Carlo.
2.1 Langevin Monte Carlo
Suppose want sample from the target distribution IT (0) exp(-U(0)) where 0 E Ed and U is known as the potential function. LMC algorithms are based on Euler-Maruyama discretization the Langevin diffusion A(t), which satisfies the stochastic differential equation
d0(t) = VU (U((())) + V2aW,
(1)

---

where (W2) is a d-dimensional standard Brownian motion. Under mild technical conditions, the Langevin diffusion (t) has its unique invariant distribution (Roberts and Tweedie 1996). With a discretization step size h, LMC updates the sample 0g by
0x+1 t 02 (VV( + 22h
(2)
where & (0, Ta).
In many applications, interested in computing the expectation f == Bg-- [f(0)] over for some z--intgrrbble function f. The LMC estimator of μ is the sample average
An
f(0x).
k=1
where n is the number of iterations.
Teh et al. 2016) provide an asymptotic bias-variance decomposition of the MSE of the weighted
average 2 (o.) and show that the optimal step size scales as 02g > k-133 leading MSE 2 he
of order O( 2/3). (n333 Here h2 is the step size used at the k-th iteration. Vollmer al. (0066) generalize this result to the non-asymptotic setting with constant step size h. They show that the MSE is of order Ohh² + n, where hi corresponds to the squared bias and x2% corresponds the variance.
2.2 Quasi-Monte Carlo
QMC an alternative to Monte Carlo for numerical integration and well-known for having much higher accuracy than Monte Carlo. QMC is primarily designed to numerically evaluate the integral = Spo a fu))d.. It estimates y by taking points u, E [0, 1] and let the estimator be
P Σ f(u)).
Unlike Monte Carlo which takes u, be identically independently distributed (i.i.d.), QMC constructs the point set { ) that aims to minimize the star discrepancy
d
1 1u E (0, - II aj
D = D (u un) sup ae[0,1]
1
The star discrepancy measures the uniformity of the point sets by comparing the fraction of points inside [0, a) and the volume I[ Qj. taking supreme over all the rectangles inside [0, 1] an- chored QMC can generate points with D; O(n Onn (log nid-1), thus QMC is also known as low-discrepancy sequence. Commonly used QMC points include Sobol' sequence (Sobol" 1967). Niederreiter's sequence (Niederreiter 1987). Halton's 's sequence, and lattice rules. For compre- hensive survey, we refer to the monograph Dick and Pillichshammer (2010). If the integrand f has bounded variation in the sense of Hardy and Krause [ then the Koksma-Hlawe inequality (see e.g. Dick and Pillichshammer 2010)) bounds the integration error by
A 11 < D || S O(n Onn-I(logn) (Iogn)d-1).
4))
While the Koksma-Hllwwke inequality shows that QMC is asymptotically better than usual Monte Carlo, docsn't provide a practical way to estimate the error. Moreover, integrands might have infinite Hardy-Krause variation.
One can apply randomization techniques to QMC address both problems. Common randomization techniques include random shifts (Cranley and Patterson, 1976) and scrambling (Owen 1995). For RQMC samples uI...., un each u; (n, 1 individually but they have the low-discrepancy property collectively with probability 1. One can estimate the error by multiple independent random replicates. For sufficiently smooth f. the scrambled Sobol' sequence has variance O(n - =(og n d-1 (Owen b).

---

2.3
Related work
The first attempt to apply quasi-random numbers to simulate stochastic differential equations was made by Hofmann and Mathé (1997). They showed that if a numerical scheme is weakly convergent with i.idd samples, then using completely uniformly distributed (CUD) sequences also leads to consistent estimation. They also demonstrated that certain owwdiicrepaancy sequences are not suitable for simulating SDEs. There have also been some efforts to apply QMC to MCMC. (Ownn and Tribble (2005) proposed to apply CUD sequences to Metropolis algorithm and showed that the method is consistent in problems with finite state spaces. Chen et al. (20TT) generalized the consistency result to continuous state spaces under the assumption that the Markov chain is contraction. More recently Dick et al. (0116); Dick and Rudolf (2014) proved that there exists constructions of the driving sequence [u J 121 such that the discrepancy between the empirical distribution of MCMC samples and the target distribution is bounded by O(n 1/2 (log
the same rate achieved by random inputs. Another line of applying QMC to Markov chains is known as array- urray-QMC proposed by Eeyyyer et al. (2008). Array-RQMC runs in parallel multiple Markov chains, and each iteration involves a complicated reordering of the states so that the low-discrepancy among the chains is maintained. Empirically, it achieves significantly smaller estimation error than usual MCMC, but theoretical guarantees remain a challenging open problem.
There has been growing interest in using QMC techniques in various machine learning tasks, such as variational inference (Buchholz et al. 2018 Liu and Owen 2021), policy learning and evaluation (Amold et al., 2022). reinforcement learning with evolution strategies (Choromanski et al. 2019] Rowland et al. 20, compression of large datasets (Dick and Feischl| [201]), example selection in
stochastic gradient descent (SGD) (Lu et al. 202I). and deep learning for solving partial differential equations (Longo et al., 2021).
Numerous efforts have been devoted improving LMC and stochastic gradient Langevin dynamics (SGLD). To overcome the instability of Euler-Maruyama discretization, various numerical schemes have been proposed, including higher-order integrators (Chen et al. 2015). underdamped LMC (Cheng et al. 2018), and stochastic Runge-Kutta diffusion (Li el al. 2019p. For SGLD, variance reduction techniques such as SAGA and SVGR (Dubey et al., 2016) and control variates (Baker et al. 2019) have been proposed. LMC also provides a useful perspective for optimization, as demonstrated by the analyses in Chen et al.] (2016); Dalalyan (2017); |Raginsky et al.] (2017); x [uu et al. (2018); EEdogdu et al. (2018). Our contribution is orthogonal to all the aforementioned work, as our algorithm only modifies the random numbers used in the algorithm. Therefore, our method can be combined with other algorithms without interference.
3 QMC for LMC
In the LMC algorithm, we can think of the Markov chain as being driven by a sequence of uniform variables uy in the unit cube [0, 1jd. For instance, the Gaussian perturbation can be represented as & = ib --1(ux). where 11 denotes the inverse Gaussian CDF applied element-wise to ug. If a stochastic gradient is employed, the randomness associated with the stochastic gradient can also be expressed as uniform variables. Therefore, we can write the transition of the Markov chain as 0x+1 = vi(A, ug+1)) ). In typical computer experiments, Ug are not really i.i.d. but are deterministic pseudo-random numbers. In this section, we will describe an alternative method of generating the pseudo-random numbers ug, which are carefully constructed and can lead to more accurate sample averages.
The idea here is to point sets that are more evenly distributed such QMC points, which can lead to significant improvement in the usual Monte Carlo estimation. However, caution is required when using QMC points to simulate an SDE like (.. This is because the correlation between successive QMC samples may introduce undesired behavior in the Markov chain, as demonstrated in (Tribble) 2007, Section 3.2). To avoid the dependence among successive values, we require that the blocks of points (v,, U,+1, " 1) for any lag d are uniformly distributed. This notion of uniformity is formally known as completely uniformly distributed (CUD, Korobov (1948)) which we define next.
vi+d-1)
We say an infinite sequence { C [0, 1]d is uniformly distributed on [o. 1]d if the star discrepancy D* ({ u, )())) goes to 0 - oc, where the star discrepancy is defined in Equation 3

---

Definition 3.1 (Completely uniformly distributed sequence (CUD)). An infinite sequence {us)i=o C [0.1] is called completely uniformly distributed, if for all positive integer the sequence {(Ubs.... Uk+d-1) S R= is uniformly distributed on [0, 1,11 trian- gular array Vn = (Un,1....., Vn,N.) is called array-CUD, if for all positive integer d, ((Inn.)..,;Pddd))
Un,d....,, (Vn,N,-d+1..... Un.N..)
oc, N₂, - oc.
In other words, the subsequent d-tuples in CUD sequence are uniformly distributed in the d- dimensional unit cube for any positive dimension d. Now we are ready to present the main algorithm.
3.1 LQMC algorithm
Let [0 be a CUD sequence. Let U1 = (Ukds.... ((8t4--d.,( E Rdd be the k-th non-overlapping d-tuple from the sequence 2 0). A CUD sequence often constructed deterministically. They can further be randomized using the Cranley-Patterson (i.c. random shift) rotation (Cranley and Patterson 1976)
+ up A mod 1,
where A ~ [n, 11 The Cranley-Patterson rotation randomly shifts each dimension us by a uniform random number separately. Then each ug is uniformly distributed on [o, 1]d we apply the inverse Gaussian CDF to each coordinate of then (ux) ~ algorithms, we will let & = $ 1-11(mu (ug) and use & as the Gaussian perturbation in the k-th iteration. Specifically, each iteration takes the form
NO, Ia). In the Langevin-type
0 k++1
0 (O) + 2h 4
- (ug+1),
2
Thus the transition map is (0, u) 0 (0) + 2h In practice, we only run finite many iterations. the following, we will describe how to construct finite CUD sequence and feed it into the LMC algorithm.
3.2 Construction of CUD sequences
finite CUD (array-CUD) sequence is often implemented by using an entire period pseudo random number generator with a small period (Tribble) 2007). There exist other constructions of CUD sequences. For further details, interested readers can refer to Levin (1999). We propose to use the finearfeeddbad shift register (LFSR) provided in Chen (2011). because has demonstrated good performance and the computational effort required comparable other commonly used PRNGs.
The binary Galois LFSR (Tausworthe generator, Tausworthe (1655) ) of order m updates the states
E {0, 1} recursively by
m-1
b, Σ
a a-bi-m-j mod 2,
m
with initial states by: bj...., bm- pre-specified. The m-tuple (b;, b. b E GF(2)" can only take 2 different values. If there is an m-tuple that is all zero, then all b, 's's this sequence must be zero. So the period of the sequence {b; h:20 is at most n = 2m 1. Moreover, the period is exactly equal 2 1 if and only if the characteristic polynomial
= + am 0m-1
+ a1 + ag
is a primitive polynomial over GP(2)
Niederreiter 1992 Lemma 9.1). Given the states h2 and 1, U/ is computed with
an offset > 0 such that ged s, 2
m
04
2
1......,,9 2m 2.
1=0
That is, for each i, we take the m-tuple (bai+) Josjsm and interpret it as the binary expansion of a For the next step, we jump s bits ahead in the sequence H2o and use the m-tuple starting from ba(t+1)- Chen (201)) provided table of the LFSR generators for 10 & < 32. They searched the offsets so that the LFSR has good equi-distribie properties. Our experiments use the LFSR generators listed there.
5

---

Given the sequence {= 1n-1 matrix
of length n, we repeat it times and arrange os 's in the following n. x d'
60
01
vd-1
ud+1
2
(5)
P(n (n--))d+1
Und-1
We run the LMC algorithm = 2 iterations. The k-th uniform vector in is the k-th row of the above matrix. The procedure summarized in Algorithm
Algorithm Langevin quasi-Monte Carlo (LQMC)
Input: Number iterations n 2 such that g2d 1, d) step size initial value A Generate LFSR sequence {t4 h200 of period 2 1.
Let 0g ()d Ukd-1) E [0, 1]6, for 1 < k: S
Apply Cranley-Patterson rotation (random shift) to ug's.
for l <- 1, do
0g < 0k-1 - AVU(O=-1) 22h0- (ux) end for Output: 0,, 0,,...,0w
If god,n, d) = 1, then each column of the matrix 5 contains no repeated values. This means that among the n = 2 iterations of the LQMC algorithm, each dimension uses one value in each sub-interval (2m k+11 at most once (0 S k S 2m 1). This perfect di-eessiiooal stratification is one of the reasons why CUD may achieve smaller estimation error than pseudo-random numbers. ged(n, d) > 1, then we take d' to be the smallest integer greater than d and co-prime with n. We then create the matrix in (5) similarly but with d' columns. In the LQMC algorithm, we take ug to be the k-thh row of the matrix but only use the first coordinates.
Algorithm| may seem to be restricted by having fixed number of iterations, n = 2 1. However, in practice, the LQMC algorithm can be started with an initial value of m. If the chain does not converge after 2 iterations, one can continue the chain with another freshly generated LFSR, possibly with larger period. This allows for flexibility in adjusting the number of iterations based on the convergence of the chain. Additionally, if burn-in period is required, one can first run the algorithm with an LFSR of small period serve as the burn-in stage and then continue with larger LFSR. Furthermore, running multiple chains with independent random shifts is embarrassingly parallel. We present the algorithm the form of the basic LMC algorithm with accurate gradient and constant learning rate. However, as we noted previously, other Langevin-type algorithms can also utilize the CUD sequence directly by substituting the pseudo-randos numbers with the LFSR sequence.
4 Theoretical guarantee
Here study the estimation error In - that is I-Lipschitz and bounded. As the first attempt to prove the convergence rate of using QMC in LMC, we impose the relatively strong conditions of smoothness and convexity.
f(2)) +(f) of LQMC for some test function
Assumption 1. The potential function U is L-smooth
(O) (U(0))la s L110 0211,. V0,0',
and M strongly convex
M 2 10'
U(0') 2
0
V 60.0'.
We will also assume a constant step size h. While LMC with vanishing step sizes converges weakly to the target distribution, practice constant step size is often used (Vollmer al. 2016 Brosse et al. 2018). With constant step size, we can derive non-asymptotic error bound for LQMC.

---

Assumption implies that if the step size h S L-M then the transition map vb is a strong contraction
2
with parameter = 11 1 h.M, i.c. [1v(0,u) u) - (00", u))l2 = 110 110 - 1 ((VU(0) vu S p|10 0
(6) See e.g. Lemma Dalalyan and Karagulyan 20199) The strong contraction implies that if we start two chains from 9 and 8' and use the same random numbers at every step, then the two chains will merge exponentially fast. other words, the state 0% largely depends on the most recent iterations
and quickly forgets about the past history. Formally, let w (t) = (Ug.... uk-e+1) denote the random
numbers used in the most recent steps. Define the g-step transition as
0 v W
and let ft(w(( (0)) denote the value of f(o)) marginalized over Ox- I,,
(0) / f ve(z, w ))(x))))
Thus fe(w), (0) only depends on the most recent / iterations. Due to the strong contraction, [w]] (0) f6x) decays exponentially fast with l. So for large E, the estimation error of n- Σズ f0x)) is close to the error of n Σx=e+1 fe(w(") ).). The latter can be viewed défdimensional numerical integration scheme based on the point set (w( ) By leveraging the discrepancy bound of the LFSR sequence and assuming that fe has bounded variation in the sense of Hardy and Krause, we can derive an error bound using the Koksma-Hlawka incquality (4). Now we state the main error bound and leave the detailed proof in the Appendix A
Theorem 4.1 (Error bound of LQMC). Let Assumption | hold. Define the step size h S L+ = 1-hh, f ((1/2) log, h1). Let 01.... O,, be the output of Algorithm [ which runs iterations with step size S L+M Assume the LFSR sequence {v.3i2o in use has period n 2m offset and gcd(m,, n)) = ged(de, n)) = 1. If fe has bounded variation in the sense of Hardy and Krause, then as + 00 we have
f(0x) x(n s Cin-1+ + C2h22/2 A5>0. 0.
k=1
Here 8 hides voly-logarithmic factors ()ogn)4, C₁ depends on d, E and [ele and C2 32/2 maxogken ||0x|| + E,, [[]@[]]
The upper bound consists of terms. The first term represents the numerical integration error, which arises from the discrepancy of the point set used in the integration scheme. By utilizing low-discrepancy CUD sequences, we reduce this numerical integration error (the first term) from the standard rate of (1/2) note that when using a constant step size h in LMC, the bias term (second term) does not vanish. This bias term includes not only the discretization error of the Langevin diffusion, but also the difference
faster rate of -1+6))) for any 8 > 0. However, it is important to
(8) Consequently, the bias term in our analysis is between f0x) and its truncated version fe (w) k larger than the bias term in Vollmer al. (2016), which employs different techniques and assumptions based on the Poisson equation.
The theorem's 's assumption of finite Hardy- Krause variation is common requirement in error bounds for QMC methods, and it can be challenging to verify in practice. |Basu and Owen (2016) provide sufficient conditions in order for f o 1/ to have finite HK variation, requiring the Ê-step transition 4 to sufficiently smooth. In the next section, we aim to assess the practical performance of the proposed LQMC algorithm through numerical experiments.
5 Numerical experiments
To comprehensively evaluate the performance of the algorithm, we will consider both convex and non-convex potentials, both low-dimensional and high-dimensional state spaces, both accurate and stochastic gradients, both smooth and discontinuous integrands, as well as different learning rate schedules. Additional numerical results can be found in the Appendix B
7

---

[x)=
1850
%
21
MC
CUD
Axl X
fix) x?
100 2,23
NSE
2
MC
CUD
Figure 2: Bayesian logistic regression with accurate gradients (top) and stochastic gradients (bottom).
5.1 Bayesian logistic regression
We first consider the Bayesian logistic model
y/ x4 ~ Bernoulli( + cxx((xx[3))1) 1 12iSN,
8 ~ Ia).
We take N 20, d 10. The features 2, are generated from (0, Σ) with Eig 2-1i- The coefficients / and the data yi 's are generated from the same model. We consider the test functions The step size h is fixed to 0.001.
f=) = , .3. (2y20) for j = 1.....d d.
We compute the MSE of the estimator based on usual LMC and the proposed LQMC with CUD sequences and report the MSE averaged over all coordinates and 20 random replicates. We do not have a closed form for the expectations EE [/]. so the ground truth is estimated using a high-accuracy estimator proposed in He al. (2033) using scrambled Sobol' sequence with very large sample size.
In Figure (top pancl), present log-log plot of the MSE against the number of iterations. Across all three test functions, we observe that LQMC reduces the MSE by a factor ranging from 4 to As the number of iterations increases, the curve corresponding LQMC reaches plateau. This behavior can be attributed to the discretization error inherent in the unadjusted LMC, which cannot be further reduced by increasing the number of iterations.
In the bottom panel of Figure 2 we increase the number of observations to = 100 and incorporate stochastic gradient estimation in the Langevin algorithm. Specifically, at each iteration, we estimate the gradient using a random subset of 10 observations. The results demonstrate that LQMC still provides big improvement when is smaller than 244, However, that the LQMC curve flattens again. It is worth noting that the improvement achieved by LQMC in this scenario is less pronounced compared the previous example, primarily due to the presence of noise in the gradient estimates.
surpasses 244, observe
5.2 Bayesian linear regression
Now we try a higher dighes--imeesiional example with Bayesian linear regression. The model is defined as 1/1 Nxx =4 IsisM N,
8 ~ I).
We take d = 100 and N = 20. We generate x; Rdd similarly in the logistic regression example. The test functions and step size are also unchanged. The posterior distribution of 3 has the closed form A (( Xx + -1 at 100 dimension, LQMÖ still brings substantial improvement over LMC in terms of MSE. In
x x XTX + T)-1). The results are shown Figure 3 We see that even

---

[X0 x
1p33
MSE
2 gar 2th 24
2
MC
CUD
Figure 3: Bayesian linear regression in 100 dimensions.
= 0.0001
h= 0.01
decreasing /
MSE
%
2 2
245
22a
Figure 4: Crossed random effect.
particular, for the integrand f(x) = zy. LQMC achieves a reduction in MSE of approximately 500-fold compared to LMC.
5.3 A hierarchical Bayesian model
We consider a hierarchical Bayesian model known as the crossed random effect model
Yi N(u ai + 1), Isist, 1 ISiSJ. 1 ~ N(O, 1), a₁ log(o?), log(o?) and 1)..
aid NM o.), by N(O,of).
The goal is to sample from the posterior distribution of (p, a, b, log(o)), log(o?)), log(o))))) which has di- We will consider the test functions f(x) Zz
mension I+ J We take 3, [} S d). The ground truth of E [z)] is estimated by Langevin dynamics with Metropolis adjustments (MALA) using large sample size.
We will compare the performance the IQMC algorithm using three different step sizes: a constant step size 10-4 4, a constant step size of 10-2, and decreasing step sizes with h4 = co(c; + k) -1/3 The choice of co and с₁ ensures that the step size decreases from 10 to 10-4 ttrooughuut the entire algorithm. The use the exponent 11/3 in the decreasing step sizes is recommended in Teh et al. (2016). The results of these comparisons
presented in Figure
In the small step size case (left panel), we observe that the errors LMC and LQMC are initially comparable for small values of n. This is because the algorithm converges slowly, and thus the error dominated by the bias. However, as n increases, the improvement of LQMC becomes evident. In the large step size case (middle panel), the MSE of LQMC is consistently smaller than that of LMC even for small values of n. This is because the algorithm converges faster to the target distribution with larger step size h. Therefore, the improvement of LQMC is more pronounced. Interestingly, in this particular example, using decreasing step sizes yields similar accuracy to using a constant step size of 10-4, It is worth noting that the MSE of LMC does not decrease at rate of n. -2/3 as Teh et al. 2016). This is because the line in the plot does not represent the accuracy against the iteration within a single training process. Instead, reflects the accuracy achieved after completing all iterations of the algorithm, considering different values of n.

---

Rx)
f00 x²
2
2222
2 22
LMC
LGMC
Figure 5: Double well potential.
5.4 Nonconvex potential
Finally we investigate a double-well potential function (xx) = 7x²² log(1 + x²) from Pagès and Panloup (2018). We know [x] and E [1 (x>00 0.5. The second moment x= is
computed by Gaussian quadrature. See the results in Figure Since the potential has separate local minimums, takes longer for the Langevin algorithm to explore the space sufficiently and converge the target distribution. Once converged, the improvement of LQMC over LMC is still significant.
Acknowledgments and Disclosure of Funding
The author thanks Prof. Art Owen for helpful conversations. This work was partially funded by the NSF grant DMS-2152780 and the Stanford Data Science Scholars program.
References
Arnold, S. M., L'Ecuyer, P., Chen, L.., Chen, Y.-F,, and Sha, F. (2022). Policy learning and evaluation with randomized quasi-Monte Carlo. International Conference on Artificial Intelligence and Statistics, pages 1041-1061. PMLR.
Baker, J., Fearnhead, P., Fox, E. B., and Nemeth, C. (2019). Control variates for stochastic gradient MCMC. Statistics and Computing, 29:599-615
Basu, K. and Owen, A. B. (2016). Transformations and Hardy-Krause variation. SIAM Journal on Numerical Analysis, (946):144619966
Brosse, N., Durmus, A., and Moulines, E. (2018). The promises and pitfalls of stochastic gradient Langevin dynamics. Advances in Neural Information Processing Systems, 31.
Buchholz, A., Wenzel, F., and Mandt, S. (2018). Quasi-Monte Carlo variational inference. In International Conference on Machine Learning, pages 668-677. PMLR.
Chen, C., Carlson, D., Gan, Z., Li, C., and Carin, L. (2016). Bridging the gap between stochastic gradient MCMC and stochastic optimization. In Artificial Intelligence and Statistics, pages 1051-1060. PMLR.
Chen, C., Ding, N., and Carin, L. (2015). On the convergence of stochastic gradient MCMC algorithms with high-order integrators. Advances in neural information processing systems, 28.
Chen, S. (2011). Consistency and convergence rate of Markov chain quasi Monte Carlo with examples. PhD thesis, Stanford University.
Chen, S., Dick. J., and Owen, A. (2011). Consistency of Markov chain quasi-Monte Carlo on continuous state spaces. The Annals of Statistics, 39(2):673-701.
Cheng, X., Chatterji, N. S., Bartlett, P. L.., and Jordan, M. (2018). Underdamped Langevin MCMC:
A non-asymptotie analysis. In Conference learning theory, pages 300-323, PMLR.
10

---

Choromanski, K. Pacchiano, A., Holder, J.. and Tang. Y. (2019). Structured Monte Carlo sampling for nonisotropic distributions via determinantal point processes. arXiv preprint arXiv:/905. 12667.
Cranley, R. and Patterson, T. N. (1976). Randomization of number theoretic methods for multiple integration. SIAM Journal on Numerical Analysis, 13(6):904-914.
Dalalyan, A. (2017). Further and stronger analogy between sampling and optimization: Langevin Monte Carlo and gradient descent. Conference on Learning Theory, pages 678-689. PMLR.
Dalalyan, A. S. and Karagulyan, A. (2019). User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient. Stochastic Processes and their Applications, 129(12):527--5111
Dick, J. and Feischl, M. (2021). A quasi-monte carlo data compression algorithm for machine learning. Journal of Complexity, 67:101587.
Dick, and Pillichshammer, (2010). Digital Sequences, Discrepancy and Quasi-Monte Carlo Integration. Cambridge University Press, Cambridge.
Dick, J. and Rudolf, D. (2014). Discrepancy estimates for variance bounding Markov chain quasi- Monte Carlo. Electron. Probab, 19(105):1-24.
Dick, J., Rudolf, D., and Zhu, H. (2016). Discrepancy bounds for uniformly ergodic Markov chain quasi-Monte Carlo. Annals of Applied Probability, 26(5):317-3322
Dubey, K. A., J Reddi, S., Williamson, S. A., Poczos, B., Smola, A. J., and Xing, E. P. (2016). Variance reduction in stochastic gradient Langevin dynamics. Advances in neural information processing systems, 29.
Erdogdu, M. A., Mackey, L., and Shamir, O. (2018). Global non-convex optimization with discretized diffusions. Advances in Neural Information Processing Systems, 31.
He, Z., Zheng, Z., and Wang, X. (2023). On the error rate of importance sampling with randomized quasi-Monte Carlo. SIAM Journal on Numerical Analysis, 61(2)5515-338
Hofmann, N. and Mathé, P. (1997). On quasi-Monte Carlo simulation of stochastic differential equations. Mathematics of computation, 66(218):573-589,
Korobov, N. (1948). On functions with uniformly distributed fractional parts. In Dokl. Akad. Nauk SSSR, volume 62, pages 21-22.
L'Ecuyer, P., Lécot, C., and Tuffin, B. (2008). A randomized guasi-Monte Carlo simulation method for Markov chains. Operations Research, 56(4):958-97)
Levin, M. B. (1999). Discrepancy estimates of completely uniformly distributed and pseudorandom number sequences. International Mathematics Research Notices, 1999-22)
Li, X., Wu, Y., Mackey, L., and Erdogdu, M. A. (2019). Stochastic Runge- Kutta accelerates Langevin Monte Carlo and beyond. Advances in neural information processing systems, 32.
Liu, and Owen, (2021). Quasi-Monte Carlo quasi- meaustneetoonn in variational Bayes. The Journal of Machine Learning Research, 1043::10033110
Longo, M., Mishra, S., Rusch, T. K., and Schwab, C. (2021). Higher-order muasi-monte carlo training of deep neural networks. SIAM Journal on Scientific Computing, 43(6):
Lu, Y., Meng, S. Y., and De Sa, C. (2021). A general analysis of example-selection for stochastic gradient descent. In International Conference on Learning Representations.
Markovic, J. and Taylor, J. (2016). Bootstrap inference after using multiple queries for model selection. arXiv preprint arXiv: 16/2.0781.
Niederreiter, H. (1987). Point sets and sequences with small discrepancy. Monatshefte für Mathematik, 104:273-337.

---

Niederreiter, H. (1992). Random Number Generation and Quasi-Mont Carlo Methods. SIAM.
Owen, A. B. (1995). Randomly permuted (t, s)-nets and (t. ssequenccs.. In Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing, pages 299-317, New York. Springer- Verlag.
Owen, A. B. (1997a). Monte Carlo variance of scrambled net quadrature. SIAM Journal of Numerical Analysis, 34(5):1884-19)) 1884-1910.
Owen, A. B. (1997b). Scrambled net variance for integrals of smooth functions. Annals of Statistics, 25(4):15411162.
Owen, A. B. and Tribble. S. D. (2005). A quasi-Monte Carlo Metropolis algorithm. Proceedings of the National Academy of Sciences, 102(25):8844 8849.
Pagès, G. and Panloup, F. (2018), Weighted multilevel langevin simulation invariant measures. Annals of Applied Probability, 28(6):3358-3417
Raginsky, M., Rakhlin, A., and Telgarsky, M. (2017). Non-convex learning via stochastic gradient
Langevin dynamics: a nonasymptotic analysis. In Conference on Learning Theory, pages 1674- 1703. PMLR.
Roberts, G. O. and Tweedie, R. L.. (1996). Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, pages 341-363.
Rowland, M., Choromanski. K. M., Chalus, F., Pacchiano, A., Sarlos, T., Turner, R. E., and Weller, A. (2018). Geometrically coupled monte carlo sampling. Advances in Neural Information Processing Systems, 31.
Shi, J., Liu, C., and Mackey, L. (2022). Sampling with mirrored stein operators. In International Conference on Learning Representations
Sobol', I. M. (1967). On the distribution of points a cube and the approximate evaluation of integrals. Zhurnal Vychislitel noi Matematiki Matematicheskoi Fiziki, 7(4):784-80) 802.
Tausworthe, R. (1965). Random numbers generated by linear recurrence modulo two. Mathematics of Computation, 29(90))201-209
Teh, Y. W. Thiery, A. H., and Vollmer, S. J. (2016). Consistency and fluctuations for stochastic gradient Langevin dynamics. Journal of Machine Learning Research, 17.
Tribble, D. (2007). Markov chain Monte Carlo algorithus using completely uniformly distributed driving sequences. PhD thesis, Citeseer.
Vollmer, S. J., Zygalakis, K. C., and Teh, Y. W. (2016). Exploration of the (non-) asymptotic bias and variance of stochastic gradient Langevin dynamics. The Journal of Machine Learning Research, 17(1):504-5548.
Welling, M. and Teh, Y. W. (2011). Bayesian learning via stochastic gradient Langevin dynamics. In Proceedings of the 28th international conference on machine learning (ICML-11), pages 681-688.
Xu, P., Chen, J., Zou, D., and Gu, Q. (2018). Global convergence of Langevin dynamics based algorithms for nonconvex optimization. Advances in Neural Information Processing Systems, 31.
12

---

A Proofs
A.1 Proof of Theorem 41
We start by decomposing the error 1, Ek=1 f(0x) (f)) into three parts
f(0)) T() - Σ f(Ox) Σ fe(w(())
(())
+ Σ fe(w - n(f) + l 211ll
k=
k=8++1
k=2+1
k=8+1
(1) + (11) + 26
||
We first upper bound (7).
Lemma (Upper bound of (I); adapted from Lemma 6.1.4 of Chen 2011)). )) If the transition map vb a contraction with parameter p and if f 1 1-Lipschitz, then
( [(0x)1 <
||0:1| + E. [1011] d².
max USiSn
Proof of Lemmal
7
Note that
(e(( ( f(0x) - 1f(06(0. w/ - f(ev((0-ee w T(d))
/ I1(eeo, wx (1e w/ (dl)
S P / 10. w (P)) (1he1 (0 (/ [1]
s r² / 110 0k -x-ell=(d())
$ rl max
|0:|1| + E. [1@11]).
0si5n
To bound (If), note that Ex=ct1 Fe( (2)) is estimating
E
[e([(())] = (he(O (( (d)))dd =: =) PP(ff).
Here, PP denote the distribution of the E-step state Og starting from Oo T. So we have the further decomposition
(TT)
< k=t+1
Fe(w(() (1)
n(5)
n
n
< <Im(s) PP(())+
Σ -(w()))=PP(ll)
TPLS))
n -1111//bb
k=t+1
L
(11) (1I)" + n
1l
The first term (ITY is due to the discretization in time. The second term (IT)" is the numerical integration error.
To bound (1TY', we use the following result.
Lemma 2 (Upper bound on discretization error (ITY), Under Assumption 7 we have for f Lipschitz.
3/2 L
(T())
TP(())
12/2d.
2 M
13

---

Proof of 2 We let f(t) be the continuous -atrrooouu---eeeeee Langevin diffusion with 0(0) = 8g sid
Wan1 - We VhEc+I, where Ex+1 ~ (O, Ia), lx = kh. So we have
0((++1)= = =((x)) - p6-1
(U(O(s))ds+ + 22h
and
0x+1 = 0k (O)) + V2h5x+1.
Combing the previous two equations gives
12t1 VU(((s)) - (U(((t)))
0(tx+1)= - -0x+1= = f((x)) - 02 - - VU((
Let D2 = 0(tx) Os. The last display reads
fx12 VU(((s)) - vU(((ts))ds..
Aut1 As (0, A+) VU(Ow)
By the contracting property 6 the main paper,
| (0x + A+) VU(O))II & pl1An1I.
Taking expectation and use L-smoothness of U, have
E [| |||
- pB | -
[1]((s) (t+x)]]] ds.
By Lemma 3 of Dalalyan and Karagulyan 2019 E (I) | - Ld. So we have E ([)(
v dE VU [ < VLd Because ((t) is a stationary process,
E [1]((s) (l+))]]
EE [[((()) (
h
l (U(()s)) v2w.l| de
E
[IIVU(((s))II] dadt + 1/' /2E [[] dt
E
00
h?
VLd+ /-¹ 2t [[] dt.
Note that
E [|| = v2 Tdd/2- - 1/2)
d+1 )1/2 = Vd+I.
S V2
I(d/2)
2
Thus,
3/2 3/2
E [1]0(s) - 0(22))ll ds - 11/2122
2
V2
3v2
x 1,3/2
1/2
2
x 3/2 1//2d
Denote r 3/2 Lh3/2d. So
E | S pE || + A+E A+
Σ p'r
2 T
h 1/2
1-p 2 M
14

---

Therefore, for any h: 2 1,
(T()) TP(( = = - [f(OA]] & & (I)((((ta)) f(xx)
E || 3/2 L 1/2d
M
If we use noisy gradient ((0k) (0() + eg where eg the noise with mean zero and bounded variance such that E(I|el] ()) S σ², then an extra term 2/a will appear in Lemma [ As o² is usually expected to be proportional to the dimension, this additional term is of the same order as the other term.
Theorem A.1 (Theorem 9.8 of Niederreiter (1922). Let vo, Ui.... be an LFSR with offset and period n 2m which satisfy ged(m, n)
Then the sequence {u: 1--) C [0,1]* with
u; (1i, 0. 1i+s-1) has, on average, star-discrepancy
O(n (nog n))+t log logn)
with an implied constant depending only on d and the average is taken over all primitive polynomials over GF(2) of degree m.
Proof of Theorem 4.1
The error on the left-hand-side is bounded by
48
(1) + (ITY + (II)" + Ilf
n
Lemmal Cermma[showw shows that (f) < (maxo |0 [1011)) I (maxosisn |1011+ [101]]) since ê = [(1/2) log, 11. Lemma shows that (IIY I 3v2 L 1/2 Denote C2 = maxosisn |00ll - M
E, [[101]] + 3v2 I. d. So (7) (rry 42/
2 M7
By Theorem A.1 and the condition that god(dl, n) 1, the discreppnyyy D' (( h221) is upper bounded by Onn (log nı) di)d+1 log log n)). Finally, by Boksma-Hlawwke incquality, have (II)" < ||fellık D ( }z21). Thus, (II)" + 14 I1fllo S Cin 1+8 where hides the poly- logarithmic terms in log n and C₁ depends on d, e, I|felluk.
Therefore, the upper bound becomes
(1) + (11) + (11)" + 48 Il S Cin-1+6 12h
B Additional numerical results
The primary contribution of this work is to improve LMC Monte Carlo sampling algorithm, not as an optimization algorithm. Therefore, our main focus is on providing better estimation of (J) for some function of interest. Downstream tasks relying on such expectations can also benefit from LQMC. For posterior prediction, solely determined by the sampling method. Even with infinite perfect samples from the posterior, the prediction error can still arise due to model misspecification. noisy data, biased sampling, etc. So the improvement achieved by LQMC might be less pronounced when assessing the prediction error.
essential to recognize that the prediction error is not
To investigate the performance LQMC in posterior prediction setting, we conducted experiments similar to those presented in Dubey et al. (2016) using three UCI datasets. Each dataset was split into a training set (70%), a validation set (10%), and a test set (20%). We performed a tuning process for the constant step size on grid using the validation set and evaluated the prediction error on the test set. Fach iteration computes the stochastic gradient using 32 data points sampled at random. Details of the datasets are in Table
15

---

Datasets
Parkinsons
Bike
Protein
(number of instances)
5875
17379
45730
p (number of features)
21
12
9
Table 1: Summary of datasets used for Bayesian posterior prediction.
The results are presented in Figure 6 The z-axes represent the total number iterations Langevin algorithm and the y-axes represent the test error. The error bars represent the variation across 10 random replicates. evident that LQMC reduced the test error, although the improvement not substantial. This aligns with our initial expectation, as the proposed method primarily enhances the accuracy of estimating the posterior mean. However, the test error often consists of other sources of error, thus the improvement achieved by the proposed method in reducing the test error might be limited.
Parkinscts
Pholdt
97225
0.00
0.18
Q11
190055
87220
390093
0:10
0.110
830089
10
214
211
24
212
212
214
3-1
Figure 6: Test error versus number of iterations for the three UCI datasets.
16


Reviews:
[{"review_id": "RIadA2LHqv", "reviewer": "Reviewer_x8C1", "summary": "The paper analyses a variant of Langevin Monte Carlo which, instead of using standard Gaussian random variables (using standard random number generation), instead uses correlated random variables following the quasi-MCMC literature. The paper demonstrates that this provably improves the efficiency of certain statistical tests, such as e.g. mean estimation.\n", "questions": "Are the assumptions of smoothness, strong convexity and bounded variation (of f) necessary? I would recommend a deeper exploration of these conditions in order to enrich the paper.\n\nAlthough stochastic gradient LMC appears in the experiments, it is not analysed. Could similar results be established in this setting?\n\nCould this approach be combined with e.g. Metropolis/HMC algorithms, or other MCMC methods based on similar walks? It seems there may be some significant theoretical barriers to this end.\n\nTypos: \nSome citations are not properly formatted, e.g. in L. 22, L. 126, L. 151 (using citet instead of citep or vice versa).\n", "limitations": "None beyond those raised in my earlier comments.", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The usage of non-random Gaussians (using an LSFR sequence instead of standard Mersenne twister) is a novelty for Langevin Monte Carlo. Consequently, the Algorithm is entirely novel, and is provably better in the sense shown in Theorem 4.1. This improvement is a significant technical innovation and is of use to both theoreticians and practitioners.\n\nThe experiments are thorough, in that many settings are considered, both in terms of convexity, dimensionality and type of gradient. The improvement for LQMC is clear and significant, which therefore lends credence to the theoretical claims.\n\nThe paper is very well written and the proofs are cleanly presented. I could not find any issues with the results.\n", "weaknesses": "The improvement of this method is weak, in the sense that the expectation of certain test functions converges . In contrast, standard Langevin MCMC results will hold in “stronger” measures of convergence, such as KL divergence, total variation, etc. This is expected since quasi-MCMC methods typically can only outperform in these specific circumstances, but is nonetheless a disadvantage of the work.\n\nThe analysis seems to be rather trivial and combines the standard inequalities in previous works on quasi-MCMC, with the standard contractivity results of LMC-type algorithms under strong convexity. \n\nOverall, while the analytic simplicity of this paper makes the result seem obvious in hindsight, I would still argue that this result is novel and impactful enough to merit publication.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\nThank you for your valuable feedback.\n\nWe appreciate your recognition regarding the differences in our analysis compared to the standard literature on Langevin Monte Carlo (LMC). Traditional results in LMC primarily focus on studying convergence through metrics like the KL divergence between the law of $x_T$ and the target distribution $\\pi$, thus investigating the distributional properties of the sample drawn at the $T$-th iteration in the LMC algorithm. In contrast, our analysis centers on the convergence of the ergodic average $\\frac{1}{T}\\sum_{k=1}^T f(x_k)$ to the target $\\pi(f)$. This perspective directly examines the behavior of the average function evaluations over the LMC samples, providing a more direct evaluation of the quality of LMC in estimating expectations. While our approach differs from the conventional ones, we believe that our result is not inherently weaker for the following reasons:\n1. Firstly, the convergence rate of the law of $x_T$ does not necessarily imply the convergence rate of the ergodic average $\\frac{1}{T}\\sum_{k=1}^T f(x_k)$. \n2. Furthermore, our result, while stated as the convergence of expectation, also implies the convergence of the empirical distribution $\\frac{1}{T}\\sum_{k=1}^T\\delta_{x_k}$ to the target distribution in the sense of the Kolmogorov-Smirnov distance (i.e., star discrepancy). In fact, this convergence of star discrepancy is the reason why we can achieve smaller errors for integrands of bounded variation in the sense of Hardy and Krause.\n3. Lastly, in practical applications, researchers often utilize all the LMC samples obtained (excluding possible burn-in and thinning) rather than focusing solely on the last single sample. Consequently, the quality of the empirical distribution $\\frac{1}{T}\\sum_{k=1}^T\\delta_{x_k}$ is of greater practical relevance than the law of $x_T$.\n\nWhile we acknowledge that our analysis introduces a nonstandard technical condition of bounded variation and establishes a different sense of convergence, we do not perceive this as a weakness. On the contrary, we believe this nonstandard approach enriches the standard literature and offers a new perspective to the understanding of LMC.\n\nBelow is a point-by-point response to your questions:\n- **Are the assumptions necessary?** Please refer to the \"global response\" for a discussion on the assumptions of smoothness, strong convexity, and bounded variation. In particular, we discussed the possibilities to relax the assumptions of smoothness and strong convexity, and provided some sufficient conditions to check bounded HK variation.\n- **Analysis for stochastic gradients:** We appreciate your suggestion. If we use a noisy gradient $\\hat g(\\theta_k)=\\nabla U(\\theta_k)+e_k$ where $e_k$ is the noise with mean zero and bounded variance such that $\\mathbb{E}(||e_k||_2^2)\\leq\\sigma^2$, then an extra term $2h\\sigma$ will appear in Lemma 2 in the proof. As $\\sigma^2$ is usually expected to be proportional to the dimension $d$, this additional term is of the same order as the other term. If the stochastic gradient $\\hat g$ is estimated by a subsample, then $\\sigma^2$ might grow with the sample size of the dataset, making this extra term dominate, as pointed out in [1]. We have included this analysis in the revision. \n- **Combining with Metropolis/HMC algorithms:** Algorithmically, incorporating a rejection step is straightforward by generating an additional uniform random variable. However, one crucial consideration that prevents us from adopting a rejection step in this work is that the rejection step introduces discontinuities, which is not QMC-friendly. In contrast, the unadjusted Langevin algorithm (ULA) is continuous in the underlying uniform random variables, making it more favorable for QMC. Beyond the issue of discontinuities, the rejection step may lead to slow mixing when the rejection rate is high. Thus, while we acknowledge the straightforward implementation of a rejection step, we agree with you that there may be theoretical barriers to this end.\n- Thank you for bringing the typos to our attention. We have addressed them in the revision.\n\n[1] Dalalyan, A. S., & Karagulyan, A. (2019). User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient. Stochastic Processes and their Applications, 129(12), 5278-5311."}, {"review_id": "ChFztsWXIy", "reviewer": "Reviewer_amBd", "summary": "In this paper, the authors proposed the Langevin quasi Monte Carlo algorithm (LQMC), that is, to use a completely uniformly distributed series (CUD quasi random number) in the Langevin Monte Carlo instead of iid pseudo random numbers. The quasi-random number is generated by the LFSR method. For a smooth and strongly convex potential, the estimation error scales as $O(n^{-1 + \\delta})$, where $\\delta$ encodes the dependence on $\\log n$ and $d$. The numerical experiment is performed for various $d$, using gradient or stochastic gradient. In all the cases, LQMC has better performance than the Langevin MC.", "questions": "Overall, I think this is already a very solid paper. However, I think it would be nice if the authors:\n* State what this algorithm could do in reality. New experiments using realistic data are welcomed. Discussions on the time scale of the realistic tasks are also interesting (for example, is there any task that was impossible and made possible by LQMC).\n* Write the error bound in the format of equation (4), that it, express the $d$ dependence explicitly.", "limitations": "My only concern regarding the limitation is that the dependence on $d$ of estimation error is not well studied. After all, $d$ is the criterion of when we should use lattice rules, quasi MC, or MC. However, I understand that the dimension of quasi MC is itself a tricky problem, and, in my humble opinion some ambiguity in this regard could be forgiven.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* Originality: It appears to me that the combination of CUD and Langevin MC is new.\n* Quality: This paper contains rigorous derivation of error bounds and numerical experiments covering a wide range of cases. Both derivation and experiment seem convincing for me.\n* Clarity: This is a well written paper in general. This paper demonstrates the backbone idea well by using examples and pseudo codes. This paper also provides a good context for the discussion. The setup for the experiments are stated clearly, and the numerical results are interpreted. I also appreciate attaching the code together with the paper.\n* Significance: This paper shows that, without making any change to the Langevin MC, changing the random number in the algorithm can lead to vast performance gain. If the performance gain shown in the paper migrates to realistic problems, the precision of a lot of computational work could be improved without a lot of software development.", "weaknesses": "* The numerical experiments are performed with synthesized data instead of real ones.\n* Unlike classic Monte Carlo, the scaling of error of quasi-Monte Carlo depends on the dimension of the system $d$. As a result, when $n$ is small and $d$ is large, quasi-Monte Carlo may perform worse than classic Monte Carlo in theory. In this paper, the error still depends on $d$ if I am not mistaken. The $d$ dependence is encoded in $\\delta$, and not shown explicitly in the main text. Although, in the numerical experiments, it is shown that the LQMC performs very well for $d=100$, I still think it would be nice to make the reader know that, in theory, the error of LQMC depends on $d$.\n* The error bound is given in a rather constraint setting.\n\nThe weaknesses and questions are addressed by the authors in their rebuttal.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\nThank you for your valuable feedback. We address your questions as follows:\n- **Numerical experiments are performed with synthesized data:** We understand the importance of real data experiments and have taken your suggestion into account. We have conducted new experiments using realistic data and tasks. Please refer to the \"global response\" and the uploaded PDF for the detailed description and results of these new experiments.\n\n- **Dimension dependence error for QMC:** We acknowledge that the standard QMC error bound of $O(n^{-1}(\\log n)^{d})$ may grow larger than $O(n^{-1/2})$ for large dimensions with insufficient $n$. In practice, however, QMC often performs better than the theoretical bounds suggest. The success of QMC in high-dimensional integration can be attributed to the concept of *effective dimension* [1]. Functions with low effective dimensions can be well approximated by a sum of low dimensional functions, which is favorable for QMC. Moreover, there exist established methods to reduce the effective dimension in the traditional QMC setting, e.g. [2-4]. How to extend these techniques to the Markov chain setting to reduce the effective dimension is an interesting future direction.\nWe appreciate your advice, and to highlight the dependence on dimension, we have updated the term $O(n^{-1+\\delta})$ to $O(n^{-1}(\\log n)^{d})$ in our paper.\n\n- **Error bound is given in a constraint setting:** We acknowledge that the assumptions of strong convexity, smoothness, and bounded variation on the integrand are restrictive. We have discussed these assumptions in the \"global response\" and mentioned possibilities to relax them. We hope that the theoretical analysis, despite its constraints, provides valuable insights into the performance of our method and its potential for more general scenarios.\n\nReferences:\n\n[1] Wang, X., \\& Fang, K. T. (2003). The effective dimension and quasi-Monte Carlo integration. Journal of Complexity, 19(2), 101-124.\\\n[2] Moskowitz, B., & Caflisch, R. E. (1996). Smoothness and dimension reduction in quasi-Monte Carlo methods. Mathematical and Computer Modelling, 23(8-9), 37-54.\\\n[3] Imai, J., \\& Tan, K. S. (2004). Minimizing effective dimension using linear transformation. In Monte Carlo and Quasi-Monte Carlo Methods 2002 (pp. 275-292).\\\n[4] Xiao, Y., & Wang, X. (2019). Enhancing quasi-Monte Carlo simulation by minimizing effective dimension for derivative pricing. Computational Economics, 54, 343-366.\n"}, {"review_id": "j32wHNSvmt", "reviewer": "Reviewer_fUNL", "summary": "For suitably smooth functions defined on a bounded support, it is well known that quasi Monte Carlo (QMC) can achieve faster rates of convergence than standard Monte Carlo when the goal is to integrate the function. This paper studies whether techniques from QMC can be beneficial for Langevin Monte Carlo (LMC), a sampling procedure based on the discretization of stochastic differential equations. The paper shows both theoretically and empirically that when the function $f$ is suitably smooth, one can augment the standard Gaussian perturbations in LMC with a specialized sequence such that the error rates converge at a faster rate than standard LMC. ", "questions": "[Q1] The main result requires $\\bar{f}_l$ to have bounded variation. Is this possible to show this for any of your examples, i.e. if f is Lipschitz and the target distribution is Gaussian? It would be great to know whether this assumption is weak or strong. \n\n[Q2] Given the main results assumptions on d, m and n, can this result always be applicable for values of n going to infinity?\n\n[Q3] The empirical examples only demonstrate the utility of LQMC on toy models. Can the authors demonstrate its usefulness in a Bayesian context on real data, say for improving the posterior predictions associated with some application?", "limitations": "It would be great if the authors mentioned the importance of doing proper inference in ML.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The paper is very easy to follow and understand. The authors do a great job of laying out the appropriate context, defining the necessary notions, and laying out their new procedure. \n\nThe idea, while relatively simple to implement, seems quite powerful. It’s simplicity seems like a feature rather than bug, and should make it relatively easy to incorporate into open source implementations.\n\nThe main result and empirical work seem sound. And while QMC is not new, incorporating it into LMC would be a nice extension of work that has previously been confined to relatively small use cases. ", "weaknesses": "The main critiques of this paper are associated with its overall applicability. On the theoretical side, the authors make some assumptions to prove their result. How strong are these assumptions? Can they apply generally as the number of samples increases?\n\nIn terms of the empirical work, all the examples are done on toy models where the final distributions are known or can be estimated for long chains. Can the authors show LQMC is useful for problems with real data? That would bolster the case for using this as a drop in replacement for vanilla LMC. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable feedback. We address your questions as follows:\n\n- **Strong assumptions/bounded variation:**  In the \"global response\", we have provided a comprehensive discussion on the smoothness, strong convexity, and bounded Hardy-Krause variation assumptions.\nIn particular, to address your question regarding the case where $f$ is Lipschitz and $\\pi$ is Gaussian, we consider integrands of the form $f\\circ \\Phi^{-1}:[0,1]^d\\to \\mathbb R$, where $\\Phi^{-1}$ is the component-wise inverse CDF of the standard Gaussian distribution. According to [1], a sufficient condition for QMC to achieve the rate $O(n^{-1+\\delta})$ for the integrand $f\\circ \\Phi^{-1}$ is that, for arbitrarily small $B_i> 0$, there exists $C$ such that\n$$\n|\\partial^u f(z)|\\leq C\\prod_{i=1}^d [1 - \\Phi(|z_i|)]^{-B_i}\n$$\nfor any $u\\subset [d]$.\nThis condition implies that as $|z_i|$ tends to infinity, the first-order mixed partial derivatives of $f(z)$ should not grow faster than $1/(1-\\Phi(|z_i|))^{B_i}$ for any $B_i>0$. Therefore, when $d>1$, Lipschitz continuity on $f$ is not sufficient nor necessary for this sufficient condition; all first-order mixed partial derivatives $|\\partial^u f(z)|$ need to grow relatively slowly.\n- **Can they apply generally as the number of samples increases?** With the increase in the number of samples, the central limit theorem (Bernstein–von Mises theorem) comes into effect, bringing the posterior distribution closer to a normal distribution, where smoothness and strong convexity assumptions hold. However, a challenge arises when utilizing a subsample to estimate the gradient, as the error of the gradient estimator usually scales with the sample size. In such cases, as highlighted by [2], the stochastic gradient's noise dominates, potentially diminishing the applicability of the paper's theorem.\n- **Can the result be applicable for $n$ going to infinity?** In our notation, $n$ is the number of iterations rather than the dataset size. However, we interpret your question as being about the dataset size approaching infinity, and this aspect has been addressed above.\n- **Real example:** We appreciate your interest in real examples and have taken your suggestion into account. We have included new experiments using realistic data and tasks. In particular, we have introduced an example of sparse regression (commonly used for Bayesian variable selection) and several prediction tasks on real datasets. We kindly direct you to the \"global response\" and the uploaded PDF for detailed descriptions and results of these new experiments.\n\nReferences\n\n[1] He, Z., Zheng, Z., & Wang, X. (2023). On the error rate of importance sampling with randomized quasi-Monte Carlo. SIAM Journal on Numerical Analysis, 61(2), 515-538.\\\n[2] Dalalyan, A. S., & Karagulyan, A. (2019). User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient. Stochastic Processes and their Applications, 129(12), 5278-5311.\n"}, {"review_id": "KhSfBWgbcz", "reviewer": "Reviewer_pgYw", "summary": "This paper analyzes the effect of using quasi-random numbers in place of the usual IID Gaussians for the driving noise of a Langevin algorithm. Assuming that the loss is strong convex and a the quasi-random numbers are completely uniformly distributed, a bound on the Monte Carlo estimation error is derived. This bound is substantially better than what is currently achieved for IID Gaussian driving noise.", "questions": "* Can the result be extended to common measures between distributions? For example, the result looks almost like a bound on the 1-Wasserstein distance between the distribution of the iterates and the stationary distribution. The main difference is that it assumed that $f$ is both bounded and Lipschitz, as opposed to just Lipschitz.\n* Are there general classes of strongly convex functions for which the you can check that the Hardy-Krause variation of $\\bar f_{\\ell}$ is finite?", "limitations": "These are adequately addressed.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The paper makes a nice contribution to the theory of sampling with Langevin dynamics, by bringing in the use of quasi-Monte Carlo methods. The proof is quite clean, and it is very clear where the power of the proposed method becomes useful. \n\nThe experiments are promising.", "weaknesses": "My only real concern is that this paper focuses on a relatively easy case, albeit an important one. It would be interesting to see how easy or difficult it is to extend the results to more complex scenarios, such as with non-convex losses or external random variables.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable feedback. We address your questions as follows:\n\n- **Extension to non-convex losses or external random variables** In the \"global response\", we provided a comprehensive discussion about the strong convexity and smoothness assumptions and discussed to what extent these assumptions might be relaxed. For instance, the strong convexity might be relaxed to functional inequalities, and the Lipschitz gradient might be relaxed to weaker Hölder-continuous conditions. In general non-convex settings, however, the improvement of LQMC might not be as significant, since the advantage of LQMC is more pronounced when the Markov chain mixes well. Nevertheless, in the optimization literature, QMC has demonstrated better exploration capability than traditional Monte Carlo. Hence, it remains an interesting open question to investigate whether LQMC exhibits better exploration than LMC in non-log-concave sampling scenarios. \\\nRegarding external variables, our current approach treats the dataset as deterministic, and the posterior distribution is treated no differently from any other target distribution. However, if the dataset is modeled as i.i.d. or a mixing Markov chain, it would be interesting to understand how the sampling of the dataset interacts with the sampling of the Gaussian variables in the LMC algorithm. In the traditional QMC setting, there are studies on combining QMC samples with external samples such as those obtained by acceptance-rejection algorithms or MCMC. How to incorporate external samples into the LMC algorithm is an interesting direction for future work.\n- **Extension to common measures between distribution** To clarify, the question asks whether we can improve the convergence of certain metric such as KL divergence between the law of $x_T$ and the target distribution $\\pi$.\\\nWe do not anticipate the distance between the law of $x_T$ and the target distribution $\\pi$ to be smaller by using LQMC as opposed to LMC. In fact, in a well-mixing scenario, the law of $x_T$ obtained by LQMC should be roughly the same as that of LMC. \nHowever, the key advantage of QMC lies in the collective behavior of the sequence of samples $\\{x_1,x_2,\\ldots,x_T\\}$. Specifically, we expect the empirical distribution $\\frac{1}{T}\\sum_{k=1}^T\\delta_{x_k}$ to be closer to the target $\\pi$ in the Kolmogorov-Smirnov distance (also called the star discrepancy). In other words, while each individual sample is not closer to the target, the ensemble of all samples collectively provides a better approximation of the target distribution. As a result, the ergodic average $\\frac{1}{T}\\sum_{k=1}^T f(X_k)$ provides a more accurate estimate of $\\pi(f)$ with LQMC. Furthermore, we believe the ergodic average is of more practical importance since researchers commonly utilize all the samples generated by the Markov chain (excluding potential burn-in or thinning) rather than relying solely on the last individual sample. Thus, the empirical distribution $\\frac{1}{T}\\sum_{k=1}^T\\delta_{x_k}$ is more closely aligned with reality than the distribution of $x_T$.\n\n- **Check bounded variation condition** In response to your question about the bounded Hardy-Krause variation condition, we provided some sufficient conditions and discussions in the \"global response\".\n\nReferences: \n\n[1] Hintz, E., Hofert, M., & Lemieux, C. (2022). Quasi-random sampling with black box or acceptance-rejection inputs. In Advances in Modeling and Simulation. Cham: Springer International Publishing.\n"}]

Author rebuttal:
# Global response
We thank the reviewers for their insightful and constructive feedback. The positive evaluations of our work as "novel, impactful, significant, promising, and important" are encouraging, and we appreciate the recognition of the clarity and cleanliness of the paper. In this global response, we address the common aspects raised by multiple reviewers regarding the strong assumptions in the theoretical analysis and the absence of real data examples. Specific response to individual reviewers are provided separately.

## Strong convexity and smoothness
- We acknowledge that these are restrictive assumptions. However, starting with these assumptions is a necessary and natural first step to study the proposed LQMC algorithm. Evaluating its performance in the simplest setting allows us to gauge its potential for more complicated scenarios. 
- The simplicity of the analysis yields valuable insights into how the driving uniform variables impact LMC. The error decomposition provides insights regarding which aspect of the algorithm might benefit from QMC. Despite its simplicity, the theoretical analysis guides the design of uniform numbers and sets expectations for the algorithm's performance, which might not be as clear in a more general setting.
- Next, we recognize the opportunity to relax these assumptions. Techniques similar to [6] can be adopted in non-strongly convex settings by introducing a quadratic penalty. A broader class of measures satisfying log-Sobolev or Poincaré inequalities also presents convergence guarantees [7]. Relaxing the Lipschitz gradient to a weaker Hölder-continuous gradient might be possible [5]. Growth conditions like dissipativity might also ensure convergence [9]. Quantifying the improvement of LQMC in these more general yet tractable settings is an interesting future direction, albeit beyond the scope of this work.
- While the improvement of LQMC might not be significant in general non-convex settings, where only first-order stationarity guarantees can be expected, the question of whether LQMC enhances exploration in non-convex settings remains open. Encouragingly, studies in optimization demonstrate QMC's superior exploration properties, such as in reinforcement learning, where QMC outperforms iid Gaussian for parameter exploration [10], and in variational inference, where QMC converges faster compared to traditional Monte Carlo [11]. These findings highlight QMC's potential benefit for exploration.

## Bounded Hardy-Krause variation
The condition of bounded HK variation is indeed a standard condition in QMC theory, ensuring the error rate of $O(n^{-1+\delta})$ for any $\delta>0$, with log terms hidden. While verifying this condition is challenging in general, there exist sufficient conditions [1,2]:
$$V_{HK}(f)\leq\sum_{u\subset1:d,u\neq\emptyset}\int_{[0,1]^{|u|}}|\partial^uf(x_{u},1_{-u})|dx_u.$$
In our analysis, we require $\bar f_{\ell}=f\circ g$ to have bounded HK, where $g=\psi_{\ell}:[0,1]^d\to G$ is the $\ell$-step transition. Consider two cases:
- **$G$ is compact:** If all the mixed partial derivatives of $g$ are square integrable, then $f\circ g$ has bounded HK for any $f\in C^d(G)$ [1]. Hence, if LMC samples are constrained within a compact set and the transition is sufficiently smooth, then $\bar f_\ell$ has bounded HK for any $f\in C^d$.
- **Otherwise:** If the function grows not too rapidly on the boundary of the unit cube, QMC can still achieve the error rate of $O(n^{-1+\delta})$. Some boundary growth conditions are introduced in [8] to ensure this rate.

Verification of these sufficient conditions is easier if such mixed partial derivatives are bounded. In other cases, one needs to examine the function's growth at singularities. It is essential to note that the bounded variation condition is a technical requirement for QMC theory. Empirically, QMC can still outperform Monte Carlo without this condition, as its performance largely depends on factors like smoothness and effective dimension.

## Real data experiments
Real data experiments were conducted in response to reviewers' suggestions. We first emphasize that the primary contribution of this work is to *improve LMC as a Monte Carlo sampling algorithm, not as an optimization algorithm*. Therefore, our main focus is on providing a better estimation of $\pi(f)$ for some function of interest. Downstream tasks relying on such expectations can also benefit from LQMC. For posterior prediction, it is essential to recognize that the prediction error is not solely determined by the sampling method. Even with infinite perfect samples from the posterior, the prediction error can still arise due to model misspecification, noisy data, biased sampling, etc. So the improvement achieved by LQMC might be less pronounced when assessing the prediction error.

The results are presented in the uploaded PDF. The additional experiments provide a more comprehensive evaluation of LQMC's effectiveness in realistic tasks.

[1] Basu et al. (2016). Transformations and HK variation. SINUM.\
[2] Owen. (2005). Multidimensional variation for QMC. Contemporary Multivariate Analysis And Design Of Experiments.\
[3] Dalalyan et al. (2012). Sparse regression learning by aggregation and LMC. Journal of Computer and System Sciences.\
[4] Dubey et al. (2016). Variance reduction in SGLD. NeurIPS.\
[5] Chatterji et al. (2020). LMC without smoothness. AISTATS.\
[6] Dalalyan et al. (2022). Bounding the error of discretized Langevin algorithms for non-strongly log-concave targets. JMLR.\
[7] Chewi et al. (2022). Analysis of LMC from Poincare to Log-Sobolev. COLT.\
[8] Owen. (2006). Halton sequences avoid the origin. SIAM Review.\
[9] Erdogdu et al. (2021). On the convergence of LMC: The interplay between tail growth and smoothness. COLT.\
[10] Choromanski et al. (2018). Structured evolution with compact architectures for scalable policy optimization. ICML.\
[11] Buchholz et al. (2018). QMC variational inference. ICML.

Meta-review: This paper considers sampling with a variation of Langevin Monte Carlo algorithm. Authors suggest using correlated random variables as in the quasi-Monte Carlo framework. Authors show that this improves the performance of vanilla LMC in which standard Gaussian random variables are used.

The paper is reviewed by 4 expert reviewers. The Rating/Confidence  scores are 6/3, 7/3, 7/3, 7/4 and reviewers all agree that the paper is well-written with significant contributions. I recommend accepting this paper to the conference program.

I recommend authors to incorporate the relevant part of the discussion period in their camera ready version. In particular, a substantial portion of the recent literature is not mentioned in the paper.
Decision: Accept (poster)