Year: 2021

Abstract:
We present a new class of Langevin based algorithms, which overcomes many of the known shortcomings of popular adaptive optimizers that are currently used for the fine tuning of deep learning models. Its underpinning theory relies on recent advances of Euler’s polygonal approximations for stochastic differential equations (SDEs) with monotone coefficients. As a result, it inherits the stability properties of tamed algorithms, while it addresses other known issues, e.g. vanishing gradients in neural networks. In particular, we provide a nonasymptotic analysis and full theoretical guarantees for the convergence properties of an algorithm of this novel class, which we named THεO POULA (or, simply, TheoPouLa). Finally, several experiments are presented with different types of deep learning models, which show the superior performance of TheoPouLa over many popular adaptive optimization algorithms.

heading: 1 INTRODUCTION
text: Modern machine learning models including deep neural networks are successfully trained when they are finely tuned via the optimization of their associated loss functions. Two aspects of such optimization tasks pose significant challenges, namely the non-convex nature of loss functions and the highly nonlinear features of many types of neural networks. Moreover, the analysis in Lovas et al. (2020) shows that the gradients of such non-convex loss functions typically grow faster than linearly and are only locally Lipschitz continuous. Naturally, stability issues are observed, which are known as the ‘exploding gradient’ phenomenon (Bengio et al., 1994; Pascanu et al., 2013), when vanilla stochastic gradient descent (SGDs) or certain types of adaptive algorithms are used for fine tuning. The sparsity of gradients of neural networks is another challenging issue, which is extensively studied in the literature. For example, momentum methods and adaptive learning rate methods such as AdaGrad (Duchi et al. (2011)), RMSProp (Tieleman & Hinton (2012)), Adam (Kingma & Ba (2015)) have been developed to tackle this problem and improve training speed by diagonally scaling the gradient by some function of the past gradients.
A family of Langevin based algorithms has been another important stream of literature on the stochastic optimization. They are built on the theoretical fact that the Langevin stochastic differential equation, (6), converges to its unique invariant measure, which concentrates on the global minimizers of the objective function as β →∞, see Hwang (1980). Since the convergence property remains true for nonconvex optimization problems, the global convergence of the stochastic gradient Langevin dynamics (SGLD) and its variants has been extensively studied in a nonconvex setting (Raginsky et al., 2017; Xu et al., 2018; Erdogdu et al., 2018; Brosse et al., 2018; Lovas et al., 2020). Moreover, it is worth noting that Langevin based algorithms have been a key element in statistics and Bayesian learning (Roberts & Tweedie, 1996; Durmus & Moulines, 2017; Dalalyan, 2017; Brosse et al., 2019; Welling & Teh, 2011; Deng et al., 2020a;b).
Motivated by the aforementioned developments in the field, we propose a new class of Langevin algorithms which is based on recent advances of Euler’s polygonal approximations for Langevin SDEs. The idea of this new form of Euler’s polygonal approximations for SDEs with monotone coefficients originates from the articles Krylov (1985) and Krylov (1990). We name this new class as polygonal unadjusted Langevin algorithms. Moreover, it is versatile enough to incorporate further
features to address other known shortcomings of adaptive optimizers. Mathematically, it is described as follows: Given an i.i.d. sequence of random variables {Xn}n≥0 of interest, which typically represent available data, the algorithm follows
θλ0 := θ0, θ λ n+1 := θ λ n − λGλ(θλn, Xn+1) + √ 2λβ−1ξn+1, n ∈ N, (1)
where θ0 is an Rd-valued random variable, λ > 0 denotes the step size of the algorithm, β > 0 is the so-called inverse temperature, (ξn)n∈N is an Rd-valued Gaussian process with i.i.d. components and Gλ : Rd × Rm → Rd satisfies the following three properties:
1. For every λ > 0, There exist constants Kλ > 0 and ρ1 ≥ 0 such that |Gλ(θ, x)| ≤ Kλ(1 + |x|)ρ1(1 + |θ|) for every θ ∈ Rd and x ∈ Rm.
2. There exist constants γ ≥ 1/2, K2 > 0 and ρ2, ρ3 ≥ 0 such that for all λ > 0, |Gλ(θ, x)−G(θ, x)| ≤ λγK2(1 + |x|)ρ2(1 + |θ|)ρ3
for every θ ∈ Rd and x ∈ Rm, whereG is the (unbiased) stochastic gradient of the objective function of the optimization problem under study.
3. There exist constants λmax and δ ∈ {1, 2} such that for any λ ≤ λmax,
lim inf |θ|→∞ E [ 〈 θ |θ|δ , Gλ(θ,X0)〉 − 2λ |θ|δ |Gλ(θ,X0)|2 ] > 0.
One obtains our new algorithm THεO POULA by considering the case where Gλ(θ, x) is the vector with entries H(i)λ,c(θ, x) as given by (8), for i ∈ {1, . . . , d}. Its name is formed from its description, namely Tamed Hybrid ε-Order POlygonal Unadjusted Langevin Algorithm and its full detailed analysis (including its convergence properties) are given in Section 3. We note that THεO POULA and TUSLA (Lovas et al. (2020)) satisfy the above three properties with δ = 2 and γ = 1/2, whereas TULA (Brosse et al. (2019)) satisfies them with δ = γ = 1 as it assumes only deterministic gradients (and thus the i.i.d. data sequence reduces to a constant).

heading: 1.1 RELATED WORK: LANGEVIN BASED ALGORITHMS AND ADAPTIVE LEARNING RATE METHODS
text: Most research on Langevin based algorithms in the literature has been focused on theoretical aspects. Raginsky et al. (2017) demonstrated the links between Langevin based algorithms and stochastic optimization in neural networks, stimulating further the development and analysis of such algorithms. Xu et al. (2018) analyzed the global convergence of GLD, SGLD and SVRG-LD. The incorporation of dependent data streams in the analysis of SGLD algorithms has been achieved in Barkhagen et al. (2021) and in Chau et al. (2019), and local conditions have been studied in Zhang et al. (2019). Recently, TUSLA of Lovas et al. (2020) has been proposed based on a new generation of tamed Euler approximations for stochastic differential equations (SDEs) with monotone coefficients in nonconvex optimization problems. See Hutzenthaler et al. (2012) and Sabanis (2013) for the rationale of taming techniques. Despite their elegant theoretical results, the use of Langevin based algorithms for training deep learning models has been limited in practice as their empirical performance lacked behind in comparison to other popular adaptive gradient methods. We refer to Appendix F.3 for the reader who is interested in recent progress on sampling and Bayesian neural networks.
Adaptive learning rate methods such as AdaGrad (Duchi et al. (2011)), RMSProp (Tieleman & Hinton (2012)) and Adam (Kingma & Ba (2015)) have been successfully applied to neural network models due to their fast training speed. Since the appearance of Adam, a large number of variants of Adam-type optimizers have been proposed to address the theoretical and practical challenges of Adam. For example, Reddi et al. (2018) provided a simple example that demonstrates the non-convergence issue of Adam and proposed a simple modification, called AMSGrad, to solve this problem. Chen et al. (2019) discussed the convergence of Adam-type optimizers in a nonconvex setting. RAdam to rectify the variance of adaptive learning rate has been proposed in Liu et al. (2020). Wilson et al. (2017) revealed that the generalization ability of adaptive learning rate methods is worse than a global learning method like SGD. AdaBound of Luo et al. (2019) attempts to overcome the drawback by employing dynamic bounds on learning rates. Recently, AdaBelief (Zhuang et al. (2020)) and AdamP (Heo et al. (2021)) demonstrated their fast convergence and good generalization via extensive experiments. Nevertheless, the convergence analysis of these (and other)
adaptive learning rate methods is still restrictive since it is only guaranteed to converge to a stationary point (which can be a local minimum or a saddle point) under strong assumptions. Namely, the stochastic gradient is globally Lipschitz continuous and bounded. Note though that none of these two assumptions hold true in a typical optimization problem involving neural networks. This is particularly evident in complex neural network architectures.

heading: 1.2 OUR CONTRIBUTIONS
text: The proposed algorithm, THεO POULA, tries to combines both advantages: namely, global convergence in Langevin based algorithms and powerful empirical performance in adaptive learning rate methods. To the best of the authors’ knowledge, our algorithm is the first Langevin based algorithm to outperform popular stochastic optimization methods such as SGD, Adam, AMSGrad, RMSProp, AdaBound and AdaBelief for deep learning tasks. The major strengths of our work over related algorithms are summarized as follows:
• (Global convergence) We provide a global convergence analysis of THεO POULA for nonconvex optimization where the stochastic gradient of the objective is locally Lipscthiz continuous. Moreover, non-asymptotic estimates for the expected excess risk are derived.
• (Stable and fast training) THεO POULA achieves a stable and fast training process using the (element-wise) taming technique, (element-wise) boosting function and averaging, which are theoretically well-designed. Furthermore, we validate the effectiveness of the taming and boosting functions through several empirical experiments.
• (Good generalization) While THεO POULA behaves like adaptive learning rate methods in the early training phase, it takes an almost global learning rate near an optimal point. That is, THεO POULA is quickly switched from adaptive methods to SGD. As a result, it inherits the good generalization ability of SGD. Our experiments support this fact by showing that THεO POULA outperforms the other optimization methods in generalization measured by test accuracy for various deep learning tasks.

heading: 2 MOTIVATING EXAMPLE
text: The local Lipschitz continuity of gradients and its effect on the performance of optimization methods are relatively under-studied. Most relevant studies assume that the stochastic gradient is global Lipscthiz continuous and bounded (Kingma & Ba, 2015; Xu et al., 2018; Brosse et al., 2018; Duchi et al., 2011; Tieleman & Hinton, 2012; Reddi et al., 2018; Chen et al., 2019; Liu et al., 2020; Luo et al., 2019; Zhuang et al., 2020) although it is not true for neural network problems. This section provides a simple, one-dimensional optimization problem that illustrates the convergence issue of popular stochastic gradient methods when the stochastic gradient is locally Lipschitz continuous, i.e., the gradient can be super-linearly growing 1.
Consider the following optimization problem:
min θ u(θ) = min θ E[U(θ,X)], (2)
where U : R× R→ R is defined as
U(θ, x) =
{ θ2 (1 + 1x≤1) + θ
30, |θ| ≤ 1, (2|θ| − 1) (1 + 1x≤1) + θ30, |θ| > 1,
and X is uniformly distributed over (−2, 2), that is, fX(x) = 141|x|≤2. Furthermore, the stochastic gradient G : R× R→ R is given by
G(θ, x) =
{ 2θ (1 + 1x≤1) + 30θ
29, |θ| ≤ 1, 2(1 + 1x≤1)sgn(θ) + 30θ 29, |θ| > 1,
where sgn(·) is the sign function. Note that the stochastic gradientG is locally Lipschitz continuous, which satisfies
|G(θ, x)−G(θ′, x)| ≤ 34(1 + |θ|+ |θ′|)28|θ − θ′| 1Lovas et al. (2020) used a similar example to show the stability of TUSLA with a different taming function.
for all x ∈ R and θ, θ′ ∈ R. Also, the optimal value is attained at θ = 0. See Appendix A for more details. Following Reddi et al. (2018), adaptive stochastic gradient methods can be generally written as follows, for n ∈ N,
mn = φn(G1, · · · , Gn), Vn = ψn(G1, · · · , Gn),
θn+1 = θn − λn mn
ε+ √ V n
(3)
where Gi := G(θi, Xi) is the stochastic gradient evaluated at the i-th iteration, λn is the step size and all operations are applied element-wise. Table 1 provides the details for some of the most popular stochastic optimization methods with corresponding averaging functions φn and ψn.
We use SGD, Adam, AMSGrad and RMSprop to solve the optimization problem with initial value θ0 = 5. For hyperparameters of optimization algorithms, we use their default settings provided in PyTorch. Figure 1(a) shows the trajectories of approximate solutions generated by each optimizer. While SGD, Adam, AMSGrad and RMSProp fail to converge to the optimal solution 0, the proposed algorithm, THεO POULA, finds the optimal solution with a reasonable step size, say, 0.01.
Intuitively, the undesirable phenomenon occurs because, in the iterating rule (3), the denominator√ V n excessively dominates the numerator mn, causing the vanishing gradient problem in the presence of the superlinear gradient. On the contrary, SGD suffers from the exploding gradient problem. Moreover, Figure 1(b) highlights that the problematic behavior cannot be simply resolved by adjusting the learning rate within the Adam-type framework, while THεO POULA perform extremely well even in the presence of such violent non-linearities.
3 NEW ALGORITHM: THεO POULA
We propose a new stochastic optimization algorithm by combining ideas from taming methods specifically designed to approximate Langevin SDEs with a hybrid approach based on recent advances of polygonal Euler approximations. The latter is achieved by identifying a suitable boosting
function (of order ε 1) to efficiently deal with the sparsity of (stochastic) gradients of neural networks. In other words, the novelty of our algorithm is to utilize a taming function and a boosting function, rather than designing a new Vn as in Adam-type optimizers.
We proceed with the necessary preliminary information, main assumptions and formal introduction of the new algorithm.

heading: 3.1 PRELIMINARIES AND ASSUMPTIONS
text: Let (Ω,F , P ) be a probability space. We denote by E[X] the expectation of a random variable X . Fix an integer k ≥ 1. For an Rk-valued random variable X , its law on B(Rk), i.e. the Borel sigma-algebra of Rk, is denoted by L(X). Scalar product is denoted by 〈·, ·〉, with | · | standing for the corresponding norm (where the dimension of the space may vary depending on the context). For any integer q ≥ 1, let P(Rq) denote the set of probability measures on B(Rq). For µ, ν ∈ P(Rk), let C(µ, ν) denote the set of probability measures ζ on B(R2k) such that its respective marginals are µ, ν. For two probability measures µ and ν, the Wasserstein distance of order p ≥ 1 is defined as
Wp(µ, ν) := inf ζ∈C(µ,ν) (∫ Rk ∫ Rk |θ − θ′|pζ(dθdθ′) )1/p , µ, ν ∈ P(Rk). (4)
Let (Xn)n∈N be a sequence of i.i.d. Rm-valued random variables generating the filtration (Gn)n∈N and (ξn)n∈N be an Rd-valued Gaussian process with independent components.
Let F : Rd×Rm → Rd be continuously differentiable function such that E[F (θ,X0)] <∞ for any θ ∈ Rd. We consider the following optimization problem
min θ u(θ) = min θ
( E[F (θ,X0)] +
η
2(r + 1) |θ|2(r+1)
) (5)
where θ ∈ Rd, η ∈ (0, 1) is the regularization parameter and r ≥ q2 +1. In the context of fine tuning of neural networks, F represents the loss function for the task at hand and θ denotes the vector of the neural network’s parameters. Note that the regularization term, η2(r+1) |θ|
2(r+1), is added in order to guarantee that the dissipativity property holds, since it is essential for the convergence analysis.
Remark 3.1. For the reader who prefers to consider the optimization problem without the regularization term, i.e. with η = 0, the dissipative condition (B.1) has to be additionally assumed as in the literature (Raginsky et al., 2017; Xu et al., 2018; Erdogdu et al., 2018). Then, the same analysis can be applied to obtain our main results without any additional effort. However, it is yet to be proven theoretically that such an assumption holds in general for neural networks and thus it becomes a case-by-case investigation. In other words, we present here the formal theoretical statement with the appropriate regularization term which covers all of these cases.
In particular, r depends on the neural network’s structure, whereas q is described in Assumption 3.1. Consequently, the stochastic gradient with the regularization term is given by
H(θ, x) := G(θ, x) + ηθ|θ|2r
where G(θ, x) := ∇θF (θ, x) for all x ∈ Rm, θ ∈ Rd and η = 0 if dissipativity holds for G. We introduce our main assumptions. The first requirement is that G is locally Lipschitz continuous.
Assumption 3.1. There exists positive constant L1, ρ and q ≥ 1 such that
|G(θ, x)−G(θ′, x)| ≤ L1(1 + |x|)ρ(1 + |θ|+ |θ′|)q−1|θ − θ′|.
for all x ∈ Rm and θ, θ′ ∈ Rd. Moreover, h(θ) := E[H(θ,X0)] for every θ ∈ Rd.
Further, conditions on the initial value θ0 and data process (Xn)n∈N are imposed as it is common to use weight initialization using the uniform or normal distribution, Assumption 3.2 is mild.
Assumption 3.2. The process (Xn)n∈N is a sequence of i.i.d. random variables with E[|X0|16ρ(2r+1)] < ∞ where ρ is given in Assumption 3.1. In addition, the initial condition is such that E[|θ0|16(2r+1)] <∞.
We refer to Appendix B for further remarks and key observations regarding the consequences of Assumptions 3.1 and 3.2. We conduct the convergence analysis of THεO POULA by employing elements of the theory of Langevin SDEs. It is shown that, under mild conditions (satisfied by Assumptions 3.1 and 3.2), the so-called (overdamped) Langevin SDE, which is given by
dZt = −h(Zt)dt+ √ 2β−1dBt, t > 0, (6)
where h = ∇u with a (possibly random) initial condition θ0 and with (Bt)t≥0 denoting a ddimensional Brownian motion, admits a unique invariant measure πβ(dz) ∝ exp(−βu(z)). Thus, for a sufficiently large β, πβ concentrates around the minimizers of (5).
3.2 MECHANISM OF THεO POULA
We introduce the mechanism of THεO POULA, which iterately updates as follows:
θλ0 := θ0, θ λ n+1 := θ λ n − λHλ,c(θλn, Xn+1) + √ 2λβ−1ξn+1, n ∈ N, (7)
where Hλ,c := (H (1) λ,c(θ, x), · · · , H (d) λ,c(θ, x)) T is given by
H (i) λ,c(θ, x) =
G(i)(θ, x)
1 + √ λ|G(i)(θ, x)|︸ ︷︷ ︸
taming function
( 1 +
√ λ
ε+ |G(i)(θ, x)|︸ ︷︷ ︸ boosting function
) + η θ(i)|θ|2r
1 + √ λ|θ|2r︸ ︷︷ ︸
regularization term
, (8)
and {ξn}n≥1 is a sequence of independent standard d-dimensional Gaussian random variables. Note that the taming and boosting functions are defined in (8).
THεO POULA has several distinct features over the existing optimization methods in the literature. We give an intuitive explanation as to how these features are complementarily harmonized to improve the performance of the algorithm, and to handle the exploding and vanishing gradient problems of neural networks. For simplicity, we omit the regularization term, that is, η = 0, and the noise term, √ 2λβ−1ξn+1, throughout the exposition. Also, we refer to λ as the learning rate and
|∆θλn| := |G(i)(θλn,Xn+1)|
1+ √ λ|G(i)(θλn,Xn+1)|
( 1 +
√ λ
ε+|G(i)(θλn,Xn+1)|
) as the stepsize by the convention in Kingma
& Ba (2015).
Firstly, the new algorithm utilizes the taming function to control the super-linearly growing gradient. In a region where the loss function is steep and narrow (the gradient is huge), it is ideal for the optimizer to take a small stepsize. This is effectively achieved since the growth of the taming function is proportional to G, but the boosting is close to one when the gradient is huge. The effectiveness of the taming function is confirmed in the motivating example in Section 2. Note that the taming function is applied element-wise to scale the effective element-wise learning rate in contrast to TUSLA of Lovas et al. (2020). This significantly improves the performance of our new algorithm in solving high-dimensional optimization problems such as the fine tuning of neural network models.
Secondly, we have designed the boosting function to accelerate training speed and prevent the vanishing gradient problem2. When the current parameter is located in a region where the loss function is flat (the gradient is small), it is desirable for the optimizer to take a large stepsize. As the gradient gets smaller, the boosting function increases the stepsize of THεO POULA by up to √ λ/ε, whereas the taming function’s contribution decreases. As a result, THεO POULA takes a larger stepsize. In other words, THεO POULA takes a desirable stepsize depending on the magnitude of the gradient. Most importantly, the taming and boosting functions do not interfere with each other in any adverse way. On the contrary, they complement each other in a harmonious way that is evident from our simulation results.
Thirdly, THεO POULA is quickly converted from adaptive learning rate methods to SGD. In the early training phase, THεO POULA certainly behaves like adaptive learning rate methods. Then,
2We provide the effectiveness of the boosting function in Appendix E.3 by comparing the performance of THεO POULA with/without the boosting function. The experiment shows that the addition of the boosting function brings a significant improvement in test accuracy across different models and data sets.
when the current position is approaching an optimal solution where |G(i)|s are close to zero, the movement of THεO POULA is similar to SGD with a learning rate (1 + √ λ/ε). Consequently, THεO POULA simultaneously attains two favorable features of fast training in adaptive learning rate methods and good generalization in SGD. The switching from adaptive learning rates to SGD has been also investigated by different strategies in Luo et al. (2019) and Keskar & Socher (2017).
Lastly, a scaled Gaussian noise, √
2λβ−1ξn+1, is added as a consequence of the discretization of the Langevin SDE. The term is essential to prove the convergence property of THεO POULA. Adding properly scaled Gaussian noise allows the new algorithm to escape local minima in a similar manner to the standard SGLD method, see Raginsky et al. (2017).

heading: 3.3 CONVERGENCE ANALYSIS
text: We present in this section the main convergence results of THεO POULA to πβ in Wasserstein-1 and Wasserstein-2 distances as defined in (4). The convergence is guaranteed when the learning rate is less than λmax, which is given by
λmax = min
{ 1
4η2 ,
1
214η2(8lC4l)2
} . (9)
where nCk is the binomial coefficient ‘n choose k’ and l = 2r + 1. Note that the learning rate restriction causes no issues a η is typically very small (η 1). Moreover, let T := 1/λ.
Theorem 3.1 and Corollary 3.1 state the non-asymptotic (upper) bounds between L ( θλn )
and πβ . An overview of the proofs of our main results can be found in Appendix C. Theorem 3.1. Let Assumptions 3.1 and 3.2 hold. Then, for every 0 < λ ≤ λmax and n ∈ N, we have
W1 ( L ( θλn ) , πβ ) ≤M1 √ λ+M2e −ċλn,
where ċ, M1 and M2 are constants independent of n and λ. The explicit form of ċ, M1 and M2 are given in Table 7. Corollary 3.1. Let Assumptions 3.1 and 3.2 hold. Then, for every 0 < λ ≤ λmax and n ∈ N, we have
W2 ( L ( θλn ) , πβ ) ≤M3 √ λ+ z2λ 1 4 +M4e −ċλn,
where ċ, z2, M3 and M4 are constants independent of n and λ. The explicit form of ċ, z2, M3 and M4 are given in Table 7.
We are now concerned with the expected excess risk of THεO POULA generated by (7), so called the optimization error of θλn, defined as
E[u(θλn)]− u(θ∗) (10) where θ∗ := arg minθ∈Rd u(θ). To derive the bound of the expected excess risk, it is again decomposed into two parts; E[u(θλn)] − E[u(θ∞)] and E[u(θ∞)] − u(θ∗). Here, θ∞ follows the invariant distribution πβ . The following theorem describes the bound of the expected excess risk of THεO POULA. Theorem 3.2. Let Assumptions 3.1 and 3.2 hold. For any n ∈ N, the expected excess risk of the n-th iterate of THεO POULA (7) is upper bounded by
E[u(θλn)]− u(θ∗) ≤M5W2(L(θλn, πβ)) + 1
β
[ d
2 log
( Ke
A
( B
d β + 1
)) + log 2 ] where W2(L(θλn), πβ) is given in Corollary 3.1 and A, B are given in Remark B.2. and M5, K are given in Table 7. All the constatns are independent of n and λ.
3.4 AVERAGED THεO POULA
One notes that Theorem 3.1 implies that THεO POULA converges, under suitable decreasing step size regime, to the invariant measure πβ and thus its performance can be further improved by averaging. It is achieved by averaging of trajectories of the parameters after a user-specified trigger Q,
1 n−Q+1 ∑n i=Q θ λ i , instead of the last updated parameter θ λ n (Polyak & Juditsky, 1992). In particular, we use a trigger strategy which starts the averaging when no improvement in the validation metric is seen for a patience number of epochs. For our experiments, we set the patience number to 5.
Our experiments show that averaged THεO POULA performs better than the other stochastic optimization methods for language modeling tasks. Moreover, while a learning rate decay, which requires additional tuning effort, has to be applied for the other optimizers to obtain their best performance, averaged THεO POULA uses a constant learning rate, which is another practical benefit of our newly proposed algorithm.

heading: 4 EMPIRICAL PERFORMANCE ON REAL DATA SETS
text: This section examines the performance of THεO POULA on real data sets by comparing it with those of other stochastic optimization algorithms including Adam (Kingma & Ba (2015)), AdaBelief (Zhuang et al. (2020)), AdamP (Heo et al. (2021)), AdaBound (Luo et al. (2019)), AMSGrad (Reddi et al. (2018)), RMSProp (Tieleman & Hinton (2012)), SWATS (Keskar & Socher (2017)), SGD (with momentum) and ASGD (Merity et al. (2018)). We conduct the following deep learning experiments: image classficiation on CIFAR-10 (Krizhevsky et al.) and CIFAR-100 (Krizhevsk (2009)) and language modeling on Penn Treebank (Marcus et al. (1999)). Each experiment is run three times to compute the mean and standard deviation of the best accuracy on the test dataset. We provide details of the experiments including learning curves and hyperparameter settings in Appendix E.
For our experiments, we consider η = 0 in (5). This is justified by the fact that some form of dissipativity may exist for specific problems such as the one considered here, although this has not been verified theoretical so far. In Appendix F.2, we perform additional experiments with η 6= 0, which show a very similar performance by THεO POULA as in Table 2 without any noticeable loss of accuracy. This demonstrates that there is no gap between theory and practice of our work.
Image classification We replicate the experiments of VGG11 (Simonya & Zisserman (2015)), ResNet34 (Ioffe & Szegedy (2016)) and DenseNet121 (Huang et al. (2017)) on CIFAR-10 and CIFAR-100 in the official implementation of Zhuang et al. (2020). They provide a reliable baseline of the experiments by comparing the performance of various stochastic optimizers with extensive hyperparameter search. We search the optimal hyperparameters for THεO POULA among λ = {1, 0.5, 0.1, 0.05, 0.01} and ε = {1, 0.1, 0.01}. β is chosen among {108, 1010, 1012} across all the experiments.
Table 2 shows the test accuracy for VGG11, ResNet34 and DenseNet121 on CIFAR-10 and CIFAR100. As reported in Table 2, our algorithm achieves the highest accuracy and significantly outperforms the other optimizers across all the experiments. In particular, even THεO POULA with the second best hyperparameter is comparable to AdaBelief and outperforms the other methods, validating that the solutions found by THεO POULA yield good generalization performance. Also, the improvement of our algorithm is increasingly prominent as the models and datasets are more complicated and large-scale.
Language modeling We perform language modeling over the Penn Treebank (PTB) with AWDLSTMs of Merity et al. (2018). It is reported that Non-monotonically Triggered ASGD (NT-ASGD) achieves state-of-the-art performance for the language modeling task with AWD-LSTMs. Motivated by this observation, we consider averaged THεO POULA for the experiment. Due to a limited computation budget, we only test ASGD and AdaBelief rather than investigating all the optimizers in this experiment 3.
For a fair comparison, the averaging scheme has also been applied to AdaBelief, but we have found that it does not improve the performance of AdaBelief. Instead, AdaBelief uses a developmentbased learning rate decay, which decreases the learning rate by a constant factor if the model does not attain a new best value for multiple epochs. For ASGD and THεO POULA, a constant learning
3Since AdaBelief significantly outperforms the other optimizers including vanilla SGD, AdaBound, Yogi (Zaheer et al. (2018)), Adam, MSVAG (Balles & Hennig (2018)), RAdam, Fromage and AdamW (Loshchilov & Hutter (2019)) in the same experiment, we believe that we do not need to explore all the optimizers.
rate is used without a learning rate decay. In order to compare with the baseline, we apply gradient clipping of 0.25 to all optimizers. See Appendix E for more information.
Table 3 shows that THεO POULA attains the lower test perplexity against the baselines for AWD-LSTM with one, two and three layers. AdaBelief shows a comparable performance with ASGD for 2-layer and 3-layer models.
Our experimental results show that THεO POULA achieves higher accuracy than AdaBelief (known as the state-of-the-art algorithm for many deep learning tasks) on image classification and language modeling tasks for various deep learning models. Furthermore, it is easier to tune parameters of THεO POULA since the number of hyperparmeters for THεO POULA is less than that of Adam-type optimizers.

heading: 5 CONCLUSION AND DISCUSSION
text: This paper begins with an example which illustrates that local Lipschitz continuous gradients can cause serious convergence issues for popular adaptive optimization methods. Such issues manifest themselves as vanishing/exploding gradient phenomena. It proceeds by proposing a novel optimization framework, which is suitable for the fine tuning of neural network models by combining elements of the theory of Langevin SDEs, tamed algorithms and carefully designed boosting functions that handle sparse and super-linearly growing gradients. Further, a detailed convergence analysis of the newly proposed algorithm THεO POULA is provided along with full theoretical guarantees for obtaining the best known convergence rates. Our experiments confirm that THεO POULA outperforms other popular stochastic optimization methods.
We believe that there is much room for improvement of our novel framework. For example, the improved performance can be further achieved by identifying more efficient taming and boosting functions, which demonstrates the potential of our framework.

heading: A DETAILS OF THE EXPERIMENT IN SECTION 2
text: This section provides the necessary theoretical details of the experiment in Section 2. We continue to consider the optimization problem (2). One calculates that
u(θ) =
{ θ30 + 74θ
2, |θ| ≤ 1, θ30 + 74 (2|θ| − 1), |θ| > 1
and
u′(θ) = { 30θ29 + 72θ, |θ| ≤ 1, 30θ29 + 72sgn(θ), |θ| > 1.
Note that u(θ) and u′(θ) are continuous since u(1) = limθ↓1 u(θ) = 114 , u(−1) = limθ↑−1 = 11 4 , u′(1) = limθ↓−1 u ′(1) = 672 and u
′(−1) = limθ↑1 u′(−1) = − 672 . Therefore, the minimum value is attained at θ = 0.
To show that G is locally Lipschitz continuous, we check that for |θ|, |θ′| > 1 and x ∈ Rd, |G(θ, x)−G(θ′, x)| ≤ (2 + 21x≤1)|sgn(θ)− sgn(θ′)|+ 30|θ29 − θ′29| ≤ 34(1 + |θ|+ |θ′|)28|θ − θ′|. For |θ|, |θ| ≤ 1, we have
|G(θ, x)−G(θ′, x)| ≤ (2 + 21x≤1)|θ − θ′|+ 30|θ29 − θ′29| ≤ 34(1 + |θ|+ |θ′|)28|θ − θ′|.
For |θ| ≤ 1, |θ| > 1, we obtain |G(θ, x)−G(θ′, x)| ≤ (2 + 21x≤1)|θ − sgn(θ′)|+ 30|θ29 − θ′29|
≤ (2 + 21x≤1)|θ − θ′|+ 30|θ29 − θ′29| ≤ 34(1 + |θ|+ |θ′|)28|θ − θ′|
where the second inequality follows from the following relations
θ − θ′ ≤ θ − 1 ≤ 0, for θ′ > 1, 0 ≤ θ + 1 ≤ θ − θ′, for θ′ < −1.

heading: B KEY OBSERVATIONS FROM ASSUMPTION 3.1 AND 3.2
text: This section introduces some useful general results, that can be obtained from Assumption 3.1 and 3.2. Note that some of the below observations can be also found in Zhang et al. (2019) and Lovas et al. (2020). However, to make our paper self-contained, we record all the results which are necessary for the convergence analysis. Remark B.1. From Assumption 3.1, one observes that for all θ ∈ Rd and x ∈ Rm
|G(θ, x)| ≤ K(x)(1 + |θ|q), where K(x) = 2q(L1(1 + |x|)ρ + |G(0, x)|). Remark B.2. From Assumptions 3.1 and 3.2, one obtains that
〈θ, h(θ)〉 = 〈θ, EG(θ,X0)〉+ 〈θ, ηθ|θ|2r〉 ≥ η|θ|2r+2 − E[K(X0)]|θ|(1 + |θ|q). Furthermore, for A = E[K(X0)] and B = (3E[K(X0)])q+2η−q−1, it holds that
〈θ, h(θ)〉 ≥ A|θ|2 −B. (B.1) Proposition B.1. (Lovas et al. (2020)) Let Assumptions 3.1 and 3.2 hold. Then, for every θ, θ′ ∈ Rd, 〈θ − θ′, h(θ)− h(θ′)〉 ≥ −a|θ − θ′|2, where a = L1E[(1 + |X0|)ρ](1 + 2|R|)q−1 and R is given by
R = max
{( 23(q−1)+1L1E[(1 + |X0|)ρ]
η
) 1 2r−1
,
( 2qL1E[(1 + |X0|)ρ]
η
) 1 2r } .
Proposition B.2. (Lovas et al. (2020)) Let Assumptions 3.1 and 3.2 hold. Then, one obtains that
|H(θ, x)−H(θ′, x)| ≤ L(1 + |x|)ρ(1 + |θ|+ |θ′|)l|θ − θ′|, for all θ ∈ Rd and x ∈ Rm, where L = L1 + 8rη and l = 2r + 1. Remark B.3. From Assumption 3.1 and the definition of H and Hλ,c, one obtains that for θ ∈ Rd, x ∈ Rm and i = 1, 2, · · · , d,
|H(i)(θ, x)−H(i)λ,c(θ, x)| ≤ ∣∣∣∣G(i)(θ, x)− G(i)(θ, x)1 +√λ|G(i)(θ, x)| ( 1 +
√ λ
ε+ |G(i)(θ, x)| )∣∣∣∣ + ∣∣∣∣ηθ(i)|θ|2r − η θ(i)|θ|2r1 +√λ|θ|2r ∣∣∣∣ ≤ |G(i)(θ, x)| √ λ|G(i)(θ, x)|
1 + √ λ|G(i)(θ, x)| +
√ λ|G(i)(θ, x)|
(1 + √ λ|G(i)(θ, x)|)(ε+ |G(i)(θ, x)|)
+ η|θ(i)||θ|2r ∣∣∣∣ √ λ|θ|2r
1 + √ λ|θ|2r ∣∣∣∣ ≤ √ λ|G(i)(θ, x)|2 + √ λ+ √ λη|θ(i)||θ|4r
which implies that
|H(θ, x)−Hλ,c(θ, x)|2 = d∑ i=1 [√ λ|G(i)(θ, x)|2 + √ λ+ √ λη|θ(i)||θ|4r ]2
≤ 3λ d∑ i=1 [ |G(i)(θ, x)|4 + 1 + η2|θ(i)|2|θ|8r ]
≤ 3λ [( d∑
i=1
|G(i)(θ, x)|2 )2 + d+ η2|θ|8r+2 ]
≤ 3λ [ |G(θ, x)|4 + d+ η2|θ|8r+2 ] ≤ 3λ [ 8|K(x)|4(1 + |θ|4q) + d+ η2|θ|8r+2 ] .
Remark B.4. From Assumption 3.1 and the definition of H and Hλ,c, one calculates that for all θ ∈ Rd and x ∈ Rm,
|H(θ, x)|2 = |G(θ, x) + ηθ|θ|2r|2 ≤ 2|G(θ, x)|2 + 2η2|θ|4r+2
≤ 4|K(x)|2(1 + |θ|2q) + 2η2|θ|4r+2
and |Hλ,c(θ, x)|2 = d∑ i=1 [ G(i)(θ, x) 1 + √ λ|G(i)(θ, x)| ( 1 +
√ λ
ε+ |G(i)(θ, x)|
) + η θ(i)|θ|2r
1 + √ λ|θ|2r
]2
≤ d∑ i=1 [ |G(i)(θ, x)| 1 + √ λ|G(i)(θ, x)| +
√ λ|G(i)(θ, x)|
(1 + √ λ|G(i)(θ, x)|)(ε+ |G(i)(θ, x)|)
+ η |θ(i)||θ|2r
1 + √ λ|θ|2r
]2
≤ d∑ i=1 ( |G(i)(θ, x)|+ √ λ+ η|θ(i)||θ|2r )2
≤ 3 d∑ i=1 ( |G(i)(θ, x)|2 + λ+ η2|θ(i)|2|θ|4r ) ≤ 3|G(θ, x)|2 + 3λd+ 3η2|θ|4r+2
≤ 6|K(x)|2(1 + |θ|2q) + 3λd+ 3η2|θ|4r+2.

heading: C OVERVIEW OF THE PROOFS
text: This section provides an overview of the proofs of our main results. We begin by introducing suitable Lyapunov functions and auxiliary processes to analyze the convergence of our newly introduced algorithm. For each m ≥ 1, define the Lyapunov function Vm by
Vm(θ) := (1 + |θ|2) m 2 , θ ∈ Rd (C.1)
and similarly vm(x) = (1+x2) m 2 for any real x ≥ 0. Both functions are continuously differentiable and lim|θ|→∞∇Vm(θ)/Vm(θ) = 0. Also, define Zλt = Zλt, which is the time-changed Lagevin dynamics governed by
dZλt = −λh(Zλt )dt+ √ 2β−1λdB̃λt (C.2)
where B̃λt = Bλt/ √ λ is a Brownian motion.
We next define the continuous-time interpolation of the new algorithm, see (7), as dθ̄λt = −λHλ ( θ̄λbtc, Xdte ) dt+ √ 2λβ−1dB̃λt (C.3)
with initial condition θ̄λ0 = θ λ 0 . Henceforth, bxc denotes the integer part of a positive real x and dxe = bxc+ 1. Remark C.1. Due to the homogeneous nature of the coefficients of the continuous-time interpolation of THεO POULA (C.3) and when one selects a version of the driving Brownian motion such that it coincides with ξn at grid points, it follows that the interpolated process (C.3) equals the process of THεO POULA (7) almost surely at grid points, i.e. θ̄λn = θ λ n (a.s), ∀n ∈ N.
Furthermore consider the continuous-time process ζs,v,λt , t ≥ s which is the solution to the SDE dζs,v,λt = −λh ( ζs,v,λt ) dt+ √ 2λβ−1dB̃λt (C.4)
with initial condition ζs,v,λs := v, v ∈ Rd. Let us also define T := 1λ , which allows us to create suitable subintervals on the positive real line in order to compare the behaviour of the aforementioned processes at each such interval.
Definition C.1. Fix k ∈ N and define ζ̄λ,kt := ζ kT,θ̄λkT ,λ t where ζ kT,θ̄λkT ,λ t is defined in (C.4).
To derive non-asymptotic (upper) bounds for W1 ( L ( θλt ) , πβ ) and W2 ( L ( θλt ) , πβ ) , the following decomposition is used in terms of the auxiliary processes θ̄λt , ζ̄ λ,n t and Z λ t as follows:
Wj ( L ( θλt ) , πβ ) ≤Wj ( L ( θ̄λt ) ,L ( ζ̄λ,nt )) +Wj ( L ( ζ̄λ,nt ) ,L ( Zλt )) +Wj ( L ( Zλt ) , πβ )
for j = 1, 2.
C.1 PRIMARY ESTIMATES We collect first the necessary estimates in order to obtain (upper) bounds for W1 ( L ( θλt ) , πβ )
and W2 ( L ( θλt ) , πβ ) . All proofs of the lemmas in this section can be found in Appendix D. The following two lemmas provide, uniform in n, moment estimates of the process (θλn)n≥1. Lemma C.1. Let Assumptions 3.1 and 3.2 hold. Then, there exist M0 > 0 and λmax, which is defined in (9), such that for any λ ∈ (0, λmax) and any n ∈ N,
E|θλn+1|2 ≤ (
1− η 2
√ λ )n E|θ0|2 + [ 5M20 + 4 √ λmaxd
η
( β−1 + 2 + 2λ2max ) + 4(1 + λmax) √ dM0
η + 4ηM20
√ λmax ] and, moreover,
sup n
E|θλn+1|2 ≤E|θ0|2 + [ 5M20 + 4 √ λmaxd
η
( β−1 + 2 + 2λ2max ) + 4(1 + λmax) √ dM0
η + 4ηM20
√ λmax ] .
Lemma C.2. Let Assumptions 3.1 and 3.2 hold. Then, there exist M0 > 0 and λmax, which is defined in (9), such that for any λ ∈ (0, λmax), n ∈ N, and p ∈ [1, 8(2r + 1)],
E|θλn+1|2p ≤ (1− η2λ)nE|θλ0 |2p + Ap η2
and sup n E|θλn+1|2p ≤ E|θλ0 |2p + Ap η2 where Ap is given in Table 7. Lemma C.3. Let Assumptions 3.1 and 3.2 hold. Then, there exist M0 > 0 and λmax, which is defined in (9), such that for any λ ∈ (0, λmax) and n ∈ N,
E[V4(θ̄λnT )] ≤ 2E|θ0|4 + 2 + 2 A2 η2
where A2, i.e. Ap for p = 2, is given in Table 7.
Proof. From the definition of the Lyapunov function and Remark C.1, we have E[V4(θ̄λnT )] = E[(1 + |θ̄λnT |2)2]
≤ 2 + 2E|θ̄λnT |4
≤ 2 + 2E|θ0|4 + 2 A2 η2 .
Moreover, the necassary moment bounds hold also for the auxiliary process {ζ̄λ,nt }t≥nT . Lemma C.4. (Lemma 3.5. of Chau et al. (2019)) Let Assumptions 3.1 and 3.2 hold. Then,
E[V2(ζ̄λ,nt )] ≤ E[V2(θ0)] + c̃(2)
c̄(2) + 2(CXη
−1 + 2M20 (2 + η) + 2d(ηβ) −1 √ λmax) + 1,
E[V4(ζ̄λ,nt )] ≤ 2E|θ0|4 + 2 + 2 A2 η2 + c̃(4) c̄(4) ,
where c̄(p), c̃(p) are given in Table 7. Let PV2 denote the subset of P(Rd) such that every µ ∈ PV2 satisfies ∫ Rd V2(θ)µ(dθ) < ∞. Moreover, let the following functional be considered
w1,2(µ, ν) := inf ζ∈C(µ,ν) ∫ Rd ∫ Rd [1 ∧ |θ − θ|′] [(1 + V2(θ) + V2 (θ′)) ζ (dθdθ′) (C.5)
where C(µ, ν) is defined in (4). The following lemma states the contraction property of the Langevin SDE (C.2) in w1,2, which yields the desired result for W1(L(Zλn), πβ). Lemma C.5. (Proposition 3.14 of Chau et al. (2019)) LetZ ′t, t ∈ R+ be the solution of the Langevin SDE (6) with initial condition Z ′0 = θ ′ 0 which is independent of G∞ and |θ′0| ∈ L2. Then,
w1,2 ( L ( Zλt ) ,L (Z ′t) ) ≤ ĉe−ċtw1,2 (L (θ0) ,L (θ′0))
where w1,2 is defined in (C.5).
The following two Lemmas combined establish the required W1(L(θ̄λt ),L(Zλt )) estimate. One recalls first that for any t > 0, there exists a unique integer m such that t ∈ [mT, (m+ 1)T ). Lemma C.6. Let Assumptions 3.1 and 3.2 hold. Then, for 0 < λ < λmax, one obtains
W2 ( L ( θ̄λt ) ,L ( ζ̄λ,mt )) ≤ √ λ √ e3a(C1 + C2 + C3)
where C1, C2, C3 are given explicitly in Table 7. Lemma C.7. Let Assumptions 3.1 and 3.2 hold. Then, for 0 < λ < λmax, one obtains
W1 ( L ( ζ̄λ,mt ) ,L ( Zλt )) ≤ √ λz1
where z1 is given explicitly in Table 7. Lemma C.8. Let Assumptions 3.1 and 3.2 hold. Then, for 0 < λ < λmax, one obtains
W2 ( L ( ζ̄λ,mt ) ,L ( Zλt )) ≤ λ 14 z2
where z2 is given explicitly in Table 7.
C.2 PROOFS OF MAIN RESULTS
It is assumed throughout the paper that the random variable θ0, G∞ := σ (∪n∈NGn) and (ξn)n∈N are independent. Proof of Theorem 3.1. Observe that W1 ( L ( θλn ) ,L ( Zλt )) is decomposed as follows:
W1 ( L ( θλn ) , πβ ) ≤W1 ( L ( θ̄λn ) ,L ( Zλn )) +W1 ( L ( Zλn ) , πβ ) .
Note that there exists a unique integer m such that n ∈ [mT, (m + 1)T ). Thus, from the results of Lemma C.6 and C.7, the first term in the right-hand side is estimated
W1 ( L ( θ̄λn ) ,L ( Zλn )) ≤ W1 ( L ( θ̄λn ) ,L ( ζ̄λ,mn )) +W1 ( L ( ζ̄λ,mn ) ,L ( Zλn ))
≤ W2 ( L ( θ̄λn ) ,L ( ζ̄λ,mn )) +W1 ( L ( ζ̄λ,mn ) ,L ( Zλn ))
≤ √ λ( √ e3a(C1 + C2 + C3) + z1).
Consequently, we derive W1 ( L ( θλn ) , πβ ) ≤ √ λ( √ e3a(C1 + C2 + C3) + z1) + w1,2 ( L ( Zλn ) , πβ )
≤ √ λ( √ e3a(C1 + C2 + C3) + z1) + ĉe −ċλnw1,2(θ0, πβ)
≤ √ λ( √ e3a(C1 + C2 + C3) + z1) + ĉe −mċ [ 1 + E[V2(θ0)] + ∫ Rd V2(θ)πβ(dθ) ] where Remark C.5 is used for the first inequality.
Proof of Corollary 3.1. Let n ∈ [mT, (m+1)T ). Then, Lemma C.6 and C.8 and Remark C.5 yield that
W2 ( L ( θλn ) , πβ ) ≤ W2 ( L ( θ̄λn ) ,L ( Zλn )) +W2 ( L ( Zλn ) , πβ )
≤ W2 ( L ( θ̄λn ) ,L ( ζ̄λ,mn )) +W2 ( L ( ζ̄λ,mn ) ,L ( Zλn )) +W2 ( L ( Zλn ) , πβ )
≤ √ e3a(C1 + C2 + C3) √ λ+ z2λ 1 4 + √ 2w1,2(L(Zλt ), πβ)
≤ √ e3a(C1 + C2 + C3) √ λ+ z2λ 1 4 + √ ĉe−ċλn/2 √ 2w1,2(θ0, πβ) ≤ √ e3a(C1 + C2 + C3) √ λ+ z2λ 1 4 + √ 2ĉe−ċm/2 ( 1 + E [V2 (θ0)] + ∫ Rd V2(θ)πβ(dθ) ) .
Proof of Theorem 3.2. We begin by decomposing expected excess risk (10) as follows:
E[u(θλn)]− u(θ∗) ≤ E[u(θλn)]− E[u(θ∞)] + E[u(θ∞)]− u(θ∗).
Let us focus on estimating the first part, E[u(θλn)]− E[u(θ∞)]. Observe that for θ ∈ Rd
|∇u(θ)| = |h(θ)| ≤ r1|θ|2r+1 + 2E[K(X0)] by separating the cases |θ| ≤ 1 and |θ| > 1 where r1 = E[K(X0)] + η due to Remark B.1. Then, we have
u(w)− u(v) = ∫ 1
0
〈∇u((1− t)v + tw), w − v〉dt
≤ ∫ 1
0
|∇u((1− t)v + tw)||w − v|dt
≤ ∫ 1
0
( a1(1− t)l|v|l + a1tl|w|l + 2E[K(X0)] ) |w − v|dt
= ( a1 l + 1 |v|l + a1 l + 1 |w|l + 2E[K(X0)] ) |w − v| (C.6)
where l = 2r+ 1 and a1 = 2lr1 = 2l(E|K(X0)|+ η). Let P denote the coupling between µ and ν that achieves W2(µ, ν) with µ = L(θλn) and ν = L(θ∞). Then, from (C.6), we obtain
Eu(θλn)− Eu(θ∞) = EP [u(θλn)− u(θ∞)] ≤ EP [(
a1 l + 1 |θλn|l + a1 l + 1
|θ∞|l + 2EP [K(X0)] ) |θλn − θ∞| ]
≤ √ EP [(
a1 l + 1 |θλn|l + a1 l + 1
|θ∞|l + 2E[K(X0)] )2]√ EP |θλn − θ∞|2
≤ (
a1 l + 1
√ E|θλn|2l +
a1 l + 1
√ E|θ∞|2l + 2E[K(X0)] ) W2 ( L(θλn), πβ ) ≤ ( a1 l + 1 √ E|θ0|2l + Al η2 + a1 l + 1 √ E|θ∞|2l + 2E[K(X0)]
) ×W2 ( L(θλn), πβ ) (C.7)
where we have used Lemma C.2 for the last inequality.
We take a similar approach in Raginsky et al. (2017) to estimate the second term. From Equation (3.18), (3.20) in Raginsky et al. (2017), we obtain
Eu(θ∞)− u(θ∗) ≤ 1
β ( − ∫ Rd e−βu(θ) Λ log e−βu(θ) Λ dθ − log Λ ) − u∗
≤ d 2β log
( 2πe(B + d/β)
Ad
) − log Λ
β − u∗ (C.8)
where Λ = ∫ Rd e −βu(θ)dθ is the normalizing constant.
Using (B.1), we obtain 〈θ∗, h(θ∗)〉 ≥ A|θ∗|2 −B
which yields
|θ∗|2 ≤ √ B
A .
Moreover, for w ∈ Rd, we have u(θ∗)− u(w) = ∫ 1
0
〈∇u(w + t(θ∗ − w)), θ∗ − w〉dt
= ∫ 1 0 〈∇u(w + t(θ∗ − w))−∇u(θ∗), θ∗ − w〉dt
= ∫ 1 0 1 t− 1 〈∇u(w + t(θ∗ − w))−∇u(θ∗), w − θ∗ + t(θ∗ − w)〉dt.
From Proposition B.2, we further obtain
−β(u(θ∗)− u(w)) = β|u(θ∗)− u(w)| ≤ β ∫ 1
0
1
t− 1 |〈h(w + t(θ∗ − w))− h(θ∗), w − θ∗ + t(θ∗ − w)〉|dt
≤ βLE(1 + |X0|)ρ ∫ 1
0
(1 + |w + t(θ∗ − w)|+ |θ∗|)l(1− t)|w − θ∗|2dt
≤ βLE(1 + |X0|)ρ ∫ 1
0
(1 + |w|+ |θ∗ − w|+ |θ∗|)l(1− t)|w − θ∗|2dt
= βLE(1 + |X0|)ρ(1 + 2|θ∗|+ 2|θ∗ − w|)l |w − θ∗|2
2
where we have used the elementary inequality 0 ≤ |w| − |θ∗| ≤ |θ∗ − w| for the last inequality.
Define R0 := max{ √ B/A, √ 2d/(βLE(1 + |X0|)ρ)} and Br(p) = {x ∈ Rd||x− p| > r}. Then, from the above inequality, one further calculates
log Λ β = −u(θ∗) + 1 β log ∫ Rd eβ(u(θ ∗)−u(w))dw
≥ −u(θ∗) + 1 β log ∫ Rd e−βLE(1+|X0|) ρ(1+2|θ∗|+2|θ∗−w|)l |w−θ ∗|2 2 dw
≥ −u(θ∗) + 1 β log ∫ BR0 (θ ∗) e−βLE(1+|X0|) ρ(1+4R0) l |w−θ∗|2 2 dw
= −u(θ∗) + 1 β log
[( 2π
βK )d/2 ∫ BR0 (θ ∗) fX(w)dw ]
≥ −u(θ∗) + 1 β log
( 1
2
( 2π
Kβ
)d/2) (C.9)
where K = LE(1 + |X0|)ρ(1 + 4R0)l and fX is the density function of a multivariate normal variable X with mean θ∗ and covariance 1Kβ Id. Here, the last inequality is obtained from the following inequality:
∫ BR0 (θ ∗) fX(w)dw = P (|X − θ∗| > R0)
= P ( |X − θ∗| > √ KβR20 d √ d Kβ ) ≤ d
KβR20
≤ 1 2(1 + 4R0)l
≤ 1 2 .
Combining (C.8) and (C.9), we derive
Eu(θ∞)− u(θ∗) ≤ d
2β log
( 2πe(B + d/β)
Ad
) − 1 β log ( 1 2 ( 2π Kβ )d/2) ≤ 1
β
[ d
2 log
( Ke
A
( B
d β + 1
)) + log 2 ] . (C.10)
Consequently, from (C.7) and (C.10), we derive
Eu(θλn)− u(θ∗) ≤ M5W2(L(θλn, πβ))
+ 1
β
[ d
2 log
( Ke
A
( B
d β + 1
)) + log 2 ]
where M5 = a1l+1 √ E|θ0|2l + Alη2 + a1 l+1 √ E|θ∞|2l + 2E[K(X0)].

heading: D PROOFS OF LEMMAS IN APPENDIX C
text: Proof of Lemma C.1. Define Ĝ(i)λ,c(θ, x) = G(i)(θ,x) 1+ √ λ|G(i)(θ,x)|
( 1 +
√ λ
ε+|G(i)(θ,x)|
) , which is part of the
adaptive gradient of H(i)λ,c(θ, x). One observes that i ∈ {1, · · · , d}
|Ĝ(i)λ,c(θ, x)| = |G(i)(θ, x)|
1 + √ λ|G(i)(θ, x)|
+ √ λ |G(i)(θ, x)| (ε+ |G(i)(θ, x)|)(1 + √ λ|G(i)(θ, x)|)
≤ 1√ λ
+ √ λ |G(i)(θ, x)|/ε
1 + |G(i)(θ, x)|/ε
≤ 1√ λ
+ √ λ (D.1)
to obtain
〈θ, Hλ,c(θ, x)〉 = d∑ i=1 θ(i) · Ĝ(i)λ,c(θ, x) + η |θ|2r+2 1 + √ λ|θ|2r
≥ d∑ i=1 |θ(i)| ( − 1√ λ − √ λ ) + η |θ|2r+2 1 + √ λ|θ|2r
≥ − (
1√ λ
+ √ λ )√ d|θ|+ η |θ| 2r+2
1 + √ λ|θ|2r
for all x ∈ Rm and θ ∈ Rd. Then, 2λE [ 〈 θ λ n
|θλn|2 , Hλ,c(θ
λ n, Xn+1)〉 ∣∣∣∣θλn] ≥ −2(√λ+ λ 32) √ d
|θλn| + 2ηλ
|θλn|2r
1 + √ λ|θλn|2r . (D.2)
On the other hand, due to (D.1),
|Hλ,c(θ, x)|2 = 〈Hλ,c(θ, x), Hλ,c(θ, x)〉 = d∑ i=1 ( Ĝ (i) λ,c(θ, x) + η θ(i)|θ|2r 1 + √ λ|θ|2r )2
≤ d∑ i=1 ( 2|Ĝ(i)λ,c(θ, x)| 2 + 2η2 |θ(i)|2|θ|4r (1 + √ λ|θ|2r)2 )
≤ 2d (
1√ λ
+ √ λ )2 + 2η2
|θ|4r+2
(1 + √ λ|θ|2r)2
≤ 4d ( 1
λ + λ
) + 2η2|θ|2 |θ| 4r
(1 + √ λ|θ|2r)2 . (D.3)
which yields that 2λE [ − λ
2|θλn|2 |Hλ,c(θλn, Xn+1)|2 ∣∣∣∣θλn] ≥ −2λ(2d (1 + λ2)|θλn|2 + η2 λ|θ| 4r (1 + √ λ|θ|2r)2 ) ≥ −4λd (1 + λ 2)
|θλn|2 − 2λη2. (D.4)
From (D.2) and (D.4), one calculates that 2λE [ 〈 θ λ n
|θλn|2 , Hλ,c(θ
λ n, Xn+1)〉 −
λ
2|θλn|2 |Hλ,c(θλn, Xn+1)|2 ∣∣∣∣θλn] ≥ −2 (√ λ+ λ 3 2 )√ d
|θλn| + 2ηλ
|θλn|2r
1 + √ λ|θλn|2r
− 4λd (1 + λ 2)
|θλn|2 − 2λη2 =: f(θλn).
Since f(θ) tends to 2η √ λ− 2λη2 as |θ| → ∞, there exists M0 > 0 such that
f(θλn) ≥ η √ λ− λη2 = η √ λ(1− √ λη)
for all |θλn| ≥M0 and λ < 1η2 . Moreover, for λ ≤ 1 4η2 , it can rewritten as there exists M0 > 0 such that
f(θλn) ≥ η √ λ− λη2 = η
√ λ 2 (D.5)
for all |θλn| ≥M0. Therefore,
E [( 2λ〈θλn, Hλ,c(θλn, Xn+1)〉 − λ2|Hλ,c(θn, Xn+1)|2 ) 1|θλn|≥M0 ∣∣∣∣θλn] ≥ η √ λ 2 |θλn|2,
implying that E [ |θλn+1|21|θλn|≥M0 ∣∣∣∣θλn] = E [( |θλn|2 − 2λ〈θλn, Hλ,c(θn, Xn+1)〉+ λ2|Hλ,c(θn, Xn+1)|2 + 2λ
β |ξn+1|2 ) 1|θλn|≥M0 ∣∣∣∣θλn] ≤ ( 1− η √ λ
2
) |θλn|2 + 2λd
β . (D.6)
Let us consider the case of |θλn| < M0. From the fact that
〈θ, Hλ,c(θ, x)〉 = d∑ i=1 θ(i) · Ĝ(i)λ,c(θ, x) + η |θ|2r+2 1 + √ λ|θ|2r
≤ d∑ i=1 |θ(i)| ( 1√ λ + √ λ ) + η |θ|2r+2 1 + √ λ|θ|2r
≤ (
1√ λ
+ √ λ )√ d|θ|+ η |θ| 2r+2
1 + √ λ|θ|2r
(D.7)
and (D.3), it can be shown that E [ |θλn+1|21|θλn|<M0 ∣∣∣∣θλn] = E [( |θλn|2 − 2λ〈θλn, Hλ,c(θλn, Xn+1)〉+ λ2|Hλ,c(θλn, Xn+1)|2 + 2λ
β |ξn+1|2 ) 1|θλn|<M0 ∣∣∣∣θλn] ≤ ( |θλn|2 + 2λd
β
) 1|θλn|<M0 + E [( 2λ|〈θλn, Hλ,c(θλn, Xn+1)〉|+ λ2|Hλ,c(θλn, Xn+1)|2 ) 1|θλn|<M0 ∣∣∣∣θλn] ≤ |θλn|2 + 2λd
β + 2
(√ λ+ λ 3 2 )√ dM0 + 2η √ λM20 + 4d(λ+ λ 3) + 2η2M20λ
≤ ( 1− η √ λ
2
) |θλn|2 + √ λ ( 5η
2 M20 +
2 √ λd
β + 2(1 + λ)
√ dM0 + 4d( √ λ+ λ 5 2 ) + 2η2M20 √ λ ) .(D.8)
Consequently, (D.6) and (D.8) yield that E [ |θλn+1|2 ∣∣∣∣θλn] ≤ (1− η √ λ
2
) |θλn|2 + √ λ ( 5η
2 M20 +
2 √ λd
β + 2(1 + λ)
√ dM0 + 4d( √ λ+ λ 5 2 ) + 2η2M20 √ λ ) .
As a result, E [ |θλn+1|2 ] ≤ ( 1− η √ λ
2
)n E|θλ0 |2
+ √ λ
( 5η
2 M20 +
2 √ λd
β + 2(1 + λ)
√ dM0 + 4d( √ λ+ λ 5 2 ) + 2η2M20 √ λ ) ∞∑ j=1 ( 1− η √ λ 2 )j ≤ ( 1− η √ λ
2
)n E|θλ0 |2 + ( 5M20 + 4 √ λd
βη + 4(1 + λ)
√ dM0η −1 + 8d( √ λ+ λ 5 2 )η−1 + 4ηM20 √ λ ) .
Proof of Lemma C.2. For any integer p > 1, |θλn+1|2p is written as
|θλn+1|2p = ( |∆n|2 + 2λ
β |ξn+1|2 + 2〈∆n,
√ 2λ
β ξn+1〉 )p where ∆n = θλn − λHλ,c(θλn, Xn+1). Then, we obtain
E[|θλn+1|2p|θλn] = E [( |∆n|2 + ∣∣∣∣ √ 2λ
β ξn+1
∣∣∣∣2 + 2〈∆n, √ 2λ
β ξn+1〉
)p∣∣∣∣θλn]
= ∑
k1+k2+k3=p
p! k1!k2!k3! E [ |∆n|2k1 ∣∣∣∣ √ 2λ β ξn+1 ∣∣∣∣2k2(2〈∆n, √ 2λ β ξn+1〉 )k3∣∣∣∣θλn]
≤ E[|∆n|2p|θλn] + 2pE [ |∆n|2p−2〈∆n, √ 2λ
β ξn+1〉
∣∣∣∣θλn]
+ 2p∑ k=2 ( 2p k ) E [ |∆n|2p−k ∣∣∣∣ √ 2λ β ξn+1 ∣∣∣∣k∣∣∣∣θλn] ≤ E[|∆n|2p|θλn] + 2(p−1)∑ l=0 ( 2p l + 2 ) E [( |∆n|2(p−1)−l ∣∣∣∣ √ 2λ β ξn+1 ∣∣∣∣q−1)2λβ |ξn+1|2 ∣∣∣∣θλn]
= E[|∆n|2p|θλn] + ( 2p 2 ) 2(p−1)∑ l=0 ( 2(p− 1) l ) E [( |∆n|2(p−1)−l ∣∣∣∣ √ 2λ β ξn+1 ∣∣∣∣l)2λβ |ξn+1|2 ∣∣∣∣θλn]
≤ E[|∆n|2p|θλn] + 22p−3p(2p− 1) ( E[|∆n|2p−2|θλn] 2λd
β +
( 2λ
β
)p E|ξn+1|2p ) . (D.9)
Define |∆n|2 = |θλn|2 + rn where rn = −2λ〈θλn, Hλ,c(θλn, Xn+1)〉+λ2|Hλ,c(θλn, Xn+1)|2 to write E [ |∆n|2p ∣∣∣∣θλn] = p∑ k=0 ( p k ) |θλn|2(p−k)E[rkn|θλn]
= |θλn|2p + p|θλn|2p−2E[rn|θλn] + p∑ k=2 ( p k ) |θλn|2(p−k)E[rkn|θλn]. (D.10)
Now, we focus on the case where |θλn| > M where
M := max { M0, 1, 2 √ λmaxd(1 + λ 2 max)
(2− √ λmaxη)η
, (1 + λmax)
√ d
η(2− η) , 22p−2p(2p− 1)d ηβ } and M0 is defined in the proof of Lemma C.1. We need the following relations to estimate the moments of rn: for all x ∈ Rd and |θ| ≥M ,
λ2|Hλ,c(θ, x)|2 ≤ 4d(λ+ λ3) + 2η2λ|θ|2 λ|θ|4r
(1 + √ λ|θ|2r)2
≤ 4dλ(1 + λ2) + 2η2λ|θ|2
≤ 4dλ(1 + λ2)|θ|+ 2η2λ|θ|2
≤ 2 √ λη
( 2 √ λmaxd(1 + λ 2 max) |θ|η + √ λmaxη ) |θ|2
≤ 2 √ λη
( 2d √ λmax(1 + λ 2 max) Mη + √ λmaxη ) |θ|2
≤ 4 √ λη|θ|2 (D.11)
where we have used the inequality (D.3), 0 ≤ η < 1 and
M > 2 √ λmaxd(1 + λ 2 max)
(2− √ λmaxη)η
⇔ ( 2d √ λmax(1 + λ 2 max) Mη + √ λmaxη ) < 2
and note that 2 √ λmaxd(1+λ 2 max)
(2− √ λmaxη)η is finite due to λmax being less than 14η2 . Moreover, from (D.7), we have the following inequality
|2λ〈θ, Hλ,c(θ, x)〉| ≤ 2( √ λ+ λ1.5) √ d|θ|+ 2η √ λ|θ|2
√ λ|θ|2r
1 + √ λ|θ|2r
≤ 2 √ λ(1 + λ) √ d|θ|+ 2η √ λ|θ|2
≤ 2 √ λη
( (1 + λmax) √ d
|θ|η + η
) |θ|2
≤ 2 √ λη
( (1 + λmax) √ d
Mη + η
) |θ|2
≤ 4 √ λη|θ|2 (D.12)
where the last inequality holds since
M > (1 + λmax)
√ d η(2− η) ⇔ ( (1 + λmax) √ d Mη + η ) ≤ 2.
Thus, rkn can be written as
E[1{|θλn|>M}|rn| k|θλn] = E [ 1{|θλn|>M} ( − 2λ〈θλn, Hλ,c(θλn, Xn+1)〉+ λ2|Hλ,c(θλn, Xn+1)|2 )k∣∣∣∣θλn] ≤ E [ 1{|θλn|>M} ( |2λ〈θλn, Hλ,c(θλn, Xn+1)〉|+ λ2|Hλ,c(θλn, Xn+1)|2
)k∣∣∣∣θλn] ≤ E [ 1{|θλn|>M}(8 √ λη|θλn|2)k
∣∣∣∣θλn] ≤ λ k2 (8η)k|θλn|2k. Moreover, (D.5) implies that
E[1{|θλn|>M}rn|θ λ n] ≤ −
η √ λ
2 |θλn|2,
equivalently,
p|θλn|2p−2E[1{|θλn|>M}rn|θ λ n] ≤ −p
η √ λ
2 |θλn|2p. (D.13)
Using (D.13), the L2p-norm of ∆n conditional on θλn > M is given by
E [ 1{|θλn|>M}|∆n| 2p ∣∣∣∣θλn] ≤ |θλn|2p + p|θλn|2p−2E[1{θλn>M}rn|θλn] + p∑ k=2 ( p k ) |θλn|2(p−k)E[1{θλn>M}|rn| k|θλn]
≤ |θλn|2p − p η √ λ
2 |θλn|2p + p∑ k=2 ( p k ) |θλn|2(p−k)λ k 2 (8η)k|θλn|2k
≤ |θλn|2p − p η √ λ
2 |θλn|2p + |θλn|2p p∑ k=2 ( p k ) λ k 2 (8η)k. (D.14)
Choose λ such that
λ ≤ 1 (27ηpCd p2 e) 2 =
1 28(8η)2pC2d p2 e ≤ 1 2 8 k−1 (8η)2(pC2d p2 e ) 2 k−1
which is equivalent to
λ k−1 2 ≤ 1
24(8η)k−1pCd p2 e
= η
2(8η)kpCd p2 e
for all integer 2 ≤ k ≤ p. Then, since the following inequality can be obtained
p∑ k=2 pCkλ k 2 (8η)k ≤ p∑ k=2 pCd p2 eλ k 2 (8η)k
≤ 1 2 p∑ k=2 √ λη = p− 2
2
√ λη,
we have
E [ 1{|θλn|>M}|∆n| 2p ∣∣∣∣θλn] ≤ (1− η√λ)|θλn|2p (D.15) and
E [ 1{|θλn|>M}|∆n| 2p−2 ∣∣∣∣θλn] ≤ (1− η√λ)|θλn|2(p−2) ≤ 1M2 (1− η√λ)|θλn|2p. (D.16)
By combining (D.9), (D.16) and (D.15), we derive
E[1{|θλn|>M}|θ λ n+1|2p|θλn] ≤ (1− η √ λ)|θλn|2p
+ 22p−2p(2p− 1)λd
M2β (1− η
√ λ)|θλn|2p + 22p−3p(2p− 1)
( 2λ
β
)p E|ξn+1|2p
≤ (1− η √ λ) ( 1 +
22p−2p(2p− 1)λd M2β ) |θλn|2p
+ 22p−3p(2p− 1) ( 2λ
β
)p E|ξn+1|2p
≤ (1− η2λ)|θλn|2p + 22p−3p(2p− 1) ( 2λ
β
)p E|ξn+1|2p (D.17)
where we used the fact that M ≥ 2 2p−2p(2p−1)d
ηβ for the last inequality.
Consider the case of |θλn| ≤M . By observing that from (D.3)
1{|θ|≤M}λ 2|Hλ,c(θ, x)|2 ≤ λ ( 4d(1 + λ2max) + 2η 2M2 )
and
1{|θ|≤M}|2λ〈θ, Hλ,c(θ, x)〉| ≤ 2λ √ |θ| √ |Hλ,c(θ, x)|
≤ 2λ √ M √ |G(θ, x)|+ d √ λ+ 2ηM2r+1
≤ 2λ √ M √ |K(x)|(1 +Mq) + d √ λ+ 2ηM2r+1,
it can be shown that E [ 1{|θλn|≤M}|rn| k ∣∣∣∣θλn] = E[1{|θλn|≤M}(|2λ〈θλn, Hλ,c(θλn, Xn+1)〉|+ λ2|Hλ,c(θλn, Xn+1)|2)k∣∣∣∣θλn] ≤ E [ 1{|θλn|≤M} ( 2λ √ M √ K(Xn+1)(1 +Mq) + d √ λmax + 2ηM2r+1
+ λ ( 4d(1 + λ2max) + 2η 2M2 ))k∣∣∣∣θλn]
≤ D̃kλk where D̃k = 2k−1 ( (2 √ M)k(E[K(X0)](1 +Mq) + d √ λmax + 2ηM 2r+1)k/2 + (4d(1 + λ2max) +
2η2M2)k ) . Hence, one calculates that
E [ 1{|θλn|≤M}|∆n| 2p ∣∣∣∣θλn] ≤ |θλn|2p + p∑ k=1 ( p k ) |θλn|2(p−k)E[1{θλn≤M}|rn| k|θλn]
≤ (1− η2λ)|θλn|2p + η2λM2p +M2pλ p∑ k=1 ( p k ) λk−1D̃k
and
E [ 1{θλn≤M}|∆n| 2p−2 ∣∣∣∣θλn] ≤ p−1∑
k=0
( p− 1 k ) |θλn|2(p−1−k)E[1{|θλn|≤M}|rn| k|θλn]
≤ M2p−2 p−1∑ k=0 ( p k ) D̃kλ k.
Consequently, we obtain
E[1{|θλn|≤M}|θ λ n+1|2p|θλn] ≤ (1− η2λ)|θλn|2p + η2λM2p + λM2p p∑ k=1 ( p k ) λk−1D̃k
+ λd
β 22p−2p(2p− 1)M2p−2 p−1∑ k=0 ( p k ) λkD̃k (D.18)
+ 22p−3p(2p− 1) ( 2λ
β
)p E|ξn+1|2p.
By defining
Ap = η 2M2p +M2p p∑ k=1 ( p k ) λk−1maxD̃k
+ 22p−3p(2p− 1) ( 2dM2p−2
β
p−1∑ k=0 ( p k ) λkD̃k + 2 β ( 2λmax β )p−1 dp(2p− 1)!! ) ,
we conclude that
E|θλn+1|2p ≤ (1− η2λ)E|θλn|2p + λAp
≤ (1− η2λ)nE|θλ0 |2p + λAp ∞∑ j=0 (1− η2λ)j
≤ (1− η2λ)nE|θλ0 |2p + Ap η2 .
Proof of Lemma C.6. We begin by observing that
E|θ̄λt − ζ̄ λ,m t |2 = −2λ ∫ t mT E〈ζ̄λ,ms − θ̄λs , h(ζ̄λ,ms )−Hλ(θ̄λbsc, Xdse)〉ds
= −2λ ∫ t mT E〈ζ̄λ,ms − θ̄λs , h(ζ̄λ,ms )− h(θ̄λs )〉ds
− 2λ ∫ t mT E〈ζ̄λ,ms − θ̄λs , h(θ̄λs )− h(θ̄λbsc)〉ds
− 2λ ∫ t mT E〈ζ̄λ,ms − θ̄λs , h(θ̄λbsc)−H(θ̄ λ bsc, Xdse)〉ds
− 2λ ∫ t mT E〈ζ̄λ,ms − θ̄λs , H(θ̄λbsc, Xdse)−Hλ(θ̄ λ bsc, Xdse)〉ds
≤ 2λa ∫ t mT E|ζ̄λ,ms − θ̄λs |2ds
+ λa
2 ∫ t mT E|ζ̄λ,ms − θ̄λs |2ds+ ∫ t mT 2λ a E|h(θ̄λs )− h(θ̄λbsc)| 2ds
+ ∫ t mT ( − 2λE〈ζ̄λ,ms − θ̄λs , h(θ̄λbsc)−H(θ̄ λ bsc, Xdse)〉 ) ds
+ λa
2 ∫ t mT E|ζ̄λ,ms − θ̄λs |2ds+ ∫ t mT 2λ a E|H(θ̄λbsc, Xdse)−Hλ,c(θ̄ λ bsc, Xdse)| 2ds
= 3λa ∫ t mT E|ζ̄λ,ms − θ̄λs |2ds+ ∫ t mT Aλ,ms +B λ,m s +D λ,m s ds (D.19)
where we have used Proposition B.1 and the Young’s inequality in the first inequality and
Aλ,mt := 2λ
a E|h(θ̄λt )− h(θ̄λbtc)| 2
Bλ,mt := −2λE〈ζ̄ λ,m t − θ̄λt , h(θ̄λbtc)−H(θ̄ λ btc, Xdte)〉
Dλ,mt := 2λ
a E|H(θ̄λbtc, Xdte)−Hλ,c(θ̄ λ btc, Xdte)| 2.
In addition, from the definition of θ̄λt and (D.3), we have
|θ̄λt − θ̄λbtc| 4 ≤
( λ ∣∣∣∣ ∫ t btc Hλ,c(θ̄ λ bsc, Xdse)ds ∣∣∣∣+√2λβ−1|B̄λt − B̄λbtc|)4 ≤ 8λ2 ( λ2|Hλ,c(θ̄λbtc, Xdte)| 4 + 4β−2|B̄λt − B̄λbtc| 4
) ≤ 8λ2 ( (4d(1 + λ2) + 2η2|θ̄λbtc| 2)2 + 4β−2|B̄λt − B̄λbtc| 4
) ≤ λ225 ( 23d2(1 + λ4) + η4|θ̄λbtc| 4 + β−2|B̄λt − B̄λbtc| 4
) which yields √
E|θ̄λt − θ̄λbtc|4 ≤ C̃1λ (D.20)
where C̃1 = 25/2 √ 8d2(1 + λ4max) + η 4(E|θ0|4 +A2/η2) + 3β2 d2.
Using Proposition B.2, Aλ,mt can be bounded as follows:
Aλ,mt ≤ 2λ
a LXE[(1 + |θ̄λt |+ |θ̄λbtc) 2l|θ̄λt − θ̄λbtc| 2]
≤ 2λ a LX √ E(1 + |θ̄λt |+ |θ̄λbtc)4l √ E|θ̄λt − θ̄λbtc|4
≤ 2λ a LX3
2l √ (1 + E|θ̄λt |4l + E|θ̄λbtc|4l) √ E|θ̄λt − θ̄λbtc|4
≤ C1λ2 (D.21) where LX = L222ρ−1(1 + E|X0|2ρ) and C1 = 2aLX9 l √
(1 + 2E|θ̄λ0 |4l + 2A2lη2 )C̃1 and (D.20) is used for the last inequality.
To estimate Bλ,mt , we observe that
Bλ,mt = −2λE〈ζ̄ λ,m t − θ̄λbtc, h(θ̄ λ btc)−H(θ̄ λ btc, Xdte)〉
− 2λE〈θ̄λbtc − θ̄ λ t , h(θ̄ λ btc)−H(θ̄ λ btc, Xdte)〉 ≤ −2λE [ E [ 〈ζ̄λ,mt − θ̄λbtc, h(θ̄ λ btc)−H(θ̄ λ btc, Xdte)〉 ∣∣∣∣ζ̄λ,mt , θ̄λbtc]] − 2λE [ 〈θ̄λbtc − θ̄ λ t , h(θ̄ λ btc)−H(θ̄ λ btc, Xdte)〉
] ≤ −2λE
[ 〈λ ∫ t btc Hλ(θ̄ λ bsc, Xdse)ds− √ 2λ β B̃λt−btc, h(θ̄ λ btc)−H(θ̄ λ btc, Xdte)〉 ] ≤ −2λ2E [ 〈Hλ(θ̄λbtc, Xdte), h(θ̄ λ btc)−H(θ̄ λ btc, Xdte)〉
] ≤ 2λ2 √ E|Hλ(θ̄λbtc, Xdte)|2 √ E|h(θ̄λbtc)−H(θ̄ λ btc, Xdte)|2
≤ 2λ2 √
6E|K(X0)|2(1 + E|θ̄λbtc|2q) + 3λd+ 3η2E|θ̄ λ btc|4r+2 √ 4E|H(θ̄λbtc, Xdte)|2
≤ 4λ2 √
6E|K(X0)|2(1 + E|θ̄λbtc|2q) + 3λd+ 3η2E|θ̄ λ btc|4r+2 × √
4E|K(X0)|2(1 + E|θ̄λbtc|2q) + 2η2E|θ̄ λ btc|4r+2
≤ C2λ2 (D.22) where
C2 = 4 √ 6E|K(X0)|2(1 + E|θ0|2q +
Aq η2 ) + 3λd+ 3η2|θ̄λbtc|4r+2
× √ 4E|K(X0)|2(1 + E|θ0|2q +
Aq η2
) + 2η2 ( E|θ̄λ0 |4r+2 +
A2r+1 η2
) .
Note that we have used the independence of θ̄λbsc and Xdse to obtain the second inequality, and used Remark B.4 and Lemma C.2 to calculate the bound of E|Hλ(θ̄λbtc, Xdte)| 2 and E|H(θ̄λbtc, Xdte)| 2.
Moreover, Dλ,mt can be estimated as follows, from Remark B.3,
Dλ,mt ≤ 6λ2
a
[ 8E|K(X0)|4(1 + E|θ̄λbtc| 4q) + d+ η2E|θ̄λbtc| 8r+2 ] ≤ C3λ2 (D.23)
where the independence of θ̄λbsc and Xdse is used and C3 is given by
C3 = 6
a
[ 8E|K(X0)|4(1 + E|θ̄λ0 |4q +A2q/η2) + d+ η2(E|θ̄λ0 |8r+2 +A4r+1/η2) ] .
Plugging (D.21), (D.22) and (D.23) into (D.19), one can derive
E|θ̄λt − ζ̄ λ,m t |2 ≤ 3λa ∫ t mT E|θ̄λs − ζ̄λ,ms |2ds+ ∫ t nT (C1 + C2 + C3)λ 2ds
≤ 3λa ∫ t mT E|θ̄λs − ζ̄λ,ms |2ds+ (C1 + C2 + C3)λ <∞
where the second inequality follows from the fact that (t −mT ) ≤ T = 1λ and the use of Grownwall’s inequality gives
E|θ̄λt − ζ̄ λ,m t |2 ≤ cλ
where c = e3a(C1 + C2 + C3).
Proof of Lemma C.7. Since Zλt = ζ̄ λ,0 t and t ∈ [mT, (m+ 1)T ), we can write
W1 ( L ( ζ̄λ,mt ) ,L ( Zλt )) ≤ m∑ k=1 W1 ( L ( ζ̄λ,kt ) ,L ( ζ̄λ,k−1t )) ≤
m∑ k=1 w1,2 ( L ( ζ̄λ,kt ) ,L ( ζ̄λ,k−1t )) where we have used the fact W1(µ, ν) ≤ w1,2(µ, ν) for µ, ν ∈ PV2 for the second inequality. Using Remark C.5 and λ(t− kT ) ≥ m− k, we further have
w1,2 ( L ( ζ̄λ,kt ) ,L ( ζ̄λ,k−1t )) ≤ ĉe−ċλ(t−kT )w1,2 ( L(θ̄λkT ),L(ζ̄ λ,k−1 kT ) ) ≤ ĉe−ċ(m−k)w1,2 ( L(θ̄λkT ),L(ζ̄ λ,k−1 kT )
) ≤ ĉe−ċ(m−k)W2 ( L(θ̄λkT ),L(ζ̄ λ,k−1 kT ) )√ E ∣∣∣∣1 + V2(θ̄λkT ) + V2(ζ̄λ,k−1kT )∣∣∣∣2
≤ ĉe−ċ(m−k)W2 ( L(θ̄λkT ),L(ζ̄ λ,k−1 kT ) ) × [ 1 + √ E[V4(θ̄λkT )] + √ E[V4(ζ̄λ,k−1kT )]
] ≤ ĉe−ċ(m−k) √ λ √ e3a(C1 + C2 + C3) [ 1 + √ 2E|θ0|4 + 2 + 2
A2 η2
+ √ 2E|θ0|4 + 2 + 2
A2 η2 + c̃(4) c̄(4)
] (D.24)
where the Cauchy-Schwarz inequality is applied to the third inequality and Lemma C.3, C.6 and C.4 are used for the last inequality. By combining the two inequalities above, we obtain
W1 ( L ( ζ̄λ,mt ) ,L ( Zλt )) ≤ ĉ √ λ √ e3a(C1 + C2 + C3) [ 1 + √ 2E|θ0|4 + 2 + 2
A2 η2
+ √ 2E|θ0|4 + 2 + 2
A2 η2 + c̃(4) c̄(4) ] m∑ k=1 e−ċ(m−k)
≤ z1 √ λ
where z1 = ĉ1−exp(−ċ) √ e3a(C1 + C2 + C3) [ 1 + √ 2E|θ0|4 + 2 + 2A2η2 +√
2E|θ0|4 + 2 + 2A2η2 + c̃(4) c̄(4)
] .
Proof of Lemma C.8. We begin by observing that
W2 ( L ( ζ̄λ,kt ) ,L ( ζ̄λ,k−1t )) ≤ √ 2w1,2 ( L ( ζ̄λ,kt ) ,L ( ζ̄λ,k−1t )) ≤ λ1/4e−ċ(m−k)/2 [ ĉ √ e3a(C1 + C2 + C3) ( 1 + √ 2E|θ0|4 + 2 + 2
A2 η2
+ √ 2E|θ0|4 + 2 + 2
A2 η2 + c̃(4) c̄(4) )]1/2 where we have used the factW2 ≤ √ 2w1,2 for the first inequality, and the second inequality follows from (D.24). Consequently, we derive
W2 ( L ( ζ̄λ,mt ) ,L ( Zλt )) ≤ m∑ k=1 W2 ( L ( ζ̄λ,kt ) ,L ( ζ̄λ,k−1t ))
≤ λ1/4 [ ĉ √ e3a(C1 + C2 + C3) ( 1 + √ 2E|θ0|4 + 2 + 2
A2 η2
+ √ 2E|θ0|4 + 2 + 2
A2 η2 + c̃(4) c̄(4) )]1/2 m∑ k=1 e−ċ(m−k)/2
≤ λ1/4z2 where
z2 =
√ ĉe3a/4(C1 + C2 + C3) 1/4
1− exp(−ċ/2)
[( 1+ √ 2E|θ0|4 + 2 + 2
A2 η2 +
√ 2E|θ0|4 + 2 + 2
A2 η2 + c̃(4) c̄(4)
)]1/2 .

heading: E DETAILS OF EXPERIMENTS
text: E.1 IMAGE CLASSIFICATION
The experiments are exactly replicated in the official implementation of Zhuang et al. (2020). More specifically, VGG11, ResNet34 and DenseNet121 are trained for 500 epochs. We apply a weight deacy of 0.0005 and decay the initial learning rate by 10 after 150 epochs to all optimizers. Batch normalization proposed in Ioffe & Szegedy (2015) is employed to prevent the models from overfitting and boost the training speed for all three models. The batch size is 128.
Regarding hyperparameter values of Adam, AdaBelief, AdamP, AdaBound, AMSGrad and RMSProp, the best hyperparameters are used across all the experiments in Luo et al. (2019) and Zhuang et al. (2020), which are λ = 0.001, β1 = 0.9, β2 = 0.999 and ε = 10−8. For SGD, we set the momentum to 0.9 for SGD. For THεO POULA, the best hyperparameters are λ = 0.1, ε = 0.1 and β = 1012.
Figure 2 shows test accuracy for VGG11, ResNet34 and DenseNet121 on CIFAR-10 and CIFAR100.
E.2 LANGUAGE MODELING
We conduct language modeling for Penn Treebank (PTB) data set. For this task, we train the AWD LSTM of Merity et al. (2018) for 750 epochs. The details of models can be found in the official implementation of AWD-LSTM 1.
For NT-ASGD and averaged THεO POULA, the constant learning rate of 30 is used for 2 and 3- layer LSTMs. For 1-layer LSTMs, we set to 10. ε = 100 and β = 1010 are set across all the experiments.
1https://github.com/salesforce/awd-lstm-lm
For AdaBelief, we used the best hyperparameters reported in Zhuang et al. (2020). Thus, we obtain the same results of AdaBelief in Zhuang et al. (2020). Specifically, we set β1 = 0.9, β2 = 0.999, ε = 10−12 and an initial learning rate of 0.01 for 2 and 3-layer LSTMs. λ = 0.001 and ε = 10−16 are used for 1-layer LSTMs.
The averaging is triggered when no improvement has been made for 5 consecutive epochs for THεO POULA. Also, AdaBelief uses a development-based learning rate decay, which decreases the learning rate by a constant factor δ if the model does not attain a new best value for k epochs. We have searched the optimal hyperparameters for the development-based learning rate decay among δ = {0.1, 0.5} and k = {5, 10, 20}. We have found δ = 0.1 and k = 5 yield the best performance. Figure 3 displays test perplexity for different AWD-LSTM models on PTB.
E.3 EFFECTIVENESS OF THE BOOSTING FUNCTION
This subsection empirically tests the effectiveness of the boosting function in our algorithm. THεO POULA without the boosting function updates the parameter as follows:
θλ0 := θ0, θ λ n+1 := θ λ n − λHλ,c(θλn, Xn+1) + √ 2λβ−1ξn+1, n ∈ N,
where Hλ,c := (H (1) λ,c(θ, x), · · · , H (d) λ,c(θ, x)) T is given by
H (i) λ,c(θ, x) =
G(i)(θ, x)
1 + √ λ|G(i)(θ, x)|
+ η θ(i)|θ|2r
1 + √ λ|θ|2r , (E.1)
and {ξn}n≥1 is a sequence of independent standard d-dimensional Gaussian random variables. Indeed, this is a special case of THεO POULA with ε = ∞. We train VGG11, ResNet34 and DenseNet121 on CIFAR-10 and CIFAR-100 using the iterating rule of (E.1). The hyperparameters are the same with the best hyperparameters of THεO POULA. As Table 4 shows, THεO POULA without the boosting function is worsen than THεO POULA. The result indicates that the addition of the boosting function leads to a meaningful increase of test accuracy, validating that the boosting function works well for the sparsity of neural networks as expected.

heading: F ADDITIONAL EXPERIMENTS
text: F.1 EFFECT OF β ON THE PERFORMANCE OF THεO POULA
This subsection examines the effect of β on the performance of THεO POULA. We conduct experiments for VGG11 and ResNet34 on CIFAR10 and CIFAR100 with different values of β ranging from 104 to 1012. As shown in Table 5, THεO POULA achieves the highest accuracy when β is 108 ∼ 1012, which is consistent with the phenomenon, so called the cold posterior effect, see Wenzel et al. (2020) and Aitchison (2021).
F.2 EXPERIMENTS WITH η 6= 0
We conduct additional experiments with η 6= 0 to demonstrate that there is no gap between theory and practice of our work. When the regularization parameter r is large (possibly, overestimated) and the dimension of θ is big, |θ|2r becomes substantially huge. As a result, the stochastic gradient of the regularization term, η θ (i)|θ|2r
1+ √ λ|θ|2r , in (8) will approximately behave like η√ λ θ(i), which is equivalent
to `2-regularization. In all the numerical experiments in Table 2, we applied a weight decay with 5× 10−4 for image classification and with 1.2× 10−6 for language modeling. That is, by choosing η = 5× 10−4 √ λ and large r, one can obtain accuracy of models with `2-regularization.
Table 6 shows that the accuracy for VGG, ResNet and DenseNet on CIFAR-10 and CIFAR-100 with r = 10 and η = 5 × 10−4 √ λ. One observes a very similar performance by THεO POULA as in Table 2 without any noticeable loss of accuracy.
F.3 SAMPLING FROM A SYNTHETIC MULTI-MODAL DISTRIBUTION
We provide a brief overview of recent progress on sampling and Bayesian neural networks to describe the potential of THεO POULA as a sampling algorithm. Deng et al. (2020b) proposed an adaptive MCMC algorithm, called CSGLD, which uses a scalable dynamic importance sampler to flatten the target distribution and reduce the energy barriers to escape local optima. Deng et al. (2020a) devloped reSGMCMC motivated by replica exchange monte carlo algorithm. In particular, reSGMCMC obtained the state-of-the art results on CIFAR-10, CIFAR-100 and SVHN in Bayesian neural networks. Zhang et al. (2020) developed cyclical stochastic gradient MCMC with a cyclical stepsize schedule.
We replicate the simulation of a synthetic multi-modal distribution in Deng et al. (2020b) to evaluate the performance of THεO POULA as a sampling method. A target distribution is π(x) ∝ e−U(x)
where U(x) = ∑2 i=1 x(i)2−10 cos(1.2πx(i)) 3 and x = (x(1), x(2)). Detail of the setting in the experiment such as hyperparameters, regularizer, training epochs can be found in Appendix D.3 of Deng et al. (2020b). THεO POULA is compared with SGLD, CSGLD (Deng et al. (2020b)), reSGLD (Deng et al. (2020a)), cycSGLD (Zhang et al. (2020)). For THεO POULA, we used λ = 0.05, ε = 1 and T = 0.3. A resampling scheme is used for CSGLD. Figure 4 illustrates that THεO POULA recovers the target multi-modal distribution successfully without the local trap issue observed in SGLD and cycSGLD.

heading: G TABLE OF CONSTANTS
text: Table 7 displays full expressions for constants which appear in the main results of this paper. In addition, Table 8 shows all main constants and their dependencies on key parameters such as d, β, the moments of K(X0) and η.

Review summaries:
The goal of this work is to bridge the gap between the theoretical qualities of Langevin based algorithms for optimisation and their empirical inability to match SOTA optimisation methods such as SGD or ADAM. To begin with, the authors propose a new iterative scheme, similar to TUSLA [1], involving an adaptive preconditioned for the gradient term only, and a temperature parameter for the noise term (equations 7-8). They show convergence of the proposed scheme, in a very similar way to what is done in [1]. Finally, they test a modified version of their algorithm, in which the returned value is the average value of past iterates, and show competitive performance compared to SOTA algorithms.  [1] A. Lovas, I. Lytas, M. Rasonyi, and S. Sabanis. Taming neural networks with tusla: Non-convex learning via adaptive stochastic gradient langevin algorithms.

The paper presents a new stochastic optimization algorithm inspired by the tamed Euler discretization of the overdamped Langevin diffusion. Convergence guarantees in the Wasserstein metric are established for this algorithm under mild assumptions, primarily local Lipschitz continuity and dissipativity conditions on the (stochastic) gradient of the loss function. In addition to a taming function, the proposed algorithm "Th$\varepsilon$O POULA", employs a boosting function to improve performance in flat regions of the loss landscape. Along with a synthetic example, the authors demonstrate the effectiveness of Th$\varepsilon$O POULA for training deep neural networks for image classification and language modelling.

The author(s) propose the TheoPouLa algorithm to overcome unstable gradient (exploding / diminishing) problems by utilizing Euler's polygonal approximation. Several experiments are presented on image classification and language modeling to justify the effectiveness of these models. The biggest contribution of this paper is to attempt to tackle stability issues generated by the fast-growing gradients and non-global Lipschitz constants.