# Langevin Monte Carlo: random coordinate descent and variance reduction

Zhiyan Ding<EMAIL_ADDRESS>
Mathematics Department
University of Wisconsin-Madison
Madison, WI 53705 USA. Qin Li<EMAIL_ADDRESS>
Mathematics Department and Wisconsin Institutes for Discovery
University of Wisconsin-Madison
Madison, WI 53705 USA

###### Abstract

Langevin Monte Carlo (LMC) is a popular Bayesian sampling method. For the log-
concave distribution function, the method converges exponentially fast, up to
a controllable discretization error. However, the method requires the
evaluation of a full gradient in each iteration, and for a problem on
$\mathbb{R}^{d}$, this amounts to $d$ times partial derivative evaluations per
iteration. The cost is high when $d\gg 1$. In this paper, we investigate how
to enhance computational efficiency through the application of RCD (random
coordinate descent) on LMC. There are two sides of the theory:

  * •

By blindly applying RCD to LMC, one surrogates the full gradient by a randomly
selected directional derivative per iteration. Although the cost is reduced
per iteration, the total number of iteration is increased to achieve a preset
error tolerance. Ultimately there is no computational gain;

  * •

We then incorporate variance reduction techniques, such as SAGA (stochastic
average gradient) and SVRG (stochastic variance reduced gradient), into RCD-
LMC. It will be proved that the cost is reduced compared with the classical
LMC, and in the underdamped case, convergence is achieved with the same number
of iterations, while each iteration requires merely one directional
derivative. This means we obtain the best possible computational cost in the
underdamped-LMC framework.

Keywords: Langevin Monte Carlo, Random coordinate descent, Variance reduction,
Bayesian inference, Wasserstein metric

##  1 Introduction

Monte Carlo Sampling is one of the core problems in Bayesian statistics, data
assimilation (Reich, 2011), and machine learning (Andrieu et al., 2003), with
wide applications in atmospheric science (Fabian, 1981), petroleum engineering
(Nagarajan et al., 2007), epidemiology (Li et al., 2020), in the form of
inverse problems, volume computation (Vempala, 2010), and bandit optimization
(Mazumdar et al., 2020; Russo et al., 2018; May et al., 2012).

Let $f(x)$ be $\mu$-strongly convex with its gradient being $L$-Lipschitz in
$\mathbb{R}^{d}$, and define the target density function $p\propto e^{-f}$.
Then $p(x)$ is a log-concave probability density function. To sample from the
distribution induced by $p(x)$ amounts to finding an $x\in\mathbb{R}^{d}$ (or
a list of samples that can be regarded as i.i.d. (independent and identically
distributed) drawn from this distribution.

The literature is very rich on sampling, and there are many approaches (Neal,
2001, 1993; Roberts and Rosenthal, 2004; Metropolis et al., 1953; Hastings,
1970; Geman and Geman, 1984; Duane et al., 1987; Neal, 2012). In this paper,
we only study variations of Langevin Monte Carlo (LMC) (Rossky et al., 1978;
Parisi, 1981; Roberts and Tweedie, 1996b; Chen et al., 2014; Ma et al., 2015).
LMC, including the classical LMC and its variations, can be regarded as one
subcategory of Markov chain Monte Carlo (MCMC). They are attractive mostly due
to the fast convergence rate. Indeed, the classical LMC and its most
variations can be regarded as discrete versions of some Langevin dynamics, a
set of stochastic differential equations (SDE) that roughly follow the
gradient of $f$, with some added Brownian motion terms. Usually, the SDEs are
designed so the distribution converges exponentially fast in time to the
target distribution (Markowich and Villani, 1999). LMC, therefore, viewed as
the discrete versions (such as the Euler-Maruyama method) of the SDEs,
converge also almost exponentially fast, up to a discretization error. The
non-asymptotic analysis for the methods is studied in detail in recent years
(Dalalyan and Riou-Durand, 2020; Dalalyan and Karagulyan, 2019; Durmus et al.,
2019; Cheng et al., 2018; Dalalyan and Riou-Durand, 2020; Eberle et al., 2019;
Dwivedi et al., 2018).

However, like many other sampling methods, the numerical cost of classical
(both overdamped and underdamped) LMC methods increases with respect to $d$,
the dimensionality of the problem. While the number of iterations to achieve a
preset error tolerance already depends on $d$, the cost per iteration may also
increase with respect to $d$. Indeed, there are examples in which, to compute
the gradients, $d$ partial derivatives need to be computed separately (see
Section 1.2), and this adds another $d$ folds of cost per iteration.

The cost of evaluating the gradient has triggered a large number of studies
centering around the “gradient-free" property. Many approaches were
investigated, including Important Sampling (Geweke, 1989; Neal, 2001; Moral et
al., 2006; Murray et al., 2016), ensemble methods (Iglesias et al., 2013;
Garbuno-Inigo et al., 2020; Reich, 2011), and random walks (Mengersen and
Tweedie, 1996; Roberts and Tweedie, 1996b, a; Dwivedi et al., 2018). These
methods shed light on eliminating the evaluation of the gradients, but to this
day we have not found a method that is proved to be more competitive with LMC.

###  1.1 Contribution

In this article, we would like to study how to enhance the numerical
efficiency of LMC, with a special eye on improving the cost’s dependence on
$d$. The strategy we take is to develop alternatives for evaluating the
gradients. More specifically, we explore how to incorporate random coordinate
descent (RCD) in LMC. RCD is a technique developed for optimization problems,
and it randomly chooses one partial derivative, as a replacement of the full
gradient in gradient descent (GD). Since only one partial derivative is
calculated in each iteration, the cost is reduced by $d$-folds per iteration.
If one carefully controls the iteration number, the overall cost can also be
reduced in certain scenarios, as proved in the optimization literature
(Nesterov, 2012).

We study if applying RCD to LMC can bring us similar benefits. The theoretical
guarantee we obtain in this paper shows that there are two sides of the story.

Firstly, we blindly apply RCD to LMC. This is to replace the full gradient in
LMC with a randomly selected partial derivative in each iteration, exactly as
RCD in optimization. This is done to both the overdamped LMC and underdamped
LMC, resulting in two algorithms: RCD-O-LMC, and RCD-U-LMC respectively. Since
only one directional derivative is calculated instead of the full gradient in
each iteration, the cost is reduced by $d$-folds per iteration. However, we
can find counterexamples that show this blind application of RCD induces a
large error term that represents the high variance produced by the random
direction selection procedure. Exactly due to this large error term, more
iterations are required. Ultimately there is no saving in the numerical cost.

Secondly, we study variance reduction techniques. Two main techniques are
explored: SAGA (stochastic average gradient) (Schmidt et al., 2017; Defazio et
al., 2014) and SVRG (stochastic variance reduced gradient) (Johnson and Zhang,
2013), hoping to reduce the large error term mentioned above, saving the
numerical cost in the end. These methods were developed as optimization
techniques to improve the convergence rate of SGD (stochastic gradient
descent) (Robbins and Monro, 1951; Kiefer and Wolfowitz, 1952; Bottou, 2010),
an algorithm that looks for the minimizer of an objective function that has
the form of $f(x)=\sum^{N}_{n=1}f_{n}(x)$. Such techniques were then applied
to improve stochastic gradient Markov chain Monte Carlo (SG-MCMC) (Welling and
Teh, 2011; Dubey et al., 2016; Durmus et al., 2016; Baker et al., 2019). The
complexities of these methods are studied in (Chatterji et al., 2018; Zou et
al., 2018a, b; Zou and Gu, 2021). A summary of these results can be found in
(Zou and Gu, 2021, Section 2.3). We investigate how to incorporate these
techniques to enhance the performance of RCD-LMC. In our case $f$ does not
necessarily have the form of $\sum^{N}_{n}f_{n}(x)$, however, viewing $\nabla
f=\sum_{k}\partial_{k}f\textbf{e}^{k}$, variance reduction techniques for SGD
can still be borrowed. The methods with variance reduction integrated are
called Randomized Coordinates Averaging Descent Overdamped/Underdamped LMC
(RCAD-O/U-LMC) and Stochastic Variance Reduced Gradient Overdamped/Underdamped
LMC (SVRG-O/U-LMC). We will show that, with either variance reduction
technique, in the underdamped setting, the new methods require the same number
of iterations as the classical U-LMC (Cheng et al., 2018) to achieve a small
preset error tolerance $\epsilon$, but the number of directional derivative
per iteration is $1$ instead of $d$. This automatically saves $d$ folds of
computation.

We summarize all the convergence results in Table 1 (assuming computing the
full gradient costs $d$ times of computing one directional derivative and
arbitrary initial distribution). We do not present the dependence on the
conditioning $\kappa$ in the table, but they are included in the later
discussions. The presented result assumes the large $\kappa$ is a secondary
concern compared to the large $d$ ($\kappa\ll d$).

Algorithm | Number of iterations | Cost
---|---|---
O-LMC | $\widetilde{O}\left(d/\epsilon\right)$ | $\widetilde{O}\left(d^{2}/\epsilon\right)$
U-LMC | $\widetilde{O}\left(d^{1/2}/\epsilon\right)$ | $\widetilde{O}\left(d^{3/2}/\epsilon\right)$
RCD-O-LMC | $\widetilde{O}\left(d^{2}/\epsilon^{2}\right)$ | $\widetilde{O}\left(d^{2}/\epsilon^{2}\right)$
RCD-U-LMC | $\widetilde{O}\left(d^{2}/\epsilon^{2}\right)$ | $\widetilde{O}\left(d^{2}/\epsilon^{2}\right)$
SVRG-O-LMC | $\widetilde{O}\left(d^{3/2}/\epsilon\right)$ | $\widetilde{O}\left(d^{3/2}/\epsilon\right)$
SVRG-U-LMC | $\widetilde{O}\left(\max\\{d^{4/3}/\epsilon^{2/3},d^{1/2}/\epsilon\\}\right)$ | $\widetilde{O}\left(\max\\{d^{4/3}/\epsilon^{2/3},d^{1/2}/\epsilon\\}\right)$
RCAD-O-LMC | $\widetilde{O}\left(d^{3/2}/\epsilon\right)$ | $\widetilde{O}\left(d^{3/2}/\epsilon\right)$
RCAD-U-LMC | $\widetilde{O}\left(\max\\{d^{4/3}/\epsilon^{2/3},d^{1/2}/\epsilon\\}\right)$ | $\widetilde{O}\left(\max\\{d^{4/3}/\epsilon^{2/3},d^{1/2}/\epsilon\\}\right)$
Table 1: Number of iterations and numerical cost to achieve
$\epsilon$-accuracy. The results for the classical O-LMC and U-LMC come from
(Dalalyan and Karagulyan, 2019) and (Cheng et al., 2018; Dalalyan and Riou-
Durand, 2020) respectively. The notation $\widetilde{O}(f)$ omits the possible
$\log$ terms. For the overdamped cases, we assume the Lipschitz continuity for
both the gradient term and the hessian term, and for the underdamped cases,
Lipschitz continuity is only assumed for the gradient term.

###  1.2 Discussions on assumptions, and relation to the literature

Throughout the paper we assume that computing one partial derivative costs
$1/d$ of computing the full gradient. This assumption may not hold in some
applications in which the full gradient can be computed efficiently. This
happens in logistic regression, support vector machines (SVM), and deep
learning where backward propagation is heavily relied on (Rumelhart et al.,
1986).

However, there are many examples where this assumption indeed holds true. For
example, given a graph with nodes ${\cal N}=\\{1,2,\dotsc,d\\}$ and directed
edges ${\cal E}\subset\\{(i,j):i,j\in{\cal N}\\}$, suppose there is a scalar
variable $x_{i}$ associated with each node $i=1,2,\dotsc,d$, and that the
function $f$ has the form

$f(x)=\sum_{(i,j)\in{\cal E}}f_{ij}(x_{i},x_{j})\,.$

Then the partial derivative of $f$ with respect to $x_{i}$ is given by

$\frac{\partial f}{\partial x_{i}}=\sum_{j:(i,j)\in{\cal E}}\frac{\partial
f_{ij}}{\partial x_{i}}(x_{i},x_{j})+\sum_{l:(l,i)\in{\cal E}}\frac{\partial
f_{li}}{\partial x_{i}}(x_{l},x_{i})\,.$

Note that the number of terms in the summations in this expression equals the
number of edges in the graph that touch node $i$, the expected value of which
is about $2/d$ times the total number of edges in the graph. Meanwhile,
evaluation of the full gradient would require evaluation of partial
derivatives with respect to both component $i$ and $j$ of the function
$f_{ij}$ for all edges in the graph, leading to a factor-of-$d$ difference in
evaluation cost. This setup is encountered in many graph-based problems, such
as graph-based label propagation in semi-supervised learning (Chapelle et al.,
2006), finding the densest $k$-subgraph (Sotirov, 2020), and finding the most
likely assignments in continuous pairwise graphical models (Rue and Held,
2005). It is also observed in large-scale SVM problems (when solved in dual
form) (Hsieh et al., 2008) as well.

The assumption also incurs in many PDE-constrained inverse problems as well.
The evaluation of the objective function $f$ typically accounts for one PDE
solve. As a comparison, each partial derivative amounts to 2 PDE solves (one
forward encoding the input data and one adjoint presenting the output test
function), and in each iteration, to have the full gradient (also termed
Fréchet derivative in that context), many partial derivatives need to be
computed (Leugering et al., 2012). This leads to a large number of PDE solves
for only one iteration.

Another assumption we make in this paper is that the conditioning of the
objective function $f$ is only moderately big, in the sense that its value is
not comparable to the largeness of $d$. This means the underperformance of
RCD-LMC and the outperformance of the variance reduced RCD-LMC hold uniformly
true for all $f$ that have moderate condition numbers. If $f$ has bad
conditioning, the Lipschitz constants for each direction are then drastically
different. When this happens, one could potentially choose coordinates at
different rates to reflect the skewness of $f$. This indeed happens in RCD
surrogating GD, where the stiffer directions get chosen more frequently. If we
drop the assumption and let $f$ be very skewed as well, a biased selection
process could potentially enhance our methods even more. This topic is
discussed in two separate contributions, see (Ding et al., 2021a) and (Ding et
al., 2021b) for overdamped and underdamped cases respectively. For the
completeness of the paper, in Remark 5 we discuss how to sample from a non-
uniform coordinate selection process, and in Remark 19 we argue why it does
not bring computational benefits when $f$ is well-conditioned.

There are strong connections between the methods proposed in the paper and the
existing literature. We first mention that most of the techniques are borrowed
from optimization. RCD is one strategy that surrogates the full gradient by a
partial derivative in GD and has been shown to be useful when the objective
function $f$ has certain structures. Similar ideas were used in SGD when the
objective function has the form of $f=\sum^{K}_{k=1}f_{k}$, and in SGD, a
$\nabla f_{i}$ is randomly selected to surrogate $\nabla f$ per iteration. It
was observed that although SGD reduces the per iteration cost by $K$, the
total number of iteration is increased due to the induced higher variance. For
this reason, multiple variance reduction techniques were developed, including
SVRG (stochastic variance reduced gradient) (Johnson and Zhang, 2013), and
SAGA (stochastic average gradient) (Schmidt et al., 2017; Defazio et al.,
2014), all adopted in the current paper.

We also would like to mention that the authors of the paper have already
presented part of the results of this paper in (Ding and Li, 2020). In that
paper, it was already observed that blindly applying RCD to LMC will not
enhance the efficiency, and some variance reduction is needed. RCAD-O/U-LMC,
built upon integrating SAGA to RCD-LMC, was proposed and studied there. That
work was our first attempt in addressing similar issues. With the insights
gained there, we greatly extended the results, and provide the general recipe
in the current paper. The main extension can be summarized as follows:

  * •

We give the non-asymptotic convergence result of RCD-O/U-LMC in Theorem 3 and
Theorem 6. These two bounds were not included in (Ding and Li, 2020).

  * •

We propose and rigorously analyze a new variance reduction method
SVRG-O/U-LMC. This method, similar to RCAD-O/U-LMC, can also be shown to
reduce the computational cost.

  * •

Moreover, we present a systematic understanding of how the variance presents
in the convergence rate, and how to integrate variance reduction techniques to
reduce the cost. The overarching results that summarize, in a general form,
the impact of involving randomness in the coordinate selection process are
presented in Theorem 15 and Theorem 23. The two theorems address the
overdamped and the underdamped cases respectively. The blind application of
RCD, and the incorporation of SAGA and SVRG, then can be viewed as three
different examples under this framework. We emphasize that though we only
study two variance reduction techniques in the current paper, Theorem 15 and
23 are general enough to treat other methods as well. Any new variance
reduction method, when cast into the random coordinate LMC framework, can be
analyzed in a similar fashion.

###  1.3 Organization

In Section 2, we discuss the essential ingredients of our methods. In Section
3, we unify the notations and assumptions used in our methods. In Section 4,
we present the convergence results of RCD-LMC and also provide a counter-
example to demonstrate the fact that RCD-LMC does not save numerical cost. In
Section 5, we introduce the methods that incorporate variance reduction
techniques: SVRG-O/U-LMC and RCAD-O/U-LMC, and present the theorems on the
convergence and the numerical cost. In Section 6, numerical evidence will be
given to demonstrate the improved efficiency of our new methods. In Section
7-8, we provide the proof of the theorems for overdamped, and underdamped
settings respectively. We note that in subsection 7.1 and subsection 8.1 we
provide the overarching results for O/U-LMC when randomness is involved in the
coordinate selection process. The later subsections are devoted to the
realization of these two theorems on specific algorithms as examples. The
discussion on the best possible numerical cost is heavily technical and thus
delayed to these two subsections (Remark 16 and Remark 24). Some technical
lemmas used in these two sections are postponed to the appendix.

##  2 Essential ingredients

In this section, we discuss the essential ingredients of our methods: the
overdamped and underdamped Langevin dynamics and the associated Monte Carlo
methods (O-LMC and U-LMC); random coordinate descent (RCD); and variance
reduction techniques, including both SVRG and SAGA.

###  2.1 Overdamped Langevin dynamics and O-LMC

The O-LMC method comes from the overdamped Langevin dynamics, a stochastic
differential equation that writes:

$\,\mathrm{d}X_{t}=-\nabla f(X_{t})\,\mathrm{d}t+\sqrt{2}\,\mathrm{d}B_{t}\,.$
(1)

This SDE characterizes the trajectory of $X_{t}$. Two forcing terms $\nabla
f(X_{t})\,\mathrm{d}t$ and $\,\mathrm{d}B_{t}$ compete: the former drives
$X_{t}$ to the minimum of $f$ and the latter provides Brownian motion and thus
small oscillations along the trajectory. The initial data $X_{0}$ is a random
variable drawn from a given distribution induced by $q_{0}(x)$. Denote
$q(x,t)$ the probability density function of $X_{t}$, it is a well-known
result that $q(x,t)$ satisfies the following Fokker-Planck equation:

$\partial_{t}q=\nabla\cdot(\nabla fq+\nabla q)\,,\quad\text{with}\quad
q(x,0)=q_{0}\,,$ (2)

and furthermore, $q(x,t)$ converges to the target density function
$p(x)\propto e^{-f}$ exponentially fast in time (Markowich and Villani, 1999).

The overdamped Langevin Monte Carlo (O-LMC), as a sampling method, is simply a
discrete-in-time version of the SDE (1). A standard Euler-Maruyama method
applied on the equation gives:

$x^{m+1}=x^{m}-\nabla f(x^{m})h+\sqrt{2h}\xi^{m}\,,$ (3)

where $h$ is the small time-step and $\xi^{m}$ is i.i.d. drawn from
$\mathcal{N}(0,I_{d})$ with $I_{d}$ being the identity matrix of size $d$.
Since (3) approximates (1), the density of $x^{m}$ becomes $p(x)$ as
$m\to\infty$, up to a discretization error. It was proved in (Dalalyan and
Karagulyan, 2019) that the convergence to $\epsilon$ is achieved within
$\widetilde{O}(d/\epsilon)$ iterations if the hessian of $f$ is Lipschitz. If
this hessian is not Lipschitz, the number of iterations is shown to be
$\widetilde{O}(d/\epsilon^{2})$. In many real applications, the gradient of
$f$ is not explicitly known and some approximation is used, introducing
another layer of numerical error. In (Dalalyan and Karagulyan, 2019), the
authors discussed the effect of such error, under the assumption that the
error term has a bounded variance.

###  2.2 Underdamped Langevin dynamics and U-LMC

The underdamped Langevin dynamics (Chen et al., 2014) is characterized by the
following SDE system:

$\left\\{\begin{aligned} &\,\mathrm{d}X_{t}=V_{t}\,\mathrm{d}t\\\
&\,\mathrm{d}V_{t}=-2V_{t}\,\mathrm{d}t-\gamma\nabla
f(X_{t})\,\mathrm{d}t+\sqrt{4\gamma}\,\mathrm{d}B_{t}\end{aligned}\right.\,,$
(4)

where $\gamma>0$ is a parameter to be tuned. Denote $q(x,v,t)$ the probability
density function of $(X_{t},V_{t})$, then $q$ satisfies the Fokker-Planck
equation

$\partial_{t}q=\nabla\cdot\left(\begin{bmatrix}-v\\\ 2v+\gamma\nabla
f\end{bmatrix}q+\begin{bmatrix}0&0\\\ 0&2\gamma\end{bmatrix}\nabla
q\right)\,,$

and under mild conditions, in $t\to\infty$, it converges to
$p_{2}(x,v)\propto\exp(-(f(x)+|v|^{2}/2\gamma))$, making the marginal density
function for $x$ the target $p(x)$.

The underdamped Langevin Monte Carlo algorithm, U-LMC, can be viewed as a
numerical solver to (4). In each step, we sample
$(x^{m+1},v^{m+1})\in\mathbb{R}^{2d}$ as a Gaussian random variable determined
by $(x^{m},v^{m})$ with the following expectation and covariance:

$\displaystyle\mathbb{E}x^{m+1}=x^{m}+\frac{1}{2}\left(1-e^{-2h}\right)v^{m}-\frac{\gamma}{2}\left(h-\frac{1}{2}\left(1-e^{-2h}\right)\right)\nabla
f(x^{m})\,,$ (5)
$\displaystyle\mathbb{E}v^{m+1}=v^{m}e^{-2h}-\frac{\gamma}{2}\left(1-e^{-2h}\right)\nabla
f(x^{m})\,,$
$\displaystyle\mathrm{Cov}\left(x^{m+1}\right)=\gamma\left[h-\frac{3}{4}-\frac{1}{4}e^{-4h}+e^{-2h}\right]\cdot
I_{d}\,,\ \mathrm{Cov}\left(v^{m+1}\right)=\gamma\left[1-e^{-4h}\right]\cdot
I_{d}\,,$
$\displaystyle\mathrm{Cov}\left(x^{m+1}\,,v^{m+1}\right)=\frac{\gamma}{2}\left[1+e^{-4h}-2e^{-2h}\right]\cdot
I_{d}\,.$

We here used the notation $\mathbb{E}$ to denote the expectation, and
$\mathrm{Cov}(a,b)$ to denote the covariance of $a$ and $b$. If $b=a$, we
abbreviate it to $\mathrm{Cov}(a)$. The scheme can be interpreted as sampling
from the following dynamics in each time interval:

$\left\\{\begin{aligned}
&\mathrm{X}_{t}=x^{m}+\int^{t}_{0}\mathrm{V}_{s}\,\mathrm{d}s\\\
&\mathrm{V}_{t}=v^{m}e^{-2t}-\frac{\gamma}{2}(1-e^{-2t})\nabla
f(x^{m})+\sqrt{4\gamma}e^{-2t}\int^{t}_{0}e^{2s}\,\mathrm{d}B_{s}\end{aligned}\right.$

with $x^{m+1}=\mathrm{X}_{h}$ and $v^{m+1}=\mathrm{V}_{h}$, and $h$ is the
time step.

The underdamped Langevin Monte Carlo demonstrates a faster convergence rate
(Cheng et al., 2018) than O-LMC. Due to the introduction of $V$ in (4), the
trajectory $X_{t}$ is smoother than that of $V_{t}$ whose smoothness is
determined by the Brownian motion $B_{t}$. As a result, a higher-order
discretization is possible for this augmented dynamics. Without the assumption
on the hessian of $f$ being Lipschitz, the number of iterations is
$\widetilde{O}(\sqrt{d}/\epsilon)$ to achieve $\epsilon$ accuracy. A high
(3rd) order discretization was discussed for (4) in (Mou et al., 2021),
further enhancing the numerical accuracy to
$\widetilde{O}(d^{1/4}/\epsilon^{1/2})$ when $f$ is smooth enough. Similar to
the discussion for O-LMC in (Dalalyan and Karagulyan, 2019), the authors in
(Dalalyan and Riou-Durand, 2020) also studied the error in estimating $\nabla
f(x^{m})$, but they also assumed the variance is bounded.

###  2.3 Random coordinate descent (RCD)

The idea of RCD is to surrogate the full gradient in Gradient Descent by a
randomly selected partial derivative in each iteration (Nesterov, 2012). Per
iteration, only one partial derivative is computed instead of $d$, so there is
some hope that total cost can be reduced. More specifically one approximates

$\nabla f\approx d\left(\nabla
f(x)\cdot\textbf{e}^{r}\right)\textbf{e}^{r}\,,$ (6)

where $\textbf{e}^{i}$ is the $i$-th unit direction and $r$ is randomly drawn
from $1,2,\cdots,d$. This approximation is consistent in the expectation sense
because

$\mathbb{E}_{r}\left(d\left(\nabla
f(x)\cdot\textbf{e}^{r}\right)\textbf{e}^{r}\right)=\nabla f(x)\,.$

It was shown in (Wright, 2015; Richtárik and Takáč, 2011; Nesterov, 2012) that
cost reduction is indeed observed, especially when the objective function $f$
is highly skewed in a high dimensional space. Similar ideas were discussed in
(Gerencsér, 1997; Kleinman et al., 1999) where SPSA, simultaneous perturbation
stochastic approximation, was proposed.

There are many different ways to compute directional derivatives. Automatic
differentiation is a strategy used often when the underlying function is
composed of many simpler operations, and such composition is explicitly known
(Baydin et al., 2018). When the function is complicated itself, one could use
the most basic finite differencing method. That is to use
$\partial_{i}f(x)\approx\frac{f(x+\eta\textbf{e}^{i})-f(x-\eta\textbf{e}^{i})}{2\eta}$
at a sacrifice of $O(\eta^{2})$ numerical error (assuming enough smoothness of
$f$). For special problems such as the PDE-based inverse problems, raising in
atmospheric science and epidemiology (Li et al., 2020), one could further
translate the derivative computation into a combination of one forward and one
adjoint PDE solves (Martin et al., 2012). Whichever strategy one employs, in
some applications, it is true that when $\nabla f$ is not explicitly known,
the computation of the full gradient costs $d$ times of computing one
directional derivative. This makes the saving on $d$ extremely important.

###  2.4 Variance reduction techniques: SVRG and SAGA

SVRG (stochastic variance reduced gradient)(Johnson and Zhang, 2013) and SAGA
(Defazio et al., 2014) (a modified version of SAG (Schmidt et al., 2017),
stochastic average gradient) are optimization techniques that are widely used
in reducing variance in randomized solvers. In particular, they are introduced
to enhance the numerical performance of SGD (stochastic gradient descent). SGD
is a stochastic version of Gradient Descent (GD) that looks for the
minimization of an objective function that has the form of
$f=\sum_{n=1}^{N}f_{n}$ with $N\gg 1$. GD requires $\nabla f=\sum_{n=1}\nabla
f_{n}$ evaluated at the current sample in each iteration, which amounts to a
computation of $N$ gradients. In SGD, one merely randomly selects one $\nabla
f_{r}$ to represent the full $\nabla f$. Per iteration, this reduces the cost
by $N$ folds. However, due to the random selection process, high variance is
induced and that brings large errors. More iterations are then needed to
achieve a preset accuracy. SVRG and SAGA are algorithms proposed to reduce
this variance, hoping to eliminate the error induced by the randomness and
save computational cost overall. We now describe the procedure of SVRG and
SAGA.

In SVRG, one has a preset iteration number $\tau$, and the full gradient
$\nabla f$ is computed only once every $\tau$ steps. Between the epochs, per
step, only one $f_{r}$ is selected at random to represent the new gradient. In
particular, call $\widetilde{x}$ the sample obtained at $k\tau$-th step with
$k$ being an integer, we have the full gradient at this point $\nabla
f(\widetilde{x})$. In the following $\tau-1$ steps, per iteration, one $f_{r}$
is chosen at random and the new gradient is approximated by

$\nabla f(x)\approx\nabla f(\widetilde{x})+N\left(\nabla f_{r}(x)-\nabla
f_{r}(\widetilde{x})\right)\,.$ (7)

After $\tau$ steps, $\widetilde{x}$ is updated and one evaluates the full
$\nabla f(\widetilde{x})$ again. This approximation (7) is unbiased if $r$ is
chosen uniformly in the sense that

$\nabla f(x)=\sum_{n}\nabla f_{n}(x)=\mathbb{E}_{r}\left[\nabla
f(\widetilde{x})+N\left(\nabla f_{r}(x)-\nabla
f_{r}(\widetilde{x})\right)\right]\,.$

SAGA is slightly different: it only requires the computation of the full
gradient at the initial step. In the later iterations, per iteration, only one
$f_{r}$ is chosen at random to update the full gradient while others are kept
untouched. Term $\\{g^{n}_{m}\\}^{N}_{i=1}$ the $m$-th step approximation to
$\nabla f_{n}$, then $g^{n}_{0}=\nabla f_{n}(x^{0})$. For the following
iterations with $m>1$, $r$ is uniformly randomly picked and $g^{r}_{m}=\nabla
f_{r}(x^{m})$ while others are kept untouched $g^{k}_{m}=g^{k}_{m-1}$. We then
approximate the full gradient using:

$\nabla f(x^{m})\approx\sum^{N}_{n=1}g^{n}_{m-1}+N(g^{r}_{m}-g^{r}_{m-1})\,.$
(8)

Similar to SVRG, the approximation (8) is unbiased in the sense that

$\nabla f(x)=\mathbb{E}_{r}\left[\sum^{N}_{n=1}g^{n}_{m-1}+N\left(\nabla
f_{r}(x)-g^{r}_{m-1}\right)\right]\,.$

Both techniques are proven to reduce variance for SGD (Johnson and Zhang,
2013; Defazio et al., 2014).

##  3 Notations and classical results

We unify notions and assumptions, and briefly summarize the classical results
in this section.

###  3.1 Assumptions and the Wasserstein distance

We make standard assumptions on $f(x)$:

###### Assumption 3.1

The function $f$ is second-order differentiable, $\mu$-strongly convex and has
an $L$-Lipschitz gradient:

  * –

Convex, meaning for any $x,x^{\prime}\in\mathbb{R}^{d}$:

$f(x)-f(x^{\prime})-\nabla
f(x^{\prime})^{\top}(x-x^{\prime})\geq(\mu/2)|x-x^{\prime}|^{2}\,.$ (9)

  * –

Gradient is Lipschitz, meaning for any $x,x^{\prime}\in\mathbb{R}^{d}$:

$|\nabla f(x)-\nabla f(x^{\prime})|\leq L|x-x^{\prime}|\,.$ (10)

These assumptions together mean $\mu{I}_{d}\preceq H(f)\preceq L{I}_{d}$ where
$H(f)$ is the hessian of $f$. We also define condition number of $f(x)$ as

$\kappa=L/\mu\geq 1\,.$ (11)

Furthermore, for O-LMC we assume Lipschitz condition of the hessian:

###### Assumption 3.2

Hessian of $f$ is H-Lipschitz, meaning for any
$x,x^{\prime}\in\mathbb{R}^{d}$:

$\|H(f)(x)-H(f)(x^{\prime})\|_{2}\leq H|x-x^{\prime}|\,.$ (12)

Throughout the analysis, we use the Wasserstein distance as the quantity to
measure the distance between two probability measures. For $p\geq 1$, define
the Wasserstein distance to be:

$W_{p}(\mu,\nu)=\left(\inf_{(X,Y)\in\Gamma(\mu,\nu)}\mathbb{E}|X-Y|^{p}\right)^{1/p}\,,$

where $\Gamma(\mu,\nu)$ is the set of distribution of
$(X,Y)\in\mathbb{R}^{2d}$ whose marginal distributions are $\mu$ and $\nu$
respectively for $X$ and $Y$. These distributions are called the couplings of
$\mu$ and $\nu$. Here $\mu$ and $\nu$ can be either probability measures
themselves or the measures induced by probability density functions $\mu$ and
$\nu$. In this paper, we only use $p=2$.

###  3.2 Classical O/U-LMC

We now review the classical results regarding O-LMC and U-LMC. The two
algorithm are summarized in the following:

Algorithm 1 Overdamped/Underdamped Langevin Monte Carlo (O/U-LMC)

Preparation:

1\. Input: $h$ (time stepsize); $\gamma$ (parameter); $d$ (dimension); $M$
(stopping index); $\nabla f(x)$.

2\. Initial: _(overdamped)_ : $x^{0}$ i.i.d. sampled from an initial
distribution induced by $q_{0}(x)$. _(underdamped)_ : $(x^{0},v^{0})$ i.i.d.
sampled from the initial distribution induced by $q_{0}(x,v)$.

Run: For $m=0\,,1\,,\cdots\,,M$ _(overdamped)_ : Draw $\xi^{m}$ from
$\mathcal{N}(0,I_{d})$:  $x^{m+1}=x^{m}-\nabla f(x^{m})h+\sqrt{2h}\xi^{m}\,.$
(13) _(underdamped)_ : Sample $(x^{m+1},v^{m+1})$ as Gaussian random variables
with expectation and covariance defined in (5).

end

Output: $\\{x^{m}\\}$.

The non-asymptotic convergence results have been thoroughly discussed in
(Dalalyan and Karagulyan, 2019; Durmus et al., 2019; Durmus and Moulines,
2019) and (Cheng et al., 2018; Dalalyan and Riou-Durand, 2020) for O-LMC and
U-LMC respectively. We merely cite the theorems.

###### Theorem 1

[(Dalalyan and Karagulyan, 2019) Theorem 5] Assume $h<\frac{2}{\mu+L}$ and $f$
satisfies Assumptions 3.1-3.2. Denote $q^{O}_{m}(x)$ the probability density
function of $x^{m}$ computed using O-LMC, and define
$W_{m}=W_{2}(q^{O}_{m},p)$, the $L_{2}$-Wasserstein distance between
$q^{O}_{m}(x)$ and $p$, then we have:

$W_{m}\leq\exp\left(-\mu
mh\right)W_{0}+\left[\frac{Hhd}{2\mu}+3\kappa^{3/2}\mu^{1/2}hd^{1/2}\right]\,.$
(14)

Note that there are two parts in the iteration formula for $W_{m}$, the
exponentially decaying term, and the remainder term. Without considering the
conditioning constants, the remainder term is of the order of $hd$. To have
the Wasserstein distance to be smaller than a preset tolerance, for example,
$W_{m}\leq\epsilon$, we need

$\exp\left(-\mu mh\right)\lesssim\epsilon\,,\quad\text{and}\quad
hd\lesssim\epsilon\,,$

which explains $m=\widetilde{O}(d/\epsilon)$ in Table 1. Since each iteration
requires the computation of the full gradient, meaning $d$ partial
derivatives, the total cost is then $\widetilde{O}(d^{2}/\epsilon)$. In Durmus
et al. (2019), the authors improved the dependence of the convergence rate on
the condition number, with the Lipschitz continuity of the hessian term
relaxed.

For the underdamped LMC, we have:

###### Theorem 2

[(Dalalyan and Riou-Durand, 2020) Theorem 2] Assume $f$ satisfies Assumption
3.1, $\gamma=\frac{1}{L}$ and $h\leq\frac{1}{8\kappa^{2}\mu}$. Denote
$q^{U}_{m}(x,v)$ the probability density function of $(x^{m},v^{m})$ computed
using U-LMC, and define $W_{m}=W_{2}(q^{U}_{m},p)$, then we have:

$W_{m}\leq\sqrt{2}\exp(-0.375mh/\kappa)W_{0}+h(\kappa d)^{1/2}\,.$ (15)

The iteration formula is once again composed of two terms: the exponential
term and the remainder, with the remainder being of order $hd^{1/2}$. Setting
both terms smaller than $\epsilon$, we have

$h=\widetilde{O}(d^{-1/2}\epsilon)\,,\quad\text{and}\quad
m=\widetilde{O}(d^{1/2}/\epsilon)\,,$

making the total number of partial derivative calculations being
$\widetilde{O}(d^{3/2}/\epsilon)$, as shown in Table 1.

##  4 Algorithms and Results of RCD-LMC

As seen in Algorithm 1, per iteration, $\nabla f$ needs to be evaluated at the
current sample, and that prompts $d$ partial derivative calculations. If $d$
is big, the computation is expensive.

RCD (random coordinate descent) is a method introduced in optimization to
reduce the number of partial derivative calculations from the classical
gradient descent. It essentially surrogates the full gradient in the gradient
descent method by one directional derivative in each iteration, and thus it
naturally reduces the computational cost by $d$ folds per iteration. If we
blindly apply this approach to O/U-LMC, we arrive at RCD-O/U-LMC. These two
methods were briefly discussed in (Ding and Li, 2020). It was also shown there
that the methods, despite easy to implement and reduce the per-iteration-cost
by $d$ folds, the overall cost is, however, not saved due to the larger number
of required iterations. This means blindly applying RCD to O/U-LMC is a bad
strategy. For the completeness of the paper, we still summarize the results
and cite the convergence rate in the following subsections.

We note that the increased numerical cost could be explained by the high
variance induced by the randomness in the coordinate selection process. The
same issue was encountered in optimization as well, see discussions in
(Richtárik and Takáč, 2011; Wright, 2015). As a consequence, a lot of efforts
have been placed on tuning the probability of the coordinate-drawing, such as
incorporating the information hidden in the directional Lipschitz constants.
The application of these strategies should be explored in sampling as well,
but it is beyond the scope of the current paper. We should also emphasize that
the counterexample is a worst-case study. There are plenty of examples where
the application of RCD to sampling already gives a numerical boost, but in
this paper, we would like to be as general as possible.

###  4.1 Algorithm

We apply RCD to both O-LMC and U-LMC. This amounts to replacing the gradient
terms in (3) and (5) using the approximation (6) . The new methods are termed
RCD-O-LMC and RCD-U-LMC respectively, and we present them in Algorithm 2.

Algorithm 2 RCD-O/U-LMC

Preparation:

1\. Input: $h$ (time stepsize); $\gamma$ (parameter); $d$ (dimension); $M$
(stopping index); $f(x)$.

2\. Initial: _(overdamped)_ : $x^{0}$ i.i.d. sampled from an initial
distribution induced by $q_{0}(x)$. _(underdamped)_ : $(x^{0},v^{0})$ i.i.d.
sampled from the initial distribution induced by $q_{0}(x,v)$.

Run: For $m=0\,,1\,,\cdots\,,M$ 1\. Prepare directional derivative: draw $r$
uniformly from $1\,,\cdots\,,d$, and compute:
$F^{m}=d\partial_{r}f(x^{m})\textbf{e}^{r}\,.$ (16) 2\. _(overdamped)_ : Draw
$\xi^{m}$ from $\mathcal{N}(0,I_{d})$:
$x^{m+1}=x^{m}-F^{m}h+\sqrt{2h}\xi^{m}\,.$ (17) _(underdamped)_ : Sample
$(x^{m+1},v^{m+1})$ as Gaussian random variables with expectation and
covariance defined in (5), replacing $\nabla f(x^{m})$ by $F^{m}$.

end

Output: $\\{x^{m}\\}$.

###  4.2 Convergence results for RCD-O/U-LMC

We discuss the convergence of Algorithm 2 in this section and compare the
results with the classical results for O-LMC (Dalalyan and Karagulyan, 2019)
and U-LMC (Cheng et al., 2018).

####  4.2.1 Convergence for RCD-O-LMC

The theorem we show is:

###### Theorem 3

Suppose $f$ satisfies Assumptions 3.1-3.2, and $h$ satisfies

$h<\min\left\\{\frac{1}{9\kappa^{2}\mu
d},\frac{2}{H^{2}/(\kappa\mu^{2})+\kappa^{2}\mu/d}\right\\}\,,$ (18)

then we have:

$W_{m}\leq\exp\left(-\mu mh/4\right)W_{0}+6d(\kappa h)^{1/2}\,,$ (19)

where $W_{m}=W_{2}(q^{O}_{m},p)$ and $q^{O}_{m}(x)$ is the probability density
function of $x^{m}$ computed by RCD-O-LMC.

We discuss the proof in Section 7.2. The statement serves as the guidance to
tune parameters and estimates the computational complexity. According to (19),
for $\epsilon$ accuracy, it suffices to choose

$6\kappa dh^{1/2}\leq\frac{\epsilon}{2},\quad\text{and}\quad\exp(-\mu
hm/4)\leq\frac{\epsilon}{2W_{0}}\,,$

which yields

$h<\frac{\epsilon^{2}}{24\kappa^{2}d^{2}}\,\quad\text{and}\quad mh\geq
4/\mu\log\left(2W_{0}/\epsilon\right)\,.$

For small enough $\epsilon$, this new $h$ restriction is stronger than (18),
and this implies $O(d^{2}/\epsilon^{2}\log(W^{O}_{0}/\epsilon))$ iterations
are needed. Since each iteration requires only one partial derivative
computation, the cost is also $\mathcal{\widetilde{O}}(d^{2}/\epsilon^{2})$.
Compared with the classical O-LMC presented in Section 3.2, this cost is in
fact $1/\epsilon$ higher.

###### Remark 4

While the exponential decay term naturally appears due to the convexity, as
was done in the classical analysis for O-LMC, the second term of (19) is at
the order of $O(h^{1/2}d)$, instead of $O(hd)$ as in O-LMC. This deterioration
comes from the extra error term induced by replacing $\nabla f(x^{m})$ by
$F^{m}$. As a consequence, the core of the proof lies in bounding

$E^{m}=\nabla f(x^{m})-F^{m}\,.$ (20)

As shown in Lemma 17, the variance of this term is:

$\mathbb{E}(|E^{m}|^{2})\lesssim O(d^{2})\,.$

When inserted in the updating formula, this contributes a term of order
$O(h^{2}d^{2})$ and it dominates all other error sources ($O(h^{3}d^{2})$ as
seen in the classical O-LMC). This eventually leads to a worse estimate. See
more detailed discussion in Remark 18.

###### Remark 5

It is a natural question to ask if uniform sampling is necessary. Indeed, in
the optimization literature, RCD outperforms GD when the objective function
$f$ is skewed, and the coordinates are chosen in a non-uniform way to reflect
such skewness. For sampling, we can also use a non-uniform sampling strategy
as well.

To this goal, we first denote the probability of choosing $i$-th direction
$\phi_{i}$ and

$\Phi:=\\{\phi_{1},\phi_{2},\dotsc,\phi_{d}\\}\,,\quad\text{with}\quad\phi_{i}>0\;(\forall
i)\,,\quad\text{and}\quad\sum_{i=1}^{d}\phi_{i}=1\,.$

Then, in the $m$-th iteration of RCD-LMC, we replace (16) by

$F^{m}=\frac{1}{\phi_{r}}\left(\nabla
f(x^{m})\cdot\textbf{e}^{r}\right)\textbf{e}^{r}\,,$ (21)

where $x^{m}$ is the $m$-th iteration sample. This is a valid calculation
because this definition is still consistent with $\nabla f$ in the expectation
sense:

$\mathbb{E}_{r}\left(\frac{1}{\phi_{r}}\left(\nabla
f(x)\cdot\textbf{e}^{r}\right)\textbf{e}^{r}\right)=\nabla f(x)\,.$

However, using a non-uniform sampling strategy will not provide better
convergence than the uniform sampling under our setting. This will be
discussed in Remark 19.

####  4.2.2 Convergence for RCD-U-LMC

In the underdamped setting, we have the following theorem:

###### Theorem 6

Suppose $f$ satisfies Assumption 3.1. Set $\gamma=1/L$. If $h$ satisfies

$h<\frac{1}{880d\kappa}\,,$ (22)

then, we have:

$W_{m}\leq 4\exp\left(-hm/(8\kappa)\right)W_{0}+Cdh^{1/2}\,,$ (23)

where $C=100(\kappa/\mu)^{1/2}$, $W_{m}=W_{2}(q^{U}_{m},p_{2})$, and
$q^{U}_{m}(x,v)$ is the probability density function of $(x^{m},v^{m})$
computed by RCD-U-LMC.

We discuss the proof in Section 8.2. To obtain $\epsilon$ accuracy, we simply
set both terms $<\epsilon/2$ in (23), and it yields, in addition to (22):

$h<\frac{\epsilon^{2}\mu}{10^{4}\kappa d^{2}}\,,\quad\text{and}\quad mh\geq
8\kappa\log\left(8W_{0}/\epsilon\right)\,.$

This leads to
$O\left(d^{2}/\epsilon^{2}\log\left(W_{0}/\epsilon\right)\right)$ number of
iterations, and thus the cost is
$\mathcal{\widetilde{O}}(d^{2}/\epsilon^{2})$. Compared with the classical
U-LMC as discussed in Section 3.2, this cost is $d^{1/2}/\epsilon$ larger.

Similar to the overdamped case, the estimate for RCD-U-LMC is worse than that
for U-LMC due to the newly induced error $E^{m}=\nabla f(x^{m})-F^{m}$. Once
again this term $\mathbb{E}|E^{m}|^{2}$ provides a larger error source than
other terms ($O(h^{3}d)$ as seen in the classical U-LMC), and is, upon the
application of the iteration formula, represented in the second term in (23).
See more details in Remark 25.

###  4.3 A counter-example

The two theoretical results above suggest that no improvement is seen when RCD
is blindly applied to O/U-LMC, if not worse. We cannot argue the sharpness of
the two results, but we do provide a counter-example to demonstrate that
including RCD to O/U-LMC may not bring computational advantage.

The example is given to show RCD-U-LMC. We assume

$q_{0}(x,v)=\frac{1}{(4\pi)^{d/2}}\exp(-|x-\textbf{u}|^{2}/2-|v|^{2}/2)\,,\quad
p_{2}(x,v)=\frac{1}{(2\pi)^{d/2}}\exp(-|x|^{2}/2-|v|^{2}/2)\,,$ (24)

where $\textbf{u}\in\mathbb{R}^{d}$ satisfies $\textbf{u}_{i}=1/8$ for all
$i$. Denote $\\{(x^{m},v^{m})\\}$ the sample computed through Algorithm 2
(underdamped) with stepsize $h$. Denote $q_{m}$ the probability density
function of $(x^{m},v^{m})$, then we can show $W_{2}(q_{m},p_{2})$ cannot
converge too fast.

###### Theorem 7

[Theorem 4.1 in Ding and Li (2020)] For the example above, choose $\gamma=1$
in Algorithm 2, and let $d,h$ satisfy

$d>1872,\quad h<\frac{1}{1440^{2}d}\,,$

then

$W_{m}\geq\exp\left(-2mh\right)\frac{\sqrt{d}}{1024}+\frac{d^{3/2}h}{2304}\,,$
(25)

where $W_{m}=W_{2}(q^{U}_{m},p_{2})$, and $q^{U}_{m}(x,v)$ is the probability
density function of $(x^{m},v^{m})$ computed by RCD-U-LMC.

We note the second term in (25) is rather big. The smallness comes from $h$,
the stepsize, and it needs to be small enough to balance out the influence
from $d^{3/2}\gg 1$. This puts a strong restriction on $h$. Indeed, to have
$\epsilon$-accuracy, $W(q^{U}_{m},p_{2})\leq\epsilon$, we need both terms
smaller than $\epsilon$, and this term suggests that
$h\leq\frac{2304\epsilon}{d^{3/2}}$ at least. And when combined with
restriction from the first term, we arrive at the conclusion that at least
$\widetilde{O}(d^{3/2}/\epsilon)$ iterations are needed, and thus
$\widetilde{O}(d^{3/2}/\epsilon)$ partial derivative calculations are
required. The $d$ dependence is $d^{3/2}$. This is exactly the same as what
U-LMC requires, meaning RCD-U-LMC brings no computational advantage over
U-LMC.

###### Remark 8

All discussions on convergence for U-LMC methods are conducted without
assuming the hessian of $f$ is Lipschitz. A natural question is: If $f$ has a
higher regularity, is the method converging faster? The example above gives a
negative answer. If we apply the classical U-LMC to this example, the number
of iterations to achieve $\epsilon$-accuracy in $W_{2}$-distance is at least
$\widetilde{O}(d^{1/2}/\epsilon)$. This coincides with the results obtained
when only $\nabla f$ is assumed to be Lipschitz. For this reason, throughout
the paper we do not require high regularity of $f$ under the U-LMC framework.

##  5 Variance reduction techniques

Blindly applying RCD to LMC does not bring any numerical savings. In the
analysis, the biggest contribution to the large error (the remainder term) is
induced by the process of randomly selecting directional derivatives. In this
section we propose two variance reduction methods inspired by SVRG (Johnson
and Zhang, 2013) and SAGA (Defazio et al., 2014) and apply them to RCD-O/U-LMC
(and thus four algorithms in total). We will prove that while the numerical
cost per iteration is reduced by $d$-folds, the number of iterations needed
for achieving the preset accuracy is mostly unchanged compared to the
classical LMC. This ultimately saves the total cost.

###  5.1 Algorithms

####  5.1.1 SVRG based variance reduction method

As presented in Section 2.4, SVRG is an optimization technique introduced to
reduce variance in SGD in its random selection process. Only one
representative gradient is computed per iteration, unless the iteration number
is an integer times $\tau$, when all gradients are evaluated. It has been
proved in (Johnson and Zhang, 2013) that the approach enhances the SGD
performance.

Inspired by this method, we propose SVRG-O/U-LMC, as summarized in Algorithm
3. In this algorithm, we compute the full gradient $\nabla f$ at $k\tau$ time-
step for all integers $k$, and during the epochs, per iteration, we merely
update one directional derivative chosen at random.

Algorithm 3 SVRG-O(U)-LMC

Preparation:

1\. Input: $h$ (time stepsize); $\tau$ (epoch length); $\gamma$ (parameter);
$M$ (stopping index); $d$ (dimension) and $f(x)$.

2\. Initial: _(overdamped)_ : $x^{0}$ i.i.d. sampled from an initial
distribution induced by $q_{0}(x)$. _(underdamped)_ : $(x^{0},v^{0})$ i.i.d.
sampled from the initial distribution induced by $q_{0}(x,v)$.

Run: For $m=0\,,1\,,\cdots\,,M$ if $m$ mod $\tau=0$ then update $\widehat{g}$
and compute flux $F^{m}$:
$F^{m}_{i}=\widehat{g}_{i}=\partial_{i}f(x^{0}),\quad\ 1\leq i\leq d\,.$ (26)
otherwise Draw a random number $r^{m}$ uniformly from $1,2,\cdots,d$ and
compute  $g^{m}_{r^{m}}=\partial_{r^{m}}f(x^{m})\,,\quad
F^{m}=\widehat{g}+d\left(g^{m}_{r^{m}}-\widehat{g}_{r^{m}}\right)\textbf{e}^{r^{m}}\,.$
(27) end if _(overdamped)_ : Draw $\xi^{m}$ from $\mathcal{N}(0,I_{d})$:
$x^{m+1}=x^{m}-F^{m}h+\sqrt{2h}\xi^{m}\,.$ (28) _(underdamped)_ : Sample
$(x^{m+1},v^{m+1})$ as Gaussian random variables with expectation and
covariance defined in (5), replacing $\nabla f(x^{m})$ by $F^{m}$.

end

Output: $\\{x^{m}\\}$.

####  5.1.2 SAGA based variance reduction method

SAGA is also a technique introduced in optimization to reduce the variance of
SGD. It starts with a full gradient, and in each step, simply updates the
gradient of one randomly selected representative function.

Inspired by the approach, we propose our Algorithm 4, in which we compute the
full gradient only in the first round of iteration, and keep updating one
randomly selected directional derivative per iteration. We term the method
Random Coordinate Averaging Descent-O/U-LMC (RCAD-O/U-LMC).

Algorithm 4 RCAD-O(U)-LMC

Preparation:

1\. Input: $h$ (time stepsize); $\gamma$ (parameter); $M$ (stopping index);
$d$ (dimension) and $f(x)$.

2\. Initial: _(overdamped)_ : $x^{0}$ i.i.d. sampled from the initial
distribution induced by $q_{0}(x)$ and calculate $g^{0}\in\mathbb{R}^{d}$:
$g^{0}_{i}=\partial_{i}f(x^{0}),\quad\ 1\leq i\leq d\,.$ (29) _(underdamped)_
: $(x^{0},v^{0})$ i.i.d. sampled from the initial distribution induced by
$q_{0}(x,v)$ and calculate $g^{0}\in\mathbb{R}^{d}$ as (29).

Run: For $m=0\,,1\,,\cdots\,,M$ 1\. Draw a random number $r^{m}$ uniformly
from $1,2,\cdots,d$. 2\. Calculate $g^{m+1}$ and flux $F^{m}\in\mathbb{R}^{d}$
by letting $g^{m+1}_{i}=g^{m}_{i}$ for $i\neq r^{m}$ and
$g^{m+1}_{r^{m}}=\partial_{r^{m}}f(x^{m})\,,\quad
F^{m}=g^{m}+d\left(g^{m+1}_{r^{m}}-g^{m}_{r^{m}}\right)\textbf{e}^{r^{m}}\,.$
(30) 3\. _(overdamped)_ : Draw $\xi^{m}$ from $\mathcal{N}(0,I_{d})$:
$x^{m+1}=x^{m}-F^{m}h+\sqrt{2h}\xi^{m}\,.$ (31) _(underdamped)_ : Sample
$(x^{m+1},v^{m+1})$ as a Gaussian random variable with expectation and
covariance defined in (5), replacing $\nabla f(x^{m})$ by $F^{m}$.

end

Output: $\\{x^{m}\\}$.

###  5.2 Convergence and numerical cost analysis

We now discuss the convergence of SVRG/RCAD-O-LMC and SVRG/RCAD-U-LMC. In the
classical papers (Durmus et al., 2019; Dalalyan and Karagulyan, 2019; Cheng et
al., 2018; Dalalyan and Riou-Durand, 2020) that discuss the convergence of
O/U-LMC, the authors indeed discussed the numerical error in approximating the
gradients, but both papers require that the variance of the error is bounded
and independent of $d$. The random coordinate selection process is rather
complicated, making their results hard to apply. One related work is
(Chatterji et al., 2018), where the authors construct a Lyapunov function to
study the convergence of SG-MCMC. Our proof for the convergence of SVRG/RCAD-
O-LMC is inspired by its technicalities. In (Cheng et al., 2018; Chatterji et
al., 2018), a contraction map is used for U-LMC, but such a map cannot be
directly applied in our situation because the variance depends on the entire
trajectory of samples. Furthermore, the history of the trajectory is reflected
in each iteration, deeming the process to be non-Markovian. We need to re-
engineer the iteration formula accordingly for tracing the error propagation.

####  5.2.1 Convergence for SVRG/RCAD-O-LMC

For SVRG-O-LMC, We have the following theorem.

###### Theorem 9

Suppose $f$ satisfies Assumptions 3.1-3.2 and $h$ satisfies

$h<\min\left\\{\frac{1}{400d\kappa^{2}\mu},\frac{1}{10\tau\max\\{\mu,1\\}}\right\\}\,.$
(32)

Denote $W_{m}=W_{2}(q^{O}_{m},p)$, $q^{O}_{m}$ the density of the sample
$x^{m}$ computed from Algorithm 3 (overdamped), and $p$ the target density
function, then we have

$W_{m}\leq\exp\left(-\frac{\mu hm}{32}\right)W_{0}+\left[h^{3/2}\tau
dC_{1}+h\tau^{1/2}dC_{2}\right]\,,$ (33)

where

$C_{1}=30\kappa^{3/2}\mu,\quad
C_{2}=50\kappa\sqrt{\mu}+5\sqrt{\kappa^{3}\mu/d+2H^{2}/\mu^{2}}\,.$

We leave the proof to Section 7.3 and simply discuss the numerical cost here.
Suppose we set $\tau=d$. And we preset the desired accuracy to be $\epsilon$,
meaning we wish to obtain $W_{2}(q^{O}_{m},p)\leq\epsilon$, then we can simply
choose to set all terms in (33) less than $\epsilon/3$. The latter two terms
give the constraints to $h$ while the first one gives the constraint on $m$,
meaning in addition to the requirement (32), we have:

$h\lesssim\min\left\\{\frac{\epsilon^{2/3}}{d^{4/3}C^{2/3}_{1}}\,,\frac{\epsilon}{d^{3/2}C_{2}}\right\\}\,,\quad\text{and}\quad
mh\gtrsim\frac{1}{\mu}\log\left(\frac{W_{0}}{\epsilon}\right)\,.$

For small $\epsilon$ and large $d$, $h<\frac{\epsilon}{d^{3/2}C_{2}}$ is the
most restrictive one, and it leads the cost
$m>\widetilde{O}(d^{3/2}/\epsilon)$. Since in most iterations (unless $k\tau$)
one only calculates one directional derivative, the total computational cost
is $\widetilde{O}(d^{3/2}/\epsilon)$.

For RCAD-O-LMC, We have the following theorem.

###### Theorem 10

[Theorem 5.1 in Ding and Li (2020)] Suppose $f$ satisfies Assumptions 3.1-3.2
and $h$ satisfies

$h<\frac{1}{3(1+9d)\kappa^{2}\mu}\,.$ (34)

Denote $W_{m}=W_{2}(q^{O}_{m},p)$, $q^{O}_{m}$ the density of the sample
$x^{m}$ computed from Algorithm 4 (overdamped), and $p$ the target density
function, then we have

$W_{m}\leq\sqrt{2}\exp(-\mu hm/4)W_{0}+2h\sqrt{d^{3}C_{1}+d^{2}C_{2}}\,.$ (35)

Here $C_{1}=77\kappa^{2}\mu$, $C_{2}=H^{2}/\mu^{2}+\kappa^{3}\mu/d$.

We omit the detailed proof for Theorem 10. Interested readers are referred to
the appendix of (Ding and Li, 2020). The sketch of the proof is presented in
Section 7.4. We here discuss the computational cost. Suppose we present the
desired accuracy to be $\epsilon$, then we can simply choose $h$ to set both
terms in (35) less than $\epsilon/2$, and it leads to, in addition to (34)

$h\lesssim\frac{\epsilon}{d^{3/2}\sqrt{C_{1}+C_{2}/d}}\,,\quad\text{and}\quad
mh\gtrsim\frac{1}{\mu}\log\left(\frac{W_{0}}{\epsilon}\right)\,.$

For small $\epsilon$, the new constraint of $h$ dominates, and this leads to
the cost $m>\widetilde{O}(d^{3/2}/\epsilon)$.

We emphasize that in Theorem 9 and Theorem 10 we require both Assumptions 3.1
and 3.2. The second assumption is on the continuity of the hessian term. If
this is relaxed, the convergence still holds true, with a degraded numerical
cost. The cost will be
$\widetilde{O}(\max\\{d^{3/2}/\epsilon,d/\epsilon^{2}\\})$. But compared to
$\widetilde{O}(d^{2}/\epsilon^{2})$ required by the classical O-LMC, the
algorithms with reduced variance still outperform. The proof is the same, and
we omit it from the paper.

###### Remark 11

We have two comments:

  * •

Compare Theorem 3, 9 and 10, the three theorems for overdamped LMC methods, it
is rather clear that the dependence on $h$ and $m$ in the exponential is the
same. Since $m$ and $h$ have the same order, then for $\epsilon$ accuracy,
$m\sim\frac{1}{h}$. The $h$ requirement, however, mainly comes from the form
of the second term. It is the quality of this term that determines the final
numerical cost.

In particular, in Theorem 3, $h$’s dependence on $d$ is at the power of $-2$,
but in Theorem 9 and 10, such dependence reduces to $-3/2$. This explains the
improvement from RCD to SVRG and RCAD in the overdamped setting, as presented
in Table 1.

  * •

As presented in Remark 4, the major error term that needs to be improved is
the second term in (19), which is mainly induced by the variance term
$\mathbb{E}|E^{m}|^{2}$ that measures the distance between $\nabla f(x^{m})$
and $F^{m}$. When a variance reduction technique is integrated, we can find a
better estimate for $\nabla f(x^{m})$, reducing the error
$\mathbb{E}|E^{m}|^{2}$, and improving the final convergence bound. In
particular, as will be shown in Lemma 20 and Lemma 21, we roughly obtain, for
SVRG-O-LMC and RCAD-O-LMC respectively,

$\mathbb{E}|E^{m}|^{2}\leq
O(h^{2}d^{2}\tau^{2}+hd^{2}\tau)\quad\text{and}\quad\mathbb{E}\left|E^{m}\right|^{2}\leq
O(h^{2}d^{4}+hd^{3})\,.$

Under the condition $\tau=d$ and $h\lesssim\frac{1}{d}$, both terms, when
inserted in the iteration, contribute an error term of order $O(h^{3}d^{3})$.
This is an improvement over RCD-O-LMC $(O(h^{2}d^{2}))$. More details can be
found at the end of Section 7.3 and 7.4.

####  5.2.2 Convergence for SVRG/RCAD-U-LMC

For SVRG-U-LMC, We have the following theorem.

###### Theorem 12

Suppose $f$ satisfies Assumptions 3.1, $\gamma=1/L$, and $h$ satisfies

$h\leq\min\left\\{\frac{1}{1648\kappa d},\frac{1}{40\tau}\right\\}\,,$ (36)

then, with $C_{1}=200\sqrt{\frac{\kappa}{\mu}}$ and $C_{2}=240/\sqrt{\mu}$:

$W_{m}<4\exp\left(-\frac{hm}{32\kappa}\right)W_{0}+hd^{1/2}C_{1}+h^{3/2}\tau
dC_{2}\,.$ (37)

Here $W_{m}=W_{2}(q^{U}_{m},p_{2})$, and $q^{U}_{m}$ denotes the density of
the sample $(x^{m},v^{m})$ derived from Algorithm 3 (underdamped), and $p_{2}$
is the target density function.

We leave the proof to Section 8.3. The theorem gives us the strategy of
choosing parameters: to achieve $\epsilon$-accuracy, meaning to have
$W_{m}\leq\epsilon$, we can set all terms in (37) less than $\epsilon/3$, and
it leads to, when $\tau=d$, in addition to (36)

$h\lesssim\min\left\\{\frac{\epsilon}{d^{1/2}C_{1}},\frac{\epsilon^{2/3}}{d^{4/3}C_{2}^{2/3}}\right\\},\quad\text{and}\quad
mh\gtrsim\kappa\log\left(\frac{4W_{0}}{\epsilon}\right)\,.$

It is hard to balance the smallness of $\epsilon$ and the largeness of $d$,
and we keep both restrictions, this means
$m>\widetilde{O}\left(\max\left\\{\frac{d^{4/3}}{\epsilon^{2/3}},\frac{d^{1/2}}{\epsilon}\right\\}\right)$.

For RCAD-U-LMC, we have the following theorem:

###### Theorem 13

[Theorem 5.2 in Ding and Li (2020)] Assume $f(x)$ satisfies Assumption 3.1,
$\gamma=1/L$, and $h$ satisfies

$h\leq\frac{1}{1648\kappa d}\,,$ (38)

then, with $C_{1}=200\sqrt{\frac{\kappa}{\mu}}$ and $C_{2}=200/\sqrt{\mu}$:

$\displaystyle W_{m}\leq
4\sqrt{2}\exp\left(-\frac{hm}{8\kappa}\right)W_{0}+hd^{1/2}C_{1}+h^{3/2}d^{2}C_{2}\,,$
(39)

where $W_{m}=W_{2}(q^{U}_{m},p_{2})$ and $q^{U}_{m}$ is the density of the
sample $(x^{m},v^{m})$ derived from Algorithm 4 (underdamped).

The proof of this theorem is referred to (Ding and Li, 2020), and we sketch
the main strategy in Section 8.4. Here we only discuss computational cost. To
achieve $\epsilon$-accuracy, meaning to have $W_{m}\leq\epsilon$, we can
choose all terms in (39) less than $\epsilon/3$. This gives, in addition to
(38).

$h\lesssim\min\left\\{\frac{\epsilon}{d^{1/2}C_{1}},\frac{\epsilon^{2/3}}{d^{4/3}C_{2}^{2/3}}\right\\},\quad\text{and}\quad
mh\gtrsim\kappa\log\left(\frac{4\sqrt{2}W_{0}}{\epsilon}\right)\,.$

For small $\epsilon$ and large $d$, the second term in the $h$ bound is
smaller than the third term, and this means
$m>\widetilde{O}\left(\max\left\\{\frac{d^{4/3}}{\epsilon^{2/3}},\frac{d^{1/2}}{\epsilon}\right\\}\right)$.

###### Remark 14

The improvement observed in the underdamped case is an analogy to that seen
for the overdamped case. By integrating variance reduction techniques, we have
a better estimate of $\nabla f(x^{m})$, and this reduces the cost of
$\mathbb{E}|E^{m}|^{2}$. The control of this term is presented in Lemma 26 and
Lemma 27 respectively. More detailed discussion can be found at the end of
Section 8.3 and 8.4.

##  6 Numerical result

In this section, we test the efficiency of variance reduction enhanced RCD
methods using the following three examples. We should mention that it is
numerically challenging to compute the Wasserstein distance on a high
dimensional space, so we present the convergence of error only in the weak
sense, by testing it on a test function.

  * •

Example 1: In this example, our target distribution is $\mathcal{N}(0,I_{d})$
with $d=1000$, meaning

$p(x)\propto\exp\left(-\sum^{d}_{i=1}\frac{|x_{i}|^{2}}{2}\right)\,.$

The initial distributions, for the overdamped and underdamped situations
respectively, are $\mathcal{N}(\textbf{0.5},I_{d})$ and
$\mathcal{N}(\textbf{0.5},I_{2d})$. We run RCD-O/U-LMC, RCAD-O/U-LMC and
SVRG-O/U-LMC with $N=5\times 10^{5}$ particles with different stepsizes $h$,
and test expectation estimation error by calculating

$\mathrm{Error}=\left|\frac{1}{N}\sum^{N}_{i=1}\phi(x^{i})-\mathbb{E}_{p}(\phi)\right|\,,$
(40)

where $\phi$ is the test function. In Figure 1, we show the error with
$\phi(x)=|x_{1}|^{2}$ using different stepsizes. In both underdamped and
overdamped cases, the improvement of convergence rate is immediate when
variance reduction techniques are added.

(a)  (b)

Figure 1: Decay of Error of LMC in overdamped (left) and underdamped (right)
settings.

  * •

Example 2: In this example, our target distribution is the summation of two
Gaussians whose density is

$p(x)\propto\exp\left(-\sum^{d}_{i=1}\frac{|x_{i}-2|^{2}}{2}\right)+\exp\left(-\sum^{d}_{i=1}\frac{|x_{i}+2|^{2}}{2}\right).$

This target distribution has two peaks and does not satisfy the assumptions in
this paper. However, the experimental results still suggest the outperformance
when variance reduction is imposed. We choose $d=1000$ and $N=10^{6}$
particles. The initial distributions, for the overdamped and underdamped
situations respectively, are $\mathcal{N}(\textbf{0},I_{d})$ and
$\mathcal{N}(\textbf{0},I_{2d})$. We run RCD-O/U-LMC, RCAD-O/U-LMC,
SVRG-O/U-LMC with different $h$, and calculate error using (40) with
$\phi(x)=|x_{1}|^{2}$. In both underdamped and overdamped settings, including
variance reduction techniques provide better accuracy, as shown in Figure 2.

(a)  (b)

Figure 2: Decay of Error of LMC in overdamped (left) and underdamped (right)
settings.

  * •

Example 3 (Generalized linear regression) : In this example, we consider the
following linear regression problem:

$b=x\cdot a+\eta\,,$ (41)

where $a,x\in\mathbb{R}^{d}$ and $\eta$, the noise is a real-valued random
variable satisfying $p_{\text{nos}}(\eta)$. Suppose $p_{\text{nos}}(\eta)$ is
known, and one is given a set of data
$\left\\{(a_{i},b_{i})\right\\}^{I}_{i=1}$ to infer $x$. According to the
Bayes’ rule, given the prior distribution of $x$, termed $p_{\text{prior}}$,
the posterior distribution, also our target distribution is:

$p(x)=p_{\text{pos}}(x)\propto
p_{\text{prior}}(x)\Pi^{I}_{i=1}p_{\text{nos}}(b_{i}-x\cdot a_{i})\,.$ (42)

Now we look for samples sampled according to $p_{\text{pos}}$. In the test, we
choose $d=100$, $x=\textbf{1}$ and generate $I=100$ data pairs. We further set
the prior to be $\mathcal{N}(\textbf{0},I_{d})$, and we evaluate the
algorithms with two different choices of $p_{\text{nos}}$. In the first
choice, this noise distribution is $\mathcal{N}(0,1)$, and in the second
choice, we set:

$p_{\text{nos}}(\eta)\propto\exp\left(-\frac{|\eta|^{2}+\cos(\eta)}{2}\right)\,.$
(43)

All six algorithms, RCD-O/U-LMC, RCAD-O/U-LMC, SVRG-O/U-LMC are run with
$N=10^{6}$. The error is calculated using (40) with
$\phi(x)=|\textsf{x}|^{2}$, where $\textsf{x}=(x_{1},x_{2},\dots,x_{10})$ is
the first ten components of $x$.

In Figure 3 and 4, $\eta$ satisfies standard Gaussian and (43) respectively.
In both cases, the variance reduction methods converge faster.

(a)  (b)

Figure 3: Decay of Error of LMC in overdamped (left) and underdamped (right)
settings.

(a)  (b)

Figure 4: Decay of Error of LMC in overdamped (left) and underdamped (right)
settings.

In all the numerical examples above, we observe, in the O-LMC framework, the
classical O-LMC gives the error saturating at $\sim h$ while the variance
reduced algorithms roughly gives saturation error $\sim h^{2}$. In the U-LMC
framework, the dependence of such error on $h$ increases from $1$ to $\sim
2.5$ when variance reduction techniques are incorporated. We should stress
that we numerically observe better convergence rates than those stated in the
theorems. This is potentially due to the fact that we are evaluating the error
in the weak sense by testing the samples on a test function. This is a weaker
criterion than the Wasserstein distance discussed in the theorems.

##  7 Proof of convergence results for O-LMC

We analyze the three O-LMC methods (RCD-O-LMC, SVRG-O-LMC, and RCAD-O-LMC)
together in this section.

Before diving into details, we quickly summarize the proving strategy. Recall
that the target distribution $p$ is merely the equilibrium of the SDE (1).
This means, if a particle prepared at the initial stage is drawn from $p$,
then following the dynamics of SDE (1), the distribution of this particle will
continue to be $p$. In the analysis below, we call the trajectory of this
particle $y_{t}$, and the sequence generated by this particle evaluated at
discrete time $y^{m}$. Essentially we evaluate how quickly $x^{m}$ converges
to $y^{m}$ as $m$ increases. In particular, we call $\Delta^{m}=x^{m}-y^{m}$
and will derive an iteration formula that shows the convergence of
$\Delta^{m}$.

In evaluating $\Delta^{m}$, two kinds of error get involved:

  * 1.

discretization error in $h$: this amounts to controlling the discretization
error of the SDE (1). To handle this part of the error we employ the estimates
in (Dalalyan and Karagulyan, 2019);

  * 2.

random coordinate selection process error: this is to measure, at each
iteration, how big can $\nabla f(x^{m})-F^{m}$ be. According to the way
$F^{m}$ is defined, it is straightforward to show that the expectation of this
error is always $0$, but the variance $\mathbb{E}|\nabla f(x^{m})-F^{m}|^{2}$
can be big. The analysis for the three different methods under O-LMC framework
mostly concentrates on giving a bound to this term.

In a way, see details in (54), we can derive the iteration formula:

$\displaystyle\Delta^{m+1}=$ $\displaystyle
y^{m+1}-x^{m+1}=\Delta^{m}+(y^{m+1}-y^{m})-(x^{m+1}-x^{m})$ $\displaystyle=$
$\displaystyle\Delta^{m}-h\left(\nabla f(y^{m})-\nabla
f(x^{m})\right)-\int^{(m+1)h}_{mh}\left(\nabla f(y_{s})-\nabla
f(y^{m})\right)\,\mathrm{d}s$ $\displaystyle-h(\nabla f(x^{m})-F^{m})\,.$

The second term on the right-hand side, by using the Lipschitz continuity,
will provide $-Lh\Delta^{m}$, and it produces desirable properties when
combined with the first $\Delta^{m}$. The third term encodes the
discretization error in $h$ and was shown to be small in (Dalalyan and
Karagulyan, 2019). The last term is the error that comes from the random
coordinate selection process, and different methods will give different
control.

It turns out that all three O-LMC methods share the same iterative formula for
$\Delta^{m}$. We will derive this formula in Theorem 15. The formula, however,
requires the boundedness of a variance term $\mathbb{E}|\nabla f-F|^{2}$, and
all three methods will provide different bounds. In subsection 7.1, we discuss
the background and derive the iteration formula. The subsequent three
subsections are dedicated to the analysis of the variance term. When combined
with Theorem 15, the theorems are immediate.

###  7.1 Iterative formula for $\Delta^{m}$

We first rewrite (17), (28), (31) into

$x^{m+1}=x^{m}-\nabla f(x^{m})h+E^{m}h+\sqrt{2h}\xi^{m}\,,$ (44)

where $E^{m}$ is defined in (20):

$E^{m}=\nabla f(x^{m})-F^{m}\,.$

It can be shown that in all three cases, if $r^{m}$ is chosen uniformly from
$1\,,\cdots\,,d$:

$\mathbb{E}_{r^{m}}(E^{m})=\textbf{0}\,.$ (45)

We simply need to recall the definition of $F^{m}$. In RCD, we have:

$F^{m}=d\partial_{r^{m}}f(x^{m})\textbf{e}^{r^{m}}\,,$ (46)

and for SVRG-O-LMC it becomes:

$F^{m}=\widehat{g}+d\left(g^{m}_{r^{m}}-\widehat{g}_{r^{m}}\right)\textbf{e}^{r^{m}}\,,\quad\text{with}\quad
g^{m}_{r^{m}}=\partial_{r^{m}}f(x^{m}),\ \widehat{g}=\nabla
f\left(x^{\left\lfloor\frac{m}{\tau}\right\rfloor}\right)\,,$ (47)

and for RCAD-O-LMC, it is:

$F^{m}=g^{m}+d\left(g^{m+1}_{r^{m}}-g^{m}_{r^{m}}\right)\textbf{e}^{r^{m}}\,,\quad\text{with}\quad
g^{m+1}_{r^{m}}=\partial_{r^{m}}f(x^{m})\,.$ (48)

Taking the expectation with respect to $r^{m}$, (45) is straightforward.

We then define the dynamics of $y_{t}$ to be:

$y_{t}=y^{0}-\int^{t}_{0}\nabla
f(y_{s})\,\mathrm{d}s+\sqrt{2}\int^{t}_{0}\,\mathrm{d}B_{s}\,,$ (49)

then, calling $y^{m}=y_{hm}$, and letting $B_{h(m+1)}-B_{hm}=\sqrt{h}\xi^{m}$,
we have:

$y^{m+1}=y^{m}-\int^{(m+1)h}_{mh}\nabla
f(y_{s})\,\mathrm{d}s+\sqrt{2h}\xi^{m}\,.$ (50)

Assume $y_{0}$ is a random vector drawn from the target distribution induced
by $p$ such that $W^{2}_{2}(q^{O}_{0},p)=\mathbb{E}|x^{0}-y^{0}|^{2}$, we have
$y_{t}$ is drawn from the distribution induced by $p$ for all $t$.

Denote

$\Delta^{m}=y^{m}-x^{m}\,,$ (51)

then

$W^{2}_{2}(q^{O}_{m},p)\leq\mathbb{E}|\Delta^{m}|^{2}=\mathbb{E}|x^{m}-y^{m}|^{2}\,,$

where $\mathbb{E}$ takes all randomness into account. We now essentially need
to give a bound to $\Delta^{m}$ by comparing (44) and (50).

In Theorem 15 we derive the iteration formula for updating $\Delta^{m}$.

###### Theorem 15

Suppose $f$ satisfies Assumptions 3.1-3.2, if $\\{x^{m}\\}$ is defined in
(44), $\\{y^{m}\\}$ is defined in (50) and $\\{\Delta^{m}\\}$ comes from (51),
then for any $a>0$ and $m\geq 0$, we have

$\mathbb{E}|\Delta^{m+1}|^{2}\leq(1+a)A\mathbb{E}|\Delta^{m}|^{2}+3(1+a)h^{2}\mathbb{E}|E^{m}|^{2}+(1+a)h^{3}B+\left(1+\frac{1}{a}\right)h^{4}C\,,$
(52)

where

$\displaystyle A=1-2\mu h+3L^{2}h^{2}\,,\quad B=2L^{2}d\,,\quad
C=\frac{1}{2}\left(H^{2}d^{2}+L^{3}d\right)\,.$ (53)

The only formula this lemma requires is (44). There is no assumption added on
$E^{m}$ besides $\mathbb{E}_{r^{m}}E^{m}=0$, and thus this result is
applicable to all three algorithms. However, as discussed before, the variance
of $E^{m}$ of these three methods are significantly different. The following
three subsections are then dedicated to analyzing this variance term in each
different method. These, together with Theorem 15, finally conclude the three
theorems for overdamped LMC.

###### Remark 16

We carefully observe the relations of the four terms in (52).

  * •

The coefficient in front of $\mathbb{E}|\Delta^{m}|^{2}$ is $(1+a)A$ where
$A\sim 1-h$. Suppose $a$ is chosen to be small enough to make $(1+a)A<1$
strictly, then we have the hope of having exponential decay. This means $a$
should be of the order of $O(h)$.

  * •

Then we need to handle the remaining three terms. The last two terms come from
translating the SDE (1) to O-LMC. This discretization error cannot be removed,
and they are of $O(h^{3}d^{2})$ (since $1/a\sim 1/h$ according to the previous
comment). In the ideal case, we would like to have $\mathbb{E}|E^{m}|^{2}$
(the variance) term to be negligible. If so, upon the iteration formula, one
$h$ drops out, and with taking the square root, the final version of the
formula should be in line with $W_{m}\leq\exp\left(-hm\right)W_{0}+dh$. This
would lead to the estimate of $m=O(\frac{d}{\epsilon})$ for an $\epsilon$
accuracy. In the overdamped setting, this would be the best possible sampling
scheme one could hope for when we assume the variance term does not contribute
much error.

Proof  We first divide $\Delta^{m+1}$ into several parts:

$\displaystyle\Delta^{m+1}=$
$\displaystyle\Delta^{m}+(y^{m+1}-y^{m})-(x^{m+1}-x^{m})$ (54)
$\displaystyle=$ $\displaystyle\Delta^{m}+\left(-\int^{(m+1)h}_{mh}\nabla
f(y_{s})\,\mathrm{d}s+\sqrt{2h}\xi_{m}\right)-\left(-\int^{(m+1)h}_{mh}F^{m}\,\mathrm{d}s+\sqrt{2h}\xi_{m}\right)$
$\displaystyle=$ $\displaystyle\Delta^{m}-\left(\int^{(m+1)h}_{mh}\left(\nabla
f(y_{s})-F^{m}\right)\,\mathrm{d}s\right)$ $\displaystyle=$
$\displaystyle\Delta^{m}-\left(\int^{(m+1)h}_{mh}\left(\nabla f(y_{s})-\nabla
f(y^{m})+\nabla f(y^{m})-\nabla f(x^{m})+\nabla
f(x^{m})-F^{m}\right)\,\mathrm{d}s\right)$ $\displaystyle=$
$\displaystyle\Delta^{m}-h\left(\nabla f(y^{m})-\nabla
f(x^{m})\right)-\int^{(m+1)h}_{mh}\left(\nabla f(y_{s})-\nabla
f(y^{m})\right)\,\mathrm{d}s$ $\displaystyle-h(\nabla f(x^{m})-F^{m})\,.$

Note that here $\nabla f(y^{m})-\nabla f(x^{m})$ can be bounded by
$\Delta^{m}$ using the Lipschitz condition, the term $\int\left(\nabla
f(y_{s})-\nabla f(y^{m})\right)\,\mathrm{d}s$ reflects the numerical error
accumulated in one-step of SDE computation. $\nabla f(x^{m})-F^{m}$ is the
$E^{m}$ term.

Denote

$\displaystyle U^{m}$ $\displaystyle=\nabla f(y^{m})-\nabla f(x^{m})\,,$
$\displaystyle V^{m}$ $\displaystyle=\int^{(m+1)h}_{mh}\left(\nabla
f(y_{s})-\nabla
f(y^{m})-\sqrt{2}\int^{s}_{mh}\nabla^{2}f(y_{r})\,\mathrm{d}B_{r}\right)\,\mathrm{d}s\,,$
$\displaystyle\Phi^{m}$
$\displaystyle=\frac{\sqrt{2}}{h}\int^{(m+1)h}_{mh}\int^{s}_{mh}\nabla^{2}f(y_{r})\,\mathrm{d}B_{r}\,\mathrm{d}s\,,$

we now have

$\displaystyle\Delta^{m+1}=$
$\displaystyle\Delta^{m}-hU^{m}-(V^{m}+h\Phi^{m})-hE^{m}$ (55)
$\displaystyle=$ $\displaystyle\Delta^{m}-V^{m}-h(U^{m}+\Phi^{m}+E^{m})\,.$

Upon getting equation (55) it is time to analyze each term and hopefully
derive an induction inequality that states
$\mathbb{E}|\Delta^{m+1}|^{2}\approx c\mathbb{E}|\Delta^{m}|^{2}+d$ with $c<1$
and $d$ being of high orders of $h$, the parameters we can tune. According to
Lemma 6 of (Dalalyan and Karagulyan, 2019), we first have

$\mathbb{E}|V^{m}|^{2}\leq\frac{h^{4}}{2}\left(H^{2}d^{2}+L^{3}d\right),\quad\mathbb{E}|\Phi^{m}|^{2}\leq\frac{2L^{2}hd}{3}\,,$
(56)

and thus we have:

$\displaystyle\mathbb{E}|U^{m}+\Phi^{m}+E^{m}|^{2}$ (57) $\displaystyle=$
$\displaystyle
3\mathbb{E}|U^{m}|^{2}+3\mathbb{E}|\Phi^{m}|^{2}+3\mathbb{E}|E^{m}|^{2}\,,$
$\displaystyle\leq$ $\displaystyle
3L^{2}\mathbb{E}|\Delta^{m}|^{2}+2L^{2}hd+3\mathbb{E}|E^{m}|^{2}$

where we used the Lipschitz continuity of $f$ for controlling $U^{m}$, and
(56) for $\Phi^{m}$.

We then handle the cross terms. For example, due to the independence, (45),
and the convexity, we have:

$\mathbb{E}\left\langle\Delta^{m},\Phi^{m}\right\rangle=0\,,\quad\mathbb{E}\left\langle\Delta^{m},E^{m}\right\rangle=0\,,\quad\left\langle\Delta^{m},U^{m}\right\rangle\geq\mu|\Delta^{m}|^{2}\,,$
(58)

this means the cross term between first and the third term in the last line
(55) leads to $-2Mh\mathbb{E}|\Delta^{m}|^{2}$. The cross term produced by the
first and the second term, however, can be hard to control, mostly because
$\mathbb{E}\left\langle\Delta^{m},V^{m}\right\rangle$ is unknown. We now
employ Young’s inequality, meaning, for any $a>0$:

$\displaystyle\mathbb{E}|\Delta^{m+1}|^{2}\leq(1+a)\mathbb{E}|\Delta^{m+1}+V^{m}|^{2}+\left(1+\frac{1}{a}\right)\mathbb{E}|V^{m}|^{2}\,.$
(59)

While the second term is already bounded in (56), the first term of (59),
according to (54) becomes:

$\begin{aligned}
\mathbb{E}|\Delta^{m+1}+V^{m}|^{2}=&\mathbb{E}|\Delta^{m}-h(U^{m}+\Phi^{m}+E^{m})|^{2}\\\
=&\mathbb{E}|\Delta^{m}|^{2}-2h\mathbb{E}\left\langle\Delta^{m},U^{m}+\Phi^{m}+E^{m}\right\rangle\\\
&+h^{2}\mathbb{E}|U^{m}+\Phi^{m}+E^{m}|^{2}\\\ \leq&(1-2\mu
h)\mathbb{E}|\Delta^{m}|^{2}+h^{2}\mathbb{E}|U^{m}+\Phi^{m}+E^{m}|^{2}\end{aligned}\,,$
(60)

where we used (58). Plug (57) into (60), using the definition of the
coefficients $A,B$ in (53):

$\displaystyle\mathbb{E}|\Delta^{m}-h(U^{m}+\Phi^{m}+E^{m})|^{2}\leq
A\mathbb{E}|\Delta^{m}|^{2}+3h^{2}\mathbb{E}|E^{m}|^{2}+h^{3}B\,.$ (61)

Finally, we plug (56) and (61) into (59) and use the definition of the
coefficients $C$ in (53). Then we conclude (52).

###  7.2 Proof of Theorem 3

As discussed before, the iteration formula for $\Delta^{m}$ in RCD-O-LMC
satisfies the formula in Theorem 15, and thus the key is to give a bound to
$\mathbb{E}|E^{m}|^{2}$ for RCD-O-LMC. For that, we use the following lemma.

###### Lemma 17

Suppose $f$ satisfies Assumption 3.1,3.2, if $\\{x^{m}\\}$ is defined in (17)
by RCD-O-LMC, $\\{y^{m}\\}$ is defined in (50) and $\\{\Delta^{m}\\}$ comes
from (51), then for all $m\geq 0$:

$\mathbb{E}(|E^{m}|^{2})<2dL^{2}\mathbb{E}|\Delta^{m}|^{2}+2d^{2}L\,.$ (62)

The proof is shown in Appendix A. Now we are ready to prove Theorem 3.

Proof [Proof of Theorem 3] First, we prove (19). Plug (62) into (52), we have

$\mathbb{E}|\Delta^{m+1}|^{2}\leq(1+a)A^{\prime}\mathbb{E}|\Delta^{m}|^{2}+(1+a)(6Ld^{2}h^{2}+h^{3}B)+\left(1+\frac{1}{a}\right)h^{4}C\,,$

where $B$ and $C$ are kept untouched as in (53) and

$A^{\prime}=1-2\mu h+3L^{2}h^{2}+6L^{2}dh^{2}\,.$

Noting $A^{\prime}<1-\mu h$ since

$h<\frac{1}{3\kappa^{2}\mu+6\kappa^{2}\mu d}=\frac{\mu}{3L^{2}+6L^{2}d}\,.$

To ensure the decay of $\mathbb{E}|\Delta^{m}|^{2}$, we need to choose $a$
such that $(1+a)A^{\prime}=1-O(h)$. Set $a=\frac{\mu h/2}{1-\mu h}$ so that

$(1+a)(1-\mu h)=1-\frac{\mu h}{2}\,,\quad\text{and}\quad 1+1/a\leq 2/\mu h\,,$

then we have

$\mathbb{E}|\Delta^{m+1}|^{2}\leq\left(1-\frac{\mu
h}{2}\right)\mathbb{E}|\Delta^{m}|^{2}+2h^{2}\left(6Ld^{2}+hB+hC/\mu\right)\,,$
(63)

where we use $1+a\leq 2$. Using (63) iteratively, we finally have

$\mathbb{E}|\Delta^{m}|^{2}\leq\left(1-\frac{\mu
h}{2}\right)^{m}\mathbb{E}|\Delta^{0}|^{2}+32h\kappa
d^{2}\leq\exp\left(-\frac{\mu
hm}{2}\right)\mathbb{E}|\Delta^{0}|^{2}+32h\kappa d^{2}\,,$ (64)

where we use

$hB/\mu\leq\kappa d^{2},\quad hC/\mu^{2}\leq\kappa d^{2}$

by
$h<\min\left\\{\frac{d}{2\kappa\mu},\frac{2}{H^{2}/(\kappa\mu^{2})+\kappa^{2}\mu/d}\right\\}$.
Since

$(W^{O}_{m})^{2}\leq\mathbb{E}|\Delta^{m}|^{2},\quad\mathbb{E}|\Delta^{0}|^{2}=(W^{O}_{0})^{2}\,,$

we have (19).

###### Remark 18

We comment that in the proof, the discretization error terms $B$ and $C$
contribute $O(h^{3}d^{2})$ but the newly induced $\mathbb{E}|E|^{2}$ term from
(62) contributes $O(h^{2}d^{2})$, which eventually becomes the dominating term
in (19) aside from the exponential decaying term. This is the term that
suggests that the number of required iterations has to be at the order of
$d^{2}$.

###### Remark 19

As mentioned in Remark 5, using non-uniform sampling to select coordinates may
not enhance the computation, when $f$ has moderate conditioning. We provide
the intuition here.

When $f$ is not skewed, $\kappa=L/\mu\sim 1$. Without loss of generality, we
assume $f$ attains its minimal at origin, then:

$\mu|x|^{2}\leq|f(x)|\leq L|x|^{2}\Rightarrow|f(x)|\sim|x|^{2}\,,$ (65)

and for every $1\leq i\leq d$

$\mu|x_{i}|^{2}\leq|\partial_{i}f(x)|^{2}\leq
L|x_{i}|^{2}\Rightarrow|\partial_{i}f(x)|^{2}\sim|x_{i}|^{2}\,.$

When $m$ is large, the distribution of $x^{m}$ is close to the target
distribution. For simplicity, we further assume $x^{m}\sim p$. The using (21),
we calculate the variance of $E^{m}$ that is defined in (20):

$\mathbb{E}_{r}\left(\left|\frac{1}{\phi_{r}}\left(\nabla
f(x^{m})\cdot\textbf{e}^{r}\right)\textbf{e}^{r}-\nabla
f(x^{m})\right|^{2}\right)=\sum^{d}_{i=1}\left(\frac{1}{\phi_{i}}-1\right)\mathbb{E}|\partial_{i}f(x^{m})|^{2}\sim\sum^{d}_{i=1}\left(\frac{1}{\phi_{i}}-1\right)\mathbb{E}_{p}(|x^{m}_{i}|^{2})\,.$

According to (65), $\mathbb{E}_{p}(|x_{i}|^{2})\sim 1$, and thus the
derivation above becomes:

$\mathbb{E}\left(\left|\frac{1}{\phi_{r}}\left(\nabla
f(x^{m})\cdot\textbf{e}^{r}\right)\textbf{e}^{r}-\nabla
f(x^{m})\right|^{2}\right)\sim\sum^{d}_{i=1}\left(\frac{1}{\phi_{i}}-1\right)\,.$
(66)

Noting $\sum_{i=1}^{d}\phi_{i}=1$, the right-hand side of (66) achieves the
smallest value when $\phi_{i}=1/d$. This recovers the uniform RCD-LMC.

This discussion suggests that the non-uniform selection of coordinates does
not bring computational efficiency to RCD-LMC when $f$ has moderate
conditioning. We note that the same scenario is observed in optimization
(Nesterov, 2012) where the outperformance of RCD over GD crucially depends on
the skewness of $f$. When $f$ is skewed, then the Lipschitz constants in
different directions are drastically different. We investigate the
conditioning’s effect on how coordinates are chosen in (Ding et al., 2021b,
a), where we show improvement can be observed if stiffer directions are chosen
more frequently.

###  7.3 Proof of Theorem 9

The strategy for showing SVRG-O-LMC is the same. We first give a bound to the
variance term, and then combine the estimate with Theorem 15 for proving the
theorem.

###### Lemma 20

Suppose $f$ satisfies Assumptions 3.1-3.2, if $\\{x^{m}\\}$ is defined in (28)
by SVRG-O-LMC, $\\{y^{m}\\}$ is defined in (50) and $\\{\Delta^{m}\\}$ comes
from (51), then for all $m\geq 0$, define

$s=\left\lfloor\frac{m}{\tau}\right\rfloor\,,$

then we have

$\mathbb{E}\left|E^{m}\right|^{2}\leq
3dL^{2}\left[\mathbb{E}\left|\Delta^{m}\right|^{2}+\mathbb{E}\left|\Delta^{s\tau}\right|^{2}+2h^{2}Ld\tau^{2}+4hd\tau\right]\,.$
(67)

We put the proof in Appendix B. Now we are ready to prove Theorem 9

Proof [Proof of Theorem 9] Recall

$s=\left\lfloor\frac{m}{\tau}\right\rfloor\,.$

Use (52), similar to (63), set $a=\frac{\mu h/2}{1-\mu h}$, use
$h<\frac{1}{48d\kappa^{2}\mu}<\frac{\mu}{3L^{2}}$, we have

$1+a=\frac{1-\mu h/2}{1-\mu h}<2,\quad 1+1/a=2/(\mu h)\,,$

then

$\displaystyle\mathbb{E}|\Delta^{m+1}|^{2}$
$\displaystyle\leq\left(1-\frac{\mu
h}{2}\right)\mathbb{E}|\Delta^{m}|^{2}+6h^{2}\mathbb{E}|E^{m}|^{2}+2h^{3}B+2h^{3}C/\mu$
(68) $\displaystyle\leq\left(1-\frac{\mu
h}{2}\right)\mathbb{E}|\Delta^{m}|^{2}+6h^{2}\mathbb{E}|E^{m}|^{2}+h^{3}D\,,$

where we set $D=2B+2C/\mu$.

Plug (67) into (68),

$\displaystyle\mathbb{E}|\Delta^{m+1}|^{2}\leq$
$\displaystyle\left(1-\frac{\mu
h}{4}\right)\mathbb{E}|\Delta^{m}|^{2}+18h^{2}dL^{2}\mathbb{E}\left|\Delta^{s\tau}\right|^{2}$
(69) $\displaystyle+36h^{4}\tau^{2}d^{2}L^{3}+72\tau
d^{2}h^{3}L^{2}+h^{3}D\,,$

where we use $h<\frac{1}{72d\kappa^{2}\mu}=\frac{\mu}{72dL^{2}}$.

Use (69) iteratively from $s\tau$ to $m$, we have

$\displaystyle\mathbb{E}|\Delta^{m}|^{2}\leq$ $\displaystyle\left(1-\frac{\mu
h}{4}\right)^{m-s\tau}\mathbb{E}|\Delta^{s\tau}|^{2}+18h^{2}(m-s\tau)dL^{2}\mathbb{E}\left|\Delta^{s\tau}\right|^{2}$
(70)
$\displaystyle+36h^{4}\tau^{3}d^{2}L^{3}+72h^{3}\tau^{2}d^{2}L^{2}+h^{3}\tau
D\,.$

Since $\mu h\tau<10$, we have $(1-\frac{\mu h}{4})^{\tau}<1-\frac{\mu
h\tau}{8}$, and with the other $h^{2}\tau dL^{2}$ term, this is still
controlled by $1-\frac{\mu h\tau}{16}$. Therefore, choose $m=(s+1)\tau$ and
use $h<\frac{\mu}{400dL^{2}}$, we have

$\displaystyle\mathbb{E}|\Delta^{(s+1)\tau}|^{2}\leq$
$\displaystyle\left(1-\frac{\mu
h\tau}{16}\right)\mathbb{E}|\Delta^{s\tau}|^{2}+36h^{4}\tau^{3}d^{2}L^{3}+72h^{3}\tau^{2}d^{2}L^{2}+h^{3}\tau
D$ (71) $\displaystyle<$ $\displaystyle\left(1-\frac{\mu
h\tau}{16}\right)^{s+1}\mathbb{E}|\Delta^{0}|^{2}+576h^{3}\tau^{2}d^{2}L^{3}/\mu+1152h^{2}\tau
d^{2}L^{3}/\mu+16h^{2}D/\mu\,.$

Plug (71) back into (70), use $18h^{2}\tau dL^{2}<1$, we have

$\displaystyle\mathbb{E}|\Delta^{m}|^{2}<$
$\displaystyle\left(\left(1-\frac{\mu h}{4}\right)^{m-s\tau}+18h^{2}\tau
dL^{2}\right)\left(1-\frac{\mu
h\tau}{16}\right)^{s}\mathbb{E}|\Delta^{0}|^{2}$
$\displaystyle+650h^{3}\tau^{2}d^{2}L^{3}/\mu+1250h^{2}\tau
d^{2}L^{2}/\mu+17h^{2}D/\mu\,.$

Use $h\tau<\frac{1}{10}$ and $18h^{2}(m-s\tau)dL^{2}<\frac{\mu
h(m-s\tau)}{16}$, we have

$\displaystyle\left(1-\frac{\mu h}{4}\right)^{m-s\tau}+18h^{2}(m-s\tau)dL^{2}$
$\displaystyle\leq 1-\frac{\mu h(m-s\tau)}{8}+18h^{2}(m-s\tau)dL^{2}\leq
1-\frac{\mu h(m-s\tau)}{16}$ $\displaystyle\leq\left(1-\frac{\mu
h}{16}\right)^{m-s\tau}\,,$

which implies

$\mathbb{E}|\Delta^{m}|^{2}<\exp\left(-\frac{\mu
hm}{16}\right)\mathbb{E}|\Delta^{0}|^{2}+650h^{3}\tau^{2}d^{2}L^{3}/\mu+1250h^{2}\tau
d^{2}L^{2}/\mu+17h^{2}D/\mu\,.$ (72)

Taking square root on both sides and use

$D=2B+2C/\mu\leq
d^{2}(4\kappa^{2}\mu^{2}/d+\kappa^{3}\mu^{2}/d+2H^{2}/\mu)\,.$

Since $\tau\geq 1,d\geq 1$, combine last two terms in (72), we have (33).

Similar to the comment at the end of Section 7.2, we here briefly review the
$d$ and $h$ dependence. As is the case in Section 7.2, the $B$ term and the
$C$ term (glued together to be called the $D$ term here) contribute $h^{3}$.
The error induced by $\mathbb{E}|E|^{2}$ is bounded by (67), when entering
(70) is also at the order of $h^{3}$. These two terms combined explain the
last term in (33). However we should note that $\tau$ is typically set at the
order of $d$, so the $d$ dependence from the newly introduced error in (71) is
in fact $d^{3}$ as a comparison of $d^{2}$ in $D$. This $d^{3}$ term
determines the $m\sim d^{3/2}$ dependence in the end.

###  7.4 Sketch proof of Theorem 10

The same strategy as was used in the previous two subsections is used here to
give an estimate to RCAD-O-LMC. We first bound the variance term
$\mathbb{E}|E^{m}|^{2}$.

We first define a few notations. Define $\beta^{m}$ using (48) with the same
$r^{m}$ but replacing $x^{m}$ with $y^{m}$:

$\beta^{0}=\nabla f(y^{0})$ (73)

and

$\beta^{m+1}_{r^{m}}=\partial_{r^{m}}f(y^{m})\quad\text{and}\quad\beta^{m+1}_{i}=\beta^{m}_{i}\quad\text{if}\quad
i\neq r^{m}\,.$ (74)

Then we have:

###### Lemma 21

Suppose $f$ satisfies Assumptions 3.1-3.2, if $\\{x^{m}\\}$ is defined in (31)
by RCAD-O-LMC, $\\{y^{m}\\}$ is defined in (50), and $\\{\Delta^{m}\\}$ comes
from (51), then for all $m\geq 0$, we have

$\displaystyle\mathbb{E}\left|E^{m}\right|^{2}\leq$ $\displaystyle
3dL^{2}\mathbb{E}|\Delta^{m}|^{2}+24hL^{2}d^{3}\left[hLd+1\right]+3d\mathbb{E}\left|\beta^{m}-g^{m}\right|^{2}\,.$
(75)

The proof is seen in (Ding and Li, 2020). One big difference between this
analysis and the counterparts in the previous two subsections is that the
bound of $\mathbb{E}|E^{m}|^{2}$ depends on $\nabla f(y^{m})-\nabla f(x^{m})$,
reflected in the $\beta^{m}-g^{m}$ term in (75). This term traces the history
of the random coordinate selection, and cannot be merely bounded by
$\Delta^{m}$, making the full iteration self-consistent. To overcome the
difficulty, we define a Lyapunov function that combines the effect of
$\Delta^{m}$ and this particular term $\beta^{m}-g^{m}$.

The Lyapunov function is defined as follows:

$T^{m}=T^{m}_{1}+c_{p}T^{m}_{2}=\mathbb{E}|\Delta^{m}|^{2}+c_{p}\mathbb{E}|g^{m}-\beta^{m}|^{2}\,,$
(76)

where $c_{p}$ will be carefully chosen later. And the relation between $T_{1}$
and $T_{2}$ are now in the following lemma:

###### Lemma 22

Under conditions of Theorem 10, if $\\{x^{m}\\}$ is defined in (31),
$\\{y^{m}\\}$ is defined in (50) and $\\{T^{m}_{1},T^{m}_{2}\\}$ comes from
(76), then for any $a>0,m\geq 0$, there are upper bounds of $T^{m}_{1}$ and
$T^{m}_{2}$:

$T^{m+1}_{1}\leq(1+a)AT^{m}_{1}+(1+a)BT^{m}_{2}+(1+a)h^{3}C+\left(1+\frac{1}{a}\right)h^{4}D\,,$
(77)

and

$T^{m+1}_{2}\leq\widetilde{A}T^{m}_{1}+\widetilde{B}T^{m}_{2}\,,$ (78)

where

$\displaystyle A=1-2\mu h+3(1+3d)L^{2}h^{2}\,,\quad B=9h^{2}d\,,$
$\displaystyle C=2L^{2}d+72L^{2}d^{3}\left[hLd+1\right]\,,\quad
D=\frac{1}{2}\left(H^{2}d^{2}+L^{3}d\right)\,,$
$\displaystyle\widetilde{A}=\frac{L^{2}}{d}\,,\quad\widetilde{B}=1-1/d\,.$

We refer interested readers to the proof in (Ding and Li, 2020), where we
present how to choose proper $a$ and $c_{p}$ in the lemma above to obtain
Theorem 10. Here, we still comment on the dependence of $h$ on $d$ again. Due
to the complicated Lyapunov function, and the involvement of
$\beta^{m}-g^{m}$, the intuition is not as straightforward as it was for
Theorem 3 and 9. However, according to Lemma 22, $\widetilde{A}$ is as small
as $1/d$ while $\widetilde{B}$ is almost $1$, meaning $T^{m}_{2}$ barely
changes. Moreover, $c_{p}$ is as small as $h^{2}$ and is neglected from the
calculation quickly. The dominating terms, in the end, are still the $C$ term
and the $D$ term in (77), both of which contribute $h^{3}$. The highest power
of $d$ in these two terms is the $d^{3}$ term that comes from the last term in
$C$, and this is introduced by $\mathbb{E}|E^{m}|^{2}$ in (75). This
$h^{3}d^{3}$ term, upon iteration, becomes the $h^{2}d^{3}$ term as shown in
the theorem, which finally explains the $d^{3/2}$ power in the numerical cost.

As shown in Theorem 3 and 9, in this theorem, we once again encounter the
situation where the error of $\mathbb{E}|E^{m}|^{2}$ dominates the term in the
final iteration formula.

##  8 Proof of convergence results for U-LMC

In this section, we prove the convergence results for RCD-U-LMC, RCAD-U-LMC,
and SVRG-U-LMC.

###  8.1 Iteration formula for U-LMC

Recall the definition $E^{m}=\nabla f(x^{m})-F^{m}$, as defined in (20).

According to Algorithm 2, 4, all three algorithms can be seen as drawing
$(x^{0},v^{0})$ from distribution induced by $q^{U}_{0}$, and update
$(x^{m},v^{m})$ using the following coupled SDEs:

$\left\\{\begin{aligned}
&\mathrm{V}_{t}=v^{m}e^{-2h}-\gamma\int^{t}_{mh}e^{-2(t-s)}\,\mathrm{d}sF^{m}+\sqrt{4\gamma}\int^{t}_{mh}e^{-2(t-s)}dB_{s}\\\
&\mathrm{X}_{t}=x^{m}+\int^{t}_{mh}\mathrm{V}_{s}ds\end{aligned}\right.\,,$
(79)

where $B_{s}$ is the Brownian motion and
$(x^{m+1},v^{m+1})=(\mathrm{X}_{(m+1)h},\mathrm{V}_{(m+1)h})$.

We then define $w^{m}=x^{m}+v^{m}$, and denote $u_{m}(x,w)$ the distribution
of $(x^{m},w^{m})$ and $u^{\ast}(x,w)$ the distribution of $(x,w)$ if
$(x,v=w-x)$ is distributed according to density function $p_{2}$. From (Cheng
et al., 2018, Lemma 8), we have:

$|x^{m}-x|^{2}+|v^{m}-v|^{2}\leq 4(|x^{m}-x|^{2}+|w^{m}-w|^{2})\leq
16(|x^{m}-x|^{2}+|v^{m}-v|^{2})\,\\\ $ (80)

and

$W^{2}_{2}(q^{U}_{m},p_{2})\leq 4W^{2}_{2}(u_{m},u^{*})\leq
16W^{2}_{2}(q^{U}_{m},p_{2})\,.$ (81)

Similar to O-LMC, define another trajectory of sampling by setting
$(\widetilde{x}^{0},\widetilde{v}^{0})$ to be drawn from distribution induced
by $p_{2}$, and let $\widetilde{x}^{m}=\widetilde{\mathrm{X}}_{mh}$,
$v^{m}=\widetilde{\mathrm{V}}_{mh}$ be samples from
$\left(\widetilde{\mathrm{X}}_{t},\widetilde{\mathrm{V}}_{t}\right)$ that
satisfy

$\left\\{\begin{aligned}
&\widetilde{\mathrm{V}}_{t}=\widetilde{v}^{0}e^{-2t}-\gamma\int^{t}_{0}e^{-2(t-s)}\nabla
f\left(\widetilde{\mathrm{X}}_{s}\right)\,\mathrm{d}s+\sqrt{4\gamma}e^{-2t}\int^{t}_{0}e^{2s}dB_{s}\\\
&\widetilde{\mathrm{X}}_{t}=\widetilde{x}^{0}+\int^{t}_{0}\widetilde{\mathrm{V}}_{s}ds\end{aligned}\right.\,$
(82)

with the same Brownian motion as before.

Clearly $\left(\widetilde{\mathrm{X}}_{t},\widetilde{\mathrm{V}}_{t}\right)$
can be seen as drawn from distribution induced by $p_{2}$ for all $t$, and
initially we can pick $(\widetilde{x}^{0},\widetilde{v}^{0})$ such that

$W^{2}_{2}(u_{0},u^{\ast})=\mathbb{E}\left(|x^{0}-\widetilde{x}^{0}|^{2}+|w^{0}-\widetilde{w}^{0}|^{2}\right)\,.$

We then also define the error $\Delta^{m}$:

$\Delta^{m}=\sqrt{|\widetilde{x}^{m}-x^{m}|^{2}+|\widetilde{w}^{m}-w^{m}|^{2}}\,.$
(83)

Similar to Theorem 15, we also have an important iteration formula for
$|\Delta^{m}|^{2}$:

###### Theorem 23

Assume $f$ satisfies Assumption 3.1, let $\\{(x^{m},v^{m})\\}$ defined in
(79), $\\{(\widetilde{x}^{m},\widetilde{v}^{m})\\}$ defined in (82) and
$\\{\Delta^{m}\\}$ defined in (83), then there exists a uniform constant $D>0$
such that if $h$ satisfies

$h<\frac{1}{100(1+D)\kappa}\,,$

then for any $m\geq 0$, we have

$\mathbb{E}|\Delta^{m+1}|^{2}\leq
A\mathbb{E}|\Delta^{m}|^{2}+B\mathbb{E}|E^{m}|^{2}+C\,,$ (84)

where

$A=1-h/(2\kappa)+480\kappa h^{3},\quad B=10\gamma^{2}h^{2}\,,\quad
C=30h^{3}d/\mu\,.$

Proof

We first divide $|\Delta^{m+1}|^{2}$ into different parts , and compare (79)
and (82) for:

$\displaystyle|\Delta^{m+1}|^{2}=$
$\displaystyle\left|(\widetilde{v}^{m}-v^{m})e^{-2h}+(\widetilde{x}^{m}-x^{m})+\int^{(m+1)h}_{mh}\widetilde{\mathrm{V}}_{s}-\mathrm{V}_{s}\,\mathrm{d}s\right.$
$\displaystyle-\gamma\int^{(m+1)h}_{mh}e^{-2((m+1)h-s)}\left[\nabla
f\left(\widetilde{\mathrm{X}}_{s}\right)-\nabla f(x^{m})\right]\,\mathrm{d}s$
$\displaystyle\left.+\gamma\int^{(m+1)h}_{mh}e^{-2((m+1)h-s)}E^{m}\,\mathrm{d}s\right|^{2}$
$\displaystyle+\left|(\widetilde{x}^{m}-x^{m})+\int^{(m+1)h}_{mh}\widetilde{\mathrm{V}}_{s}-\mathrm{V}_{s}\,\mathrm{d}s\right|^{2}$
$\displaystyle=$
$\displaystyle\left|\mathrm{J}^{m}_{1}\right|^{2}+\left|\mathrm{J}^{m}_{2}\right|^{2}=\left|\mathrm{J}^{r,m}_{1}+\mathrm{J}^{E,m}_{1}\right|^{2}+\left|\mathrm{J}^{m}_{2}\right|^{2}\,,$

where we denote

$\displaystyle\mathrm{J}^{r,m}_{1}$
$\displaystyle=(\widetilde{v}^{m}-v^{m})e^{-2h}+(\widetilde{x}^{m}-x^{m})+\int^{(m+1)h}_{mh}\widetilde{\mathrm{V}}_{s}-\mathrm{V}_{s}\,\mathrm{d}s$
(85) $\displaystyle-\gamma\int^{(m+1)h}_{mh}e^{-2((m+1)h-s)}\left[\nabla
f\left(\widetilde{\mathrm{X}}_{s}\right)-\nabla
f(x^{m})\right]\,\mathrm{d}s\,$

and

$\displaystyle\mathrm{J}^{E,m}_{1}=\gamma\int^{(m+1)h}_{mh}e^{-2((m+1)h-s)}E^{m}\,\mathrm{d}s\,,\quad\mathrm{J}^{m}_{2}=(\widetilde{x}^{m}-x^{m})+\int^{(m+1)h}_{mh}\widetilde{\mathrm{V}}_{s}-\mathrm{V}_{s}\,\mathrm{d}s\,.$
(86)

To control $\mathrm{J}^{m}_{1}$, we realize all terms in
$\mathrm{J}^{r,m}_{1}$, except $\mathrm{V}_{s}$, are independent of $r^{m}$,
and thus,

$\mathbb{E}\langle\mathrm{J}^{r,m}_{1},\mathrm{J}^{E,m}_{1}\rangle=-\mathbb{E}\left(\mathbb{E}_{r^{m}}\left\langle\int^{(m+1)h}_{mh}\mathrm{V}_{s}\,\mathrm{d}s,\mathrm{J}^{E,m}_{1}\right\rangle\right)\,,$
(87)

where all other terms are eliminated during the process of taking
$\mathbb{E}_{r^{m}}$. Therefore, we can bound
$\mathbb{E}\left|\mathrm{J}^{m}_{1}\right|^{2}$ by

$\displaystyle\mathbb{E}\left|\mathrm{J}^{m}_{1}\right|^{2}=\mathbb{E}\left|\mathrm{J}^{r,m}_{1}+\mathrm{J}^{E,m}_{1}\right|^{2}$
$\displaystyle=$
$\displaystyle\mathbb{E}\left|\mathrm{J}^{r,m}_{1}\right|^{2}+\mathbb{E}\left|\mathrm{J}^{E,m}_{1}\right|^{2}+2\mathbb{E}\left\langle\mathrm{J}^{r,m}_{1},\mathrm{J}^{E,m}_{1}\right\rangle$
$\displaystyle\leq$
$\displaystyle\mathbb{E}\left|\mathrm{J}^{r,m}_{1}\right|^{2}+\gamma^{2}h^{2}\mathbb{E}\left|E^{m}\right|^{2}+2\gamma^{2}h^{3}\mathbb{E}|E^{m}|^{2}\,,$

where we used

$\mathbb{E}\left|\mathrm{J}^{E,m}_{1}\right|^{2}\leq\gamma^{2}h^{2}\mathbb{E}\left|E^{m}\right|^{2}\,$

and

$\begin{aligned}
&2\mathbb{E}\langle\mathrm{J}^{r,m}_{1},\mathrm{J}^{E,m}_{1}\rangle\\\
=&-2\mathbb{E}\left\langle\int^{(m+1)h}_{mh}\mathrm{V}_{s}\,\mathrm{d}s,\gamma\int^{(m+1)h}_{mh}e^{-2((m+1)h-s)}\,\mathrm{d}sE^{m}\right\rangle\\\
=&2\mathbb{E}\left\langle\gamma\int^{(m+1)h}_{mh}\int^{s}_{mh}e^{-2(s-t)}\,\mathrm{d}t\,\mathrm{d}sE^{m},\gamma\int^{(m+1)h}_{mh}e^{-2((m+1)h-s)}\,\mathrm{d}sE^{m}\right\rangle\\\
\leq&2\gamma^{2}h^{3}\mathbb{E}|E^{m}|^{2}\end{aligned}\,,$

In conclusion, we have

$\displaystyle\mathbb{E}\left|\Delta^{m+1}\right|^{2}\leq$
$\displaystyle\mathbb{E}\left|\mathrm{J}^{r,m}_{1}\right|^{2}+\left|\mathrm{J}^{m}_{2}\right|^{2}+\left(\gamma^{2}h^{2}\mathbb{E}|E^{m}|^{2}+2\gamma^{2}h^{3}\mathbb{E}|E^{m}|^{2}\right)\,,$
(88)

which implies

$\displaystyle\mathbb{E}\left|\Delta^{m+1}\right|^{2}\leq$
$\displaystyle\mathbb{E}\left|\mathrm{J}^{r,m}_{1}\right|^{2}+\left|\mathrm{J}^{m}_{2}\right|^{2}+2\gamma^{2}(h^{2}+3h^{3})\mathbb{E}|E^{m}|^{2}\,.$
(89)

According to Proposition 29 (114), for any $a>0$, we have

$\mathbb{E}\left|\mathrm{J}^{r,m}_{1}\right|^{2}+\left|\mathrm{J}^{m}_{2}\right|^{2}\leq
C_{1}\mathbb{E}|\Delta^{m}|^{2}+5(1+1/a)\gamma^{2}h^{4}\mathbb{E}(|E^{m}|^{2})+5(1+1/a)\gamma
h^{4}d\,,$

where

$C_{1}=(1+a)[1-h/\kappa+200h^{2}]+80(1+1/a)h^{4}\,.$

Plug into (89), we have

$\displaystyle\mathbb{E}|\Delta^{m}|^{2}\leq$ $\displaystyle
C_{1}\mathbb{E}|\Delta^{m}|^{2}+C_{2}\mathbb{E}|E^{m}|^{2}+C_{3}\,,$ (90)

where we use $\gamma L=1,h<1$ for:

$C_{2}=\gamma^{2}\left[5(1+1/a)h^{4}+8h^{2}\right]\,,\quad
C_{3}=5(1+1/a)\gamma h^{4}d\,.$

Considering $h\leq\frac{1}{600\kappa}$, we have

$1-h/\kappa+200h^{2}\leq 1-2h/(3\kappa)\,,$

and thus by setting $a$ so that

$1+a=\frac{1-h/(2\kappa)}{1-2h/(3\kappa)}\,,\quad\text{which leads}\quad
1+1/a\leq 6\kappa/h\,,$

we have

$C_{1}\leq 1-h/(2\kappa)+480\kappa h^{3}\,,\quad C_{2}\leq 30\gamma^{2}\kappa
h^{3}+8\gamma^{2}h^{2}\,,\quad C_{3}\leq 30\kappa\gamma h^{3}d\,.$

Noting $h<\frac{1}{600\kappa}$, we have $30\kappa h^{3}<2h^{2}$, making:

$C_{1}\leq A\,,\quad C_{2}\leq B=10\gamma^{2}h^{2}\,,\quad C_{3}\leq
C=30h^{3}d/\mu\,,$

finishing the proof.

Similar to before we need to make a list of comments:

###### Remark 24

Several comments are in order.

  * •

The proof is completely independent of $\mathbb{E}|E^{m}|^{2}$, meaning the
lemma holds true for all U-LMC. The differences between different sampling
methods come in through the $\mathbb{E}|E^{m}|^{2}$ term that will be analyzed
in the following subsections respectively.

  * •

Similar to O-LMC, the last term of (84) corresponds to the discretization
error of U-LMC. In the ideal case, let us suppose $\mathbb{E}|E^{m}|^{2}$ is
negligible, then besides the iteration term ($A$ term), the remainder term $C$
is at the order of $h^{3}d$. Upon iterations, this will finally lead to a
formula similar to $W_{m}<\exp(-hm)W_{0}+hd^{1/2}$. With the same argument as
before, we have $h<\frac{\epsilon}{d^{1/2}}$, and $m>1/h$ for an $\epsilon$
accuracy. This is the optimal result one could hope for in the underdamped LMC
case. As will be shown below, the $\mathbb{E}|E^{m}|^{2}$ is a dominating term
for RCD-U-LMC, and such optimal rate is not achieved. However, in the other
two cases, $\mathbb{E}|E^{m}|^{2}$ is not a dominating term if $\epsilon$ is
small enough, this means a rate the same as U-LMC is possible for small
$\epsilon$.

###  8.2 Proof of Theorem 6

To prove Theorem 6, we note that the definition of $E$ is the same as in (20),
and the results in Lemma 17 can be recycled:

Proof [Proof of Theorem 6] Recall Lemma 17, we have:

$\mathbb{E}(|E^{m}|^{2})<2dL^{2}\mathbb{E}|\Delta^{m}|^{2}+2L^{2}d^{2}/\mu\,.$
(91)

Plug (91) into (84), we have

$\displaystyle\mathbb{E}|\Delta^{m+1}|^{2}$
$\displaystyle\leq(A+2dL^{2}B)\mathbb{E}|\Delta^{m}|^{2}+2L^{2}d^{2}B/\mu+C$
(92) $\displaystyle\leq C_{1}\mathbb{E}|\Delta^{m}|^{2}+2L^{2}d^{2}B/\mu+C\,.$

Recalling

$A=1-h/(2\kappa)+480\kappa h^{3},\quad B=10\gamma^{2}h^{2}\,,\quad
C=30h^{3}d/\mu\,,$

we have, due to $\gamma L=1$:

$C_{1}=1-h/(2\kappa)+480\kappa h^{3}+20dh^{2}\,.$

For $h$ small enough $\left(h<\frac{1}{880d\kappa}\right)$, it is clear that

$C_{1}<1-h/(4\kappa)\,.$ (93)

Plug (93) into (92), we have

$\mathbb{E}|\Delta^{m+1}|^{2}\leq\left(1-h/(4\kappa)\right)\mathbb{E}|\Delta^{m}|^{2}+20d^{2}h^{2}/\mu+30h^{3}d/\mu\,,$
(94)

which implies

$\displaystyle\mathbb{E}|\Delta^{m}|^{2}$
$\displaystyle\leq\exp\left(-hm/(4\kappa)\right)\mathbb{E}|\Delta^{0}|^{2}+80\kappa
d^{2}h/\mu+120\kappa h^{2}d/\mu$
$\displaystyle\leq\exp\left(-hm/(4\kappa)\right)\mathbb{E}|\Delta^{0}|^{2}+100\kappa
d^{2}h/\mu$

after iteration. Taking square root on each term and using (81), we prove
(23).

###### Remark 25

The result comes with no surprise. The discretization error term $C$
contributes $O(h^{3}d)$ but the newly induced $\mathbb{E}|E|^{2}$ term from
(62) contributes $O(h^{2}d^{2})$ . This $h^{2}d^{2}$ term, upon iteration,
drops down to $hd^{2}$, which explains the $d^{2}$ dependence in the end.

###  8.3 Proof of Theorem 12

We first have the lemma bounding $\mathbb{E}|E^{m}|^{2}$:

###### Lemma 26

Assume $f$ satisfies Assumption 3.1, and let $\\{(x^{m},v^{m})\\}$ be defined
in algorithm 3 (underdamped), $\\{(\widetilde{x}^{m},\widetilde{v}^{m})\\}$ be
defined in (82), then for any $m\geq 0$, call

$s=\left\lfloor\frac{m}{\tau}\right\rfloor\,,$

we have

$\mathbb{E}\left|E^{m}\right|^{2}\leq
3dL^{2}\left[\mathbb{E}\left|\Delta^{m}\right|^{2}+\mathbb{E}\left|\Delta^{s\tau}\right|^{2}+2\tau^{2}h^{2}\gamma
d\right]\,.$ (95)

We put the proof in Appendix C. Now, we are ready to prove the Theorem 12.

Proof [Proof of Theorem 12] Use Theorem 23 (84), we first bound
$\mathbb{E}|E^{m}|^{2}$ by plugging in Lemma 26 equation (95),

$\mathbb{E}|\Delta^{m+1}|^{2}\leq(A+3dL^{2}B)\mathbb{E}|\Delta^{m}|^{2}+3dL^{2}B\mathbb{E}|\Delta^{s\tau}|^{2}+6Ld^{2}\tau^{2}h^{2}B+C\,.$

Recalling $A=1-h/(2\kappa)+480\kappa h^{3}$, $B=10\gamma^{2}h^{2}$ and
$C=30h^{3}d/\mu$, we have

$A+3dL^{2}B\leq 1-h/(4\kappa)\,,\quad\text{when}\quad h<\frac{1}{1648d\kappa}$

and thus

$\displaystyle\mathbb{E}|\Delta^{m+1}|^{2}$
$\displaystyle\leq\left(1-\frac{h}{4\kappa}\right)\mathbb{E}|\Delta^{m}|^{2}+30h^{2}d\mathbb{E}|\Delta^{s\tau}|^{2}+30h^{3}d/\mu+60\gamma
h^{4}\tau^{2}d^{2}\,.$ (96)

Use $\gamma=\frac{1}{L}=\frac{1}{\kappa\mu}$, we have:

$\displaystyle\mathbb{E}|\Delta^{m}|^{2}<$
$\displaystyle\left(1-\frac{h}{4\kappa}\right)^{m-s\tau}\mathbb{E}|\Delta^{s\tau}|^{2}+30h^{2}(m-s\tau)d\mathbb{E}|\Delta^{s\tau}|^{2}+30h^{3}\tau
d/\mu+60h^{4}\tau^{2}d^{2}/(\kappa\mu)\,.$ (97)

Similar to (71), by setting $m=(s+1)\tau$, we also have:

$\displaystyle\mathbb{E}|\Delta^{(s+1)\tau}|^{2}<$
$\displaystyle\left(1-\frac{h\tau}{8\kappa}\right)\mathbb{E}|\Delta^{s\tau}|^{2}+30h^{2}\tau
d\mathbb{E}|\Delta^{s\tau}|^{2}+30h^{3}\tau
d/\mu+60h^{4}\tau^{3}d^{2}/(\kappa\mu)$ (98) $\displaystyle<$
$\displaystyle\left(1-\frac{h\tau}{16\kappa}\right)\mathbb{E}|\Delta^{s\tau}|^{2}+30h^{3}\tau
d/\mu+60h^{4}\tau^{3}d^{2}/(\kappa\mu)\,,$

where we absorb the $30h^{2}\tau d$ term into $\frac{h\tau}{16\kappa}$ for
small enough $h$ (for example when $h<\frac{1}{1648d\kappa}$) in the second
inequality. This finally gives:

$\mathbb{E}|\Delta^{s\tau}|^{2}<\left(1-\frac{h\tau}{16\kappa}\right)^{s}\mathbb{E}|\Delta^{0}|^{2}+480h^{2}d\kappa/\mu+960h^{3}\tau^{2}d^{2}/\mu\,.$
(99)

Plug (99) back into (97), we have:

$\displaystyle\mathbb{E}|\Delta^{m}|^{2}<$
$\displaystyle\left(\left(1-\frac{h}{4\kappa}\right)^{m-s\tau}+30h^{2}(m-s\tau)d\right)\left(1-\frac{h\tau}{16\kappa}\right)^{s}\mathbb{E}|\Delta^{0}|^{2}+510h^{2}d\kappa/\mu+1020h^{3}\tau^{2}d^{2}/\mu$
$\displaystyle<$
$\displaystyle\left(1-\frac{h}{16\kappa}\right)^{m}\mathbb{E}|\Delta^{0}|^{2}+510h^{2}d\kappa/\mu+1020h^{3}\tau^{2}d^{2}/\mu$
$\displaystyle<$
$\displaystyle\exp\left(-\frac{hm}{16\kappa}\right)\mathbb{E}|\Delta^{0}|^{2}+510h^{2}d\kappa/\mu+1020h^{3}\tau^{2}d^{2}/\mu\,,$

where in the first inequality we used $h\tau<1\leq\kappa$, and in the second
inequality we used $h(m-s\tau)\leq h\tau<\frac{1}{40}$ and
$30h^{2}(m-s\tau)d<1-\frac{(m-s\tau)h}{16\kappa}$ to obtain

$\displaystyle\left(1-\frac{h}{4\kappa}\right)^{m-s\tau}+30h^{2}(m-s\tau)d$
$\displaystyle<1-\frac{(m-s\tau)h}{8\kappa}+30h^{2}(m-s\tau)d<1-\frac{(m-s\tau)h}{16\kappa}$
$\displaystyle\leq\left(1-\frac{h}{16\kappa}\right)^{m-s\tau}\,.$

Taking square root on both sides gives the conclusion (37).

Again, we discuss the newly introduced error by $\mathbb{E}|E^{m}|^{2}$. Take
$\tau=d$ in Lemma 26, this part of error is at the order of $h^{4}d^{4}$ when
entering the iteration formula with an extra $h^{2}$ term from $B$. By
induction, one $h$ drops out and upon taking square root, this part of error
is at the order of $h^{3/2}d^{2}$. Making it smaller than $\epsilon$ gives
$h<\frac{\epsilon^{2/3}}{d^{4/3}}$. When $\epsilon$ is small enough, this is a
more relaxed constraint than the original discretization constraint
$h<\frac{\epsilon^{1/2}}{d}$.

###  8.4 Sketch proof of Theorem 13

Similar to section 7.4, we have the following estimation for
$\mathbb{E}\left|E^{m}\right|^{2}$:

###### Lemma 27

Assume $f$ satisfies Assumption 3.1, if $\\{(x^{m},v^{m})\\}$ is defined in
algorithm 4 (underdamped), $\\{(\widetilde{x}^{m},\widetilde{v}^{m})\\}$ is
defined in (82), then for any $m\geq 0$, we have

$\mathbb{E}(|E^{m}|^{2})\leq
3dL^{2}\mathbb{E}|\Delta^{m}|^{2}+24Lh^{2}d^{4}+3d\mathbb{E}\left|\beta^{m}-g^{m}\right|^{2}\,.$
(100)

We leave out the proof and refer interested readers to (Ding and Li, 2020).
Here, $\beta^{m}$ and $g$ are defined same as (74) and (48) ($y$ change to
$\widetilde{x}$) and we will show the decay of the following Lyapunov
function:

$T^{m}\triangleq
T^{m}_{1}+c_{p}T^{m}_{2}=\mathbb{E}\left(|\widetilde{x}^{m}-x^{m}|^{2}+|\widetilde{w}^{m}-w^{m}|^{2}\right)+c_{p}\mathbb{E}|g^{m}-\beta^{m}|^{2}\,,$
(101)

where $c_{p}$ will be carefully chosen later. And the relation between $T_{1}$
and $T_{2}$ are now in the following lemma:

###### Lemma 28

Under conditions of Theorem 13, if $\\{(x^{m},v^{m})\\}$ is defined in
algorithm 4 (underdamped), $\\{(\widetilde{x}^{m},\widetilde{v}^{m})\\}$ is
defined in (82) and $\\{(T^{m}_{1},T^{m}_{2})\\}$ comes from (101), then for
any $m\geq 0$, we have

$\displaystyle T^{m+1}_{1}<$ $\displaystyle
D_{1}T^{m}_{1}+D_{2}T^{m}_{2}+D_{3}\,,$ (102)
$T^{m+1}_{2}\leq\frac{L^{2}}{d}T^{m}_{1}+\left(1-\frac{1}{d}\right)T^{m}_{2}\,,$
(103)

where

$D_{1}=1-h/(2\kappa)+40h^{2}d,\ D_{2}=30\gamma^{2}h^{2}d,\ D_{3}=240\gamma
h^{4}d^{4}+30h^{3}d/\mu\,.$

The proof is found in (Ding and Li, 2020). As in the case of an overdamped
algorithm, one needs to tune $a$ and $c_{p}$ in the lemma above to obtain
Theorem 13. The result here is also similar to the previous ones. $T^{m}_{2}$
term enters in the Lyapunov function with an $h^{2}$ rate and is negligible
and the iteration formula is dominated by the $D_{3}$ term in (102). There are
two main terms here, the $h^{4}d^{4}$ mainly comes from
$\mathbb{E}|E^{m}|^{2}$ while the $h^{3}d$ is from the original iteration
formula (84). It is hard to compare which term dominates. Indeed with smaller
$\epsilon$, the requirement from the $h^{3}d$ term is more restrictive. This
gives the last two terms in (39).

Acknowledgments

Q.L. acknowledges support from Vilas Early Career award. The research is
supported in part by NSF via grant DMS-1750488, TRIPODS 1740707 and 2023239
and Office of the Vice Chancellor for Research and Graduate Education at the
University of Wisconsin Madison with funding from the Wisconsin Alumni
Research Foundation. Both authors appreciate valuable suggestions from the two
anonymous reviewers and discussions with Stephen J. Wright.

##  A Proof of Lemma 17

We prove (62) in this section. Recall the definition of $E^{m}$ to be

$E^{m}=\nabla f(x^{m})-F^{m}\,.$

To bound variance of $E^{m}$, we consider each component

$\mathbb{E}(|E^{m}|^{2})=\mathbb{E}\left(\mathbb{E}_{r}(|\partial_{i}f(x^{m})-d\partial_{r}f(x^{m})\textbf{e}^{r}_{i}|^{2})\right)=(d-1)\mathbb{E}|\partial_{i}f(x^{m})|^{2}\,,$

where we use $r$ is uniformly chosen from $1,2,\dots,d$, this implies

$\mathbb{E}(|E^{m}_{i}|^{2})=\sum^{d}_{i=1}\mathbb{E}(|E^{m}_{i}|^{2})<\mathbb{E}|\nabla
f(x^{m})|^{2}d\,.$ (104)

Since gradient of $f$ is $L$-Lipschitz, we have

$\displaystyle|\nabla f(x^{m})|^{2}$ $\displaystyle\leq 2|\nabla
f(x^{m})-\nabla f(y^{m})|^{2}+2|\nabla f(y^{m})|^{2}$ (105) $\displaystyle\leq
2L^{2}|x^{m}-y^{m}|^{2}+2|\nabla f(y^{m})|^{2}\,,$

and

$\mathbb{E}|\nabla f(y^{m})|^{2}\leq dL\,,$ (106)

where we use $y_{t}\sim p$ for any $t$ and Lemma 3 in (Dalalyan and
Karagulyan, 2019) in the last inequality. Plug (105),(106) into (104), we have

$\mathbb{E}(|E^{m}|^{2})<2dL^{2}\mathbb{E}|\Delta^{m}|^{2}+2d^{2}L\,,$ (107)

which proves (62).

##  B Proof of Lemma 20

Let

$\widetilde{x}=x^{s\tau},\quad\widetilde{y}=y^{s\tau}$

For $m=s\tau$, according to (27),

$\mathbb{E}\left|E^{m}\right|^{2}=(d-1)\mathbb{E}\left|\nabla f(x^{m})-\nabla
f(\widetilde{x})\right|^{2}=0\,.$

Therefore, we only need to consider $m>s\tau$, according to (27), we directly
have

$\mathbb{E}\left|E^{m}\right|^{2}=(d-1)\mathbb{E}\left|\nabla f(x^{m})-\nabla
f(\widetilde{x})\right|^{2}\leq
dL^{2}\mathbb{E}\left|x^{m}-\widetilde{x}\right|^{2}\,,$ (108)

where the last inequality comes from Hölder’s inequality. Then, use Hölder’s
inequality again, we have

$\mathbb{E}\left|x^{m}-\widetilde{x}\right|^{2}\leq
3\mathbb{E}\left|x^{m}-y^{m}\right|^{2}+3\mathbb{E}\left|y^{m}-\widetilde{y}\right|^{2}+3\mathbb{E}\left|\widetilde{y}-\widetilde{x}\right|^{2}\,,$
(109)

The second term can be further bounded by (50)

$\mathbb{E}\left|y^{m}-\widetilde{y}\right|^{2}=\mathbb{E}\left|\int^{mh}_{s\tau
h}\nabla
f(y_{s})\,\mathrm{d}s-\sqrt{2h}\sum^{m-1}_{j=s\tau}\xi^{j}\right|^{2}\leq
2h^{2}\tau^{2}\mathbb{E}_{p}|\nabla f(y)|^{2}+4h\tau d\,,$ (110)

where we use Hölder’s inequality and $y_{t}\sim p$ for any $t$ in the last
inequality. Finally, according to Lemma 3 in (Dalalyan and Karagulyan, 2019)
to obtain

$\mathbb{E}_{p}|\nabla f(y)|^{2}\leq Ld$

Plug this into (110), combine (108),(109), we obtain (67).

##  C Proof of Lemma 26

The proof is similar to Appendix B

For $m=s\tau$, according to (27),

$\mathbb{E}\left|E^{m}\right|^{2}=(d-1)\mathbb{E}\left|\nabla f(x^{m})-\nabla
f(x^{s\tau})\right|^{2}=0\,.$

Therefore, we only need to consider $m>s\tau$, according to (27), we directly
have

$\mathbb{E}\left|E^{m}\right|^{2}=(d-1)\mathbb{E}\left|\nabla f(x^{m})-\nabla
f(x^{s\tau})\right|^{2}\leq
dL^{2}\mathbb{E}\left|x^{m}-x^{s\tau}\right|^{2}\,,$ (111)

where the last inequality comes from Hölder’s inequality. Then, use Hölder’s
inequality again, we have

$\mathbb{E}\left|x^{m}-x^{s\tau}\right|^{2}\leq
3\mathbb{E}\left|x^{m}-\widetilde{x}^{m}\right|^{2}+3\mathbb{E}\left|\widetilde{x}^{m}-\widetilde{x}^{s\tau}\right|^{2}+3\mathbb{E}\left|\widetilde{x}^{s\tau}-x^{s\tau}\right|^{2}\,,$
(112)

The second term can be further bounded by (82)

$\mathbb{E}\left|\widetilde{x}^{m}-\widetilde{x}^{s\tau}\right|^{2}=\mathbb{E}\left|\int^{mh}_{s\tau
h}\widetilde{\mathrm{V}}_{s}\,\mathrm{d}s\right|^{2}\leq
2h^{2}\tau^{2}\mathbb{E}_{p_{2}}|\widetilde{\mathrm{V}}|^{2}\,,$ (113)

where we use Hölder’s inequality and $\widetilde{\mathrm{V}}_{t}\sim p$ for
any $t$ in the last inequality. Finally, since we have

$\mathbb{E}_{p}|\widetilde{\mathrm{V}}|^{2}=\gamma d\,,$

plug this into (110), combine (108),(109), we obtain (95).

##  D Key lemmas in the proof of Section 8

In this section, we always assume $h<1$ and $f$ satisfies assumption 3.1, and
all notations come from Section 8.

###### Proposition 29

$\mathrm{J}^{r,m}_{1},\mathrm{J}^{m}_{2}$ are defined in (85) and (86), then
for any $a>0$, we have

$\mathbb{E}\left|\mathrm{J}^{r,m}_{1}\right|^{2}+\left|\mathrm{J}^{m}_{2}\right|^{2}\leq
C_{1}\mathbb{E}|\Delta^{m}|^{2}+5(1+1/a)\gamma^{2}h^{4}\mathbb{E}(|E^{m}|^{2})+5(1+1/a)\gamma
h^{4}d\,,$ (114)

where

$C_{1}=(1+a)[1-h/\kappa+200h^{2}]+80(1+1/a)h^{4}\,.$

We first introduce four new terms to prove Proposition 29. Denote

$\displaystyle A^{m}=$
$\displaystyle(\widetilde{v}^{m}-v^{m})(h+e^{-2h})+(\widetilde{x}^{m}-x^{m})$
(115) $\displaystyle-\gamma\int^{(m+1)h}_{mh}e^{-2((m+1)h-s)}\left[\nabla
f(\widetilde{x}^{m})-\nabla f(x^{m})\right]\,\mathrm{d}s\,,$ $\begin{aligned}
B^{m}=&\int^{(m+1)h}_{mh}\widetilde{\mathrm{V}}_{s}-\mathrm{V}_{s}-(\widetilde{v}^{m}-v^{m})\,\mathrm{d}s\\\
&-\gamma\int^{(m+1)h}_{mh}e^{-2((m+1)h-s)}\left[\nabla
f\left(\widetilde{\mathrm{X}}_{s}\right)-\nabla
f(\widetilde{x}^{m})\right]\,\mathrm{d}s\end{aligned}\,,$ (116)
$C^{m}=(\widetilde{x}^{m}-x^{m})+\int^{(m+1)h}_{mh}\widetilde{v}^{m}-v^{m}\,\mathrm{d}s=(\widetilde{x}^{m}-x^{m})+h(\widetilde{v}^{m}-v^{m})\,,$
(117)
$D^{m}=\int^{(m+1)h}_{mh}\widetilde{\mathrm{V}}_{s}-\mathrm{V}_{s}-(\widetilde{v}^{m}-v^{m})\,\mathrm{d}s\,,$
(118)

Then, we have the following three lemmas:

###### Lemma 30

For any $m\geq 0$

$\mathbb{E}\int^{(m+1)h}_{mh}\left|\widetilde{\mathrm{X}}_{t}-\widetilde{x}^{m}\right|^{2}\,\mathrm{d}t\leq\frac{h^{3}\gamma
d}{3}$ (119)

and

$\displaystyle\mathbb{E}\int^{(m+1)h}_{mh}\left|\left(\widetilde{\mathrm{V}}_{t}-\mathrm{V}_{t}\right)-\left(\widetilde{v}^{m}-v^{m}\right)\right|^{2}\,\mathrm{d}t\leq$
$\displaystyle
16h^{3}\mathbb{E}|\Delta^{m}|^{2}+\gamma^{2}h^{3}\mathbb{E}(|E^{m}|^{2})+0.4\gamma
h^{5}d\,,$ (120)

###### Lemma 31

$B^{m},D^{m}$ are defined in (116),(118), we have for any $m\geq 0$

$\mathbb{E}|B^{m}|^{2}\leq
32h^{4}\mathbb{E}|\Delta^{m}|^{2}+2\gamma^{2}h^{4}\mathbb{E}(|E^{m}|^{2})+2\gamma
h^{4}d$ (121) $\mathbb{E}|D^{m}|^{2}\leq
16h^{4}\mathbb{E}|\Delta^{m}|^{2}+\gamma^{2}h^{4}\mathbb{E}(|E^{m}|^{2})+0.4\gamma
h^{6}d$ (122)

###### Lemma 32

$A^{m},C^{m}$ defined in (115),(117), there exists a uniform constant $D$ such
that for any $m\geq 0$

$\mathbb{E}(|A^{m}|^{2}+|C^{m}|^{2})\leq\left[1-h/\kappa+200h^{2}\right]\mathbb{E}|\Delta^{m}|^{2}\,$
(123)

where $\kappa=L/\mu$ is the condition number of $f$.

Now, we are ready to prove Proposition 29:

Proof [Proof of Proposition 29] First, we notice that

$\mathrm{J}^{r,m}_{1}=A^{m}+B^{m},\quad\mathrm{J}^{m}_{2}=C^{m}+D^{m}\,.$

By Young’s inequality, we have

$\displaystyle\mathbb{E}|\mathrm{J}^{r,m}_{1}|^{2}+\mathbb{E}|\mathrm{J}^{m}_{2}|^{2}=$
$\displaystyle\mathbb{E}|A^{m}+B^{m}|^{2}+\mathbb{E}|C^{m}+D^{m}|^{2}$ (124)
$\displaystyle\leq$
$\displaystyle(1+a)\left(\mathbb{E}|A^{m}|^{2}+\mathbb{E}|C^{m}|^{2}\right)$
$\displaystyle+(1+1/a)(\mathbb{E}|B^{m}|^{2}+\mathbb{E}|D^{m}|^{2})\,,$

where $a>0$ will be carefully chosen later. Now, the first term of (124) only
contains information from previous step, using $f$ is strongly convex, we can
bound it using $\left|\Delta^{m}\right|^{2}$ (shown in Appendix D Lemma 32).
To bound the second term, we need to consider difference between $x,v$ at
$t_{m+1}$ and $t_{m}$, which can be bounded by $|\Delta^{m}|^{2}$ and
$|E^{m}|^{2}$ (shown in Appendix D Lemma 31).

According to Lemma 31-32, we first have

$\displaystyle\mathbb{E}|\mathrm{J}^{r,m}_{1}|^{2}+\mathbb{E}|\mathrm{J}^{m}_{2}|^{2}$
$\displaystyle\leq$
$\displaystyle(1+a)\left[1-h/\kappa+200h^{2}\right]\mathbb{E}|\Delta^{m}|^{2}$
$\displaystyle+(1+1/a)\left[80h^{4}\mathbb{E}|\Delta^{m}|^{2}+5\gamma^{2}h^{4}\mathbb{E}(|E^{m}|^{2})+5\gamma
h^{4}d\right]$ $\displaystyle=$ $\displaystyle
C_{1}\mathbb{E}|\Delta^{m}|^{2}+5(1+1/a)\gamma^{2}h^{4}\mathbb{E}(|E^{m}|^{2})+5(1+1/a)\gamma
h^{4}d\,,$

where in the first inequality we use $1+h^{2}<2$ and

$C_{1}=(1+a)[1-h/\kappa+200h^{2}]+80(1+1/a)h^{4}\,.$

To complete the proof, we prove three Lemmas one by one:

Proof [Proof of Lemma 30] First we prove (119). According to (82), we have

$\displaystyle\mathbb{E}\int^{(m+1)h}_{mh}\left|\widetilde{\mathrm{X}}_{t}-\widetilde{x}^{m}\right|^{2}dt$
$\displaystyle=\mathbb{E}\int^{(m+1)h}_{mh}\left|\int^{t}_{mh}\widetilde{\mathrm{V}}_{s}ds\right|^{2}dt$
(125)
$\displaystyle\leq\int^{(m+1)h}_{mh}(t-mh)\int^{t}_{mh}\mathbb{E}\left|\widetilde{\mathrm{V}}_{s}\right|^{2}dsdt$
$\displaystyle=\int|v|^{2}p_{2}(x,v)\,\mathrm{d}x\,\mathrm{d}v\int^{(m+1)h}_{mh}(t-mh)^{2}dt=\frac{h^{3}\gamma
d}{3}\,,$

where in the first inequality we use Hölder’s inequality, and for the second
equality we use $p_{2}$ is a stationary distribution so that
$\left(\widetilde{\mathrm{X}}_{t},\widetilde{\mathrm{V}}_{t}\right)\sim p_{2}$
and $\widetilde{\mathrm{V}}_{t}\sim\exp(-|v|^{2}/(2\gamma))$ for any $t$.

Second, to prove (120), using (79),(82), we first rewrite
$\left(\widetilde{\mathrm{V}}_{t}-\mathrm{V}_{t}\right)-\left(\widetilde{v}^{m}-v^{m}\right)$
as

$\displaystyle\left(\widetilde{\mathrm{V}}_{t}-\mathrm{V}_{t}\right)-\left(\widetilde{v}^{m}-v^{m}\right)=$
$\displaystyle\left(\widetilde{v}^{m}-v^{m}\right)(e^{-2(t-mh)}-1)$ (126)
$\displaystyle-\gamma\int^{t}_{mh}e^{-2(t-s)}\left[\nabla
f(\widetilde{\mathrm{X}}_{s})-\nabla f(x^{m})\right]\,\mathrm{d}s$
$\displaystyle+\gamma\int^{t}_{mh}e^{-2(t-s)}\,\mathrm{d}sE^{m}$
$\displaystyle=$
$\displaystyle\mathrm{I}(t)+\mathrm{II}(t)+\mathrm{III}(t)\,.$

for $mh\leq t\leq(m+1)h$. Then we bound each term separately:

  * •

$\displaystyle\mathbb{E}\int^{(m+1)h}_{mh}\left|\mathrm{I}(t)\right|^{2}\,\mathrm{d}t$
$\displaystyle\leq
h\mathbb{E}\int^{(m+1)h}_{mh}\left|\left(\widetilde{v}^{m}-v^{m}\right)(e^{-2(t-mh)}-1)\right|^{2}\,\mathrm{d}t$
(127) $\displaystyle\leq
h\int^{(m+1)h}_{mh}(2(t-mh))^{2}\mathbb{E}\left|\widetilde{v}^{m}-v^{m}\right|^{2}\,\mathrm{d}t$
$\displaystyle\leq\frac{4h^{3}}{3}\mathbb{E}\left|\widetilde{v}^{m}-v^{m}\right|^{2}\,,$

where we use Hölder’s inequality in the first inequality and $1-e^{-x}<x$ in
the second inequality.

  * •

$\displaystyle\mathbb{E}\int^{(m+1)h}_{mh}\left|\mathrm{II}(t)\right|^{2}\,\mathrm{d}t\leq\gamma^{2}\mathbb{E}\int^{(m+1)h}_{mh}\left|\int^{t}_{mh}e^{-2(t-s)}\left[\nabla
f(\widetilde{\mathrm{X}}_{s})-\nabla
f(x^{m})\right]\,\mathrm{d}s\right|^{2}\,\mathrm{d}t$ (128)
$\displaystyle\leq$ $\displaystyle
2\gamma^{2}\mathbb{E}\int^{(m+1)h}_{mh}\left|\int^{t}_{mh}e^{-2(t-s)}\left[\nabla
f(\widetilde{\mathrm{X}}_{s})-\nabla
f(\widetilde{x}^{m})\right]\,\mathrm{d}s\right|^{2}\,\mathrm{d}t$
$\displaystyle+2\gamma^{2}\mathbb{E}\int^{(m+1)h}_{mh}\left|\int^{t}_{mh}e^{-2(t-s)}\left[\nabla
f(\widetilde{x}^{m})-\nabla
f(x^{m})\right]\,\mathrm{d}s\right|^{2}\,\mathrm{d}t$ $\displaystyle\leq$
$\displaystyle
2\gamma^{2}\int^{(m+1)h}_{mh}(t-mh)\mathbb{E}\int^{t}_{mh}\left|\nabla
f(\widetilde{\mathrm{X}}_{s})-\nabla
f(\widetilde{x}^{m})\right|^{2}\,\mathrm{d}s\,\mathrm{d}t$
$\displaystyle+2\gamma^{2}\int^{(m+1)h}_{mh}(t-mh)\mathbb{E}\int^{t}_{mh}\left|\nabla
f(\widetilde{x}^{m})-\nabla f(x^{m})\right|^{2}\,\mathrm{d}s\,\mathrm{d}t$
$\displaystyle\leq$ $\displaystyle
2\gamma^{2}L^{2}\int^{(m+1)h}_{mh}(t-mh)\mathbb{E}\int^{t}_{mh}\left|\widetilde{\mathrm{X}}_{s}-\widetilde{x}^{m}\right|^{2}\,\mathrm{d}s\,\mathrm{d}t$
$\displaystyle+2\gamma^{2}L^{2}\int^{(m+1)h}_{mh}(t-mh)\mathbb{E}\int^{t}_{mh}\left|\widetilde{x}^{m}-x^{m}\right|^{2}\,\mathrm{d}s\,\mathrm{d}t$
$\displaystyle\leq$ $\displaystyle
2\gamma^{3}L^{2}d\int^{(m+1)h}_{mh}\frac{(t-mh)^{4}}{3}\,\mathrm{d}t+2\gamma^{2}L^{2}\int^{(m+1)h}_{mh}(t-mh)^{2}\,\mathrm{d}t\mathbb{E}\left|\widetilde{x}^{m}-x^{m}\right|^{2}$
$\displaystyle\leq$
$\displaystyle\frac{2\gamma^{3}L^{2}h^{5}d}{15}+\frac{2\gamma^{2}L^{2}h^{3}}{3}\mathbb{E}\left|\widetilde{x}^{m}-x^{m}\right|^{2}\,,$

where in the third inequality we use the gradient of $f$ is $L$-Lipschitz
function and we use (119) in the fourth inequality.

  * •

$\displaystyle\mathbb{E}\int^{(m+1)h}_{mh}\left|\mathrm{III}(t)\right|^{2}\,\mathrm{d}t$
$\displaystyle=\gamma^{2}\mathbb{E}\int^{(m+1)h}_{mh}\left|\int^{t}_{mh}e^{-2(t-s)}\,\mathrm{d}sE^{m}\right|^{2}\,\mathrm{d}t$
(129)
$\displaystyle\leq\gamma^{2}\int^{(m+1)h}_{mh}(t-mh)^{2}\,\mathrm{d}t\mathbb{E}(|E^{m}|^{2})$
$\displaystyle\leq\frac{\gamma^{2}h^{3}}{3}\mathbb{E}(|E^{m}|^{2})\,,$

Plug (127),(128),(129) into (126) and using $\gamma L=1$, we have

$\displaystyle\mathbb{E}\int^{(m+1)h}_{mh}\left|\left(\widetilde{\mathrm{V}}_{t}-\mathrm{V}_{t}\right)-\left(\widetilde{v}^{m}-v^{m}\right)\right|^{2}\,\mathrm{d}t$
$\displaystyle\leq$ $\displaystyle
3\left(\mathbb{E}\int^{(m+1)h}_{mh}\left|\mathrm{I}(t)\right|^{2}\,\mathrm{d}t+\mathbb{E}\int^{(m+1)h}_{mh}\left|\mathrm{II}(t)\right|^{2}\,\mathrm{d}t+\mathbb{E}\int^{(m+1)h}_{mh}\left|\mathrm{III}(t)\right|^{2}\,\mathrm{d}t\right)$
$\displaystyle\leq$ $\displaystyle
4h^{3}\left(\mathbb{E}\left|\widetilde{x}^{m}-x^{m}\right|^{2}+\mathbb{E}\left|\widetilde{v}^{m}-v^{m}\right|^{2}\right)+\gamma^{2}h^{3}\mathbb{E}(|E^{m}|^{2})+0.4\gamma
h^{5}d\,,$

using (80), we get the desired result.

Proof [Proof of Lemma 31] First, we separate $B^{m}$ into two parts:

$\displaystyle\mathbb{E}|B^{m}|^{2}\leq$ $\displaystyle
2\mathbb{E}\left|\int^{(m+1)h}_{mh}\left(\widetilde{\mathrm{V}}_{t}-\mathrm{V}_{t}\right)-\left(\widetilde{v}^{m}-v^{m}\right)\,\mathrm{d}t\right|^{2}$
$\displaystyle+2\mathbb{E}\left|\gamma\int^{(m+1)h}_{mh}e^{-2((m+1)h-t)}\left[\nabla
f(\widetilde{\mathrm{X}}_{t})-\nabla
f(\widetilde{x}^{m})\right]\,\mathrm{d}t\right|^{2}\,.$

And each terms can be bounded:

  * •

$\displaystyle\mathbb{E}\left|\int^{(m+1)h}_{mh}\left(\widetilde{\mathrm{V}}_{t}-\mathrm{V}_{t}\right)-\left(\widetilde{v}^{m}-v^{m}\right)\,\mathrm{d}t\right|^{2}$
(130) $\displaystyle\leq$ $\displaystyle
h\mathbb{E}\int^{(m+1)h}_{mh}\left|\left(\widetilde{\mathrm{V}}_{t}-\mathrm{V}_{t}\right)-\left(\widetilde{v}^{m}-v^{m}\right)\right|^{2}\,\mathrm{d}t$
$\displaystyle\leq$ $\displaystyle
16h^{4}\mathbb{E}|\Delta^{m}|^{2}+\gamma^{2}h^{4}\mathbb{E}(|E^{m}|^{2})+0.4\gamma
h^{6}d\,,$

where we use Lemma 30 (120) in the second inequality.

  * •

$\displaystyle\mathbb{E}\left|\gamma\int^{(m+1)h}_{mh}e^{-2((m+1)h-t)}\left[\nabla
f(\widetilde{\mathrm{X}}_{t})-\nabla
f(\widetilde{x}^{m})\right]\,\mathrm{d}t\right|^{2}$ (131) $\displaystyle\leq$
$\displaystyle
h\gamma^{2}\mathbb{E}\int^{(m+1)h}_{mh}\left|e^{-2((m+1)h-t)}\left[\nabla
f(\widetilde{\mathrm{X}}_{t})-\nabla
f(\widetilde{x}^{m})\right]\right|^{2}\,\mathrm{d}t$ $\displaystyle\leq$
$\displaystyle
h\gamma^{2}L^{2}\mathbb{E}\int^{(m+1)h}_{mh}\left|\widetilde{\mathrm{X}}_{t}-\widetilde{x}^{m}\right|^{2}\,\mathrm{d}t$
$\displaystyle\leq$
$\displaystyle\frac{h^{4}\gamma^{3}L^{2}d}{3}\leq\frac{h^{4}\gamma d}{3}\,,$

where we use Lemma 30 (119) and $\gamma L=1$ in the last two inequalities.

Combine (130),(131) together, we finally have

$\mathbb{E}|B|^{2}\leq
32h^{4}\mathbb{E}|\Delta^{m}|^{2}+2\gamma^{2}h^{4}\mathbb{E}(|E^{m}|^{2})+0.8h^{6}\gamma
d+2h^{4}\gamma d/3\,,$

which implies (121) if we further use $h<1$.

Next, estimation of $\left(\mathbb{E}|D|^{2}\right)^{1/2}$ is a direct result
of (130).

Proof [Proof of Lemma 32] Let $\widetilde{x}^{m}-x^{m}=a$ and
$\widetilde{w}^{m}-w^{m}=b$. First, by the mean-value theorem, there exists a
matrix $H$ such that $\mu{I}_{d}\preceq H\preceq L{I}_{d}$ and

$\nabla f(\widetilde{x}^{m})-\nabla f(x^{m})=Ha\,.$

By calculation,
$\int^{(m+1)h}_{mh}e^{-2((m+1)h-t)}\,\mathrm{d}t=\frac{1-e^{-2h}}{2}$ and

$\begin{aligned}
A^{m}&=(h+e^{-2h})(\widetilde{v}^{m}-v^{m})+\left(I_{d}-\frac{(1-e^{-2h})}{2}\gamma
H\right)(\widetilde{x}^{m}-x^{m})\\\
&=\left(\left(1-h-e^{-2h}\right)I_{d}-\frac{(1-e^{-2h})}{2}\gamma
H\right)a+(h+e^{-2h})b\end{aligned}\,.$ $C^{m}=(1-h)a+hb\,.$

Since $\|\gamma H\|_{2}\leq 1$ and we also have the following calculation

$|1-e^{-2h}-2h|\leq 5h^{2}\,,$

where we use $h<\frac{1}{1648}.$

If we further define matrix $\mathcal{M}_{A}$ and $\mathcal{M}_{C}$ such that

$|A^{m}|^{2}=\left(a,b\right)^{\top}\mathcal{M}_{A}\left(a,b\right)\,,\quad|C^{m}|^{2}=\left(a,b\right)^{\top}\mathcal{M}_{C}\left(a,b\right)\,,$

then, we have

$\left\|\mathcal{M}_{A}-\begin{bmatrix}0&hI_{d}-\gamma hH\\\ hI_{d}-\gamma
hH&(1-2h)I_{d}\end{bmatrix}\right\|_{2}\leq 100h^{2}\,,$

and

$\left\|\mathcal{M}_{B}-\begin{bmatrix}(1-2h)I_{d}&hI_{d}\\\
hI_{d}&0\end{bmatrix}\right\|_{2}\leq 100h^{2}\,,$

where we use $h<1/1648$ by (38) and $\|\gamma H\|_{2}\leq 1$. This further
implies

$|A^{m}|^{2}+|C^{m}|^{2}=\left(a,b\right)^{\top}\begin{bmatrix}(1-2h)I_{d}&2hI_{d}-\gamma
hH\\\ 2hI_{d}-\gamma
hH&(1-2h)I_{d}\end{bmatrix}\left(a,b\right)+h^{2}\left(a,b\right)^{\top}Q\left(a,b\right)$

where $\|Q\|_{2}\leq 200$. Calculate the eigenvalue of the dominating matrix
(first term), we need to solve

$\mathrm{det}\left\\{(1-2h-\lambda)^{2}I_{d}-(2hI_{d}-\gamma
hH)^{2}\right\\}=0\,,$

which implies eigenvalues $\\{\lambda_{j}\\}^{d}_{j=1}$ solve

$(1-2h-\lambda_{j})^{2}-(2h-\gamma h\Lambda_{j})^{2}=0\,,$

where $\Lambda_{j}$ is $j$-th eigenvalue of $H$. Since
$\gamma\Lambda_{j}\leq\gamma L=1$ and $h<1$, we have

$\lambda_{j}\leq 1-\gamma\Lambda_{j}h\leq 1-\mu h\gamma=1-h/\kappa$

for each $j=1,\dots,d$. This implies

$\left\|\begin{bmatrix}(1-2h)I_{d}&2hI_{d}-\gamma hH\\\ 2hI_{d}-\gamma
hH&(1-2h)I_{d}\end{bmatrix}\right\|_{2}\leq 1-h/\kappa\,,$

and

$|A^{m}|^{2}+|C^{m}|^{2}\leq(1-h/\kappa+200h^{2})(|a|^{2}+|b|^{2})\,.$

Taking expectation on both sides, we obtain (123).

## References

  * Andrieu et al. (2003) Christophe Andrieu, Nando de Freitas, Arnaud Doucet, and Michael I. Jordan.  An introduction to MCMC for machine learning.  _Machine Learning_ , 50:5–43, 2003.
  * Baker et al. (2019) Jack Baker, Paul Fearnhead, Emily Fox, and Christopher Nemeth.  Control variates for stochastic gradient MCMC.  _Statistics and Computing_ , 29, 05 2019.
  * Baydin et al. (2018) Atilim Gunes Baydin, Barak A. Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind.  Automatic differentiation in machine learning: a survey.  _Journal of Machine Learning Research_ , 18(153):1–43, 2018.
  * Bottou (2010) Léon Bottou.  Large-scale machine learning with stochastic gradient descent.  In _Proceedings of COMPSTAT’2010_ , pages 177–186, 2010.
  * Chapelle et al. (2006) Olivier Chapelle, Bernhard Schölkopf, and Alexander Zien.  _Semi-Supervised Learning_.  The MIT Press, Cambridge, MA, USA, 2006.
  * Chatterji et al. (2018) Niladri S. Chatterji, Nicolas Flammarion, Yi-An Ma, Peter L. Bartlett, and Michael I. Jordan.  On the theory of variance reduction for stochastic gradient Monte Carlo.  In _Proceedings of the 35th International Conference on Machine Learning_ , pages 764–773, 2018.
  * Chen et al. (2014) Tianqi Chen, Emily B. Fox, and Carlos Guestrin.  Stochastic gradient Hamiltonian Monte Carlo.  In _Proceedings of the 31st International Conference on International Conference on Machine Learning_ , pages 1683–1691, 2014.
  * Cheng et al. (2018) Xiang Cheng, Niladri S. Chatterji, Peter L. Bartlett, and Michael I. Jordan.  Underdamped Langevin MCMC: A non-asymptotic analysis.  In _Proceedings of the 31st Conference On Learning Theory_ , pages 300–323, 2018.
  * Dalalyan and Karagulyan (2019) 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.
  * Dalalyan and Riou-Durand (2020) Arnak S. Dalalyan and Lionel Riou-Durand.  On sampling from a log-concave density using kinetic Langevin diffusions.  _Bernoulli_ , 26(3):1956–1988, 2020.
  * Defazio et al. (2014) Aaron Defazio, Francis Bach, and Simon Lacoste-Julien.  SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives.  In _Proceedings of the 27th International Conference on Neural Information Processing Systems_ , pages 1646–1654, 2014.
  * Ding and Li (2020) Zhiyan Ding and Qin Li.  Variance reduction for random coordinate descent-Langevin Monte Carlo.  In _Proceedings of the 34th International Conference on Neural Information Processing Systems_ , pages 3748–3760, 2020.
  * Ding et al. (2021a) Zhiyan Ding, Qin Li, Jianfeng Lu, and Stephen J. Wright.  Random coordinate underdamped langevin monte carlo.  In _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , pages 2701–2709, 2021a.
  * Ding et al. (2021b) Zhiyan Ding, Qin Li, Jianfeng Lu, and Stephen J. Wright.  Random coordinate langevin monte carlo.  In _Proceedings of 34th Conference on Learning Theory_ , pages 1683–1710, 2021b.
  * Duane et al. (1987) Simon Duane, A.D. Kennedy, Brian J.Pendleton, and Duncan Roweth.  Hybrid Monte Carlo.  _Physics Letters B_ , 195(2):216–222, 1987.
  * Dubey et al. (2016) Kumar Avinava Dubey, Sashank J. Reddi, Sinead A Williamson, Barnabas Poczos, Alexander J Smola, and Eric P Xing.  Variance reduction in stochastic gradient Langevin dynamics.  In _Advances in Neural Information Processing Systems_ , volume 29, 2016.
  * Durmus and Moulines (2019) Alain Durmus and Eric Moulines.  High-dimensional Bayesian inference via the unadjusted Langevinalgorithm.  _Bernoulli_ , 25(4A):2854–2882, 2019.
  * Durmus et al. (2016) Alain Durmus, Umut Simsekli, Eric Moulines, Roland Badeau, and Gaël RICHARD.  Stochastic gradient richardson-romberg Markov Chain Monte Carlo.  In _Advances in Neural Information Processing Systems_ , volume 29, 2016.
  * Durmus et al. (2019) Alain Durmus, Szymon Majewski, and Blażej Miasojedow.  Analysis of Langevin Monte Carlo via convex optimization.  _Journal of Machine Learning Research_ , 20(73):1–46, 2019.
  * Dwivedi et al. (2018) Raaz Dwivedi, Yuansi Chen, M. Wainwright, and Bin Yu.  Log-concave sampling: Metropolis-Hastings algorithms are fast!  _Journal of Machine Learning Research_ , 20(183):1–42, 2018.
  * Eberle et al. (2019) Andreas Eberle, Arnaud Guillin, and Raphael Zimmer.  Couplings and quantitative contraction rates for Langevin dynamics.  _Annals of Probability_ , 47(4):1982–2010, 2019\.
  * Fabian (1981) Peter Fabian.  Atmospheric sampling.  _Advances in Space Research_ , 1(11):17–27, 1981\.
  * Garbuno-Inigo et al. (2020) Alfredo Garbuno-Inigo, Franca Hoffmann, Wuchen Li, and Andrew M. Stuart.  Interacting Langevin diffusions: gradient structure and ensemble Kalman sampler.  _SIAM Journal on Applied Dynamical Systems_ , 19(1):412–441, 2020.
  * Geman and Geman (1984) Stuart Geman and Donald Geman.  Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 6(6):721–741, 1984.
  * Gerencsér (1997) László Gerencsér.  Rate of convergence of moments of Spall’s SPSA method.  In _1997 European Control Conference (ECC)_ , pages 2192–2197, 1997\.
  * Geweke (1989) John Geweke.  Bayesian inference in econometric models using Monte Carlo integration.  _Econometrica_ , 57(6):1317–1339, 1989.
  * Hastings (1970) Wilfred Keith Hastings.  Monte Carlo sampling methods using Markov chains and their applications.  _Biometrika_ , 57(1):97–109, 1970.
  * Hsieh et al. (2008) Cho-Jui Hsieh, Kai-Wei Chang, Chih-Jen Lin, S. Sathiya Keerthi, and S. Sundararajan.  A dual coordinate descent method for large-scale linear svm.  In _Proceedings of the 25th International Conference on Machine Learning_ , pages 408–415, 2008.
  * Iglesias et al. (2013) Marco A. Iglesias, Kody J. H. Law, and Andrew M. Stuart.  Ensemble Kalman methods for inverse problems.  _Inverse Problems_ , 29(4):045001, 2013.