Acceleration of stochastic gradient descent with momentum by averaging: finite-sample rates and asymptotic normality

Kejie Tang ¹¹1School of Mathematical Sciences, Shanghai Jiao Tong University Weidong Liu ²²2School of Mathematical Sciences, MoE Key Lab of Artificial Intelligence, Shanghai Jiao Tong University Yichen Zhang³³3Daniels School of Business, Purdue University Xi Chen⁴⁴4Stern School of Business, New York University

Abstract

Stochastic gradient descent with momentum (SGDM) has been widely used in many machine learning and statistical applications. Despite the observed empirical benefits of SGDM over traditional SGD, the theoretical understanding of the role of momentum for different learning rates in the optimization process remains widely open. We analyze the finite-sample convergence rate of SGDM under the strongly convex settings and show that, with a large batch size, the mini-batch SGDM converges faster than the mini-batch SGD to a neighborhood of the optimal value. Additionally, our findings, supported by theoretical analysis and numerical experiments, indicate that SGDM permits broader choices of learning rates. Furthermore, we analyze the Polyak-averaging version of the SGDM estimator, establish its asymptotic normality, and justify its asymptotic equivalence to the averaged SGD. The asymptotic distribution of the averaged SGDM enables uncertainty quantification of the algorithm output and statistical inference of the model parameters.

1 Introduction

In this paper, we are interested in solving the following stochastic optimization problem:

\displaystyle x^{*}=\text{argmin}_{x\in\mathbb{R}^{d}}\mathbb{E}[f_{\xi}(x)],

(1.1)

where $f_{\xi}(x)$ is a stochastic convex loss function, $\xi$ is a random element, and $\mathbb{E}[\cdot]$ is the expectation respected to $\xi$ . Stochastic optimization plays an important role in many statistics and machine learning problems. For a wide range of applications, the objective function is strongly convex, and stochastic gradient descent (SGD) is often preferred over gradient descent (GD) (see, e.g., Nesterov, 2003; Nocedal and Wright, 2006) due to its computational advantage. SGD is a first-order optimization algorithm that approximates the expected loss by averaging the loss function over a mini-batch of training examples. At each iteration, the algorithm updates the model parameters in the direction of the negative gradient of the mini-batch loss, scaled by a learning rate parameter.

While SGD is simple and easy to implement, it may suffer from slow convergence rates or oscillations in high-dimensional optimization problems, particularly when the loss function is noisy and ill-conditioned. Momentum-based methods enhance SGD by introducing an exponentially weighted moving average of the past gradients to the update rule, which serves to dampen oscillations and accelerate convergence, allowing the algorithm to maintain a more consistent direction of movement even in the presence of noisy gradients. SGD with momentum (SGDM) has become increasingly popular in modern applications, e.g., large-scale deep neural networks. Evident numerical studies have shown that the use of momentum-based optimization methods improves the convergence rate and the generalization performance, as well as reduces the sensitivity to the choice of hyperparameters. However, the theoretical analysis of SGDM is still an active area of research.

In this paper, we focus on the theoretical analysis from an optimization perspective under strong convexity that ensures the existence of a unique and well-defined global minimizer. Our goal is to establish the convergence properties of SGDM to provide insights into the role of momentum and other hyperparameters in the optimization process. For solving (1.1), SGDM updates the target estimator by a weighted combination with historical gradients

	$\displaystyle m_{t+1}$	$\displaystyle=$	$\displaystyle\gamma m_{t}+(1-\gamma)\triangledown f_{\xi_{t}}(x_{t}),$		(1.2)
	$\displaystyle x_{t+1}$	$\displaystyle=$	$\displaystyle x_{t}-\alpha m_{t+1},$		(1.3)

where $\xi_{t}$ , $t\geq 1$ are i.i.d. samplings from $\xi$ , $\gamma$ is the momentum weight and $\alpha$ is the learning rate. The classical SGD is a special case of SGDM with $\gamma=0$ and $m_{t+1}=\nabla f_{\xi_{t}}(x_{t})$ .

In practice, a large momentum weight is often placed to accelerate the algorithm, for example, $\gamma=0.9$ . That being said, it is still an open problem to investigate the theoretical properties of SGDM with a general specification of $\gamma$ . The theoretical analysis in existing studies has not given an affirmative answer to the open problem by Liu et al. (2020) and the assertion by Kidambi et al. (2018). In this paper, we consider the mini-batch SGDM with batch size $B$ . We establish finite sample rates for SGDM and the averaging of its trajectories under smooth and strongly convex loss functions. By our results, we can give some partial answers to these open questions and some related contributions in several aspects.

1.

Drawing upon the assertion presented in Kidambi et al. (2018), this study addresses the open question posited by Liu et al. (2020) by demonstrating that mini-batch SGDM, with appropriate momentum weights, converges to a local neighborhood of the minimum with a faster rate than SGD. We rigorously establish the finite-sample convergence rate, and we further provide an adaptive choice of the momentum weight which theoretically attains the optimal convergence rate. Our experiments support the theoretical results of convergence and the optimal momentum weight.
2.

We provide a non-asymptotic analysis of the Polyak-Ruppert averaging version of SGDM. The averaging of SGDM trajectories $\bar{x}_{t}=\frac{1}{t-n_{0}}\sum_{i=n_{0}+1}^{t}x_{i}$ can accelerate SGDM with a wide range of learning rates $\alpha$ to the rate $\mathcal{O}(\frac{\sigma}{\sqrt{t}})$ , where $\sigma^{2}$ is the variance of stochastic gradient. Furthermore, we show theoretically that averaged SGDM converges faster than averaged SGD in the early phases of the iteration process. Moreover, through both theoretical analysis and numerical experiments, we demonstrate that SGDM is less sensitive to the choice of learning rates, and in addition to this, the averaged SGDM is less sensitive to the start of the averaging iteration $n_{0}$ .
3.

We further establish the asymptotic normality for the averaged $\bar{x}_{t}$ as $t\rightarrow\infty$ with a decaying learning rate $\alpha$ or a diverging batch size $B$ . Particularly, as $\alpha=\Theta(t^{-\epsilon})$ for $\epsilon\in(\frac{1}{2},1)$ , the asymptotic normality holds for averaged SGDM under strongly convex loss functions. The asymptotic covariance matrix depends only on the Hessian and Gram matrix of the stochastic gradients at $x^{*}$ , and the batch size $B$ . Interestingly, mini-batch averaged SGDM is asymptotically equivalent to mini-batch averaged SGD as $t\rightarrow\infty$ . We further demonstrate that the optimal learning rates of averaged SGDM that correspond to asymptotic normality are $\alpha=\Theta(t^{-1/2})$ for quadratic losses and $\alpha=\Theta(t^{-2/3})$ for general strongly convex losses. To our best knowledge, this is the first work to analyze the convergence of averaged SGDM to asymptotic normality under mini-batching. The results enable us to perform uncertainty quantification for the algorithm outputs of the averaged SGDM algorithm $\bar{x}_{t}$ , and statistical inference for model parameters $x^{*}$ .

In addition to the convergence rate under pre-specified learning rates, the selection of appropriate learning rates is a critical aspect when evaluating the performance of optimizers, as either an excessively high or low learning rate can detrimentally affect convergence. We study the sensitivity of the convergence over different learning rates. Extensive research has been conducted to tackle this challenge and accelerate convergence, such as Zeiler (2012); Kingma and Ba (2014). Recent results (Paquette and Paquette, 2021; Bollapragada et al., 2022) demonstrated the acceleration of SGDM on quadratic forms, however, under a shrinking range of allowable learning rates for larger momentum. This finding contradicts the advantages of SGDM and indicates an increased sensitivity to the choice of learning rates in theoretical analysis.

1.1 Related Works

Since the seminal work by Robbins and Monro (1951), the convergence properties of SGD have been extensively studied in the literature (See, e.g., Moulines and Bach, 2011; Bottou et al., 2018; Nguyen et al., 2018). On the contrary, the theory for SGDM has rarely been explored until recently, although it is very popular in training many modern machine learning models such as neural networks to improve the training speed and accuracy of various models.

Kidambi et al. (2018) showed that (1.2) cannot achieve any improvement over SGD for a specially constructed linear regression problem. Together with some numerical experiments, they asserted that the only reason for the superiority of stochastic momentum methods in practice is mini-batching. In a recent paper, Liu et al. (2020) proved SGDM can be as fast as SGD. They established the identical convergence bound as SGD. They also posed an open problem of whether it is possible to show that SGDM converges faster than SGD for special objectives such as quadratic ones. Other studies on the theoretical analysis of SGDM can be found in Loizou and Richtárik (2017), Loizou and Richtárik (2020), Gitman et al. (2019) for linear system and quadratic loss; Sebbouh et al. (2021) for almost sure convergence rates under smooth and convex loss functions with a time-varying momentum weight. Mai and Johansson (2020) studied a class of general convex loss and obtained convergence rates of time averages regardless of the momentum weight. We refer the readers to the latter paper and Liu et al. (2020) for a few more papers on the convergence rate for SGDM. Richtárik and Takác (2020) considered the problem of solving a consistent linear system $Ax=b$ on least square regression. SGDM is also called the stochastic heavy-ball method, originated from Polyak’s heavy-ball method for deterministic optimization (Polyak, 1964). Some works (Loizou and Richtárik, 2017; Kidambi et al., 2018; Loizou and Richtárik, 2020; Paquette and Paquette, 2021; Bollapragada et al., 2022; Lee et al., 2022) based on the stochastic heavy ball method (SHB) achieved the linear convergence rate. Yang et al. (2016); Defazio (2020); Jin et al. (2022); Li et al. (2022) investigated the impact of momentum on the convergence properties of non-convex optimization problems and provided insights into its practical use. Gitman et al. (2019) focused on the noise reduction properties of momentum.

The averaged SGDM is not widely analyzed in the literature, but averaging tools have been studied for SGD and its variants since Ruppert (1988); Polyak and Juditsky (1992); Moulines and Bach (2011) to accelerate the convergence and establish the asymptotic normality. Building statistical inference and uncertainty quantification of SGD iterates has been an emerging topic, and an extensive list of literature follows, including Chen et al. (2020); Su and Zhu (2023); Zhu et al. (2023) that studied averaged SGD and provided inference procedures based on plug-in, batch-means, and tree-based construction of the confidence intervals, respectively. Zhu and Dong (2021); Lee et al. (2022) applied process-level function central limit theorem and utilized it to construct confidence regions seamlessly. Toulis et al. (2021); Chen et al. (2023) studied the uncertainty quantification of the implicit and gradient-free variants of the SGD procedures.

2 Preliminaries

In this section, we present the mini-batch SGDM settings considered in this paper. Particularly, we define a mini-batch stochastic loss

\displaystyle g_{\eta_{t}}(x)=\frac{1}{B}\sum_{i=1}^{B}f_{\xi_{ti}}(x),

(2.4)

where $\eta_{t}=\{\xi_{t1},...,\xi_{tB}\}$ is a mini-batch of size $B$ and elements $\xi_{ti}$ are i.i.d. sampled from the distribution of $\xi$ . The update of $m_{t}$ in (1.2) uses the mini-batch stochastic gradient $\triangledown g_{\eta_{t}}(x_{t})$ instead of the individual stochastic gradient $\triangledown f_{\xi_{ti}}(x_{t})$ .

We first introduce some regularity assumptions and discuss their use in the theoretical results. We denote the $\ell_{2}$ -norm of a vector as $\left\|\cdot\right\|$ and the operator norm as $\left\|A\right\|=\max_{\left\|x\right\|=1}\left\|Ax\right\|$ .

(A1). Assume that the loss function $f_{\xi}(x)$ is twice differentiable. Let $\kappa_{1}\leq\cdots\leq\kappa_{d}$ be the eigenvalues of Hessian matrix $\Sigma=\mathbb{E}[\triangledown^{2}f_{\xi}(x^{*})]$ . Assume that $\mu:=\kappa_{1}>0$ , and $L:=\kappa_{d}<\infty.$

(A2). There exists a constant $\overline{L}\geq 0$ such that the Hessian matrix of the loss $f(x):=\mathbb{E}[f_{\xi}(x)]$ , $\Sigma(x)=\triangledown^{2}f(x)$ satisfies $\left\|\Sigma(x)-\Sigma(x^{*})\right\|\leq\overline{L}\left\|x-x^{*}\right\|$ , $\forall x$ .

(A3). Define $L_{\xi}=\sup_{x}\frac{\left\|\nabla f_{\xi}(x)-\nabla f_{\xi}(x^{*})\right\|}{% \left\|x-x^{*}\right\|}$ . Assume that $\mathbb{E}[\left\|\triangledown f_{\xi}(x^{*})\right\|^{2}]\leq\sigma^{2}$ and $\mathbb{E}[L_{\xi}^{2}]\leq L_{f}^{2}$ .

(A3’). Assume that $\triangledown f_{\xi}(x^{*})$ and $L_{\xi}$ satisfy

\displaystyle\sup_{\left\|v\right\|=1}\mathbb{E}\left[\exp\left(\left|v^{\top}% \triangledown f_{\xi}(x^{*})\right|/\sigma\right)\right]\leq 2,\quad\mathbb{E}% \left[\exp\left(L_{\xi}/L_{f}\right)\right]\leq 2.

A few discussions follow concerning the assumptions above. Assumption (A1) is a regularity condition on the smoothness and strong convexity of the loss function $f_{\xi}(x)$ at $x^{*}$ . Beyond that, (A2) assumes the Lipschitz condition on the Hessian matrix of the loss function. Notably, the $\Sigma$ defined in (A1) is a brief notation for $\mathbb{E}[\Sigma_{\xi}(x^{*})]$ defined in (A2). For quadratic losses, the Hessian matrix $\Sigma(x)$ is identical for different $x$ , and therefore (A2) holds with $\overline{L}=0$ . For general strongly convex losses, we illustrate (A2) under a logistic regression setting in Example 1 below. In the following, we consider two scenarios, quadratic losses, and general (strongly convex) losses.

Assumptions (A3) and (A3’) are two separate conditions on the smoothness of quadratic and general stochastic loss functions, respectively. Assumption (A3) is weaker than (A3’) and indeed weaker than many in the literature, such as those in Yan et al. (2018); Liu et al. (2020) that $\mathbb{E}[\left\|\nabla f_{\xi}(x)-\nabla f(x)\right\|^{2}]\leq\sigma^{2}$ for all $x$ , which does not hold for linear regression problems with unbounded domain of $x$ . On the other hand, Bottou et al. (2018); Wang and Johansson (2022) assumed that $\mathbb{E}[\left\|\nabla f_{\xi}(x)-\nabla f(x)\right\|^{2}]\leq\rho\left\|% \nabla f(x)\right\|^{2}+\sigma^{2}$ for $\rho>0$ . Meanwhile, our assumptions lead to $\mathbb{E}[\left\|\nabla f_{\xi}(x)-\nabla f(x)\right\|^{2}]\leq\frac{3(L_{f}^% {2}+L^{2})}{\mu^{2}}\left\|\nabla f(x)\right\|^{2}+3\sigma^{2},$ due to the triangle inequality and strongly convexity. For any bounded domain $\{\left\|x\right\|\leq r\}$ , (A3) and (A3’) are easily satisfied for bounded gradients and smoothness. For unbounded problems, they are usually satisfied given certain design properties in many popular statistical models. For instance, we illustrate the assumptions under logistic regression in the following example.

Example 1

Consider an $\ell_{2}$ -regularized logistic regression model with samples $\{a_{\xi},b_{\xi}\}$ such that $a_{\xi}\in\mathbb{R}^{d}$ and $b_{\xi}\in\{0,1\}$ are generated by $b_{\xi}=1$ with probability $p_{x}(a_{\xi})$ and $b_{\xi}=0$ otherwise, where $p_{x}(a)=1/(1+\exp(-x^{\top}a))$ . We consider the loss function defined as:

\displaystyle f_{\xi}(x)=-b_{\xi}\log(p_{x}(a_{\xi}))-(1-b_{\xi})\log(1-p_{x}(% a_{\xi}))+\frac{\nu}{2}\left\|x\right\|^{2}.

The gradient and the Hessian matrix are

\displaystyle\nabla f_{\xi}(x)=(p_{x}(a_{\xi})-b_{\xi})a_{\xi}+\nu x,\quad% \Sigma(x)=\mathbb{E}[p_{x}(a_{\xi})(1-p_{x}(a_{\xi}))a_{\xi}a_{\xi}^{\top}]+% \nu I_{d}.

By the fact that $\left|p_{x}(a)-p_{y}(a)\right|\leq\frac{\sqrt{3}}{18}\left\|a\right\|\left\|x-% y\right\|$ , (A2) is satisfied with $\overline{L}=\frac{\sqrt{3}}{6}\mathbb{E}\left\|a_{\xi}\right\|^{3}+\nu.$ Suppose $\{a_{\xi}\}$ in logistic regression satisfy $\mathbb{E}[a_{\xi}]=0_{d}$ and $\sup_{\left\|v\right\|=1}\mathbb{E}\big{[}\exp\big{(}\left|v^{\top}a_{\xi}% \right|^{2}/\sigma^{2}\big{)}\big{]}\leq 2$ , we have that $\left|p_{x}(a)-b\right|\leq 2$ , $\left\|(p_{x}(a)-p_{y}(a))a\right\|\leq\frac{\sqrt{3}}{18}\left\|a\right\|^{2}% \left\|x-y\right\|$ , and therefore,

\displaystyle\sup_{\left\|v\right\|=1}\mathbb{E}\left[\exp\left(\left|v^{\top}% \triangledown f_{\xi}(x^{*})\right|/(c\sigma)\right)\right]\leq 2,~{}\mathbb{E% }\left[\exp\left(L_{\xi}/(cL_{f})\right)\right]\leq 2.

for some absolute constant $c>0$ and $L_{f}=\mathbb{E}\left\|a_{\xi}\right\|^{2}+\nu$ .

3 Finite-sample Convergence Rates for SGDM

In this section, we first present the finite-sample convergence results for SGDM with general momentum weight $\gamma$ . We consider the two cases separately: $\overline{L}=0$ corresponds to the quadratic losses, and $\overline{L}>0$ corresponds to general strongly convex losses.

3.1 The finite-sample rates for SGDM on quadratic losses

We first establish the finite-sample rates under quadratic losses, where

\displaystyle f_{\xi}(x)=\frac{1}{2}x^{\top}A_{\xi}x-b^{\top}_{\xi}x+c.

(3.5)

The Hessian matrix of the loss function $\Sigma=\mathbb{E}[A_{\xi}]$ . The stochastic loss is $g_{\eta_{t}}(x)=\frac{1}{B}\sum_{i=1}^{B}f_{\xi_{ti}}(x)$ in the mini-batch setting. The following theorem shows the convergence rate of the last iterate $x_{t+1}$ .

Theorem 1

Under (A1)-(A3) and $\overline{L}=0$ , for any momentum $\gamma\in[0,1)$ and fixed $\delta\in(0,1]$ , assume the learning rate $\alpha>0$ satisfies $\alpha L<2(1+\gamma)/(1-\gamma)$ and $16M^{2}\alpha^{2}L_{f}^{2}\leq B\delta\lambda^{2(1-\delta)}(1-\lambda)$ , where

\displaystyle M=\frac{4}{\sqrt{\Delta}}\left(2(1-\gamma)(1+\alpha L+L)+3\alpha% \gamma\right),~{}\Delta=\min_{k}\left\{\left|\left(\gamma+1-\alpha(1-\gamma)% \kappa_{k}\right)^{2}-4\gamma\right|\right\}>0,

(3.6)

and $\lambda$ is the spectral radius of the matrix

\displaystyle\Gamma=\left(\begin{array}[]{cc}\gamma I,&(1-\gamma)\Sigma\\ -\alpha\gamma I,&I-\alpha(1-\gamma)\Sigma\\ \end{array}\right).

(3.9)

Let $\widetilde{m}_{t+1}=(1-\gamma)\sum_{j=1}^{t}\gamma^{t-j}\Sigma(x_{j}-x^{*})$ , we have for $t\geq 1$ ,

\displaystyle\mathbb{E}[\left\|\widetilde{m}_{t+1}\right\|^{2}+\left\|{x}_{t+1% }-x^{*}\right\|^{2}]\leq 2M^{2}\left(\frac{4}{B(1-\lambda)}\alpha^{2}\sigma^{2% }+\left\|x_{1}-x^{*}\right\|^{2}\lambda^{2(1-\delta)t}\right).

(3.10)

Theorem 1 provides a finite-sample bound simultaneously for the error of the last iterate $x_{t+1}-x^{*}$ , and a weighted average of the stochastic gradients $\widetilde{m}_{t+1}$ . The first term in the error bound (3.10), $\frac{8M^{2}}{B(1-\lambda)}\alpha^{2}\sigma^{2}$ , corresponds to a non-decaying bias, due to the noisy observation in the stochastic gradient. The bias term is proportional to the squared learning rate $\alpha^{2}$ and the variance of the stochastic gradient $\sigma^{2}$ , which is the same as the one in SGD. The second term is exponentially decaying when $\lambda<1$ , and establishes the convergence of the SGDM algorithm from any initialization $x_{1}$ to the true solution $x^{*}$ . The linear convergence convergence is up to a neighborhood of $x^{*}$ with size $\sqrt{\frac{8M^{2}}{B(1-\lambda)}}\alpha\sigma$ . With a larger batch size $B$ or a smaller $\sigma$ , the size of the neighborhood will be smaller, which aligns with the experimental findings reported in Kidambi et al. (2018) that the superiority of momentum methods is mainly due to mini-batching. Meanwhile, the second term in (3.10), $\left\|x_{1}-x^{*}\right\|^{2}\lambda^{2(1-\delta)t}$ , remains important to determine the convergence rate in the initial stage. Subsequently, we will demonstrate that this term, particularly the quantity $\lambda$ in SGDM, is improved compared to the $\lambda$ in SGD.

3.1.1 Linear convergence to a local neighborhood

The second term in (3.10) determines a linear convergence of SGDM to a local neighborhood of $x^{*}$ determined by its first term. The rate of this linear convergence is determined by $\gamma$ , the spectral radius of the matrix (3.9). We first provide some intuition in deriving $\lambda$ . Considering the noiseless setting where the stochastic gradient is the same as the true gradient, i.e., $\nabla f_{\xi}(x)=\nabla f(x)$ for any $x$ and $\xi$ . We have $\nabla g(x)=\Sigma(x-x^{*})$ and the SGDM updating rule (1.2) can be rewritten as

\displaystyle\left(\begin{array}[]{c}m_{t+1}\\ x_{t+1}-x^{*}\end{array}\right)=\Gamma\left(\begin{array}[]{c}m_{t}\\ x_{t}-x^{*}\end{array}\right)=\Gamma^{t}\left(\begin{array}[]{c}m_{1}\\ x_{1}-x^{*}\end{array}\right),

where $\Gamma$ is the matrix in (3.9). When the spectral radius of $\Gamma$ is less than 1, the full-batch gradient descent with momentum enjoys linear convergence. In Theorem 1, the matrix $\Gamma$ satisfies $\left\|\Gamma^{t}\right\|\leq M\lambda^{t}$ with $\lambda<1$ , which is proved in the appendix.

Remark 2

The assumption that $\Delta>0$ in (3.6) is placed for the diagonalization of the matrix $\Gamma$ , which is required in our analysis to provide a last-iterate convergence analysis and can be relaxed if we aim for a time-average convergence analysis of the sum $\sum_{t}\mathbb{E}[\left\|\widetilde{m}_{t}\right\|^{2}+\left\|x_{t}-x^{*}% \right\|^{2}]$ . More particularly, a time-average convergence analysis computes $\left\|\sum_{t=0}^{\infty}\Gamma^{t}\right\|=\left\|(I-\Gamma)^{-1}\right\|$ and requires only that the spectral radius $\lambda<1$ .

In Theorem 1, the convergence rate mainly depends on the spectral radius $\lambda$ of the matrix (3.9). We characterize its explicit form in the following theorem.

Theorem 3

For any momentum weight $\gamma\in[0,1)$ , let $\alpha,\gamma$ satisfy $\alpha L<2(1+\gamma)/(1-\gamma)$ and define $\phi=\min\left\{\alpha\mu,2(1+\gamma)/(1-\gamma)-\alpha L\right\}.$ If the momentum $\gamma<(1-\phi)^{2}/(1+\phi)^{2}$ , we have that the spectral radius $\lambda$ of $\Gamma$ defined by (3.9) satisfies

\displaystyle\lambda=\frac{\gamma+1-(1-\gamma)\phi+\sqrt{\left(\gamma+1-(1-% \gamma)\phi\right)^{2}-4\gamma}}{2}.

(3.12)

On the other hand, if the momentum $\gamma\geq(1-\phi)^{2}/(1+\phi)^{2}$ , we have

\lambda=\sqrt{\gamma}.

Theorem 3 reveals that the behavior of the convergence rate $\lambda$ is essentially different in two ranges of momentum weights $\gamma$ , exhibiting a phase transition at $\gamma=(1-\phi)^{2}/(1+\phi)^{2}$ . With Theorem 3, the following remark sheds light on the optimal choice of $\alpha$ to achieve the fastest convergence rate.

Remark 4

When $\gamma$ increases from $0$ to $1$ , the spectral radius $\lambda$ first decreases and then increases, and the minimal spectral radius is achieved under the condition $\gamma=(1-\phi)^{2}/(1+\phi)^{2}$ . Particularly, for $\gamma\leq(1-\phi)^{2}/(1+\phi)^{2}$ , the quantity $\phi$ in Theorem 3 is non-decreasing, and $\lambda$ decreases as $\gamma$ increases. For $\gamma\geq(1-\phi)^{2}/(1+\phi)^{2}$ , the spectral radius $\lambda=\sqrt{\gamma}$ increases as $\gamma$ increases. Moreover, the minimal spectral radius is

\displaystyle\lambda^{*}=\frac{\sqrt{L}-\sqrt{\mu}}{\sqrt{L}+\sqrt{\mu}},

if we specify $\alpha=1/\sqrt{\mu L}$ and $\gamma=(\sqrt{L}-\sqrt{\mu})^{2}/(\sqrt{L}+\sqrt{\mu})^{2}$ , and it follows that $\phi=\sqrt{\mu/L}$ .

Figure 1 illustrates the spectral radium $\lambda$ presented in Theorem 3 with respect to different $\gamma$ and $\alpha$ when $L/\mu$ is set to $5$ . A brighter color corresponds to smaller $\lambda$ so that the convergence is faster. Figure 1 verifies that the fastest convergence rate is achieved when $\alpha$ is approximately $1/\sqrt{\mu L}$ , $\gamma$ is approximately $(\sqrt{5}-1)^{2}/(\sqrt{5}+1)^{2}\approx 0.146$ and the optimal $\lambda^{*}$ is near $\frac{\sqrt{5}-\sqrt{1}}{\sqrt{5}+\sqrt{1}}\approx 0.382$ .

Remark 5

Figure 1 shows that SGDM with large batch sizes converges faster than SGD to a local neighborhood of $x^{*}$ . This is reflected by the observation that, in the figure, the minimal radium $\lambda$ for SGDM (with best $\gamma\approx 0.146$ ) is smaller than the minimal $\lambda$ of SGD ( $\gamma=0$ ). A smaller $\lambda$ for SGDM corresponds to a smaller second term in (3.10) of Theorem 1, and thus implies that SGDM will converge faster to enter a local neighborhood of $x^{*}$ determined by the first term in (3.10). This indicates a faster convergence rate of SGDM in the initial stage with mini-batching.

Moreover, compared to SGD ( $\gamma=0$ ), SGDM with a larger $\gamma$ permits more flexible choices of the learning rate $\alpha$ , i.e., the colored region in the figure is larger for larger $\gamma$ . This means that the convergence of SGDM is less sensitive to learning rates. The conclusion is validated by numerical experiments in the subsequent section. Particularly, Figure 6 in Section 5.2 shows that SGDM permits a wider range of learning rates to achieve convergence.

Refer to caption — Figure 1: The spectral radius of $\Gamma$ with respect to $\gamma$ and $\alpha$ in Theorem 3, where $L/\mu=5/1$ .

3.1.2 Explicit convergence rates for SGDM with specified $\gamma$ and $\alpha$

In the following, we formalize the conclusions drawn from Theorem 1 and Remark 4 in two corollaries on the explicit convergence rates of SGDM within the two phases of specifications of $\gamma$ .

Corollary 6 (Small momentum weight $\gamma$ )

Under the assumptions in Theorem 1, for $0\leq\gamma\leq(1-\alpha L)/(4+4\alpha^{2}L^{2})$ , we have that $\Delta\geq(1-\alpha L)^{2}/2$ , $\alpha L^{2}L_{f}^{2}=\mathcal{O}(B\mu)$ and

\displaystyle\mathbb{E}[\left\|\widetilde{m}_{t+1}\right\|^{2}+\left\|{x}_{t+1% }-x^{*}\right\|^{2}]=\mathcal{O}\left(\frac{L^{2}}{B\mu}\alpha\sigma^{2}+L^{2}% \lambda^{2(1-\delta)t}\right).

Corollary 6 establishes the explicit finite-sample rate of convergence of SGDM under small $\gamma$ , and shows that, the spectral radius $\lambda$ decreases slowly as $\gamma$ increases from $0$ to $(1-\alpha L)/(4+4\alpha^{2}L^{2})$ .

Corollary 7 (Large momentum weight $\gamma$ )

Under the assumptions in Theorem 1, for $0<1-\gamma\leq 2\alpha\mu/(1+\alpha\mu)^{2}$ , we have that $\Delta\geq 2(1-\gamma)\alpha\mu$ , $\alpha L_{f}^{2}\left(L^{2}+\frac{\alpha^{2}}{(1-\gamma)^{2}}\right)=\mathcal{% O}\left(B\mu\right),$ and

\displaystyle\mathbb{E}[\left\|\widetilde{m}_{t+1}\right\|^{2}+\left\|{x}_{t+1% }-x^{*}\right\|^{2}]=\mathcal{O}\left(\frac{1}{B\mu}\left(L^{2}+\frac{\alpha^{% 2}}{(1-\gamma)^{2}}\right)\alpha\sigma^{2}+\left(L^{2}+\frac{\alpha}{(1-\gamma% )\mu}\right)\lambda^{2(1-\delta)t}\right),

where $\lambda$ is defined in (3.12) above. Furthermore, if $\alpha<\frac{2(1+\gamma)}{(\mu+L)(1-\gamma)}$ , we have $\lambda=\sqrt{\gamma}$ .

Corollary 7 demonstrates that the spectral radius $\lambda$ increases with $\gamma$ when it approaches $1$ , exhibiting an opposite behavior compared to the small $\gamma$ settings that in Corollary 6. As we demonstrated in Remark 5 above, compared to mini-batch SGD ( $\gamma=0$ ), mini-batch SGDM with an appropriate $\gamma$ converges faster to the local neighborhood of $x^{*}$ by achieving a smaller second term in their convergence rates. The following Remark 8 implements an optimal choice of the momentum weight $\gamma$ in SGDM, to provide explicit convergence rates of $\widetilde{m}_{t+1}$ and ${x}_{t+1}$ with respect to $L$ , $\mu$ and $B$ .

Remark 8

As one increases $\gamma$ from $0$ to $1$ , the corresponding spectral radius $\lambda$ first decreases and then increases. Particularly, if one specifies $\gamma$ and $\alpha$ as

\displaystyle\gamma=\frac{(1-c\alpha\mu)^{2}}{(1+c\alpha\mu)^{2}},\quad\text{% and}\quad\alpha\leq\sqrt{\frac{1}{\mu L}},

for $c<1$ sufficiently close to 1, we have that $\Delta=\mathcal{O}(\alpha^{2}\mu^{2})$ and

\displaystyle\mathbb{E}[\left\|\widetilde{m}_{t+1}\right\|^{2}+\left\|{x}_{t+1% }-x^{*}\right\|^{2}]=\mathcal{O}\left(\frac{(L+1/\mu)^{2}}{B\mu}\alpha\sigma^{% 2}+(L+1/\mu)^{2}\gamma^{(1-\delta)t}\right).

We now compare our rate of convergence with several results in the existing literature.

Comparison to SGD.

Theorem 4.6 in Bottou et al. (2018) proved the convergence rate of SGD under strong convexity $\mathcal{O}(\frac{L}{\mu}\alpha\sigma^{2}+(1-\alpha\mu)^{t})$ under the assumption $\mathbb{E}[\left\|\nabla f_{\xi}(x)-\nabla f(x)\right\|^{2}]\leq\sigma^{2}$ for all $x$ . Comparing with their result, our analysis in Theorem 3 can apply with $\gamma=0$ and $\alpha<2/(\mu+L)$ such that $\phi=\alpha\mu$ , which is consistent with Bottou et al. (2018). Our Remark 4 also shows that SGDM with large momentum weight can accelerate the convergence rate in Bottou et al. (2018).

Comparison to existing results for SGDM.

Liu et al. (2020) proved that the convergence rate for strongly convex loss is $\mathcal{O}(\frac{L}{\mu}\alpha\sigma^{2}+\max\{1-\alpha\mu,\gamma\}^{t})$ for $t\geq t_{0}=\mathcal{O}(-1/\log(\gamma))$ . Comparing that with Corollary 6 above, our result improves the rate when the momentum weight $\gamma$ is small. Moreover, Liu et al. (2020) required learning rate $\alpha=\mathcal{O}(1-\gamma)$ , which is considerably more restrictive than the range permitted in our Theorems that $\alpha L<2(1+\gamma)/(1-\gamma)$ . This requirement is in opposition to our finding that SGDM permits a wider range of learning rates. The stochastic heavy ball (SHB) method is equivalent to SGDM for $\alpha^{\prime}=\alpha(1-\gamma)$ and $\gamma^{\prime}=\gamma$ . Under the noiseless setting, for $\alpha^{\prime}=1/L$ and $\gamma^{\prime}=(\sqrt{L/\mu}-1)/(\sqrt{L/\mu}+1)$ , Nesterov (1983) gave the $\mathcal{O}((1-\sqrt{\mu/L})^{t})$ convergence rate. For $\alpha^{\prime}=\mathcal{O}(1/L)$ and $\alpha=\mathcal{O}(\sqrt{1/\mu L})$ in Remark 8, our convergence factor is $\gamma=\mathcal{O}(1-\sqrt{\mu/L})$ , which improves the factor $\mathcal{O}(1-\mu/L)$ in SGD and is consistent with Nesterov (1983). Theorem 3.5 in Bollapragada et al. (2022) proved that for $\sqrt{\alpha^{\prime}}=\frac{2}{\sqrt{L+\epsilon}+\sqrt{\mu-\epsilon}}$ and $\sqrt{\gamma^{\prime}}=\frac{\sqrt{(L+\epsilon)/(\mu-\epsilon)}-1}{\sqrt{(L+% \epsilon)/(\mu-\epsilon)}+1}$ , $\epsilon\in(0,\mu)$ , the convergence rate of the SHB method on a consistent linear system $Ax=b$ is $\mathcal{O}(\frac{L^{2}}{\epsilon}t^{2}(\gamma^{\prime})^{t}+\frac{L^{4}}{% \epsilon\mu^{2}}\sigma^{2})$ . Their convergence rate is the same as ours in Remark 8 under a different setting, while their results are more restrictive on the choice of learning rates and momentum weights.

3.2 The finite-sample rates for SGDM on general losses

For the general strongly convex loss where $\overline{L}>0$ , we provide the following convergence rate.

Theorem 9

Under (A1), (A2), and (A3’) and $\overline{L}>0$ , for any momentum $\gamma\in[0,1)$ , assume the learning rate $\alpha>0$ satisfies $\alpha L<2(1+\gamma)/(1-\gamma)$ and

\displaystyle 6\sqrt{2}c(2\log T+4d)M\left(\frac{2L_{f}}{\lambda^{(1-\delta)}% \sqrt{\delta}(1-\lambda)^{1/2}}+\frac{3\overline{L}M\alpha\sigma}{(1-\lambda)^% {3/2}}\right)\alpha\leq\sqrt{B},

for fixed $\delta\in(0,1]$ and an absolute constant $c>0$ , where $M$ is defined in (3.6). In addition, the initialization $x_{1}$ satisfies $\left\|x_{1}-x^{*}\right\|\leq\frac{\delta(1-\lambda)\lambda^{1-\delta}}{36M^{% 2}\alpha\overline{L}}$ . Then with probability $1-2T^{-1}$ ,

\displaystyle\sqrt{\left\|\widetilde{m}_{t+1}\right\|^{2}+\left\|x_{t+1}-x^{*}% \right\|^{2}}\leq\frac{3\sqrt{2}c\left(2\log T+4d\right)M}{\sqrt{B(1-\lambda)}% }\alpha\sigma+3M\left\|x_{1}-x^{*}\right\|\lambda^{(1-\delta)t},~{}\text{for }% 1\leq t\leq T.

Compared to Theorem 1, the convergence rate in Theorem 9 is worse but with only up to a logarithm term, which indicates that SGDM for general losses achieves almost the same rate as that for quadratic ones. This convergence bound is presented with high probability for general losses, in contrast to the convergence in expectation outlined for the quadratic losses previously. Nonetheless, it remains uniformly applicable across all iterations $T$ .

Mini-batch SGDM for general convex loss still converges faster than mini-batch SGD with high probability, as discussed in Theorem 3. In contrast to the quadratic settings $\overline{L}=0$ , the general convex setting requires additionally that the initialization $x_{1}$ is close to $x^{*}$ , i.e., $\left\|x_{1}-x^{*}\right\|$ is small. This assumption is expected since the loss function is not guaranteed to be convex everywhere under the weak assumptions in (A1)–(A2). Particularly, for $\gamma$ close to 0, we assume $\left\|x_{1}-x^{*}\right\|=\mathcal{O}\left(\mu/\overline{L}L^{2}\right)$ . For $\gamma$ close to $1$ , we assume $\left\|x_{1}-x^{*}\right\|=\mathcal{O}\left(\mu/\overline{L}\left(L^{2}+\alpha% ^{2}/(1-\gamma)^{2}\right)\right)$ .

Based on Theorem 9, we have the following corollaries of the convergence rates for small and large specifications of $\gamma$ under general losses, similar to Corollaries 6–7 for quadratic losses.

Corollary 10 (Small momentum weight $\gamma$ )

Following Theorem 9, for $0\leq\gamma\leq(1-\alpha L)/(4+4\alpha^{2}L^{2})$ , we have that $\Delta\geq(1-\alpha L)^{2}/2$ , $\sqrt{\alpha}L_{f}L\left(\mu+\overline{L}L\right)=\mathcal{O}\left(B^{1/2}\mu^% {3/2}\right)$ and with high probability,

\displaystyle\sqrt{\left\|\widetilde{m}_{t+1}\right\|^{2}+\left\|x_{t+1}-x^{*}% \right\|^{2}}=\mathcal{O}\left(\frac{L\log T}{\sqrt{B\mu}}\sqrt{\alpha}\sigma+% L\lambda^{(1-\delta)t}\right),

where $\lambda$ is defined in (3.12) above.

Corollary 11 (Large momentum weight $\gamma$ )

Following Theorem 9, for $0<1-\gamma\leq 2\alpha\mu/(1+\alpha\mu)^{2}$ , we have that $\Delta\geq 2(1-\gamma)\alpha\mu$ ,

\sqrt{\alpha}L_{f}\left(L+\frac{\alpha}{1-\gamma}\right)\left\{1+\overline{L}% \left(L+\frac{\alpha}{1-\gamma}\right)\sqrt{\frac{\alpha}{(1-\gamma)\mu}}% \right\}=\mathcal{O}\left(B\mu\right),

and with high probability

\displaystyle\sqrt{\left\|\widetilde{m}_{t+1}\right\|^{2}+\left\|x_{t+1}-x^{*}% \right\|^{2}}=\mathcal{O}\left(\frac{\log T}{\sqrt{B\mu}}\left(L+\frac{\alpha}% {1-\gamma}\right)\sqrt{\alpha}\sigma+\left(L+\frac{\sqrt{\alpha}}{\sqrt{(1-% \gamma)\mu}}\right)\lambda^{(1-\delta)t}\right),

where $\lambda$ is defined in (3.12) above. Furthermore, if $\alpha<\frac{2(1+\gamma)}{(\mu+L)(1-\gamma)}$ , we have $\lambda=\sqrt{\gamma}$ .

In the following remark, we choose the optimal momentum weights $\gamma$ to show the explicit convergence results with respect to $L$ , $\mu$ , and $B$ .

Remark 12

If one specifies $\gamma$ and $\alpha$ as

\displaystyle\gamma=\frac{(1-c\alpha\mu)^{2}}{(1+c\alpha\mu)^{2}},\quad\text{% and}\quad\alpha\leq\sqrt{\frac{1}{\mu L}},

for $c<1$ sufficiently close to 1, we have that $\Delta=\mathcal{O}(\alpha^{2}\mu^{2})$ and with high probability

\displaystyle\left\|\widetilde{m}_{t+1}\right\|^{2}+\left\|x_{t+1}-x^{*}\right% \|^{2}=\mathcal{O}\left(\frac{(L+1/\mu)^{2}(\log T)^{2}}{B\mu}\alpha\sigma^{2}% +(L+1/\mu)^{2}\gamma^{(1-\delta)t}\right).

We arrive at the same conclusions as we had in quadratic settings that SGDM can accelerate convergence over SGD under general strongly convex losses when $\bar{L}>0$ .

4 Acceleration by Averaging and Asymptotic Normality

In this section, we study the Polyak-averaging SGDM (referred to as averaged SGDM). Particularly, we aim at building the convergence result of $\frac{1}{n-n_{0}}\sum_{t=n_{0}+1}^{n}x_{t}$ , which is an average of all iterates $x_{t}$ starting from a constant period $n_{0}$ . We show that the averaging leads to acceleration for SGDM, compared to the last-iterate bounds established in the previous section.

4.1 Averaged SGDM under quadratic losses

For the mini-batch model (2.4) under quadratic losses $(\overline{L}=0)$ , we now characterize a decomposition of the error of the averaged SGDM $\frac{1}{n-n_{0}}\sum_{t=n_{0}+1}^{n}(x_{t}-x^{*})$ , which leads to an acceleration over the last iterate of SGDM $x_{n}$ .

Theorem 13

Under (A1)-(A3) and $\overline{L}=0$ , suppose that the conditions in Theorem 1 hold, and the $t$ -th iteration $(\widetilde{m}_{t},x_{t})$ satisfies

\displaystyle\mathbb{E}[\left\|\widetilde{m}_{t}\right\|^{2}+\left\|x_{t}-x^{*% }\right\|^{2}]\leq\frac{C_{1}}{B(1-\lambda)}\alpha^{2}\sigma^{2}+C_{2}\left\|x% _{1}-x^{*}\right\|^{2}\lambda^{2(1-\delta)(t-1)},\quad t\geq 1.

(4.13)

Then for $n\geq 2n_{0}$ where $\lambda^{2n_{0}}=B^{-1}(1-\lambda)$ , we have

\displaystyle\frac{\sum_{t=n_{0}+1}^{n}(x_{t}-x^{*})}{n-n_{0}}=-\frac{\sum_{t=% n_{0}+1}^{n}\sum_{i=1}^{B}\Sigma^{-1}(A_{\xi_{ti}}x^{*}-b_{\xi_{ti}})}{B(n-n_{% 0})}+R_{n},

(4.14)

where $\mathbb{E}\left\|R_{n}\right\|^{2}\leq\frac{\tilde{C}_{1}}{Bn^{2}}+\frac{% \tilde{C}_{2}}{B^{2}n}$ , and

	$\displaystyle\tilde{C}_{1}=40\alpha^{2}L_{f}^{2}M^{2}\frac{C_{2}\left\\|x_{1}-x% ^{}\right\\|^{2}}{(1-\delta)(1-\lambda)^{3}}+120\alpha^{2}\sigma^{2}M^{2}\frac% {1}{(1-\lambda)^{3}}+\frac{4M^{2}\left\\|x_{1}-x^{}\right\\|^{2}}{1-\lambda},$
	$\displaystyle\tilde{C}_{2}=30\alpha^{2}\sigma^{2}L_{f}^{2}M^{2}\frac{C_{1}% \alpha^{2}}{(1-\lambda)^{3}}.$

and $M$ is defined in (3.6).

In the above theorem, we obtain the convergence rate of averaged SGDM in (4.14). While (4.13) is assumed for the purpose of explicitly establishing the dependence of $\tilde{C}_{1}$ , $\tilde{C}_{2}$ to the constants $C_{1}$ , $C_{2}$ in (4.13), it is nothing but a rephrasing of the last-iterate bounds established in (3.10). Comparing the convergence rates in Theorem 13, the averaged SGDM converges to $x^{*}$ without a bias as $n$ increases. The leading term $-\frac{\sum_{t=n_{0}+1}^{n}\Sigma^{-1}(A_{\xi_{ti}}x^{*}-b_{\xi_{ti}})}{n-n_{0}}$ is an average of $(n-n_{0})$ i.i.d. vectors when $\xi_{ti}$ are i.i.d. sampled. Therefore the leading term is bounded by $\mathcal{O}(\frac{1}{Bn})$ in squared expectation, and the remainder term $R_{n}$ in (4.14) is bounded by $\frac{\tilde{C}_{1}}{Bn^{2}}+\frac{\tilde{C}_{2}}{B^{2}n}$ . We demonstrate that the application of the averaging technique enhances the convergence rate of SGDM from the biased expression presented in (4.13) to an unbiased rate of $\mathcal{O}(\frac{1}{n})$ , which asymptotically approaches zero as $n\rightarrow\infty$ . It is noteworthy to mention that, for averaged SGDM, the initialization error $\|x_{1}-x^{*}\|$ is forgotten at the rate of $n^{-2}$ , slower than the exponential initialization-forgetting for the last-iterate convergence established in Theorem 1.

Remark 14

Comparing the convergence rates of SGDM and averaged SGDM in Theorem 1 and Theorem 13, respectively, we establish that the leading term in the averaged SGDM does not depend on the learning rate $\alpha$ . On the other hand, the averaging process starts from $n_{0}$ , effectively excluding the initial iterations that may have larger deviations. Practically, the starting point $n_{0}$ is manually selected due to $\mu$ and $L$ . With a small $\lambda$ , averaged SGDM affords the selection of a smaller $n_{0}$ to fulfill the condition $\lambda^{2n_{0}}=B^{-1}(1-\lambda)$ , thereby demonstrating reduced sensitivity to $n_{0}$ . In addition, the leading term in (4.14) is independent of momentum weight $\gamma$ and learning rate $\alpha$ , which indeed indicates that the averaged SGDM has the same rate of convergence as the averaged SGD.

The following Corollary 15 establishes the asymptotic distribution of averaged SGDM, which is indeed the distribution of the leading term in (4.14), and the remainder term $R_{n}$ in (4.14) establishes the convergence of the averaged SGDM algorithm to the asymptotic normal distribution.

Corollary 15

Following Theorem 13, we have

\displaystyle\frac{\sqrt{B}}{\sigma}\frac{\sum_{t=n_{0}+1}^{n}({x}_{t}-x^{*})}% {\sqrt{n-n_{0}}}\Rightarrow\mathcal{N}(0,\Sigma^{-1}\Omega\Sigma^{-1}),\quad% \text{as}~{}\left(\alpha n,\frac{B}{\alpha}\right)\to\infty,

where

\displaystyle\Omega=\frac{1}{\sigma^{2}}\mathbb{E}[(A_{\xi}x^{*}-b_{\xi})(A_{% \xi}x^{*}-b_{\xi})^{\top}].

(4.15)

The asymptotic distribution of the averaged SGDM is the same as that of the averaged SGD in the existing literature (Polyak and Juditsky, 1992), with different remainder terms $R_{n}$ . As the averaged SGDM is close to $x^{*}$ , the limiting distribution is established with covariance matrix $\Sigma^{-1}\Omega\Sigma^{-1}$ where $\Sigma$ and $\Omega$ are respectively the Hessian and Gram matrix of $A_{\xi}x^{*}-b_{\xi}$ .

Remark 16

Asymptotic normality in Corollary 15 hold under asymptotics $\alpha n,B/\alpha\rightarrow\infty$ . We specify three common scenarios of the learning rate $\alpha$ , batch size $B$ , and sample size $n$ .

•

For a constant learning rate $\alpha$ , Corollary 15 holds as $n$ and $B$ tend to infinity. As the batch size $B$ increases, the bias term in (4.13) approaches zero, reducing the difference between the stochastic gradient at $x_{t}$ and $x^{*}$ .
•

For decaying learning rate $\alpha$ , the asymptotic normality holds with a fixed batch size $B$ as long as $\alpha n$ diverges. The idea is intuitive: as $\alpha$ decreases, the random effect reduces. For instance, Corollary 15 holds when the batch size $B$ is fixed and $\alpha=\Theta(n^{-\epsilon})$ with $\epsilon\in(0,1)$ as $n\to\infty$ .
•

By minimizing the remainder term $R_{n}$ in (4.13), the optimal learning rate is $\alpha=\Theta((n/B)^{-1/2})$ . With such a specified learning rate, the averaged SGDM correspondingly converges to asymptotic normality with rate $\mathcal{O}((nB)^{-1/2})$ .

In addition, Corollary 15 provides a rigorous foundation for constructing asymptotically valid confidence intervals based on the asymptotic normality of averaged SGDM. By estimating the covariance matrix, uncertainty quantification and statistical inference can be performed based on SGDM methods. Based on the asymptotic normality and the given covariance matrix, the $95\%$ confidence region of $x^{*}$ can be constructed as

\displaystyle I_{n}=\left\{x\in\mathbb{R}^{d}:\frac{B(n-n_{0})}{\sigma^{2}}% \left(\frac{\sum_{t=n_{0}+1}^{n}x_{t}}{n-n_{0}}-x\right)^{\top}\Sigma\Omega^{-% 1}\Sigma\left(\frac{\sum_{t=n_{0}+1}^{n}x_{t}}{n-n_{0}}-x\right)\leq\chi^{2}_{% d,0.05}\right\},

where $\chi^{2}_{d,0.05}$ is the $0.95$ -quantile of the chi-squared distribution with $d$ degrees of freedom. Notably, for $\omega\in\mathbb{R}^{d}$ with $\left\|\omega\right\|=1$ as a test vector, we can define the random variable

\displaystyle Z=\frac{\sqrt{B}}{\sigma\sqrt{\omega^{\top}\Sigma^{-1}\Omega% \Sigma^{-1}\omega}}\frac{\sum_{t=n_{0}+1}^{n}\omega^{\top}(x_{t}-x^{*})}{\sqrt% {n-n_{0}}}\Rightarrow\mathcal{N}(0,1),\quad\text{as}\quad n\to\infty,

(4.16)

based on which, one can construct a one-dimensional asymptotic exact confidence interval for $\omega^{\top}x^{*}$ ,

I_{n}^{\omega}=\left[\frac{\omega^{\top}\sum_{t=n_{0}+1}^{n}x_{t}}{n-n_{0}}-z_% {0.025}\frac{\sigma\sqrt{\omega^{\top}\Sigma^{-1}\Omega\Sigma^{-1}\omega}}{% \sqrt{B(n-n_{0})}},\frac{\omega^{\top}\sum_{t=n_{0}+1}^{n}x_{t}}{n-n_{0}}+z_{0% .025}\frac{\sigma\sqrt{\omega^{\top}\Sigma^{-1}\Omega\Sigma^{-1}\omega}}{\sqrt% {B(n-n_{0})}}\right],

where $z_{0.025}$ is the $0.975$ -quantile of the standard normal distribution. Particularly, $\Pr(\omega^{\top}x^{*}\in I_{n}^{\omega})\rightarrow 0.95$ , as $n\rightarrow\infty$ . In Section 5.2 below, we conduct a simulation experiment on constructing the confidence intervals and report the outcomes in Figure 9. The coverage of the proposed construction is persuasive for averaged SGDM, and its performance benefits from the less sensitivity to the learning rates. The confidence interval can also serve as an uncertainty quantification of the averaged SGDM estimator $\frac{\omega^{\top}\sum_{t=n_{0}+1}^{n}x_{t}}{n-n_{0}}$ , and one may consider the length of the confidence interval as a criterion for stopping the algorithm or other decision-making purposes.

4.2 Averaged SGDM under general losses

For the general strongly convex loss function, we provide the convergence rate of the averaged SGDM in the following theorem.

Theorem 17

For $\overline{L}>0$ , suppose that the conditions in Theorem 9 hold, and the $t$ -th iteration $(\widetilde{m}_{t},x_{t})$ satisfies for fixed $\delta\in(0,\frac{1}{2}]$ ,

\displaystyle\sqrt{\left\|\widetilde{m}_{t}\right\|^{2}+\left\|{x}_{t}-x^{*}% \right\|^{2}}\leq\frac{C_{1}}{\sqrt{B(1-\lambda)}}\alpha\sigma+C_{2}\left\|x_{% 1}-x^{*}\right\|\lambda^{(1-\delta)(t-1)},\quad 1\leq t\leq T.

(4.17)

Then for $2n_{0}\leq n\leq T$ , where $\lambda^{n_{0}}=B^{-1/2}(1-\lambda)^{1/2}$ , with probability of at least $1-4T^{-1}$ , we have

\displaystyle\frac{\sum_{t=n_{0}+1}^{n}(x_{t}-x^{*})}{n-n_{0}}=-\frac{\sum_{t=% n_{0}+1}^{n}\sum_{i=1}^{B}\Sigma^{-1}\triangledown f_{\xi_{ti}}(x^{*})}{B(n-n_% {0})}+R_{n},

(4.18)

where

\displaystyle\left\|R_{n}\right\|\leq\frac{\tilde{C}_{1}+\tilde{C}_{3}}{\sqrt{% B}n}+\frac{\tilde{C}_{2}}{B\sqrt{n}}+\frac{\tilde{C}_{4}}{B},

(4.19)

	$\displaystyle\tilde{C}_{1}=\alpha L_{f}M\frac{16c(2\log T+4d)C_{2}\left\\|x_{1}% -x^{}\right\\|}{(1-\lambda)^{3/2}}+\alpha\sigma M\frac{4\sqrt{3}c(2\log T+4d)}% {(1-\lambda)^{3/2}}+\frac{2M\left\\|x_{1}-x^{}\right\\|}{(1-\lambda)^{1/2}},$
	$\displaystyle\tilde{C}_{2}=\alpha L_{f}\sigma M\frac{8\sqrt{2}c(2\log T+4d)C_{% 1}\alpha}{(1-\lambda)^{3/2}},~{}\tilde{C}_{3}=8\alpha M\overline{L}\frac{C_{2}% ^{2}\left\\|x_{1}-x^{*}\right\\|^{2}\left(n_{0}(1-\lambda)+1\right)}{(1-\lambda)% ^{3/2}},~{}\tilde{C}_{4}=4\alpha\sigma^{2}M\overline{L}\frac{C_{1}^{2}\alpha^{% 2}}{(1-\lambda)^{2}}.$

Here $M$ is defined in (3.6).

For non-quadratic settings of $\overline{L}>0$ , there are two additional terms $\tilde{C}_{3}/(\sqrt{B}n)$ and $\tilde{C}_{4}/B$ in (4.19) compared to the case of $\overline{L}=0$ . As $n\rightarrow\infty$ , the term $\tilde{C}_{4}/B$ appears as a non-vanishing bias, if the batch size $B$ and learning rate $\alpha$ both stays fixed, due to which the asymptotic normality does not hold. Nonetheless, if one decays $\alpha$ or increases $B$ appropriately as $n$ increases, the asymptotic normality result remains to hold in the following corollary.

Corollary 18

Following Theorem 17, we have

\displaystyle\frac{\sqrt{B}}{\sigma}\frac{\sum_{t=n_{0}+1}^{n}({x}_{t}-x^{*})}% {\sqrt{n-n_{0}}}\Rightarrow\mathcal{N}(0,\Sigma^{-1}\Omega\Sigma^{-1}),\quad% \text{as}~{}\left(\frac{\sqrt{\alpha n}}{\log(B/\alpha)},\frac{\sqrt{B}}{% \alpha\sqrt{n}}\right)\to\infty,

where

\displaystyle\Omega=\frac{1}{\sigma^{2}}\mathbb{E}[\triangledown f_{\xi}(x^{*}% )\triangledown f_{\xi}(x^{*})^{\top}].

(4.20)

Remark 19

The conditions on $n$ , $B$ , and $\alpha$ in this corollary are more restrictive than those in Corollary 15 due to the presence of the approximation errors $\overline{L}\left\|x-x^{*}\right\|^{2}$ . For a decaying learning rate $\alpha$ , the condition on the batch size is much relaxed since a small learning rate reduces the error caused by randomness, as shown in the $C_{1}$ term in (4.17). Particularly, when the batch size $B$ is fixed, the condition of Corollary 18 is met for $\alpha=\Theta(n^{-\epsilon})$ with $\epsilon\in(\frac{1}{2},1)$ as $n\to\infty$ . Specifically, with diverging $n$ and $B$ , the nearly optimal learning rate is $\alpha=\Theta((n/\sqrt{B})^{-2/3})$ , and the corresponding rate of convergence to asymptotic normality is $\mathcal{O}((nB)^{-1/3}\log^{2}(nB))$ .

5 Experiments

In this section, we support our theoretical results with simulations in a quadratic example and a logistic loss example of the general strongly convex loss, and real-data experiments of a multinomial logistic regression on MNIST hand-written digit classification.

In the numerical experiments, we verify the convergence results built in the paper on a training set of size $N$ . Setting $\xi$ to be a discrete uniform random variable with values in $\{1,2,\dots,N\}$ , we consider to minimize $\min_{x}\frac{1}{N}\sum_{i=1}^{N}f_{i}(x)$ , as an example of the stochastic optimization model (1.1), where the objective is either the empirical risk $\mathbb{E}[f_{\xi}(x)]=\frac{1}{N}\sum_{i=1}^{N}f_{i}(x)$ , or more specifically in a maximum likelihood estimation, the negative log-likelihood. In all the simulation results below, we repeat the experiment 200 times. For the real-data analysis on MNIST, we repeat the experiment in a replicable setting while we fix the random seeds to 1, 2, and 3 in training.

5.1 Simulation: quadratic loss

In a quadratic loss model (3.5), we first conduct a simulation study with a sample of size $N=20,000$ and dimension $d=10$ . We generate the parameters $\{(A_{i},b_{i})\}_{i=1}^{N}$ in (3.5) i.i.d and $b_{i}\sim\mathcal{N}(0_{d},I_{d})$ , and the positive definitive matrix $A_{i}=\rho V_{i}^{\top}V_{i}+10I_{d}$ . Here $V_{i}$ is a matrix in $\mathbb{R}^{d\times d}$ and each row of $V_{i}$ generates from the normal distribution $\mathcal{N}(0_{d},I_{d})$ , where $\rho>0$ affect the conditional number $L/\mu$ . We choose $\rho=1$ and the average conditional number $L/\mu$ is $35/10$ . The deterministic minimizer $x^{*}$ is computed by $x^{*}=\left(\sum_{i=1}^{N}A_{i}\right)^{-1}\left(\sum_{i=1}^{N}b_{i}\right)$ .

We consider the mini-batch SGDM with batch size $B=0.2N$ with replacement. The learning rate is fixed at $\alpha=0.001$ . According to Theorem 3, we choose the momentum weight following

\displaystyle\gamma=\left(\frac{1-\mu\alpha}{1+\mu\alpha}\right)^{2}.

(5.21)

We refer to the above $\gamma$ as the adaptive momentum weight in SGDM (SGDM-adap), while we also compare it with other fixed $\gamma$ , as well as SGD as a special case of SGDM with $\gamma=0$ in Figure 2.

Figure 2 supports the theoretical results in Theorems 1 and 3. Both SGD and SGDM enjoy linear convergence at the beginning and have the same order of error at the end. The adaptive momentum SGDM-adap with $\gamma=(1-\phi)^{2}/(1+\phi)^{2}$ in (5.21) converges fastest, where the average value over 200 experiments is $0.95$ . The acceleration of convergence with small momentum weights is not significant, while SGDM with a very large $\gamma$ becomes slower and unstable. Particularly, SGDM with $\gamma=0.5$ converges as fast as SGD, and SGDM with momentum weight 0.9 converges faster than SGD. For large momentum weight $\gamma=0.99$ , the linear convergence factor is $\sqrt{\gamma}$ , and we can see from the experiment that its convergence is the slowest in Figure 2.

We can see that appropriate momentum weights lead to faster convergence as shown in Figure 2(a). However, when the momentum weight is exceedingly large, SGDM may cause the $\ell_{2}$ error to oscillate, as shown in Figure 2(b). This is because momentum causes the algorithm to continue moving in the direction of past gradients, even if the current gradient is opposite to the momentum. As a result, SGDM oscillates back and forth near the minimum of the loss function instead of converging steadily towards it. Moreover, it also leads to a decrease in the convergence speed. The observation matches the finding in Theorem 3.

We further compare the averaged SGD and SGDM. We set $n_{0}$ in Theorem 13 as $n_{0}=200$ and $500$ , and report the convergence of the averaged SGD and SGDM in

\displaystyle\bar{x}_{n_{0}+t}=\frac{\sum_{j=n_{0}+1}^{n_{0}+t}x_{j}}{t}.

(5.22)

In Figures 3–4, averaged SGDM with adaptive momentum weight defined as (5.21) converges faster than the others. Specifically, it converges significantly faster than the averaged SGD. When $n_{0}=500$ , Figure 4 illustrates that the averaging technique reduces the order of the bias from $10^{-8}$ , as seen in Figure 2, to $10^{-9}$ .

We further conduct a simulation study to verify the asymptotic normality result of Corollary 15. In this experiment, we evaluate the one-dimensional projection statistic $Z$ in (4.16) with $n_{0}=1000$ , $n=2000$ and $\Omega=\left(N\sigma^{2}\right)^{-1}\sum_{i=1}^{N}(A_{i}x^{*}-b_{i})(A_{i}x^{*% }-b_{i})^{\top}$ . We generate the trajectory $\{x_{1},\cdots,x_{n}\}$ 1000 times and get the replications $Z^{(1)},Z^{(2)},\cdots,Z^{(1000)}$ . Figure 5 shows the frequency of $\{Z^{(1)},Z^{(2)},\cdots,Z^{(1000)}\}$ . From Figure 5 we can see that both averaged SGD and SGDM well approximate the normal distribution. We can further see that the frequency of averaged SGDM is slightly closer to the normal distribution than that of averaged SGD under finite rounds $n=2000$ . These observations reflect the theoretical results we build in Theorem 13.

5.2 Sensitivity to learning rates

In this section, we will show that the performance of SGDM and averaged SGDM have a wider range of tunable learning rates compared to SGD and averaged SGD. In the quadratic loss model with $N=20,000$ and $d=10$ , we generate the parameters $\{(A_{i},b_{i})\}_{i=1}^{N}$ in (3.5) i.i.d and $b_{i}\sim\mathcal{N}(0_{d},I_{d})$ , and the positive definitive matrix $A_{i}=V_{i}^{\top}V_{i}+I_{d}$ , where each row of $V_{i}$ generates from the normal distribution $\mathcal{N}(0_{d},I_{d})$ . The average conditional number $L/\mu$ is $26/1$ . The learning rates we set are chosen from $\{2^{1},2^{0},2^{-1},2^{-2},\cdots\}$ . The batch size is $B=0.2N$ .

Figure 6 shows the convergence behaviors of SGD and SGDM for a wide range of learning rates $\alpha=2^{-1},2^{-3},2^{-5},2^{-7},2^{-9}$ . Notably, SGD fails to converge with learning rates of $\alpha\geq 2^{-3}$ . In contrast, SGDM with $\gamma=0.8$ exhibits robust convergence properties, successfully converging even with a larger learning rate $\alpha=2^{-1}$ . This comparison emphasizes that SGDM is less sensitive to the learning rate. Figure 7 further illustrates the finite-sample errors of SGD and SGDM across different learning rates, given a fixed number of iterations at $T=500$ . Intuitively, smaller learning rates necessitate more iterations to converge and are associated with an increased bias for each method. Notably, SGDM is capable of converging with a relatively large learning rate. The maximum learning rate ensuring convergence is $\alpha=2^{-4}$ for SGD, $\alpha=2^{-1}$ for SGDM with $\gamma=0.8$ and $\alpha=2^{0}$ for SGDM with $\gamma=0.9$ . Reflecting upon the learning rate condition $\alpha L<(1+\gamma)/(1-\gamma)$ , it can be confirmed that $(1+\gamma)/(1-\gamma)\approx 2^{-1}/2^{-4}$ for $\gamma=0.8$ and $(1+\gamma)/(1-\gamma)\approx 2^{0}/2^{-4}$ for $\gamma=0.9$ . This alignment validates that the figure is in agreement with the theoretical analysis.

Moreover, we examine the performance of averaged SGD and averaged SGDM across a range of learning rates. For both algorithms, we set $n_{0}=500$ and define $\bar{x}_{n_{0}+t}$ as (5.22). Figure 8 illustrates that all algorithms exhibit a sublinear rate of convergence. SGDM with learning rates ranging from $\alpha=2^{-1}$ to $2^{-7}$ , as well as SGD with learning rates $\alpha=2^{-5},2^{-7}$ , achieve convergence to an similar error level. Compared to Figure 6, we can see that the averaging technique improves convergence and is less dependent on learning rates. However, it is noteworthy that algorithms with learning rate $\alpha=2^{-9}$ require a substantially larger initial iteration count $n_{0}$ to mitigate the effects of iterations with significant deviations.

Additionally, we investigate the asymptotic normality as outlined in Corollary 15. Figure 9 illustrates the asymptotic behavior of the statistic $Z$ in (4.16) with $n=1000$ and $n_{0}=500$ . This illustration confirms its convergence to $\mathcal{N}(0,1)$ , under a range of learning rates. The $y$ -axis displays the empirical probability $P(|Z|<1.96)$ , which indicates the frequency with which the statistic $Z$ falls within the critical range $(-1.96,1.96)$ across $1000$ trials, corresponding to a $95\%$ confidence interval. The $x$ -axis specifies the learning rates employed in the numerical experiments. Figure 9 supports the assertions in the corollary, demonstrating the robustness of the asymptotic normality with respect to different learning rates. The averaged SGDM with large momentum permits a broader selection of learning rates and exhibits reduced sensitivity to their variation.

5.3 Simulation: logistic regression

In Example 1, we generate the data $a_{i}\in\mathbb{R}^{d}$ i.i.d from $\mathcal{N}(0_{d},I_{d})$ , and $b_{i}\in\{0,1\}$ where generated by $b_{i}=1$ with probability $p_{x}(a_{i})$ and $b_{i}=0$ otherwise. Here $p_{x}(a)=1/(1+\exp(-x^{\top}a))$ , and we use $x=\frac{1}{\sqrt{d}}(1,1,\cdots,1)^{\top}$ . We fix the sample size $N=20,000$ and the dimension $d=10$ . We set the regularization parameter $\nu$ to zero. Before each simulation, we run the full-batch gradient descent to get the minimizer $x^{*}$ , then compute the Hessian matrix at $x^{*}$ ,

\displaystyle\Sigma=\frac{1}{N}\sum_{i=1}^{N}p_{x^{*}}(a_{i})(1-p_{x^{*}}(a_{i% }))a_{i}a_{i}^{\top}.

We specify the batch size $B=0.2N$ and the learning rate is $\alpha=0.5$ . We consider SGDM with fixed momentum weights $\gamma=0.3,0.5,0.7,0.8,0.9$ , as well as the adaptive momentum weight as (5.21). Under the data generation procedure, the average of the adaptive momentum weight is $0.75$ .

Figure 10 illustrates that for SGDM with momentum weights $\gamma$ less than 0.5, the convergence rate of SGDM is slightly faster than SGD, while for $\gamma=0.7$ and adaptive weights $\gamma$ specified as in (5.21), the convergence is much faster than SGD. For $\gamma=0.8$ , the convergence is similar to that of the adaptive weight but less stable. For $\gamma=0.9$ , the convergence is much slower and more unstable.

To further compare the convergence of averaged SGD and SGDM, we set $n_{0}$ in Theorem 13 as $n_{0}=10$ and $40$ . We plot the comparison of the $\ell_{2}$ error of the averaged SGDM and SGD. Figure 11(a) shows that for momentum weights smaller and equal to the adaptive momentum weight in (5.21), the averaged SGDM converges faster than the averaged SGD. In Figure 11(b), we can see the averaged SGDM with $\gamma=0.8$ converges even faster than the averaged SGDM with adaptive $\gamma$ , while for $\gamma=0.9$ , the averaged SGDM converges faster at the beginning but fluctuates, which cause the increase of $\ell_{2}$ error for $\gamma=0.9$ . In comparison with Figure 2(b) for the quadratic loss, we can see that although large momentum weight causes fluctuation in both cases, the convergence is more flattened for the the logistic loss in Figure 10(b). This may also explain why the averaged SGDM performs well with $\gamma$ slightly larger than the adaptive momentum weight in Figure 11(b). Figure 12 demonstrates that with $n_{0}=40$ , there is an improvement in the bias order from $10^{-3}$ , as seen in Figure 10, to $10^{-4}$ . This improvement reflects a sublinear convergence that is in agreement with the behavior observed in Figure 4 for quadratic loss. Consequently, for averaged SGDM, this suggests reduced sensitivity to the choice of a smaller $n_{0}$ .

5.4 Real data: MNIST classification

In this section, we consider the multinomial logistic regression on the MNIST dataset, which consists of $N=60,000$ images of handwritten digits with size $28\times 28$ . We reshape the images to vectors of size $784\times 1$ . For the samples $(a_{i},b_{i})$ , where $a_{i}\in\mathbb{R}^{784}$ are the vectorized images and $b_{i}\in\mathbb{R}^{10}$ are one-hot indicators corresponding the digits $0,1,\cdots,9$ , the loss function is $f_{i}(X)=H(X^{\top}a_{i},b_{i}),$ where $X\in\mathbb{R}^{784\times 10}$ and $H$ is the multi-class cross entropy. We specify the batch size $B=256$ and the learning rate $\alpha=1.0$ . Due to the difficulty of computing the condition number $L/\mu$ , we fix the momentum weight $\gamma=0.1,0.3,0.5,0.7,0.9,0.99$ in training.

Figure 13 shows the convergence of the training loss over iterations $t$ . For the purpose of clear representation, the reported training loss for each iteration $t$ is based on an average of mini-batch stochastic losses $g_{\eta_{t}}$ of the past $\lfloor N/B\rfloor$ batches, with the batch size $B$ .

From Figure 13, we see that SGDM with momentum weight $\gamma$ ranging from $0.1$ to $0.9$ greatly outperforms SGD. SGDM with $\gamma=0.5$ converges fastest in the earlier iterations, while the training losses are almost identical for different momentum weights in the later iterations, except for $\gamma=0.99$ . The convergence rate is not very sensitive to the momentum weights and a wide range of $\gamma\in[0.3,0.9]$ leads to similar performance.

6 Conclusion

Our study sheds light on the performance of mini-batch SGDM in solving optimization problems under strongly convex loss functions. The convergence rate of mini-batch SGDM is influenced by several factors, including the batch size, the momentum weight, and the learning rate. Our analysis rigorously shows that certain choices of momentum weight with a reasonably large batch size can lead to faster convergence compared to SGD. This finding is consistent with previous numerical studies on SGDM. Additionally, our findings, supported by theoretical analysis and numerical experiments, indicate that SGDM permits a broader selection of learning rates.

We further establish the asymptotic normality of averaged SGDM in the quadratic settings and reveal the non-vanishing bias of that in the general settings. Our investigation reveals that averaged SGDM is asymptotically equivalent to averaged SGD. By minimizing the remainder term, we give the optimal learning rate and the corresponding rate of convergence to asymptotic normality. In addition, we present the asymptotic covariance matrix for the averaged SGDM, enabling the uncertainty quantification of the algorithm outputs and statistical inference of the true model parameters based on SGDM, as opposed to SGD.

In summary, our study contributes to the theoretical understanding of mini-batch SGDM and has practical implications for developing efficient optimization algorithms in machine learning. It is noteworthy to mention that our study assumed that the loss function was smooth and strongly convex. In future work, our work can be extended in several ways. It would be interesting to investigate the performance of mini-batch SGDM under non-convex loss functions, since this would have important implications for the application of SGDM in deep learning. It is also of interest to extend the analysis to the case when the learning rate decays over time.

References

Bollapragada et al. (2022) Bollapragada, R., T. Chen, and R. Ward (2022). On the fast convergence of minibatch heavy ball momentum. arXiv preprint arXiv:2206.07553.
Bottou et al. (2018) Bottou, L., F. E. Curtis, and J. Nocedal (2018). Optimization methods for large-scale machine learning. SIAM Review 60(2), 223–311.
Chen et al. (2023) Chen, X., Z. Lai, H. Li, and Y. Zhang (2023). Online statistical inference for stochastic optimization via Kiefer-Wolfowitz methods. Journal of the American Statistical Association (To appear).
Chen et al. (2020) Chen, X., J. D. Lee, X. T. Tong, and Y. Zhang (2020). Statistical inference for model parameters in stochastic gradient descent. The Annals of Statistics 48(1), 251 – 273.
Defazio (2020) Defazio, A. (2020). Understanding the role of momentum in non-convex optimization: Practical insights from a lyapunov analysis. arXiv preprint arXiv:2010.00406.
Gitman et al. (2019) Gitman, I., H. Lang, P. Zhang, and L. Xiao (2019). Understanding the role of momentum in stochastic gradient methods. Advances in Neural Information Processing Systems 32.
Jin et al. (2022) Jin, R., Y. Xing, and X. He (2022). On the convergence of mSGD and AdaGrad for stochastic optimization. arXiv preprint arXiv:2201.11204.
Kidambi et al. (2018) Kidambi, R., P. Netrapalli, P. Jain, and S. Kakade (2018). On the insufficiency of existing momentum schemes for stochastic optimization. In Information Theory and Applications Workshop, pp. 1–9.
Kingma and Ba (2014) Kingma, D. P. and J. Ba (2014). Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.
Lee et al. (2022) Lee, K., A. Cheng, E. Paquette, and C. Paquette (2022). Trajectory of mini-batch momentum: Batch size saturation and convergence in high dimensions. Advances in Neural Information Processing Systems 35, 36944–36957.
Lee et al. (2022) Lee, S., Y. Liao, M. H. Seo, and Y. Shin (2022). Fast and robust online inference with stochastic gradient descent via random scaling. In Proceedings of the AAAI Conference on Artificial Intelligence, Volume 36, pp. 7381–7389.
Li et al. (2022) Li, X., M. Liu, and F. Orabona (2022). On the last iterate convergence of momentum methods. In International Conference on Algorithmic Learning Theory, pp. 699–717.
Liu et al. (2020) Liu, Y., Y. Gao, and W. Yin (2020). An improved analysis of stochastic gradient descent with momentum. Advances in Neural Information Processing Systems 33, 18261–18271.
Loizou and Richtárik (2017) Loizou, N. and P. Richtárik (2017). Linearly convergent stochastic heavy ball method for minimizing generalization error. arXiv preprint arXiv:1710.10737.
Loizou and Richtárik (2020) Loizou, N. and P. Richtárik (2020). Momentum and stochastic momentum for stochastic gradient, newton, proximal point and subspace descent methods. Computational Optimization and Applications 77(3), 653–710.
Mai and Johansson (2020) Mai, V. and M. Johansson (2020). Convergence of a stochastic gradient method with momentum for non-smooth non-convex optimization. In International Conference on Machine Learning, pp. 6630–6639.
Moulines and Bach (2011) Moulines, E. and F. Bach (2011). Non-asymptotic analysis of stochastic approximation algorithms for machine learning. Advances in Neural Information Processing Systems 24.
Nesterov (1983) Nesterov, Y. (1983). A method for unconstrained convex minimization problem with the rate of convergence $o(1/k^{2})$ . In Doklady AN USSR, Volume 269, pp. 543–547.
Nesterov (2003) Nesterov, Y. (2003). Introductory lectures on convex optimization: A basic course, Volume 87. Springer.
Nguyen et al. (2018) Nguyen, L., P. H. Nguyen, M. Dijk, P. Richtárik, K. Scheinberg, and M. Takác (2018). SGD and Hogwild! convergence without the bounded gradients assumption. In International Conference on Machine Learning, pp. 3750–3758.
Nocedal and Wright (2006) Nocedal, J. and S. J. Wright (2006). Numerical optimization. springer series in operations research. SIAM J Optimization.
Paquette and Paquette (2021) Paquette, C. and E. Paquette (2021). Dynamics of stochastic momentum methods on large-scale, quadratic models. Advances in Neural Information Processing Systems 34, 9229–9240.
Polyak (1964) Polyak, B. T. (1964). Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics 4(5), 1–17.
Polyak and Juditsky (1992) Polyak, B. T. and A. B. Juditsky (1992). Acceleration of stochastic approximation by averaging. SIAM Journal on Control and Optimization 30(4), 838–855.
Richtárik and Takác (2020) Richtárik, P. and M. Takác (2020). Stochastic reformulations of linear systems: algorithms and convergence theory. SIAM Journal on Matrix Analysis and Applications 41(2), 487–524.
Robbins and Monro (1951) Robbins, H. and S. Monro (1951). A stochastic approximation method. The Annals of Mathematical Statistics, 400–407.
Ruppert (1988) Ruppert, D. (1988). Efficient estimations from a slowly convergent Robbins-Monro process. Technical report, Cornell University Operations Research and Industrial Engineering.
Sebbouh et al. (2021) Sebbouh, O., R. M. Gower, and A. Defazio (2021). Almost sure convergence rates for stochastic gradient descent and stochastic heavy ball. In Conference on Learning Theory, pp. 3935–3971. PMLR.
Su and Zhu (2023) Su, W. J. and Y. Zhu (2023). Higrad: Uncertainty quantification for online learning and stochastic approximation. Journal of Machine Learning Research 24(124), 1–53.
Toulis et al. (2021) Toulis, P., T. Horel, and E. M. Airoldi (2021). The proximal Robbins-Monro method. Journal of the Royal Statistical Society Series B: Statistical Methodology 83(1), 188–212.
Wang and Johansson (2022) Wang, X. and M. Johansson (2022). On uniform boundedness properties of sgd and its momentum variants. arXiv preprint arXiv:2201.10245.
Yan et al. (2018) Yan, Y., T. Yang, Z. Li, Q. Lin, and Y. Yang (2018). A unified analysis of stochastic momentum methods for deep learning. arXiv preprint arXiv:1808.10396.
Yang et al. (2016) Yang, T., Q. Lin, and Z. Li (2016). Unified convergence analysis of stochastic momentum methods for convex and non-convex optimization. arXiv preprint arXiv:1604.03257.
Zeiler (2012) Zeiler, M. D. (2012). Adadelta: an adaptive learning rate method. arXiv preprint arXiv:1212.5701.
Zhu et al. (2023) Zhu, W., X. Chen, and W. B. Wu (2023). Online covariance matrix estimation in stochastic gradient descent. Journal of the American Statistical Association 118(541), 393–404.
Zhu and Dong (2021) Zhu, Y. and J. Dong (2021). On constructing confidence region for model parameters in stochastic gradient descent via batch means. In 2021 Winter Simulation Conference (WSC), pp. 1–12.

Appendix A Proofs of finite-time convergence rates of SGDM

A.1 Proof of Theorem 1

Recall that

\displaystyle m_{t+1}=\gamma m_{t}+(1-\gamma)\triangledown g_{\eta_{t}}(x_{t})% ,\quad x_{t+1}=x_{t}-\alpha m_{t+1}.

Put

\displaystyle\Gamma=\left(\begin{array}[]{ccc}\gamma I,&(1-\gamma)\Sigma\\ -\alpha\gamma I,&I-\alpha(1-\gamma)\Sigma\\ \end{array}\right).

Let $\widetilde{x}_{t}=x_{t}-x^{*}$ and $\widetilde{m}_{t+1}=(1-\gamma)\sum_{j=1}^{t}\gamma^{t-j}\Sigma\widetilde{x}_{j}$ , we can write

	$\displaystyle\left(\begin{array}[]{ccc}\widetilde{m}_{t+1}\\ \widetilde{x}_{t+1}\\ \end{array}\right)$		$\displaystyle=\Gamma\left(\begin{array}[]{ccc}\widetilde{m}_{t}\\ \widetilde{x}_{t}\\ \end{array}\right)-\alpha(1-\gamma)\left(\begin{array}[]{ccc}0\\ \sum_{j=1}^{t}\gamma^{t-j}(\triangledown g_{\eta_{j}}(x_{j})-\Sigma\widetilde{% x}_{j})\\ \end{array}\right)$		(A.35)
			$\displaystyle=\Gamma^{t}\left(\begin{array}[]{ccc}\widetilde{m}_{1}\\ \widetilde{x}_{1}\\ \end{array}\right)-\alpha(1-\gamma)\sum_{j=1}^{t}\Gamma^{t-j}\left(\begin{% array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}(\triangledown g_{\eta_{k}}(x_{k})-\Sigma\widetilde{% x}_{k})\\ \end{array}\right).$		(A.35)

Define

\displaystyle q_{t}:=\mathbb{E}[\left\|\widetilde{m}_{t}\right\|^{2}+\left\|% \widetilde{x}_{t}\right\|^{2}].

Recall for $\overline{L}=0$ , there hold that $\triangledown g(x)=\mathbb{E}[\triangledown g_{\eta}(x)]=\Sigma(x-x^{*})$ and

\displaystyle\triangledown g_{\eta_{k}}(x_{k})-\Sigma\widetilde{x}_{k}=\left(% \triangledown g_{\eta_{k}}(x_{k})-\triangledown g_{\eta_{k}}(x^{*})+% \triangledown g(x^{*})-\triangledown g(x_{k})\right)+\triangledown g_{\eta_{k}% }(x^{*}).

By the $L_{f}$ -smooth of the individual function $f_{\xi}(x)$ in (A3) and the independence of $L_{\xi}$ and $\widetilde{x}_{k}$ , we have

\displaystyle\mathbb{E}\left\|\triangledown f_{\xi}(x_{k})-\triangledown f_{% \xi}(x^{*})+\triangledown f(x^{*})-\triangledown f(x_{k})\right\|^{2}\leq% \mathbb{E}\left\|\triangledown f_{\xi}(x_{k})-\triangledown f_{\xi}(x^{*})% \right\|^{2}\leq L_{f}^{2}\mathbb{E}\left\|\widetilde{x}_{k}\right\|^{2}.

Then from the iteration (A.35), we have

	$\displaystyle q_{t+1}$		$\displaystyle\leq\left\\|\Gamma^{t}\right\\|^{2}q_{1}+2\alpha^{2}(1-\gamma)^{2}% \frac{\sigma^{2}}{B}\sum_{k=1}^{t}\left(\sum_{j=k}^{t}\left\\|\Gamma^{t-j}% \right\\|\gamma^{j-k}\right)^{2}$		(A.36)
			$\displaystyle+2\alpha^{2}(1-\gamma)^{2}\frac{L_{f}^{2}}{B}\sum_{k=1}^{t}\left(% \sum_{j=k}^{t}\left\\|\Gamma^{t-j}\right\\|\gamma^{j-k}\right)^{2}q_{k}.$		(A.36)

In the remaining part of proof, we use the inductive method to derive the bound of $q_{t}$ . Let $C_{1}\geq 1$ and $C_{2}\geq 1$ be some constants which will be specified later. We will prove that, for any $t\geq 1$ , the inequalities

\displaystyle q_{j}\leq\frac{C_{1}}{B(1-\lambda)}\alpha^{2}\sigma^{2}+C_{2}q_{% 1}\lambda^{2(1-\delta)(j-1)},\quad 1\leq j\leq t

(A.37)

with fixed $\delta\in(0,1]$ , imply that

\displaystyle q_{t+1}\leq\frac{C_{1}}{B(1-\lambda)}\alpha^{2}\sigma^{2}+C_{2}q% _{1}\lambda^{2(1-\delta)t}.

(A.38)

By Lemma 20 below, we have $\left\|\Gamma^{j}\right\|\leq M\lambda^{j}$ for all $j\geq 0$ . Theorem 3 provides that $\lambda\geq\sqrt{\gamma}$ . By applying Lemma 27, we are equipped to manage the summation presented in (A.36). Consequently, from (A.36) and (A.37), we deduce that

	$\displaystyle q_{t+1}$		$\displaystyle\leq M^{2}\lambda^{2t}q_{1}+2\alpha^{2}\frac{\sigma^{2}}{B}M^{2}% \frac{2}{1-\lambda}$
			$\displaystyle+2\alpha^{2}(1-\gamma)^{2}\frac{L_{f}^{2}}{B}M^{2}\left(\frac{2C_% {1}\alpha^{2}\sigma^{2}}{B(1-\gamma)^{2}(1-\lambda)^{2}}+C_{2}q_{1}\lambda^{2(% 1-\delta)t}\frac{1}{\lambda^{2(1-\delta)}}\frac{4}{\delta(1-\gamma)^{2}(1-% \lambda)}\right)$
			$\displaystyle\leq\frac{C_{1}}{B(1-\lambda)}\alpha^{2}\sigma^{2}M^{2}\left(% \frac{4}{C_{1}}+\frac{4\alpha^{2}L_{f}^{2}}{B(1-\lambda)}\right)+C_{2}q_{1}% \lambda^{2(1-\delta)t}M^{2}\left(\frac{\lambda^{2\delta t}}{C_{2}}+\frac{8% \alpha^{2}L_{f}^{2}}{B\lambda^{2(1-\delta)}\delta(1-\lambda)}\right).$

Let the learning rate $\alpha$ and the batch size $B$ satisfy

\displaystyle 8M^{2}\frac{1}{\lambda^{2(1-\delta)}}\frac{1}{\delta(1-\lambda)}% \alpha^{2}\frac{L_{f}^{2}}{B}\leq\frac{1}{2},

we can take

\displaystyle C_{1}=8M^{2},\quad C_{2}=2M^{2}.

It is straightforward to verify that (A.38) is satisfied, and the induction is completed.

A.2 Proof of Theorem 3

Put

\displaystyle\Gamma=\left(\begin{array}[]{ccc}U,&0_{d}\\ 0_{d},&U\\ \end{array}\right)\left(\begin{array}[]{ccc}\gamma I,&(1-\gamma)A\\ -\alpha\gamma I,&I-\alpha(1-\gamma)A\\ \end{array}\right)\left(\begin{array}[]{ccc}U^{\top},&0_{d}\\ 0_{d},&U^{\top}\\ \end{array}\right)

where $\Sigma=Udiag(\kappa_{1},...,\kappa_{d})U^{T}=:UAU^{T}$ , $\mu=\kappa_{1}\leq\cdots\leq\kappa_{d}=L$ , $U$ is an orthogonal matrix. Now we can define

\displaystyle\lambda^{\pm}_{k}=\frac{\gamma+1-\alpha(1-\gamma)\kappa_{k}\pm% \sqrt{(\alpha(1-\gamma)\kappa_{k}-\gamma-1)^{2}-4\gamma}}{2},

for $1\leq k\leq d$ , and the diagonal matrix

\displaystyle\Lambda=\left(\begin{array}[]{ccc}diag(\lambda^{+}_{1},...,% \lambda^{+}_{d}),&0_{d}\\ 0_{d},&diag(\lambda^{-}_{1},...,\lambda^{-}_{d})\\ \end{array}\right).

(A.42)

Notice that if $(\alpha(1-\gamma)\kappa_{k}-\gamma-1)^{2}<4\gamma$ , $\lambda^{\pm}_{k}$ is complex and the definition is

\displaystyle\lambda^{\pm}_{k}=\frac{\gamma+1-\alpha(1-\gamma)\kappa_{k}\pm% \sqrt{-1}\sqrt{4\gamma-(\alpha(1-\gamma)\kappa_{k}-\gamma-1)^{2}}}{2}.

Prior to establishing Theorem 3, we first demonstrate the diagonalization of $\Gamma$ , with $\Lambda$ representing the resultant diagonal matrix. Subsequent to this, we proceed to delineate the bounds of the spectral radius of $\Lambda$ .

Lemma 20

For $\alpha,\gamma,L$ satisfying that $\alpha L<2(1+\gamma)/(1-\gamma)$ and $\left|\gamma+1-\alpha(1-\gamma)\kappa_{k}\right|\neq 2\sqrt{\gamma}$ for all $1\leq k\leq d$ , we have

\displaystyle\left(\begin{array}[]{ccc}\gamma I,&(1-\gamma)A\\ -\alpha\gamma I,&I-\alpha(1-\gamma)A\\ \end{array}\right)=P\Lambda P^{-1},

for some invertible matrix $P$ that satisfies

\displaystyle\left\|P\right\|\leq 2,\quad\left\|P^{-1}\right\|\leq\frac{2}{% \sqrt{\Delta}}\left(2(1-\gamma)(1+\alpha L+L)+3\alpha\gamma\right),

where $\Delta=\min_{k}\left\{\left|\left(\gamma+1-\alpha(1-\gamma)\kappa_{k}\right)^{% 2}-4\gamma\right|\right\}$ . Therefore, $\left\|\Gamma^{j}\right\|\leq M\lambda^{j}$ for any $j\geq 1$ , where

\displaystyle M=\frac{4}{\sqrt{\Delta}}\left(2(1-\gamma)(1+\alpha L+L)+3\alpha% \gamma\right),

and $\lambda$ is the spectral radius of $\Lambda$ .

Proof. The eigenvalues satisfy

	$\displaystyle(\lambda^{\pm}_{k})^{2}+(-\gamma-1+\alpha(1-\gamma)\kappa_{k})% \lambda^{\pm}_{k}+\gamma=0,$
	$\displaystyle\lambda^{+}_{k}+\lambda^{-}_{k}=\gamma+1-\alpha(1-\gamma)\kappa_{% k},\quad\lambda^{+}_{k}\lambda^{-}_{k}=\gamma.$

Given that $\alpha L<2(1+\gamma)/(1-\gamma)$ and considering $\mu\leq\kappa_{k}\leq L$ , it follows that $|\gamma+1-\alpha(1-\gamma)\kappa_{k}|<1+\gamma$ for all $k$ . Considering the product of the eigenvalues $\lambda^{+}_{k}\lambda^{-}_{k}=\gamma<1$ , subsequent calculations reveal that $\max_{k}\{|\lambda^{\pm}_{k}|\}<1$ , which holds even in the case of complex eigenvalues. Moreover, since $|\gamma+1-\alpha(1-\gamma)\kappa_{k}|\neq 2\sqrt{\gamma}$ for any $k$ , it follows that $\lambda^{+}_{k}\neq\lambda^{-}_{k}$ , confirming that the matrix is indeed diagonalizable.

Let $e_{k}$ be the unit vector which has 1 in the $k$ -th coordinate and others zero. By the definition of $A$ , we have

\displaystyle(\alpha\gamma(1-\gamma)+\alpha(\lambda^{\pm}_{k}-\gamma)(1-\gamma% ))Ae_{k}=(\lambda^{\pm}_{k}-\gamma)(1-\lambda^{\pm}_{k})e_{k},

which is equivalent to

\displaystyle(1-\gamma)Ae_{k}=\frac{\lambda^{\pm}_{k}-\gamma}{\alpha\gamma}% \Big{\{}(1-\lambda^{\pm}_{k})e_{k}-\alpha(1-\gamma)Ae_{k}\Big{\}}.

(A.44)

Now define $\boldsymbol{z}^{\pm}_{k}$ by

\displaystyle(\lambda^{\pm}_{k}-\gamma)\boldsymbol{z}^{\pm}_{k}=(1-\gamma)Ae_{% k},

which can be written as $\gamma\boldsymbol{z}_{k}^{\pm}+(1-\gamma)Ae_{k}=\lambda^{\pm}_{k}\boldsymbol{z% }_{k}^{\pm}$ . Combining this equation with (A.44), we have

\displaystyle(\lambda^{\pm}_{k}-\gamma)\boldsymbol{z}^{\pm}_{k}=\frac{\lambda^% {\pm}_{k}-\gamma}{\alpha\gamma}\Big{\{}(1-\lambda^{\pm}_{k})e_{k}-\alpha(1-% \gamma)Ae_{k}\Big{\}},

which yields that

\displaystyle-\alpha\gamma\boldsymbol{z}^{\pm}_{k}+e_{k}-\alpha(1-\gamma)Ae_{k% }=\lambda^{\pm}_{k}e_{k}.

Thus

\displaystyle\left(\begin{array}[]{ccc}\gamma I,&(1-\gamma)A\\ -\alpha\gamma I,&I-\alpha(1-\gamma)A\\ \end{array}\right)\left(\begin{array}[]{c}\boldsymbol{z}^{\pm}_{k}\\ e_{k}\\ \end{array}\right)=\lambda^{\pm}_{k}\left(\begin{array}[]{c}\boldsymbol{z}^{% \pm}_{k}\\ e_{k}\\ \end{array}\right)=\left(\begin{array}[]{c}\boldsymbol{z}^{\pm}_{k}\\ e_{k}\\ \end{array}\right)\lambda^{\pm}_{k}.

(A.53)

Define

\displaystyle\boldsymbol{z}^{+}=(\boldsymbol{z}^{+}_{1},...,\boldsymbol{z}^{+}% _{d}),\quad\boldsymbol{z}^{-}=(\boldsymbol{z}^{-}_{1},...,\boldsymbol{z}^{-}_{% d}),

and

	$\displaystyle P$		$\displaystyle=\left(\begin{array}[]{cc}\boldsymbol{z}^{+}\left(\left\|% \boldsymbol{z}^{+}\right\|^{2}+I\right)^{-\frac{1}{2}},&\boldsymbol{z}^{-}\left% (\left\|\boldsymbol{z}^{-}\right\|^{2}+I\right)^{-\frac{1}{2}}\\ I\left(\left\|\boldsymbol{z}^{+}\right\|^{2}+I\right)^{-\frac{1}{2}},&I\left(% \left\|\boldsymbol{z}^{-}\right\|^{2}+I\right)^{-\frac{1}{2}}\end{array}\right),$
	$\displaystyle P^{-1}$		$\displaystyle=\left(\begin{array}[]{cc}\left(\left\|\boldsymbol{z}^{+}\right\|^{% 2}+I\right)^{\frac{1}{2}}(\boldsymbol{z}^{+}-\boldsymbol{z}^{-})^{-1},&-\left(% \left\|\boldsymbol{z}^{+}\right\|^{2}+I\right)^{\frac{1}{2}}\boldsymbol{z}^{-}(% \boldsymbol{z}^{+}-\boldsymbol{z}^{-})^{-1}\\ -\left(\left\|\boldsymbol{z}^{-}\right\|^{2}+I\right)^{\frac{1}{2}}(\boldsymbol{% z}^{+}-\boldsymbol{z}^{-})^{-1},&\left(\left\|\boldsymbol{z}^{-}\right\|^{2}+I% \right)^{\frac{1}{2}}\boldsymbol{z}^{+}(\boldsymbol{z}^{+}-\boldsymbol{z}^{-})% ^{-1}\end{array}\right).$

Note that $\boldsymbol{z}^{+}$ and $\boldsymbol{z}^{-}$ are diagonal matrices. Then by (A.53),

\displaystyle\left(\begin{array}[]{ccc}\gamma I,&(1-\gamma)A\\ -\alpha\gamma I,&I-\alpha(1-\gamma)A\\ \end{array}\right)P=P\Lambda,

where $\Lambda$ is defined by (A.42). By the definition of $\boldsymbol{z}^{\pm}$ , we have

	$\displaystyle\boldsymbol{z}^{\pm}$		$\displaystyle=diag\left(\frac{\lambda_{k}^{\mp}-\gamma}{\alpha\gamma}\right),$
	$\displaystyle\boldsymbol{z}^{+}\boldsymbol{z}^{-}$		$\displaystyle=diag\left(\frac{(1-\gamma)\kappa_{k}}{\alpha\gamma}\right),$
	$\displaystyle(\boldsymbol{z}^{+}-\boldsymbol{z}^{-})^{-1}$		$\displaystyle=diag\left(\frac{-\alpha\gamma}{\sqrt{(\gamma+1-\alpha(1-\gamma)% \kappa_{k})^{2}-4\gamma}}\right).$

For $\left|\boldsymbol{z}_{k}^{+}\right|+\left|\boldsymbol{z}_{k}^{-}\right|$ , we have

\displaystyle\left|\boldsymbol{z}_{k}^{+}\right|+\left|\boldsymbol{z}_{k}^{-}% \right|\leq 2\sqrt{\left|\boldsymbol{z}_{k}^{+}\right|^{2}+\left|\boldsymbol{z% }_{k}^{-}\right|^{2}}=\left\{\begin{array}[]{cc}\frac{2}{\alpha\gamma}\sqrt{2(% \lambda_{k}^{+}-\gamma)(\lambda_{k}^{-}-\gamma)},&(\gamma+1-\alpha(1-\gamma)% \kappa_{k})^{2}<4\gamma,\\ \frac{2}{\alpha\gamma}\sqrt{(\lambda_{k}^{+}-\gamma)^{2}+(\lambda_{k}^{-}-% \gamma)^{2}},&(\gamma+1-\alpha(1-\gamma)\kappa_{k})^{2}>4\gamma,\end{array}\right.

and

	$\displaystyle\frac{2}{\alpha\gamma}\sqrt{2(\lambda_{k}^{+}-\gamma)(\lambda_{k}% ^{-}-\gamma)}$		$\displaystyle=\frac{2}{\alpha\gamma}\sqrt{2\alpha\gamma(1-\gamma)\kappa_{k}},$
	$\displaystyle\frac{2}{\alpha\gamma}\sqrt{(\lambda_{k}^{+}-\gamma)^{2}+(\lambda% _{k}^{-}-\gamma)^{2}}$		$\displaystyle=\frac{2}{\alpha\gamma}\sqrt{(\lambda_{k}^{+}+\lambda_{k}^{-}-% \gamma)^{2}-2\lambda_{k}^{+}\lambda_{k}^{-}+\gamma^{2}}$
			$\displaystyle=\frac{2}{\alpha\gamma}\sqrt{(1-\gamma)^{2}-2\alpha(1-\gamma)% \kappa_{k}+\alpha^{2}(1-\gamma)^{2}\kappa_{k}^{2}}\leq\frac{2}{\alpha\gamma}(1% -\gamma)\left(1+\alpha\kappa_{k}\right).$

So it holds that

	$\displaystyle\left\|\boldsymbol{z}_{k}^{+}\right\|+\left\|\boldsymbol{z}_{k}^{-}\right\|$		$\displaystyle\leq\frac{2}{\alpha\gamma}\left(\sqrt{2\alpha\gamma(1-\gamma)% \kappa_{k}}+(1-\gamma)(1+\alpha\kappa_{k})\right)$
			$\displaystyle\leq\frac{1}{\alpha\gamma}\left(2\alpha\gamma+(1-\gamma)(2+2% \alpha\kappa_{k}+\kappa_{k})\right).$

Then for $P$ and $P^{-1}$ , it holds that

	$\displaystyle\left\\|P\right\\|$		$\displaystyle\leq\sqrt{4}=2,$
	$\displaystyle\left\\|P^{-1}\right\\|$		$\displaystyle\leq 2\max_{k}\left\|\sqrt{\left\|\boldsymbol{z}^{+}_{k}\right\|^{2}% +1}\sqrt{\left\|\boldsymbol{z}^{-}_{k}\right\|^{2}+1}(\boldsymbol{z}^{+}-% \boldsymbol{z}^{-})_{k}^{-1}\right\|$
			$\displaystyle\leq 2\max_{k}\left(\left\|(\boldsymbol{z}^{+}\boldsymbol{z}^{-})_% {k}\right\|+\left\|\boldsymbol{z}^{+}_{k}\right\|+\left\|\boldsymbol{z}^{-}_{k}% \right\|+1\right)\left\|(\boldsymbol{z}^{+}-\boldsymbol{z}^{-})_{k}^{-1}\right\|$
			$\displaystyle\leq 2\max_{k}\left(\frac{(1-\gamma)\kappa_{k}}{\alpha\gamma}+% \frac{1}{\alpha\gamma}\left(2\alpha\gamma+(1-\gamma)(2+2\alpha\kappa_{k}+% \kappa_{k})\right)+1\right)\frac{\alpha\gamma}{\sqrt{\Delta}}$
			$\displaystyle\leq\frac{2}{\sqrt{\Delta}}\left(2(1-\gamma)(1+\alpha L+L)+3% \alpha\gamma\right),$

where $\Delta=\min_{k}\left\{\left|\left(\gamma+1-\alpha(1-\gamma)\kappa_{k}\right)^{% 2}-4\gamma\right|\right\}$ .

Put

\displaystyle\tilde{P}=\left(\begin{array}[]{cc}U,&0_{d}\\ 0_{d},&U\end{array}\right)P.

Then we have

\displaystyle\left\|\Gamma^{j}\right\|=\left\|\tilde{P}\Lambda^{j}\tilde{P}^{-% 1}\right\|\leq\left\|P\right\|\left\|P^{-1}\right\|\lambda^{j}\leq\frac{4}{% \sqrt{\Delta}}\left(2(1-\gamma)(1+\alpha L+L)+3\alpha\gamma\right)\lambda^{j}:% =M\lambda^{j}.

After diagonalizing $\Gamma$ , we can obtain an upper bound for the spectral radius $\lambda$ by evaluating the maximal absolute eigenvalue $\max_{k}\left|\lambda^{\pm}_{k}\right|$ in the diagonal matrix $\Lambda$ , as specified in definition (A.42).

Proof. [Theorem 3] From Lemma 20, the eigenvalues satisfy

\displaystyle\lambda^{+}_{k}+\lambda^{-}_{k}=\gamma+1-\alpha(1-\gamma)\kappa_{% k}\in[\gamma+1-\alpha(1-\gamma)L,\gamma+1-\alpha(1-\gamma)\mu],

For real eigenvalues, there exists $k$ such that $\left|\gamma+1-\alpha(1-\gamma)\kappa_{k}\right|>2\sqrt{\gamma}$ . By the fact that $\mu\leq\kappa_{k}\leq L$ , we have

$\displaystyle\lambda=\max_{k}\left\|\lambda^{\pm}_{k}\right\|$	$\displaystyle=\max\left\{\frac{\gamma+1-\alpha(1-\gamma)\mu+\sqrt{\left(\gamma% +1-\alpha(1-\gamma)\mu\right)^{2}-4\gamma}}{2},\right.$	(A.59)
	$\displaystyle\quad\left.\frac{\alpha(1-\gamma)L-\gamma-1+\sqrt{\left(\alpha(1-% \gamma)L-\gamma-1\right)^{2}-4\gamma}}{2}\right\}$
	$\displaystyle=\frac{\gamma+1-(1-\gamma)\phi+\sqrt{\left(\gamma+1-(1-\gamma)% \phi\right)^{2}-4\gamma}}{2}:=h(\gamma),$

where

\displaystyle\phi=\min\left\{\alpha\mu,\frac{2(1+\gamma)}{1-\gamma}-\alpha L% \right\}.

In this condition, there holds that $\gamma+1-(1-\gamma)\phi>2\sqrt{\gamma}$ . Define $\eta=\left(\gamma+1-(1-\gamma)\phi\right)^{2}-4\gamma>0$ , then

\displaystyle h^{\prime}(\gamma)=\frac{(1+\phi)h(\gamma)-1}{\sqrt{\eta}},\quad h% ^{\prime\prime}(\gamma)=\frac{h^{\prime}(\gamma)\left((1+\phi)\sqrt{\eta}+2-(1% +\phi)(\gamma+1-(1-\gamma)\phi)\right)}{2\eta}.

By the fact that $h(0)=1-\phi,h^{\prime}(0)<0$ and $2-(1+\phi)(\gamma+1-(1-\gamma)\phi)\geq 0$ , where $0\leq\gamma<(1-\phi)^{2}/(1+\phi)^{2}$ , we have that $h$ is concave and

\displaystyle 1-\phi-\frac{\phi(1+\phi)}{1-\phi}\gamma\leq\lambda\leq 1-\phi-% \frac{\phi^{2}}{1-\phi}\gamma.

For complex eigenvalues, which means that $\left|\gamma+1-\alpha(1-\gamma)\kappa_{k}\right|\leq 2\sqrt{\gamma}$ for all $k$ . Then $\gamma\geq(1-\phi)^{2}/(1+\phi)^{2}$ and

\displaystyle\lambda=\max_{k}\left|\lambda^{\pm}_{k}\right|=\max_{k}\left\{% \left|\frac{\gamma+1-\alpha(1-\gamma)\kappa_{k}\pm\sqrt{-1}\sqrt{4\gamma-(% \alpha(1-\gamma)\kappa_{k}-\gamma-1)^{2}}}{2}\right|\right\}=\sqrt{\gamma}.

A.3 Proof of Theorem 9

Recall the iteration (A.35) that

\displaystyle\left(\begin{array}[]{ccc}\widetilde{m}_{t+1}\\ \widetilde{x}_{t+1}\\ \end{array}\right)=\Gamma^{t}\left(\begin{array}[]{ccc}\widetilde{m}_{1}\\ \widetilde{x}_{1}\\ \end{array}\right)-\alpha(1-\gamma)\sum_{j=1}^{t}\Gamma^{t-j}\left(\begin{% array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}(\triangledown g_{\eta_{k}}(x_{k})-\Sigma\widetilde{% x}_{k})\\ \end{array}\right).

(A.66)

Define the events

\displaystyle E_{j}=\left\{\left\|\left(\begin{array}[]{ccc}\widetilde{m}_{j}% \\ \widetilde{x}_{j}\\ \end{array}\right)\right\|\leq\frac{C_{1}}{\sqrt{B(1-\lambda)}}\alpha\sigma+C_% {2}q_{1}^{1/2}\lambda^{(1-\delta)(j-1)}\right\},\quad 1\leq j\leq t.

For certain constants $C_{1},C_{2}\geq 1$ , which will be defined subsequently, the probability of event $E_{1}$ is $\mathbb{P}(E_{1})=1$ . In the remaining part of proof, we use the inductive method to demonstrate that, for any $t\geq 1$ , on the events $\cap_{j=1}^{t}E_{j}$ , $t\geq 1$ , the event $E_{t+1}$ occurs with at least $1-2/T^{2}$ probability.

For $\overline{L}>0$ , it holds that

	$\displaystyle\triangledown g_{\eta_{k}}(x_{k})-\Sigma\widetilde{x}_{k}$		$\displaystyle=\left(\triangledown g_{\eta_{k}}(x_{k})-\triangledown g_{\eta_{k% }}(x^{})+\triangledown g(x^{})-\triangledown g(x_{k})\right)+\triangledown g% _{\eta_{k}}(x^{*})+\triangledown g(x_{k})-\Sigma\widetilde{x}_{k}$		(A.68)
			$\displaystyle=\left(\triangledown g_{\eta_{k}}(x_{k})-\triangledown g_{\eta_{k% }}(x^{})+\triangledown g(x^{})-\triangledown g(x_{k})\right)+\triangledown g% _{\eta_{k}}(x^{*})+V_{k}\widetilde{x}_{k},$		(A.68)

where $\left\|V_{k}\right\|=\left\|\int_{0}^{1}\Sigma(x^{*}+y\widetilde{x}_{k})-% \Sigma(x^{*})dy\right\|\leq\overline{L}\left\|\widetilde{x}_{k}\right\|$ . By applying Lemmas 21, 22 and 23, we establish bounds for each of the three components in (A.68). Then from the iteration (A.35), with at least $1-2/T^{2}$ probability, we have that on the events $\cap_{j=1}^{t}E_{j}$

	$\displaystyle\left\\|\left(\begin{array}[]{ccc}\widetilde{m}_{t+1}\\ \widetilde{x}_{t+1}\\ \end{array}\right)\right\\|\leq M\lambda^{t}q_{1}^{1/2}+\alpha\frac{\sigma}{% \sqrt{B}}M\frac{\sqrt{2}c\left(2\log T+4d\right)}{(1-\lambda)^{1/2}}+\alpha M% \overline{L}\left(\frac{2C_{1}^{2}\alpha^{2}\sigma^{2}}{B(1-\lambda)^{2}}+% \frac{4C_{2}^{2}q_{1}\lambda^{(1-\delta)(t-1)}}{\delta(1-\lambda)}\right)$
	$\displaystyle+\alpha\frac{L_{f}}{\sqrt{B}}Mc(2\log T+4d)\left(\frac{4C_{1}% \alpha\sigma}{\sqrt{B}(1-\lambda)}+\frac{4\sqrt{2}C_{2}q_{1}^{1/2}\lambda^{(1-% \delta)(t-1)}}{\sqrt{\delta}(1-\lambda)^{1/2}}\right)$
	$\displaystyle\leq\frac{C_{1}}{\sqrt{B(1-\lambda)}}\alpha\sigma M\left(\frac{% \sqrt{2}c\left(2\log T+4d\right)}{C_{1}}+c(2\log T+4d)\frac{4\alpha L_{f}}{% \sqrt{B}(1-\lambda)^{1/2}}+\frac{2\overline{L}C_{1}\alpha^{2}\sigma}{\sqrt{B}(% 1-\lambda)^{3/2}}\right)$
	$\displaystyle+C_{2}q_{1}^{1/2}\lambda^{(1-\delta)t}M\left(\frac{\lambda^{% \delta t}}{C_{2}}+\frac{4\sqrt{2}c(2\log T+4d)\alpha L_{f}}{\sqrt{B}\lambda^{(% 1-\delta)}\sqrt{\delta}(1-\lambda)^{1/2}}+\frac{4\overline{L}\alpha C_{2}q_{1}% ^{1/2}}{\lambda^{(1-\delta)}\delta(1-\lambda)}\right),$

where $c>0$ is an absolute constant. To ensure the occurrence of event $E_{t+1}$ , the learning rate $\alpha$ and the batch size $B$ must fulfill the condition

\displaystyle M\left(\frac{4\sqrt{2}c(2\log T+4d)}{\lambda^{(1-\delta)}\sqrt{% \delta}(1-\lambda)^{1/2}}\frac{\alpha L_{f}}{\sqrt{B}}+\frac{2\overline{L}C_{1% }\alpha^{2}\sigma}{\sqrt{B}(1-\lambda)^{3/2}}\right)\leq\frac{1}{3},

where

\displaystyle C_{1}=3\sqrt{2}c\left(2\log T+4d\right)M.

Furthermore, we require that the initialization $q_{1}=\left\|\widetilde{m}_{1}\right\|^{2}+\left\|\widetilde{x}_{1}\right\|^{2% }=\left\|\widetilde{x}_{1}\right\|^{2}$ satisfy

\displaystyle 4M\overline{L}\alpha C_{2}q_{1}^{1/2}\frac{1}{\lambda^{1-\delta}% \delta(1-\lambda)}\leq\frac{1}{3},

where

\displaystyle C_{2}=3M.

It is straightforward to verify that

\displaystyle\mathbb{P}\left(E_{t+1}\mid\cap_{j=1}^{t}E_{j}\right)=\mathbb{P}% \left(\left\|\left(\begin{array}[]{ccc}\widetilde{m}_{t+1}\\ \widetilde{x}_{t+1}\\ \end{array}\right)\right\|\leq\frac{C_{1}}{\sqrt{B(1-\lambda)}}\alpha\sigma+C_% {2}q_{1}^{1/2}\lambda^{(1-\delta)t}\mid\cap_{j=1}^{t}E_{j}\right)\geq 1-\frac{% 2}{T^{2}},

and the induction is completed. Consequently, we deduce that the probability of the intersection of events from $E_{1}$ through $E_{T}$ is at least $\mathbb{P}\left(\bigcap_{j=1}^{T}E_{j}\right)\geq 1-\frac{2}{T}$ .

A.4 Supporting Lemmas for Bounds in Theorem 9 Proof

Lemma 21

Under (A1)-(A2) and (A3’), we have

\displaystyle\mathbb{P}\left(\left\|\sum_{j=1}^{t}\Gamma^{t-j}\left(\begin{% array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}\triangledown g_{\eta_{k}}(x^{*})\\ \end{array}\right)\right\|\geq M\frac{\sigma}{\sqrt{B}}\frac{\sqrt{2}c\left(2% \log T+4d\right)}{(1-\lambda)^{1/2}(1-\gamma)}\right)\leq\frac{1}{T^{2}},

for $1\leq t\leq T$ , where $c>0$ is an absolute constant.

Proof. Define

\displaystyle Y_{t,k}=\sum_{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\left(\begin{array% }[]{ccc}0\\ \triangledown g_{\eta_{k}}(x^{*})\\ \end{array}\right).

By the definition of the sub-exponential, $Y_{t,k}$ is $\left\|\sum_{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\right\|\frac{c\sigma}{\sqrt{B}}$ -sub-exponential random vector for some absolute constant $c>0$ , and $\{Y_{t,k}\}_{k=1}^{t}$ are independent. Then by Lemmas 27 and 20, we have

\displaystyle\sum_{k=1}^{t}\left\|\sum_{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\right% \|^{2}\leq M^{2}\frac{2}{(1-\lambda)(1-\gamma)^{2}}.

By applying Lemma 28 on $\{Y_{t,k}\}_{k=1}^{t}$ , with probability $1-1/T^{2}$ , we can obtain

\displaystyle\left\|\sum_{k=1}^{t}Y_{t,k}\right\|\leq M\frac{\sigma}{\sqrt{B}}% \frac{\sqrt{2}c\left(2\log T+4d\right)}{(1-\lambda)^{1/2}(1-\gamma)}.

Lemma 22

Under (A1)-(A2) and (A3’), suppose that the $t$ -th iteration $(\widetilde{m}_{t},\widetilde{x}_{t})$ satisfies

\displaystyle\sqrt{\left\|\widetilde{m}_{j}\right\|^{2}+\left\|\widetilde{x}_{% j}\right\|^{2}}\leq\frac{C_{1}}{\sqrt{B(1-\lambda)}}\alpha\sigma+C_{2}q_{1}^{1% /2}\lambda^{(1-\delta)(j-1)},

for $1\leq j\leq t$ , we have

\displaystyle\mathbb{P}\left(\left\|\sum_{j=1}^{t}\Gamma^{t-j}\left(\begin{% array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}\left(\triangledown g_{\eta_{k}}(x_{k})-% \triangledown g_{\eta_{k}}(x^{*})+\triangledown g(x^{*})-\triangledown g(x_{k}% )\right)\\ \end{array}\right)\right\|\geq C\right)\leq\frac{1}{T^{2}},

where

\displaystyle C=c(2\log T+4d)M\frac{L_{f}}{\sqrt{B}}\left(\frac{4C_{1}\alpha% \sigma}{\sqrt{B}(1-\lambda)(1-\gamma)}+\frac{4\sqrt{2}C_{2}q_{1}^{1/2}\lambda^% {(1-\delta)(t-1)}}{\sqrt{\delta}(1-\lambda)^{1/2}(1-\gamma)}\right),

and $c>0$ is an absolute constant.

Proof. Define

\displaystyle Y_{t,k}=\sum_{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\left(\begin{array% }[]{ccc}0\\ \triangledown g_{\eta_{k}}(x_{k})-\triangledown g_{\eta_{k}}(x^{*})+% \triangledown g(x^{*})-\triangledown g(x_{k})\\ \end{array}\right).

According to the iteration of SGDM, we have $\xi_{k}$ and $x_{k}$ are independent and $\{Y_{t,k}\}_{k=1}^{t}$ are martingale difference. Due to the $L_{f}$ -smooth of the individual gradient $\triangledown f_{\xi}(x)$ , it holds that

\displaystyle\max_{\left\|v\right\|=1}\mathbb{E}\left[\exp\left(\left|v^{\top}% Y_{t,k}\right|/\sigma_{k}\right)\right]\leq 2,

where

\displaystyle\sigma_{k}=c\frac{L_{f}}{\sqrt{B}}\left\|\sum_{j=k}^{t}\Gamma^{t-% j}\gamma^{j-k}\right\|\left\|\widetilde{x}_{k}\right\|,

and $c>0$ is an absolute constant. Then according to Lemmas 27 and 20, we have

\displaystyle\sum_{k=1}^{t}\left\|\sum_{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\right% \|^{2}\left\|\widetilde{x}_{k}\right\|^{2}\leq M^{2}\left(2C_{1}^{2}\alpha^{2}% \frac{\sigma^{2}}{B}\frac{2}{(1-\gamma)^{2}(1-\lambda)^{2}}+2C_{2}^{2}q_{1}% \frac{4\lambda^{2(1-\delta)(t-1)}}{\delta(1-\gamma)^{2}(1-\lambda)}\right).

By applying Lemma 28 on $\{Y_{t,k}\}_{k=1}^{t}$ , with probability $1-1/T^{2}$ , we can obtain the following result:

\displaystyle\left\|\sum_{j=1}^{t}Y_{t,j}\right\|\leq c(2\log T+4d)M\frac{L_{f% }}{\sqrt{B}}\left(\frac{4C_{1}\alpha\sigma}{\sqrt{B}(1-\lambda)(1-\gamma)}+% \frac{4\sqrt{2}C_{2}q_{1}^{1/2}\lambda^{(1-\delta)(t-1)}}{\sqrt{\delta}(1-% \lambda)^{1/2}(1-\gamma)}\right).

Lemma 23

Under (A1)-(A2) and (A3’), suppose that the $t$ -th iteration $(\widetilde{m}_{t},\widetilde{x}_{t})$ satisfies

\displaystyle\sqrt{\left\|\widetilde{m}_{j}\right\|^{2}+\left\|\widetilde{x}_{% j}\right\|^{2}}\leq\frac{C_{1}}{\sqrt{B(1-\lambda)}}\alpha\sigma+C_{2}q_{1}^{1% /2}\lambda^{(1-\delta)(j-1)},

for $1\leq j\leq t$ , for $\overline{L}>0$ , we have

\displaystyle\left\|\sum_{j=1}^{t}\Gamma^{t-j}\left(\begin{array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}V_{k}\widetilde{x}_{k}\\ \end{array}\right)\right\|\leq M\overline{L}\left(\frac{2C_{1}^{2}\alpha^{2}% \sigma^{2}}{B(1-\lambda)^{2}(1-\gamma)}+\frac{4C_{2}^{2}q_{1}\lambda^{(1-% \delta)(t-1)}}{\delta(1-\lambda)(1-\gamma)}\right).

Proof. By the fact that $\sum_{j=1}^{t}\sum_{k=1}^{j}\lambda^{t-j}\gamma^{j-k}\leq(1-\lambda)^{-1}(1-% \gamma)^{-1}$ and Lemma 27, we have

	$\displaystyle\left\\|\sum_{j=1}^{t}\Gamma^{t-j}\left(\begin{array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}V_{k}\widetilde{x}_{k}\\ \end{array}\right)\right\\|\leq\sum_{j=1}^{t}M\lambda^{t-j}\sum_{k=1}^{j}\gamma% ^{j-k}\overline{L}\left\\|\widetilde{x}_{k}\right\\|^{2}$
	$\displaystyle\leq M\overline{L}\sum_{j=1}^{t}\lambda^{t-j}\sum_{k=1}^{j}\gamma% ^{j-k}\left(\frac{2C_{1}^{2}}{B(1-\lambda)}\alpha^{2}\sigma^{2}+2C_{2}^{2}q_{1% }\lambda^{2(1-\delta)(k-1)}\right)$
	$\displaystyle\leq M\overline{L}\left(\frac{2C_{1}^{2}\alpha^{2}\sigma^{2}}{B(1% -\lambda)^{2}(1-\gamma)}+2C_{2}^{2}q_{1}\sum_{k=1}^{t}\sum_{j=k}^{t}\lambda^{t% -j}\gamma^{j-k}\lambda^{(1-\delta)(k-1)}\right)$
	$\displaystyle\leq M\overline{L}\left(\frac{2C_{1}^{2}\alpha^{2}\sigma^{2}}{B(1% -\lambda)^{2}(1-\gamma)}+2C_{2}^{2}q_{1}\frac{2\lambda^{(1-\delta)(t-1)}}{% \delta(1-\gamma)(1-\lambda)}\right).$

Appendix B Proofs of convergence rates of averaged SGDM

B.1 Proof of Theorem 13

By the definition of $\Gamma$ , one can easily check that

\displaystyle\alpha(I-\Gamma)^{-1}\left(\begin{array}[]{ccc}0\\ \triangledown g_{\eta_{t}}(x^{*})\\ \end{array}\right)=\alpha\left(\begin{array}[]{cc}I,&\frac{1}{\alpha}I\\ \frac{\gamma}{\gamma-1}\Sigma^{-1},&\frac{1}{\alpha}\Sigma^{-1}\end{array}% \right)\left(\begin{array}[]{ccc}0\\ \triangledown g_{\eta_{t}}(x^{*})\\ \end{array}\right)=\left(\begin{array}[]{ccc}\triangledown g_{\eta_{t}}(x^{*})% \\ \Sigma^{-1}\triangledown g_{\eta_{t}}(x^{*})\\ \end{array}\right),

and we have

	$\displaystyle(I-\Gamma)^{-1}\sum_{t=1}^{n}\left(\begin{array}[]{ccc}0\\ \triangledown g_{\eta_{t}}(x^{})\\ \end{array}\right)-(1-\gamma)\sum_{t=1}^{n}\sum_{j=1}^{t}\Gamma^{t-j}\left(% \begin{array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}\triangledown g_{\eta_{k}}(x^{})\\ \end{array}\right)$
	$\displaystyle=\sum_{k=1}^{n}\left((I-\Gamma)^{-1}-(1-\gamma)\sum_{t=k}^{n}\sum% _{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\right)\left(\begin{array}[]{ccc}0\\ \triangledown g_{\eta_{k}}(x^{*})\\ \end{array}\right).$

Then summing (A.35) over $t$ from $n_{0}+1$ to $n$ , for $\overline{L}=0$ , we obtain

	$\displaystyle\sum_{t=n_{0}+1}^{n}\left(\begin{array}[]{ccc}\widetilde{m}_{t+1}% \\ \widetilde{x}_{t+1}\\ \end{array}\right)+\sum_{t=n_{0}+1}^{n}\left(\begin{array}[]{ccc}\triangledown g% _{\eta_{t}}(x^{})\\ \Sigma^{-1}\triangledown g_{\eta_{t}}(x^{})\\ \end{array}\right)=\sum_{t=n_{0}+1}^{n}\left(\begin{array}[]{ccc}\widetilde{m}% _{t+1}\\ \widetilde{x}_{t+1}\\ \end{array}\right)+\alpha\sum_{t=n_{0}+1}^{n}(I-\Gamma)^{-1}\left(\begin{array% }[]{ccc}0\\ \triangledown g_{\eta_{t}}(x^{*})\\ \end{array}\right)$		(B.88)
	$\displaystyle=-\alpha(1-\gamma)\sum_{t=1}^{n}\sum_{j=1}^{t}\Gamma^{t-j}\left(% \begin{array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}\left(\triangledown g_{\eta_{k}}(x_{k})-% \triangledown g_{\eta_{k}}(x^{})+\triangledown g(x^{})-\triangledown g(x_{k}% )\right)\\ \end{array}\right)$		(B.91)
	$\displaystyle+\alpha(1-\gamma)\sum_{t=1}^{n_{0}}\sum_{j=1}^{t}\Gamma^{t-j}% \left(\begin{array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}\left(\triangledown g_{\eta_{k}}(x_{k})-% \triangledown g_{\eta_{k}}(x^{})+\triangledown g(x^{})-\triangledown g(x_{k}% )\right)\\ \end{array}\right)$		(B.94)
	$\displaystyle+\alpha\sum_{k=1}^{n}\left((I-\Gamma)^{-1}-(1-\gamma)\sum_{t=k}^{% n}\sum_{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\right)\left(\begin{array}[]{ccc}0\\ \triangledown g_{\eta_{k}}(x^{*})\\ \end{array}\right)$		(B.97)
	$\displaystyle-\alpha\sum_{k=1}^{n_{0}}\left((I-\Gamma)^{-1}-(1-\gamma)\sum_{t=% k}^{n_{0}}\sum_{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\right)\left(\begin{array}[]{% ccc}0\\ \triangledown g_{\eta_{k}}(x^{*})\\ \end{array}\right)+\sum_{t=n_{0}+1}^{n}\Gamma^{t}\left(\begin{array}[]{ccc}% \widetilde{m}_{1}\\ \widetilde{x}_{1}\\ \end{array}\right).$		(B.102)

The subsequent proof aims to calculate the expected values of various components within the equation (B.102). For the $t$ -th iteration $(\widetilde{m}_{t},\widetilde{x}_{t})$ satisfying

\displaystyle\mathbb{E}[\left\|\widetilde{m}_{t}\right\|^{2}+\left\|\widetilde% {x}_{t}\right\|^{2}]\leq\frac{C_{1}}{B(1-\lambda)}\alpha^{2}\sigma^{2}+C_{2}q_% {1}\lambda^{2(1-\delta)(t-1)},\quad t\geq 1,

by the fact that $\sum_{t=k}^{n}\sum_{j=k}^{t}\lambda^{t-j}\gamma^{j-k}\leq(1-\lambda)^{-1}(1-% \gamma)^{-1}$ and $\frac{1-x}{1-x^{\delta}}\leq\delta^{-1}$ for $\delta\in(0,1]$ , we have

	$\displaystyle\mathbb{E}\left\\|\sum_{t=1}^{n}\alpha(1-\gamma)\sum_{j=1}^{t}% \Gamma^{t-j}\left(\begin{array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}\left(\triangledown g_{\eta_{k}}(x_{k})-% \triangledown g_{\eta_{k}}(x^{})+\triangledown g(x^{})-\triangledown g(x_{k}% )\right)\\ \end{array}\right)\right\\|^{2}$		(B.105)
	$\displaystyle\leq\alpha^{2}(1-\gamma)^{2}\frac{L_{f}^{2}}{B}M^{2}\sum_{k=1}^{n% }\left(\sum_{t=k}^{n}\sum_{j=k}^{t}\lambda^{t-j}\gamma^{j-k}\right)^{2}\left(% \frac{C_{1}}{B(1-\lambda)}\alpha^{2}\sigma^{2}+C_{2}q_{1}\lambda^{2(1-\delta)(% k-1)}\right)$
	$\displaystyle\leq\alpha^{2}\frac{L_{f}^{2}}{B}M^{2}\left(\frac{C_{1}\alpha^{2}% \sigma^{2}}{B(1-\lambda)^{3}}n+\frac{C_{2}q_{1}}{(1-\lambda)^{2}(1-\lambda^{2(% 1-\delta)})}\right)\leq\alpha^{2}\frac{L_{f}^{2}}{B}M^{2}\left(\frac{C_{1}% \alpha^{2}\sigma^{2}}{B(1-\lambda)^{3}}n+\frac{C_{2}q_{1}}{(1-\delta)(1-% \lambda)^{3}}\right).$		(B.106)

Incorporating Lemmas 27 and 29, it is established that

	$\displaystyle\mathbb{E}\left\\|\alpha\sum_{k=1}^{n}\left((I-\Gamma)^{-1}-(1-% \gamma)\sum_{t=k}^{n}\sum_{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\right)\left(\begin% {array}[]{ccc}0\\ \triangledown g_{\eta_{k}}(x^{*})\\ \end{array}\right)\right\\|^{2}$		(B.109)
	$\displaystyle\leq\alpha^{2}\frac{\sigma^{2}}{B}M^{2}\sum_{k=1}^{n}\left(2% \gamma^{2}\left(\sum_{j=k}^{n}\lambda^{n-j}\gamma^{j-k}\right)^{2}+\frac{2% \lambda^{2(n-k+1)}}{(1-\lambda)^{2}}\right)$
	$\displaystyle\leq 2\alpha^{2}\frac{\sigma^{2}}{B}M^{2}\left(\frac{2\gamma^{2}}% {(1-\gamma)^{2}(1-\lambda)}+\frac{\lambda^{2}}{(1-\lambda)^{3}}\right)\leq 6% \alpha^{2}\frac{\sigma^{2}}{B}M^{2}\frac{1}{(1-\lambda)^{3}},$		(B.110)

where the last inequality is due to the fact that $\sqrt{\gamma}\leq\lambda<1$ in Theorem 3. Furthermore, according to Lemma 20, it holds that

\displaystyle\left\|\sum_{t=n_{0}+1}^{n}\Gamma^{t}\left(\begin{array}[]{ccc}% \widetilde{m}_{1}\\ \widetilde{x}_{1}\\ \end{array}\right)\right\|^{2}\leq M^{2}\frac{\lambda^{2n_{0}}}{(1-\lambda)^{2% }}q_{1}.

(B.113)

Combining the above inequalities (B.106)(B.110)(B.113), we deduce that

	$\displaystyle(n-n_{0})^{2}\mathbb{E}\left\\|\frac{\sum_{t=n_{0}+1}^{n}% \widetilde{x}_{t}}{n-n_{0}}+\frac{\sum_{k=n_{0}+1}^{n}\Sigma^{-1}\triangledown g% _{\eta_{k}}(x^{*})}{n-n_{0}}\right\\|^{2}$
	$\displaystyle\leq 5\alpha^{2}\frac{L_{f}^{2}}{B}M^{2}\left(\frac{C_{1}\alpha^{% 2}\sigma^{2}}{B(1-\lambda)^{3}}(n+n_{0})+\frac{2C_{2}q_{1}}{(1-\delta)(1-% \lambda)^{3}}\right)$
	$\displaystyle+30\alpha^{2}\frac{\sigma^{2}}{B}M^{2}\frac{1}{(1-\lambda)^{3}}+M% ^{2}\frac{\lambda^{2n_{0}}}{(1-\lambda)^{2}}q_{1}.$

Take $n_{0}$ such that $\lambda^{2n_{0}}=B^{-1}(1-\lambda)$ and $n\geq 2n_{0}$ , it holds that

\displaystyle\mathbb{E}\left\|\frac{\sum_{t=n_{0}+1}^{n}\widetilde{x}_{t}}{n-n% _{0}}+\frac{\sum_{k=n_{0}+1}^{n}\Sigma^{-1}\triangledown g_{\eta_{k}}(x^{*})}{% n-n_{0}}\right\|^{2}\leq\frac{\tilde{C}_{1}}{Bn^{2}}+\frac{\tilde{C}_{2}}{B^{2% }n},

where

	$\displaystyle\tilde{C_{1}}$		$\displaystyle=40\alpha^{2}L_{f}^{2}M^{2}\frac{C_{2}q_{1}}{(1-\delta)(1-\lambda% )^{3}}+120\alpha^{2}\sigma^{2}M^{2}\frac{1}{(1-\lambda)^{3}}+\frac{4M^{2}q_{1}% }{1-\lambda},$
	$\displaystyle\tilde{C_{2}}$		$\displaystyle=30\alpha^{2}\sigma^{2}L_{f}^{2}M^{2}\frac{C_{1}\alpha^{2}}{(1-% \lambda)^{3}}.$

B.2 Proof of Theorem 17

From (A.35), for $\overline{L}>0$ , we have

	$\displaystyle\sum_{t=n_{0}+1}^{n}\left(\begin{array}[]{ccc}\widetilde{m}_{t+1}% \\ \widetilde{x}_{t+1}\\ \end{array}\right)+\sum_{t=n_{0}+1}^{n}\left(\begin{array}[]{ccc}\triangledown g% _{\eta_{t}}(x^{})\\ \Sigma^{-1}\triangledown g_{\eta_{t}}(x^{})\\ \end{array}\right)$
	$\displaystyle=-\alpha(1-\gamma)\sum_{t=1}^{n}\sum_{j=1}^{t}\Gamma^{t-j}\left(% \begin{array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}\left(\triangledown g_{\eta_{k}}(x_{k})-% \triangledown g_{\eta_{k}}(x^{})+\triangledown g(x^{})-\triangledown g(x_{k}% )\right)\\ \end{array}\right)$
	$\displaystyle+\alpha(1-\gamma)\sum_{t=1}^{n_{0}}\sum_{j=1}^{t}\Gamma^{t-j}% \left(\begin{array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}\left(\triangledown g_{\eta_{k}}(x_{k})-% \triangledown g_{\eta_{k}}(x^{})+\triangledown g(x^{})-\triangledown g(x_{k}% )\right)\\ \end{array}\right)$
	$\displaystyle+\alpha\sum_{k=1}^{n}\left((I-\Gamma)^{-1}-(1-\gamma)\sum_{t=k}^{% n}\sum_{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\right)\left(\begin{array}[]{ccc}0\\ \triangledown g_{\eta_{k}}(x^{*})\\ \end{array}\right)$
	$\displaystyle-\alpha\sum_{k=1}^{n_{0}}\left((I-\Gamma)^{-1}-(1-\gamma)\sum_{t=% k}^{n_{0}}\sum_{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\right)\left(\begin{array}[]{% ccc}0\\ \triangledown g_{\eta_{k}}(x^{*})\\ \end{array}\right)$
	$\displaystyle-\alpha(1-\gamma)\sum_{t=n_{0}+1}^{n}\sum_{j=1}^{t}\Gamma^{t-j}% \left(\begin{array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}V_{k}\widetilde{x}_{k}\\ \end{array}\right)+\sum_{t=n_{0}+1}^{n}\Gamma^{t}\left(\begin{array}[]{ccc}% \widetilde{m}_{1}\\ \widetilde{x}_{1}\\ \end{array}\right).$

By applying latter Lemmas 24, 25 and 26, we establish bounds for each respective component. Consequently, with at least probability $1-4/T$ , for fixed $\delta\in(0,\frac{1}{2}]$ , we have

	$\displaystyle(n-n_{0})\left\\|\frac{\sum_{t=n_{0}+1}^{n}\widetilde{x}_{t}}{n-n_% {0}}+\frac{\sum_{k=n_{0}+1}^{n}\Sigma^{-1}\triangledown g_{\eta_{k}}(x^{*})}{n% -n_{0}}\right\\|$
	$\displaystyle\leq c(2\log T+4d)\alpha\frac{L_{f}}{\sqrt{B}}M\left(\frac{2\sqrt% {2}C_{1}\alpha\sigma}{\sqrt{B}(1-\lambda)^{3/2}}(\sqrt{n}+\sqrt{n_{0}})+\frac{% 8C_{2}q_{1}^{1/2}}{(1-\lambda)^{3/2}}\right)$
	$\displaystyle+2\alpha\frac{\sigma}{\sqrt{B}}M\frac{\sqrt{3}c(2\log T+4d)}{(1-% \lambda)^{3/2}}+\alpha M\overline{L}\left(\frac{2C_{1}^{2}\alpha^{2}\sigma^{2}% }{B(1-\lambda)^{2}}n+\frac{4C_{2}^{2}q_{1}\left(n_{0}(1-\lambda)+1\right)% \lambda^{n_{0}}}{(1-\lambda)^{2}}\right)+M\frac{\lambda^{n_{0}}}{1-\lambda}q_{% 1}^{1/2},$

where $c>0$ is an absolute constant. Take $n_{0}$ such that $\lambda^{n_{0}}=B^{-1/2}(1-\lambda)^{1/2}$ and $n\geq 2n_{0}$ , then we have

\displaystyle\left\|\frac{\sum_{t=n_{0}+1}^{n}\widetilde{x}_{t}}{n-n_{0}}+% \frac{\sum_{t=n_{0}+1}^{n}\Sigma^{-1}\triangledown g_{\eta_{t}}(x^{*})}{n-n_{0% }}\right\|\leq\frac{\tilde{C}_{1}+\tilde{C}_{3}}{\sqrt{B}n}+\frac{\tilde{C}_{2% }}{B\sqrt{n}}+\frac{\tilde{C}_{4}}{B},

(B.120)

where

	$\displaystyle\tilde{C}_{1}$		$\displaystyle=\alpha L_{f}M\frac{16c(2\log T+4d)C_{2}q_{1}^{1/2}}{(1-\lambda)^% {3/2}}+\alpha\sigma M\frac{4\sqrt{3}c(2\log T+4d)}{(1-\lambda)^{3/2}}+\frac{2% Mq_{1}^{1/2}}{(1-\lambda)^{1/2}},$
	$\displaystyle\tilde{C}_{2}$		$\displaystyle=\alpha L_{f}\sigma M\frac{8\sqrt{2}c(2\log T+4d)C_{1}\alpha}{(1-% \lambda)^{3/2}},$
	$\displaystyle\tilde{C}_{3}$		$\displaystyle=8\alpha M\overline{L}\frac{C_{2}^{2}q_{1}(n_{0}(1-\lambda)+1)}{(% 1-\lambda)^{3/2}},\quad\tilde{C}_{4}=4\alpha\sigma^{2}M\overline{L}\frac{C_{1}% ^{2}\alpha^{2}}{(1-\lambda)^{2}}.$

So (B.120) holds for $2n_{0}\leq n\leq T$ with probability $1-4/T$ .

B.3 Supporting Lemmas for Bounds in Theorem 17 Proof

Lemma 24

Under (A1)-(A2) and (A3’), suppose that the $t$ -th iteration $(\widetilde{m}_{t},\widetilde{x}_{t})$ satisfies

\displaystyle\sqrt{\left\|\widetilde{m}_{j}\right\|^{2}+\left\|\widetilde{x}_{% j}\right\|^{2}}\leq\frac{C_{1}}{\sqrt{B(1-\lambda)}}\alpha\sigma+C_{2}q_{1}^{1% /2}\lambda^{(1-\delta)(j-1)},

for $1\leq j\leq t$ , we have

\displaystyle\mathbb{P}\left(\left\|\alpha(1-\gamma)\sum_{t=1}^{n}\sum_{j=1}^{% t}\Gamma^{t-j}\left(\begin{array}[]{ccc}0\\ \sum_{k=1}^{j}\gamma^{j-k}\left(\triangledown g_{\eta_{k}}(x_{k})-% \triangledown g_{\eta_{k}}(x^{*})+\triangledown g(x^{*})-\triangledown g(x_{k}% )\right)\\ \end{array}\right)\right\|\geq C\right)\leq\frac{1}{T^{2}},

where

\displaystyle C

\displaystyle=c(2\log T+4d)\alpha\frac{L_{f}}{\sqrt{B}}M\left(\frac{2\sqrt{2}C% _{1}\alpha\sigma}{\sqrt{B}(1-\lambda)^{3/2}}\sqrt{n}+\frac{2\sqrt{2}C_{2}q_{1}% ^{1/2}}{(1-\delta)^{1/2}(1-\lambda)^{3/2}}\right),

and $c>0$ is an absolute constant.

Proof. Define

\displaystyle Y_{n,k}=(1-\gamma)\sum_{t=k}^{n}\sum_{j=k}^{t}\Gamma^{t-j}\gamma% ^{j-k}\left(\begin{array}[]{ccc}0\\ \triangledown g_{\eta_{k}}(x_{k})-\triangledown g_{\eta_{k}}(x^{*})+% \triangledown g(x^{*})-\triangledown g(x_{k})\\ \end{array}\right).

\displaystyle\max_{\left\|v\right\|=1}\mathbb{E}\left[\exp\left(\left|v^{\top}% Y_{t,k}\right|/\sigma_{k}\right)\right]\leq 2,

where

\displaystyle\sigma_{k}=c\frac{L_{f}}{\sqrt{B}}(1-\gamma)\left\|\sum_{t=k}^{n}% \sum_{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\right\|\left\|\widetilde{x}_{k}\right\|,

and $c>0$ is an absolute constant. By the fact that $\sum_{t=k}^{n}\sum_{j=k}^{t}\lambda^{t-j}\gamma^{j-k}\leq(1-\gamma)^{-1}(1-% \lambda)^{-1}$ and $\frac{1-x}{1-x^{\delta}}\leq\delta^{-1}$ for $\delta\in(0,1]$ , it holds that

	$\displaystyle\sum_{k=1}^{n}(1-\gamma)^{2}\left\\|\sum_{t=k}^{n}\sum_{j=k}^{t}% \Gamma^{t-j}\gamma^{j-k}\right\\|^{2}\left\\|\widetilde{x}_{k}\right\\|^{2}$
	$\displaystyle\leq M^{2}(1-\gamma)^{2}\sum_{k=1}^{n}\left(\sum_{t=k}^{n}\sum_{j% =k}^{t}\lambda^{t-j}\gamma^{j-k}\right)^{2}\left(\frac{2C_{1}^{2}}{B(1-\lambda% )}\alpha^{2}\sigma^{2}+2C_{2}^{2}q_{1}\lambda^{2(1-\delta)(k-1)}\right)$
	$\displaystyle\leq M^{2}\left(\frac{2C_{1}^{2}\alpha^{2}\sigma^{2}}{B(1-\lambda% )^{3}}n+\frac{2C_{2}^{2}q_{1}}{(1-\delta)(1-\lambda)^{3}}\right).$

By applying Lemma 28 on $\{Y_{n,k}\}_{k=1}^{n}$ , with probability $1-1/T^{2}$ , we have

\displaystyle\left\|\sum_{k=1}^{n}Y_{n,k}\right\|\leq c(2\log T+4d)\frac{L_{f}% }{\sqrt{B}}M\left(\frac{2\sqrt{2}C_{1}\alpha\sigma}{\sqrt{B}(1-\lambda)^{3/2}}% \sqrt{n}+\frac{2\sqrt{2}C_{2}q_{1}^{1/2}}{(1-\delta)^{1/2}(1-\lambda)^{3/2}}% \right).

Lemma 25

Under (A1)-(A2) and (A3’), we have

\displaystyle\mathbb{P}\left(\left\|\alpha\sum_{k=1}^{n}\left((I-\Gamma)^{-1}-% (1-\gamma)\sum_{t=k}^{n}\sum_{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\right)\left(% \begin{array}[]{ccc}0\\ \triangledown g_{\eta_{k}}(x^{*})\\ \end{array}\right)\right\|\geq C\right)\leq\frac{1}{T^{2}},

where

\displaystyle C

\displaystyle=\alpha\frac{\sigma}{\sqrt{B}}M\frac{\sqrt{3}c(2\log T+4d)}{(1-% \lambda)^{3/2}},

and $c>0$ is an absolute constant.

Proof. Define

\displaystyle Y_{n,k}=\left((I-\Gamma)^{-1}-(1-\gamma)\sum_{t=k}^{n}\sum_{j=k}% ^{t}\Gamma^{t-j}\gamma^{j-k}\right)\left(\begin{array}[]{ccc}0\\ \triangledown g_{\eta_{k}}(x^{*})\\ \end{array}\right).

Then $\{Y_{n,k}\}_{k=1}^{n}$ are independent and sub-exponential. By Lemma 29 and the fact $\sqrt{\gamma}\leq\lambda<1$ , it holds that

\displaystyle\sum_{k=1}^{n}\left\|(I-\Gamma)^{-1}-(1-\gamma)\sum_{t=k}^{n}\sum% _{j=k}^{t}\Gamma^{t-j}\gamma^{j-k}\right\|^{2}

\displaystyle\leq M^{2}\left(\frac{2\gamma^{2}}{(1-\gamma)^{2}(1-\lambda)}+% \frac{\lambda^{2}}{(1-\lambda)^{3}}\right)\leq\frac{3M^{2}}{(1-\lambda)^{3}}.

By applying Lemma 28 on $\{Y_{n,k}\}_{k=1}^{n}$ , with probability $1-1/T^{2}$ , we have

\displaystyle\left\|\sum_{k=1}^{n}Y_{n,k}\right\|\leq\frac{\sigma}{\sqrt{B}}M% \frac{\sqrt{3}c(2\log T+4d)}{(1-\lambda)^{3/2}}.

Lemma 26

Under (A1)-(A2) and (A3’), suppose that the $t$ -th iteration $(\widetilde{m}_{t},\widetilde{x}_{t})$ satisfies

\displaystyle\sqrt{\left\|\widetilde{m}_{j}\right\|^{2}+\left\|\widetilde{x}_{% j}\right\|^{2}}\leq\frac{C_{1}}{\sqrt{B(1-\lambda)}}\alpha\sigma+C_{2}q_{1}^{1% /2}\lambda^{(1-\delta)(j-1)},

for $1\leq j\leq t$ and fixed $\delta\in(0,\frac{1}{2}]$ , then for $\overline{L}>0$ , we have

	$\displaystyle\left\\|\alpha(1-\gamma)\sum_{t=n_{0}+1}^{n}\sum_{j=1}^{t}\Gamma^{% t-j}\left(\begin{array}[]{ccc}0\\ -\sum_{k=1}^{j}\gamma^{j-k}V_{k}\widetilde{x}_{k}\\ \end{array}\right)\right\\|$
	$\displaystyle\leq\alpha M\overline{L}\left(\frac{2C_{1}^{2}\alpha^{2}\sigma^{2% }}{B(1-\lambda)^{2}}n+4C_{2}^{2}q_{1}\frac{\left(n_{0}(1-\lambda)+1\right)% \lambda^{n_{0}}}{(1-\lambda)^{2}}\right).$

Proof. By the fact that $\sum_{t=k}^{n}\sum_{j=k}^{t}\lambda^{t-j}\gamma^{j-k}\leq(1-\gamma)^{-1}(1-% \lambda)^{-1}$ , we have

	$\displaystyle\left\\|\alpha(1-\gamma)\sum_{t=n_{0}+1}^{n}\sum_{j=1}^{t}\Gamma^{% t-j}\left(\begin{array}[]{ccc}0\\ -\sum_{k=1}^{j}\gamma^{j-k}V_{k}\widetilde{x}_{k}\\ \end{array}\right)\right\\|$
	$\displaystyle\leq\alpha(1-\gamma)M\overline{L}\sum_{k=1}^{n}\left(\sum_{t=\max% \{n_{0}+1,k\}}^{n}\sum_{j=k}^{t}\lambda^{t-j}\gamma^{j-k}\right)\left(\frac{2C% _{1}^{2}}{B(1-\lambda)}\alpha^{2}\sigma^{2}+2C_{2}^{2}q_{1}\lambda^{2(1-\delta% )(k-1)}\right)$
	$\displaystyle\leq\alpha M\overline{L}\frac{2C_{1}^{2}\alpha^{2}\sigma^{2}}{B(1% -\lambda)^{2}}n+2\alpha(1-\gamma)M\overline{L}C_{2}^{2}q_{1}\sum_{k=1}^{n_{0}}% \left(\sum_{t=n_{0}+1}^{n}\sum_{j=k}^{t}\lambda^{t-j}\gamma^{j-k}\right)% \lambda^{2(1-\delta)(k-1)}$
	$\displaystyle+2\alpha(1-\gamma)M\overline{L}C_{2}^{2}q_{1}\sum_{k=n_{0}+1}^{n}% \left(\sum_{t=k}^{n}\sum_{j=k}^{t}\lambda^{t-j}\gamma^{j-k}\right)\lambda^{2(1% -\delta)(k-1)}.$

For the second term, due to that $\lambda\geq\sqrt{\gamma}$ and $\delta\in(0,\frac{1}{2}]$ , we have

	$\displaystyle\sum_{k=1}^{n_{0}}\left(\sum_{t=n_{0}+1}^{n}\sum_{j=k}^{t}\lambda% ^{t-j}\gamma^{j-k}\right)\lambda^{2(1-\delta)(k-1)}\leq\sum_{k=1}^{n_{0}}\left% (\lambda^{n_{0}+1-k}\sum_{t=n_{0}+1}^{n}\lambda^{t-n_{0}-1}\sum_{j=k}^{t}(% \gamma/\lambda)^{j-k}\right)\lambda^{2(1-\delta)(k-1)}$
	$\displaystyle\leq\sum_{k=1}^{n_{0}}\frac{2}{(1-\lambda)(1-\gamma)}\lambda^{n_{% 0}+1-k+2(1-\delta)(k-1)}\leq\frac{2n_{0}\lambda^{n_{0}}}{(1-\lambda)(1-\gamma)}.$

For the third term, we have

	$\displaystyle\sum_{k=n_{0}+1}^{n}\left(\sum_{t=k}^{n}\sum_{j=k}^{t}\lambda^{t-% j}\gamma^{j-k}\right)\lambda^{2(1-\delta)(k-1)}$		$\displaystyle\leq\sum_{k=n_{0}+1}^{n}\frac{1}{(1-\lambda)(1-\gamma)}\lambda^{2% (1-\delta)(k-1)}$
			$\displaystyle\leq\frac{2\lambda^{n_{0}}}{(1-\lambda)^{2}(1-\gamma)}.$

Thus we draw the conclusion.

Appendix C Several Useful Lemmas

Lemma 27

For $0\leq\lambda,\gamma<1$ , we have

\displaystyle\sum_{k=1}^{t}\left(\sum_{j=k}^{t}\lambda^{t-j}\gamma^{j-k}\right% )^{2}\leq\frac{2}{(1-\gamma)^{2}(1-\lambda)}.

In addition, if $\lambda\geq\sqrt{\gamma}$ , for fixed $\delta\in(0,1]$ , we have

\displaystyle\sum_{k=1}^{t}\left(\sum_{j=k}^{t}\lambda^{t-j}\gamma^{j-k}\right% )^{2}\lambda^{2(1-\delta)(k-1)}\leq\frac{4}{\delta(1-\gamma)^{2}(1-\lambda)}% \lambda^{2(1-\delta)(t-1)},

and

\displaystyle\sum_{k=1}^{t}\sum_{j=k}^{t}\lambda^{t-j}\gamma^{j-k}\lambda^{(1-% \delta)(k-1)}\leq\frac{2}{\delta(1-\gamma)(1-\lambda)}\lambda^{(1-\delta)(t-1)}.

Proof.

	$\displaystyle\sum_{k=1}^{t}\left(\sum_{j=k}^{t}\lambda^{t-j}\gamma^{j-k}\right% )^{2}$		$\displaystyle=\sum_{k=1}^{t}\sum_{i=k}^{t}\sum_{j=k}^{t}\lambda^{2t-i-j}\gamma% ^{i+j-2k}=\sum_{i=1}^{t}\sum_{j=1}^{t}\lambda^{2t-i-j}\sum_{k=1}^{\min\{i,j\}}% \gamma^{i+j-2k}$
			$\displaystyle\leq\frac{1}{1-\gamma}\sum_{i=1}^{t}\sum_{j=1}^{t}\lambda^{2t-i-j% }\gamma^{\left\|i-j\right\|}\leq\frac{2}{1-\gamma}\sum_{k=0}^{t-1}\sum_{j=1}^{t-% k}\lambda^{2t-2j-k}\gamma^{k}$
			$\displaystyle\leq\frac{2}{1-\gamma}\sum_{k=0}^{t-1}\frac{\lambda^{k}}{1-% \lambda^{2}}\gamma^{k}\leq\frac{2}{(1-\gamma)(1-\lambda)(1-\lambda\gamma)}\leq% \frac{2}{(1-\gamma)^{2}(1-\lambda)}.$

If $\lambda\geq\sqrt{\gamma}$ , we have

	$\displaystyle\lambda^{-2(1-\delta)(t-1)}\sum_{k=1}^{t}\left(\sum_{j=k}^{t}% \lambda^{t-j}\gamma^{j-k}\right)^{2}\lambda^{2(1-\delta)(k-1)}=\sum_{k=1}^{t}% \sum_{i=k}^{t}\sum_{j=k}^{t}\lambda^{2k-i-j}\gamma^{i+j-2k}\lambda^{2\delta(t-% k)}$
	$\displaystyle\leq\sum_{i=1}^{t}\sum_{j=1}^{t}\lambda^{2\delta(t-\min\{i,j\})}% \sum_{k=1}^{\min\{i,j\}}\left(\frac{\gamma}{\lambda}\right)^{i+j-2k}\leq\frac{% 1}{1-\gamma^{2}/\lambda^{2}}\sum_{i=1}^{t}\sum_{j=1}^{t}\lambda^{2\delta(t-% \min\{i,j\})}\left(\sqrt{\gamma}\right)^{\left\|i-j\right\|}$
	$\displaystyle\leq\frac{2}{1-\gamma}\sum_{i=1}^{t}\sum_{k=0}^{t-i}\lambda^{2% \delta(t-i)}\left(\sqrt{\gamma}\right)^{k}\leq\frac{4}{(1-\gamma)^{2}(1-% \lambda^{2\delta})}\leq\frac{4}{\delta(1-\gamma)^{2}(1-\lambda)},$

where the last inequality is due to the fact that $\frac{1-x}{1-x^{\delta}}\leq\frac{1}{\delta}$ for $\delta\in(0,1]$ . Similarly, we have

	$\displaystyle\lambda^{-(1-\delta)(t-1)}\sum_{k=1}^{t}\sum_{j=k}^{t}\lambda^{t-% j}\gamma^{j-k}\lambda^{(1-\delta)(k-1)}=\sum_{k=1}^{t}\sum_{j=k}^{t}\lambda^{% \delta(t-k)}\left(\frac{\gamma}{\lambda}\right)^{j-k}$
	$\displaystyle\leq\sum_{j=1}^{t}\lambda^{\delta(t-j)}\sum_{k=1}^{j}(\sqrt{% \gamma})^{j-k}\leq\frac{2}{1-\gamma}\frac{1}{1-\lambda^{\delta}}\leq\frac{2}{% \delta(1-\gamma)(1-\lambda)}.$

Lemma 28

Suppose $x_{k}\in\mathbb{R}^{d}$ are independent, zero-mean $\sigma_{k}$ -sub-exponential random vector, we have

\displaystyle\mathbb{P}\left(\left\|\sum_{k=1}^{t}x_{k}\right\|>c\left(\log(1/% \epsilon)+2d\right)\sqrt{\sum_{k=1}^{n}\sigma_{k}^{2}}\right)\leq\epsilon,

where $c>0$ is an absolute constant.

Proof. Let $\mathcal{N}$ be the 1/2-net of the unit ball $\{z\in\mathbb{R}^{d}:\left\|z\right\|\leq 1\}$ with respect to the Euclidean norm that satisfies $\left|\mathcal{N}\right|\leq 6^{d}$ . By the inequality that

\displaystyle\max_{\left\|v\right\|=1}v^{\top}x\leq\max_{z\in\mathcal{N}}z^{% \top}x+\max_{\left\|v\right\|=\frac{1}{2}}v^{\top}x=\max_{z\in\mathcal{N}}z^{% \top}x+\frac{1}{2}\max_{\left\|v\right\|=1}v^{\top}x,

we have

\displaystyle\mathbb{P}[\max_{\left\|v\right\|=1}v^{\top}\sum_{k=1}^{n}x_{k}>% \delta]\leq\mathbb{P}[2\max_{z\in\mathcal{N}}z^{\top}\sum_{k=1}^{n}x_{k}>% \delta]\leq 6^{d}\sup_{\left\|z\right\|=1}\mathbb{P}[2z^{\top}\sum_{k=1}^{n}x_% {k}>\delta].

According to Bernstein’s inequality, there holds

\displaystyle\mathbb{P}[2z^{\top}\sum_{k=1}^{n}x_{k}>\delta]\leq\exp\left(-c_{% 0}\min\left\{\frac{\delta^{2}}{4\sum_{k=1}^{n}\sigma_{k}^{2}},\frac{\delta}{2% \max_{k}\sigma_{k}}\right\}\right),

for some absolute constant $c_{0}$ . To bound $\mathbb{P}[\max_{\left\|v\right\|=1}v^{\top}\sum_{k=1}^{n}x_{k}>\delta]\leq\epsilon$ , we find $\delta$ such that

\displaystyle\delta

\displaystyle=2\max\left\{\sqrt{\sum_{k=1}^{n}\sigma_{k}^{2}}\frac{1}{\sqrt{c_% {0}}}\sqrt{\log(1/\epsilon)+2d},\max_{k}\sigma_{k}\frac{1}{c_{0}}\left(\log(1/% \epsilon)+2d\right)\right\}\leq c\sqrt{\sum_{k=1}^{n}\sigma_{k}^{2}}\left(\log% (1/\epsilon)+2d\right),

where $c=2\max\left\{\frac{1}{c_{0}},\frac{1}{\sqrt{c_{0}}}\right\}.$

Lemma 29

For $\Gamma$ satisfying the conditions in Lemma 20, we have

\displaystyle\left\|(I-\Gamma)^{-1}-(1-\gamma)\sum_{t=k}^{n}\sum_{j=k}^{t}% \Gamma^{t-j}\gamma^{j-k}\right\|

\displaystyle\leq M\left(\gamma\sum_{j=k}^{n}\lambda^{n-j}\gamma^{j-k}+\frac{% \lambda^{n-k+1}}{1-\lambda}\right),

and

\displaystyle\left\|(I-\Gamma)^{-1}-\sum_{t=0}^{n}\Gamma^{t}\right\|

\displaystyle\leq M\frac{\lambda^{n+1}}{1-\lambda}.

Proof. From Lemma 20, we have $\Gamma=\tilde{P}\Lambda\tilde{P}^{-1}$ . Then

	$\displaystyle\left\\|(I-\Gamma)^{-1}-(1-\gamma)\sum_{t=k}^{n}\sum_{j=k}^{t}% \Gamma^{t-j}\gamma^{j-k}\right\\|\leq\left\\|P\right\\|\left\\|P^{-1}\right\\|\left% \\|(I-\Lambda)^{-1}-(1-\gamma)\sum_{t=k}^{n}\sum_{j=k}^{t}\Lambda^{t-j}\gamma^{% j-k}\right\\|$
	$\displaystyle\leq M\left\\|\sum_{j=0}^{\infty}\Lambda^{j}-(1-\gamma)\sum_{j=0}^% {n-k}\Lambda^{j}\sum_{l=0}^{n-k-j}\gamma^{l}\right\\|=M\left\\|\sum_{j=0}^{n-k}% \Lambda^{j}(1-(1-\gamma)\sum_{l=0}^{n-k-j}\gamma^{l})+\sum_{j=n-k+1}^{\infty}% \Lambda^{j}\right\\|$
	$\displaystyle\leq M\left(\gamma\sum_{j=k}^{n}\lambda^{n-j}\gamma^{j-k}+\frac{% \lambda^{n-k+1}}{1-\lambda}\right).$

Also we have

\displaystyle\left\|(I-\Gamma)^{-1}-\sum_{t=0}^{n}\Gamma^{t}\right\|\leq M\sum% _{j=n+1}^{\infty}\lambda^{j}\leq M\frac{\lambda^{n+1}}{1-\lambda}.

	$\displaystyle\left\\|P\right\\|$		$\displaystyle\leq\sqrt{4}=2,$
	$\displaystyle\left\\|P^{-1}\right\\|$		$\displaystyle\leq 2\max_{k}\left\|\sqrt{\left\|\boldsymbol{z}^{+}_{k}\right\|^{2}% +1}\sqrt{\left\|\boldsymbol{z}^{-}_{k}\right\|^{2}+1}(\boldsymbol{z}^{+}-% \boldsymbol{z}^{-})_{k}^{-1}\right\|$
			$\displaystyle\leq 2\max_{k}\left(\left\|(\boldsymbol{z}^{+}\boldsymbol{z}^{-})_% {k}\right\|+\left\|\boldsymbol{z}^{+}_{k}\right\|+\left\|\boldsymbol{z}^{-}_{k}% \right\|+1\right)\left\|(\boldsymbol{z}^{+}-\boldsymbol{z}^{-})_{k}^{-1}\right\|$
			$\displaystyle\leq 2\max_{k}\left(\frac{(1-\gamma)\kappa_{k}}{\alpha\gamma}+% \frac{1}{\alpha\gamma}\left(2\alpha\gamma+(1-\gamma)(2+2\alpha\kappa_{k}+% \kappa_{k})\right)+1\right)\frac{\alpha\gamma}{\sqrt{\Delta}}$
			$\displaystyle\leq\frac{2}{\sqrt{\Delta}}\left(2(1-\gamma)(1+\alpha L+L)+3% \alpha\gamma\right),$

Acceleration of stochastic gradient descent with momentum by averaging: finite-sample rates and asymptotic normality

Abstract

1 Introduction

1.1 Related Works

2 Preliminaries

Example 1

3 Finite-sample Convergence Rates for SGDM

3.1 The finite-sample rates for SGDM on quadratic losses

Theorem 1

3.1.1 Linear convergence to a local neighborhood

Remark 2

Theorem 3

Remark 4

Remark 5

3.1.2 Explicit convergence rates for SGDM with specified γ𝛾\gammaitalic_γ and α𝛼\alphaitalic_α

Corollary 6 (Small momentum weight γ𝛾\gammaitalic_γ)

Corollary 7 (Large momentum weight γ𝛾\gammaitalic_γ)

Remark 8

Comparison to SGD.

Comparison to existing results for SGDM.

3.2 The finite-sample rates for SGDM on general losses

Theorem 9

Corollary 10 (Small momentum weight γ𝛾\gammaitalic_γ)

Corollary 11 (Large momentum weight γ𝛾\gammaitalic_γ)

Remark 12

4 Acceleration by Averaging and Asymptotic Normality

4.1 Averaged SGDM under quadratic losses

Theorem 13

Remark 14

Corollary 15

Remark 16

4.2 Averaged SGDM under general losses

Theorem 17

Corollary 18

Remark 19

5 Experiments

5.1 Simulation: quadratic loss

5.2 Sensitivity to learning rates

5.3 Simulation: logistic regression

5.4 Real data: MNIST classification

6 Conclusion

References

Appendix A Proofs of finite-time convergence rates of SGDM

A.1 Proof of Theorem 1

A.2 Proof of Theorem 3

Lemma 20

A.3 Proof of Theorem 9

A.4 Supporting Lemmas for Bounds in Theorem 9 Proof

Lemma 21

Lemma 22

Lemma 23

Appendix B Proofs of convergence rates of averaged SGDM

B.1 Proof of Theorem 13

B.2 Proof of Theorem 17

B.3 Supporting Lemmas for Bounds in Theorem 17 Proof

Lemma 24

Lemma 25

Lemma 26

Appendix C Several Useful Lemmas

Lemma 27

Lemma 28

Lemma 29

3.1.2 Explicit convergence rates for SGDM with specified $\gamma$ and $\alpha$

Corollary 6 (Small momentum weight $\gamma$ )

Corollary 7 (Large momentum weight $\gamma$ )

Corollary 10 (Small momentum weight $\gamma$ )

Corollary 11 (Large momentum weight $\gamma$ )