Accelerated Linear Convergence of Stochastic Momentum Methods in Wasserstein Distances

Introduction

Many key problems in machine learning can be formulated as convex optimization problems. Prominent examples in supervised learning include linear and non-linear regression problems, support vector machines, logistic regression or more generally risk minimization problems [Vap13]. Accelerated first-order optimization methods based on momentum averaging and their stochastic and proximal variants have been of significant interest in the machine learning community due to their scalability to large-scale problems and good performance in practice both in convex and non-convex settings, including deep learning (see e.g. [SMDH13, Nit14, HPK09, Xia10]).

Accelerated optimization methods for unconstrained problems based on momentum averaging techniques go back to Polyak who proposed the heavy ball (HB) method [Pol64] and are closely related to Tschebyshev acceleration, conjugate gradient and under-relaxation methods from numerical linear algebra [Var09, KV17]. Another popular momentum-based method is the Nesterov’s accelerated gradient (AG) method [Nes04]. For deterministic strongly convex problems, with access to the gradients of the objective, there is a well-established convergence theory for momentum methods. In particular, for minimizing strongly convex smooth objectives with Lipschitz gradients AG method requires $O(\sqrt{\kappa}\log(1/\varepsilon))$ iterations to find an $\varepsilon$-optimal solution where $\kappa$ is the condition number, this improves significantly over the $O(\kappa\log(1/\varepsilon))$ complexity of the gradient descent (GD) method. HB method also achieves a similar accelerated rate asymptotically in a local neighborhood around the global minimum. Also, for the special case of quadratic objectives, HB method can achieve the accelerated linear rate globally. In the absence of strong convexity, for convex functions, AG has an iteration complexity of $O(1/\sqrt{\varepsilon})$ in function values which accelerates the standard $O(1/\varepsilon)$ convergence rate of GD. In particular, it can be argued that AG method achieves an optimal convergence rate among all the methods that has access to only first-order information [Nes04]. For constrained problems, a variant of AG, the accelerated projected gradient (APG) method [OC15] can also achieve similar accelerated rates [Nes04, FRMP17].

On the other hand, in many applications, the true gradient of the objective function $\nabla f(x)$ is not available but we have access to a noisy but unbiased estimated gradient $\hat{\nabla}f(x)$ of the true gradient instead. The common choice of the noise that arises frequently in (stochastic oracle) models is the centered, statistically independent noise with a finite variance where for every $x\in\mathcal{X}$,

It is well recognized that momentum-based accelerated methods are quite sensitive to gradient noise [Har14, DGN14, FB15, DGN13], and need higher accuracy of the gradients to perform well [d’A08, DGN14] compared to standard methods like GD. In fact, with the standard choice of their stepsize and momentum parameter, numerical experiments show that they lose their superiority over a simple method like GD in the noisy setting [Har14], yet alone they can diverge [FB15]. On the other hand, numerical studies have also shown that carefully tuned constant stepsize and momentum parameters can lead to good practical performance for both HB and AG under noisy gradients in deep learning [SMDH13]. Overall, there has been a growing interest for obtaining convergence guarantees for stochastic momentum methods, i.e. momentum methods subject to noise in the gradients.

Several works provided sublinear convergence rates for stochastic momentum methods. [Lan12, GL12] developed the AC-SA method which is an adaptation of the AG method to the stochastic composite convex and strongly convex optimization problems and obtained an optimal $O(1/\sqrt{k})$ for the convex case. In a follow-up paper, [GL13] obtained an optimal $O(1/k)$ convergence bound for the constrained strongly convex optimization employing a domain shrinking procedure. However, these results do not apply to stochastic HB (SHB). [YLL16] provided a uniform analysis of SHB and accelerated stochastic gradient (ASG) showing $O(1/\sqrt{k})$ convergence rate for weakly convex stochastic optimization. [GPS18] obtained a number of sublinear convergence guarantees for SHB, showing that with decaying stepsize $\alpha_{k}=O(1/k^{\theta})$ for some $\theta\in(0,1]$, SHB method converges with rate $O(1/k^{\theta})$. Several other works focused on proper averaging for reducing the variance of the gradient error in the iterates for strongly convex linear regression problems [JKK+17, FB15, DFB17] and obtained a $O(1/k)$ convergence rate that achieves the minimax estimation rate. Recently, [LR17] studied the SHB algorithm for optimizing the least squares problems arising in the solution of consistent linear systems where the gradient noise comes from sampling the rows of the associated linear system and therefore the gradient errors have a multiplicative form vanishing at the optimum (see [LR17, Sec 2.5]), in which case SGD enjoys linear rates to the optimum with constant stepsize. The authors show that using a constant stepsize the expected SHB iterates converge linearly to a global minimizer with the accelerated rate and provide a first linear (but not an accelerated linear) rate for the expected suboptimality in function values, however the rate provided is not better than the linear rate of SGD and does not reflect the acceleration behavior compared to SGD. We note however that the results of this paper do not apply to our setting as our noise assumptions (H1)–(H2) are more general. In our setting, due to the persistence of the noise, it is not possible for the iterates of stochastic momentum methods converge to a global minimum, but rather converge to a stationary distribution around the global minimum. To our knowledge, a linear convergence result for momentum-based methods has never been established under this setting. For SGD, [DDB17] showed that when $f$ is strongly convex, the distribution of the SGD iterates with constant stepsize converges linearly to a unique stationary distribution $\pi_{\alpha}$ in the 2-Wasserstein distance requiring $O(\kappa\log(1/\varepsilon))$ iterations to be $\varepsilon$ close to the stationary distribution when $\alpha=1/L$ which is similar to the iteration complexity of (deterministic) gradient descent. A natural question is whether stochastic momentum methods admit a stationary distribution, if so whether the convergence to this distribution can happen faster compared to SGD. As the momentum methods are quite sensitive to gradient noise [Har14, CDO18] in terms of performance; a precise characterization of how much noise can be tolerated to achieve accelerated convergence rates under stochastic momentum methods remains understudied.

Contributions: We obtain a number of accelerated convergence guarantees for the SHB, ASG and accelerated stochastic projected gradient (ASPG) methods on both (weakly) convex and strongly convex smooth problems. We note that existing convergence bounds obtained for finite-sum problems that approximate stochastic optimization problems [Nit14] do not apply to our setting as our noise is more general, allowing us to deal directly with the stochastic optimization problem itself.

Third, we focus on the accelerated stochastic projected gradient (ASPG) algorithm for constrained stochastic strongly convex optimization on a bounded domain. We obtain fast accelerated convergence rate to a stationary distribution in the $p$-Wasserstein distance for any $p\geq 1$. Finally, we extend our results to the weakly convex setting where we show an accelerated $O(\frac{1}{\sqrt{\varepsilon}}\log(1/\varepsilon))$ convergence rate as long as the noise level is smaller than explicit bounds we provide. To our knowledge, accelerated rates in the presence of non-zero noise was not reported in the literature before.

Preliminaries

2 AG method

For $f\in\mathcal{S}_{\mu,L}$, the deterministic AG method consists of the iterations

and rewrite the AG iterations in terms of $\xi_{k}$. To simplify the presentation and the analysis, we build on the representation of optimization algorithms as a dynamical system from [HL17] and rewrite the AG iterations as

and $w_{k}:=\nabla f\left((1+\beta)x_{k}-\beta x_{k-1}\right)$. The standard analysis of deterministic AG is based on the following Lyapunov function that combines the state vector and function values:

In particular, Theorem 1 can recover existing convergence rate results for deterministic AG. For example, for the particular choice of

and $(\alpha,\beta)=(\alpha_{AG},\beta_{AG})$ with

in Theorem 1, we obtain the accelerated convergence rate of

However, as outlined in the introduction, in a variety of applications in machine learning and stochastic optimization, we do not have access to the true gradient $\nabla f(y_{k})$ as in the deterministic AG iterations but we have access to a (noisy) stochastic version $\hat{\nabla}f(y_{k})=\nabla f(y_{k})+\varepsilon_{k+1}$, where $\varepsilon_{k+1}$ is the random gradient noise. AG algorithm with stochastic gradients has the form

which is called the accelerated stochastic gradient (ASG) method (see e.g. [JKK+17]). We note that due to the existence of noise, the standard Lyapunov analysis from the literature (see e.g. [WRJ16, SBC14]) does not apply directly. We make the assumption that the random gradient errors are centered, statistically independent from the past iterates and have a finite second moment following the literature [CDO18, Har14, NVL+15, AFGO18, FB15]. The following assumption is a more formal statement of (H1)–(H2) adapting to the iterations $\xi_{k}$.

Under Assumption 2, the iterations $\xi_{k}$ forms a time-homogeneous Markov chain which we will study further in Sections 3 and 4.

3 HB method

For $f\in\mathcal{S}_{\mu,L}$, the HB method was proposed by [Pol64]. It consists of the iterations

where $\alpha>0$ is the step size and $\beta$ is the momentum parameter. The following asymptotic convergence rate result for HB is well known.

Then, $\|x_{k}-x_{*}\|\leq(\rho_{HB}+\delta_{k})^{k}\cdot\|\xi_{0}-\xi_{\ast}\|,$ where $\delta_{k}$ is a non-negative sequence that goes to zero and

Furthermore, $f(x_{k})-f(x_{*})\leq\frac{L}{2}(\rho_{HB}+\delta_{k})^{2k}\cdot\|\xi_{0}-\xi_{\ast}\|^{2}.$

This result has an asymptotic nature as the sequence $\delta_{k}$ is not explicit. There exist non-asymptotic linear convergence results for HB, but to our knowledge, known linear rate guarantees are slower than the accelerated rate $\rho_{HB}$; with a rate similar to the rate of gradient descent [GFJ14]. In Section 3.2, we will derive a new non-asymptotic version of this theorem that can guarantee suboptimality for finite $k$ with explicit constants and the accelerated rate $\rho_{HB}$. Note that the asymptotic rate $\rho_{HB}$ of HB in (17) on quadratic problems is strictly (smaller) faster than the rate $\rho_{AG}$ of AG from (13) in general (except in the particular special case of $\kappa=1$, we have $\rho_{AG}=\rho_{HB}=0$). However, for strongly convex functions, HB iterates given by (15) is not globally convergent with parameters $\alpha_{HB}$ and $\beta_{HB}$ [LRP16], but if the iterates are started in a small enough neighborhood around the global minimum of a strongly convex function, this rate can be achieved asymptotically [Pol87]. Since known guarantees for deterministic AG is stronger than deterministic HB on non-quadratic strongly convex functions, we will focus on the AG method for non-quadratic objectives in our paper.

We will analyze the HB method under noisy gradients:

where the noise satisfies Assumption 2. This method is called the stochastic HB method [GPS18, LR18, Flå04].

In the next section, we show that stochastic momentum methods admit an invariant distribution towards which they converge linearly in a sense we make precise. For illustrative purposes, we first analyze the special case when the objective is a quadratic function, and then move on to the more general case when $f$ is smooth and strongly convex. Also, for quadratic functions we can obtain stronger guarantees exploiting the linearity properties of the gradients.

Special case: strongly convex quadratics

First, we assume that the objective $f\in\mathcal{S}_{\mu,L}$ and is a quadratic function of the form

where $\mathcal{W}_{2,S_{\alpha,\beta}}$ is the $2$-Wasserstein distance (1) equipped with the $\|\cdot\|_{S_{\alpha,\beta}}$ norm. In particular, with $(\alpha,\beta)=(\alpha_{AG},\beta_{AG})$ and $P=P_{AG}$ with $P_{AG}$ defined in (11), we obtain the optimal accelerated linear rate of convergence:

with $\rho_{AG}=1-\frac{1}{\sqrt{\kappa}}$ as in (13).

For the AG method, the choice of $(\alpha,\beta)=(\alpha_{AG},\beta_{AG})$ is popular in practice, however a faster rate can be achieved asymptotically if

so that the asymptotic linear convergence rate in distance to the optimality becomes $\rho_{AG}^{*}:=1-\frac{2}{\sqrt{3\kappa+1}}$, which translates into the rate $(\rho_{AG}^{*})^{2}$ in function values that is (smaller) faster than $\rho_{AG}$ [LRP16]; improving the iteration complexity by a factor of $4/\sqrt{3}\approx 2.3$ when $\kappa$ is large. However, these results are asymptotic. Below we provide a first non-asymptotic bound with the faster rate $\rho_{AG}^{*}$.

where $\rho_{AG}^{*}=1-\frac{2}{\sqrt{3\kappa+1}}$ and

where $\{\lambda_{i}\}_{i=1}^{d}$ are the eigenvalues of the Hessian $Q$.

The constants $C_{k}^{*}$ grows linearly with $k$ in Theorem 5 and this dependency is tight in the sense that there are examples achieving it (see the proof in the supplementary file). Our bounds improves the existing results that provide a slower rate $\rho_{AG}$ with bounded constants in front of the linear rate [Nes04, Bub14], if $k$ is large enough (larger than a constant that can be made explicit).

Building on this non-asymptotic convergence result for the deterministic AG method, we obtain similar non-asymptotic convergence guarantees for the ASG method in $p$-Wasserstein distances towards convergence to a stationary distribution.

where $\rho_{AG}^{*}=1-\frac{2}{\sqrt{3\kappa+1}}$, $C_{k}^{*}$ is defined in (25) and $\mathcal{W}_{p}$ is the standard the $p$-Wasserstein distance.

With the same assumptions as in Theorem 7,

where $\rho_{AG}^{*}=1-\frac{2}{\sqrt{3\kappa+1}}$, $C_{k}^{*}$ is defined in (25), $X_{AG}^{*}$ is the covariance matrix of $\xi_{\infty}-\xi_{\ast}$ and $V_{AG}^{*}(\xi_{0})$ is a constant depending on any initial state $\xi_{0}$ and both $X$ and $V_{AG}^{*}(\xi_{0})$ will be spelled out in explicit form in the supplementary file.

2 Accelerated linear convergence of HB and SHB

We first give a non-asymptotic convergence result for the deterministic HB method with explicit constants, which also implies a bound on the suboptimality $f(x_{k})-f(x_{\ast})$. This refines the asymptotic results in the literature (Theorem 3).

with $\bar{C}:=\max_{i:\mu<\lambda_{i}<L}\frac{\mu+L}{2\sqrt{(\lambda_{i}-\mu)(L-\lambda_{i})}}$, where $\{\lambda_{i}\}_{i=1}^{d}$ are the eigenvalues of the Hessian matrix of $f$.

It is clear from the definition of $C_{k}$ in Theorem 9 that the leading coefficient $C_{k}$ grows at most linearly in the number of iterates $k$ and this dependency cannot be removed in the sense that there are some examples achieving our upper bounds in terms of $k$ dependency (see the supplementary file).

Building on this non-asymptotic convergence result for the deterministic HB method, we obtain similar non-asymptotic convergence guarantees for the SHB method in Wasserstein distances towards convergence to a stationary distribution.

where $\rho_{HB}=1-\frac{2}{\sqrt{k}+1}$ as defined in (17), $C_{k}$ is defined in (29) and $\mathcal{W}_{p}$ is the standard the $p$-Wasserstein distance.

With the same assumptions as in Theorem 11,

where $\rho_{HB}=1-\frac{2}{\sqrt{\kappa}+1}$ as in (17), $C_{k}$ is defined in (29), $X_{HB}$ is the covariance matrix of $\xi_{\infty}-\xi_{\ast}$, $V_{HB}(\xi_{0})$ is a constant depending on any initial state $\xi_{0}$ and both $X$ and $V_{HB}(\xi_{0})$ will be spelled out in explicit form in the supplementary file.

Strongly convex smooth optimization

for some explicit constant $c_{0}$ (to be given in the supplementary file), where $\mathcal{W}_{1}$ is the standard $1$-Wasserstein distance.

Given any $\eta\in(0,1)$ and $M>0$ so that $\int_{\|x-x_{\ast}\|\leq M}p(\xi_{\ast},x)dx\geq\sqrt{\eta},$ and any $R>0$ so that

Then there is a unique stationary distribution $\pi_{\alpha,\beta}$ so that

where $\mathcal{W}_{1}$ is the standard $1$-Wasserstein distance and $\psi:=\frac{\eta}{2K_{\alpha,\beta}}$ and

Next, we obtain the optimal convergence rate and provide a bound on the expected suboptimality by choosing $(\alpha,\beta)=(\alpha_{AG},\beta_{AG})$.

Given $(\alpha,\beta)=(\alpha_{AG},\beta_{AG})$. Define $M$ and $R$ as in Theorem 13 with $\eta=1/\kappa^{1/2}$. Also assume that the noise has small variance, i.e. $\sigma^{2}\leq RL/(4\sqrt{\kappa}).$ Then, with $\psi:=\frac{L}{2\sqrt{\kappa}\sigma^{2}}$, we have

where $\mathcal{W}_{1}$ is the standard $1$-Wasserstein distance and for any initial state $\xi_{0}$,

The bound (33) is similar in spirit to Corollary 4.7. in [AFGO18] but with a different assumption on noise. We can see that the expected value of the objective with respect to the $k$-th iterate is close to the true minimum of the objective if $k$ is large, and the variance of the noise $\sigma^{2}$ is small. In the special case when the noise are i.i.d. Gaussian, one can compute the constants in closed-form.

If the noise $\varepsilon_{k}$ are i.i.d. Gaussian $\mathcal{N}(0,\Sigma)$, where $\Sigma\prec L^{2}I_{d}$. Then, Proposition 14 holds with

If we take $\mu=\Theta(1)$, then $L=\Theta(\kappa)$ and it follows that we have $M=O(\kappa^{-1/8})$ and $R=O\left(\kappa^{-13/4}\log^{2}(\kappa)\right)$.

We note that Proposition 14 and Corollary 15 provide explicit bounds on the admissable noise level $\sigma^{2}$ to ensure accelerated convergence with respect to Wasserstein distances and expected suboptimality after $k$ iterations.

ASPG and the weakly convex setting

Weakly convex functions. If the objective is (weakly) convex but not strongly convex and the constraint set is bounded, our analysis for the strongly convex case can be adapted with minor modifications. Following standard regularization techniques (see e.g. [LRP16, Bub14]), that allow to approximate a weakly convex function with a strongly convex function, we provide explicit bounds on the noise level to obtain the accelerated $O(\varepsilon^{-1/2})$ rate up to a log factor on $\varepsilon$ in expected suboptimality in function values (see the supplementary file).

Conclusion

We have studied accelerated convergence guarantees for a number of stochastic momentum methods (SHB, ASG, ASPG) for strongly and (weakly) convex smooth problems. First, we studied the special case when the objective is quadratic and the gradient noise is additive and i.i.d. with a finite second moment. Non-asymptotic guarantees for accelerated linear convergence are obtained for the deterministic and stochastic AG and HB methods for any $p$-Wasserstein distance ($p\geq 1$), and also for the ASG method in the weighted $2$-Wasserstein distance, which builds on the dissipativity theory from the deterministic setting. Our analysis for HB and AG also leads to improved non-asymptotic convergence bounds in suboptimality after $k$ iterations for both deterministic and stochastic settings which is of independent interest. Second, we studied the (non-quadratic) strongly convex optimization under the stochastic oracle model (H1)–(H2). Accelerated linear convergence rate is obtained for the ASG method in the $1$-Wasserstein distance. Third, we studied the ASPG method for constrained stochastic strongly convex optimization on a bounded domain. Accelerated linear convergence rate is obtained in any $p$-Wasserstein distance ($p\geq 1$), and extension to the (weakly) convex setting will be discussed in the supplementary file. Our results provide performance bounds for stochastic momentum methods in expected suboptimality and in Wasserstein distances. Finally, the proofs of all the results in our paper will be given in the supplementary file.

Acknowledgements

Mert Gürbüzbalaban and Bugra Can acknowledge support from the grants NSF DMS-1723085 and NSF CCF-1814888. Lingjiong Zhu is grateful to the support from the grant NSF DMS-1613164.

References

Appendix A Constrained Optimization and ASPG

where $\varepsilon_{k}$ is the random gradient error satisfying Assumption 2, $\alpha,\beta>0$ are the stepsize and momentum parameter and $\mathcal{P}_{\mathcal{C}}(x)$ denotes the projection of a point $x$ to the compact set $\mathcal{C}$. For constrained problems, algorithms based on projection steps that restricts the iterates to the constraint set are more natural compared to the standard AG algorithm primarily designed for the unconstrained optimization [Bub14]. Accelerated projected gradient methods can also be viewed as a special case of the accelerated proximal gradient methods as the proximal operator reduces to a projection in a special case (see e.g. [PB+14]).

We will show in Proposition 28 that the metric $d_{\psi}$ implies the standard $p$-Wasserstein metric in the sense that for any two probability measures $\mu_{1},\mu_{2}$ on the product space $\mathcal{C}^{2}:=\mathcal{C}\times\mathcal{C}$,

where $\mathcal{D}_{\mathcal{C}^{2}}=\sqrt{2}D_{C}$ is the diameter of $\mathcal{C}^{2}$.

Given any $\eta\in(0,1)$ and $R>0$ so that

where $\mathcal{W}_{p}$ is the standard $p$-Wasserstein metric ($p\geq 1$) and

We can see from (37) that the expected value of the objective with respect to the $k$-th iterate is close to the true minimum of the objective if $k$ is large, and the stepsize $\alpha$ or the variance of the noise $\sigma^{2}$ is small. By choosing $(\alpha,\beta)=(\alpha_{AG},\beta_{AG})$, we obtain the optimal convergence in the next theorem.

Given $(\alpha,\beta)=(\alpha_{AG},\beta_{AG})$. Define $R$ as in Theorem 16 with $\eta=1/\kappa^{1/2}$. Also assume that the noise has small variance, i.e.

where $a_{1}:=\frac{1}{L^{2}}\left(\frac{\mu}{2}((1-\sqrt{\kappa})^{2}+\kappa)+\frac{L}{2}\right)$ and $b_{1}:=\frac{1}{L}\left(\mathcal{D}_{\mathcal{C}}\mu((1-\sqrt{\kappa})^{2}+\kappa)+G_{M}\right)$. Then, we have

where $\mathcal{W}_{p}$ is the standard $p$-Wasserstein metric ($p\geq 1$) and

Appendix B Weakly Convex Constrained Optimization

In this section, we extend the constrained optimization for the accelerated stochastic projected gradient method (ASPG) from the strongly convex objectives studied in Section A to the (weakly) convex objectives.

This shows that if the noise is small is enough, it suffices to have

many iterations to sample an $\varepsilon$-optimal point in expectation.

Appendix C Proofs of Results in Section 3

In this section, we prove the results for Section 3, in which the objective is quadratic: $f(x)=\frac{1}{2}x^{T}Qx+a^{T}x+b$ and $f\in\mathcal{S}_{\mu,L}$, which satisfies the inequalities:

Before we proceed to the proofs of the results in Section 3.1, we first show that the matrix $S_{\alpha,\beta}$ defined in (21) is positive definite so that the weighted 2-Wasserstein metric $\mathcal{W}_{2,S_{\alpha,\beta}}$ given in (1) is well-defined.

Next, before we proceed to the proofs of the results in Section 3.1, let us first recall that throughout Section 3, the noise $\varepsilon_{k}$ are assumed to be i.i.d. Let us define the coupling

Before we proceed, let us recall the following lemma from [HL17].

The proof of Theorem 4 relies on the following lemma.

with $j=1,2$. Assume that $f$ is quadratic and $f(x)=\frac{1}{2}x^{T}Qx+a^{T}x+b$, where $Q$ is positive definite.

Let $\rho=\rho_{\alpha,\beta}\in(0,1)$ that can depend on $\alpha$ and $\beta$ so that there exists some $P=P_{\alpha,\beta}$ symmetric and positive semi-definite that can depend on $\alpha$ and $\beta$ such that

Let us first consider the simpler case $f(x)=\frac{1}{2}x^{T}Qx$. Since $f$ is quadratic, $\nabla f$ is linear. Applying (52) and the linearity of $\nabla f$, we get

Applying (53) and the linearity of $\nabla f$, we get

By Lemma 19 and the definition of $\rho_{\alpha,\beta}$, $P_{\alpha,\beta}$ the inequality (51) holds. Thus

Since $f$ is quadratic, and we assumed that $f(x)=\frac{1}{2}x^{T}Qx$, where $Q$ is positive definite, we get

Previously, we assumed $f(x)=\frac{1}{2}x^{T}Qx$, so that $\nabla f(x-y)=\nabla f(x)-\nabla f(y)$. In general, the quadratic function takes the form

Hence, by Lemma 19 and the definition of $\rho_{\alpha,\beta}$, $P_{\alpha,\beta}$ so that (51) holds, we get the same result as before:

By taking $\alpha=\alpha_{AG}$, $\beta=\beta_{AG}$, $\rho=\rho_{AG}$ and $P_{AG}$ in definition (11), we recall the following result from [HL17].

We immediately obtain the following result.

Assume the coupling (49)-(50). Assume that $f$ is quadratic and $f(x)=\frac{1}{2}x^{T}Qx+a^{T}x+b$, where $Q$ is positive definite. Then, we have

Now, we are ready to state the proof of Theorem 4.

Recall the iterates $\xi_{k}=(x_{k}^{T},x_{k-1}^{T})^{T}$, the Markov kernel $\mathcal{P}_{\alpha,\beta}$ and the definition of the weighted $2$-Wasserstein distance (1) with the weighted norm (20)-(21) and $P=P_{\alpha,\beta}$. Then showing Theorem 4 is equivalent to show

Let $(((x_{k}^{(i)})^{T},(x_{k-1}^{(i)})^{T})^{T})_{k=0}^{\infty}$, $i=1,2$ be a coupling of $((x_{k}^{T},x_{k-1}^{T})^{T})_{k=0}^{\infty}$ defined as before. We have shown before that for every $k$,

By taking expectation and since $\frac{1}{2}x^{T}Qx\geq 0$ for any $x$, we get

By taking $\lambda_{2}=\mathcal{P}_{\alpha,\beta}\lambda_{1}$, we get

Hence $\mathcal{P}_{\alpha,\beta}^{k}\lambda_{1}$ is a Cauchy sequence and converges to a limit $\pi_{\alpha,\beta}^{\lambda_{1}}$:

Next, let us show that $\pi_{\alpha,\beta}^{\lambda_{1}}$ does not depend on $\lambda_{1}$. Assume that there exists $\pi_{\alpha,\beta}^{\lambda_{2}}$ so that $\lim_{k\rightarrow\infty}\mathcal{W}_{2,S_{\alpha,\beta}}(\mathcal{P}_{\alpha,\beta}^{k}\lambda_{2},\pi_{\alpha,\beta}^{\lambda_{2}})=0$. Since $\mathcal{W}_{2,S_{\alpha,\beta}}$ is a metric, by the triangle inequality,

which goes to zero as $k\rightarrow\infty$. Hence, $\pi_{\alpha,\beta}^{\lambda_{1}}=\pi_{\alpha,\beta}^{\lambda_{2}}$. The limit is therefore the same for any initial distributions and we can denote it by $\pi_{\alpha,\beta}$. Indeed,

which goes to zero as $k\rightarrow\infty$. Hence $\mathcal{P}_{\alpha,\beta}\pi_{\alpha,\beta}=\pi_{\alpha,\beta}$ gives the invariant distribution. We can also show similarly as before that it is unique. ∎

If $\alpha\in(0,1/L]$ and $\beta=\frac{1-\sqrt{\alpha\mu}}{1+\sqrt{\alpha\mu}}$, then we can take the matrix $P_{\alpha,\beta}$ appearing in Theorem 4 according to the $P_{\alpha}$ matrix defined in [AFGO19, Theorem 2.3] to obtain $\rho(\alpha,\beta)=1-\sqrt{\alpha\mu}$. For $\alpha=\frac{\log^{2}(k)}{\mu k^{2}}$, then this leads to $\mathcal{W}_{2,S_{\alpha,\beta}}\left(\nu_{k,\alpha,\beta},\pi_{\alpha,\beta}\right)\leq\frac{1}{k}\mathcal{W}_{2,S_{\alpha,\beta}}(\nu_{0,\alpha,\beta},\pi_{\alpha,\beta})$ and it can be shown with an analysis similar to that of [AFGO19] that the second moment of $\pi_{\alpha,\beta}$ is also $O(1/k)$; ignoring some logarithmic factors in $k$. Therefore, our results do not violate (and are in agreement with) the $\Omega(1/k)$ lower bounds studied in [CDLZ16, RR11, AWBR09] for strongly convex stochastic optimization.

where $\alpha>0$ is the step size and $\beta$ is the momentum parameter. In the case when $f$ is quadratic and $f(x)=\frac{1}{2}x^{T}Qx+a^{T}x+b$, we can compute that

and we aim to provide an upper bound to the 2-norm of the matrix, that is:

Let us assume that $Q$ has the decomposition

where $D$ is diagonal consisting of eigenvalues $\lambda_{i}$, $1\leq i\leq d$ in increasing order:

which has the same eigenvalues as the matrix:

are $2\times 2$ matrices with eigenvalues:

Next, we upper bound $\|T_{i}^{k}\|$. We recall the choice:

Therefore $\Delta_{i}=0$ if and only if $\lambda_{i}=\mu$ or $\lambda_{i}=\frac{3L+\mu}{4}$, and moreover $\Delta_{i}<0$ for $\mu<\lambda_{i}<\frac{3L+\mu}{4}$ and $\Delta_{i}>0$ for $\lambda_{i}>\frac{3L+\mu}{4}$.

(1) Consider the case $\mu<\lambda_{i}<\frac{3L+\mu}{4}$. Then $\Delta_{i}<0$. It is known that the $k$-th power of a $2\times 2$ matrix $A$ with distinct eigenvalues $\mu_{\pm}$ is given by

where $I$ is the $2\times 2$ identity matrix [Wil92]. In our context, $A=T_{i}$ and $\mu_{\pm}=\mu_{i,\pm}$, we get

Hence, it follows from (63), (66), (67), (68) and (69) that

(2) Consider the case $\frac{3L+\mu}{4}<\lambda_{i}<L$. Then, $\Delta_{i}>0$. As before, we have

Hence, it follows from (70), (71), (72), (73) and (74) that

(3) Consider the case $\lambda_{i}=\mu$. Then $\Delta_{i}=0$. It is known that the $k$-th power of a $2\times 2$ matrix $A$ with two equal eigenvalues $\mu_{+}=\mu_{-}=\mu$ is given by

where $I$ is the $2\times 2$ identity matrix [Wil92]. In our context, $A=T_{i}$ and

Therefore, with $\lambda_{i}=\mu$, we have

Furthermore, we see that the sequence $T_{i}^{k}/k$ converges to a non-zero matrix. Therefore, $\|T_{i}^{k}\|\geq ck$ for some constant $c$ for every $k$. This means that the linear dependency to $k$ of our upper bound in (78) is tight. This behavior is expected due to the fact that $T_{i}^{k}$ has double roots.

(4) Consider the case $\lambda_{i}=\frac{3L+\mu}{4}$. Then $\Delta_{i}=0$. We can compute that

Finally, combining the three cases (1) $\mu<\lambda_{i}<\frac{3L+\mu}{4}$; (2) $\lambda_{i}>\frac{3L+\mu}{4}$; (3) $\lambda_{i}=\mu$; (4) $\lambda_{i}=\frac{3L+\mu}{4}$, and recall (60), we get

where $\alpha>0$ is the step size and $\beta$ is the momentum parameter. In the case when $f$ is quadratic and $f(x)=\frac{1}{2}x^{T}Qx+a^{T}x+b$, we can compute that

so that with two couplings $x_{k}^{(1)},x_{k}^{(2)}$:

Following from the proof of Theorem 4, we can show by constructing a Cauchy sequence that there exists a unique stationary distribution $\pi_{\alpha,\beta}$. Finally, we assume that $(x_{0}^{(1)},x_{-1}^{(1)})$ starts from the given $(x_{0},x_{-1})$ distributed as $\nu_{0,\alpha,\beta}$ and $(x_{0}^{(2)},x_{-1}^{(2)})$ starts from the stationary distribution $\pi_{\alpha,\beta}$ so that their $L_{p}$ distance is exactly the $\mathcal{W}_{p}$ distance. Then we get

and the proof is complete by taking the power $1/p$ in the above equation. ∎

Before we state the proof of Theorem 8, let us spell out $X$ and $V_{AG}^{*}(\xi_{0})$ in the statement of Theorem 8 explicitly here. We will show that Theorem 8 holds with $V_{AG}^{*}(\xi_{0})$ given by

In the special case $\Sigma=c^{2}I_{d}$ for some constant $c\geq 0$, it follows from [AFGO18] that

where $\{\lambda_{i}\}_{i=1}^{d}$ are the eigenvalues of $Q$.

where we consider the quadratic objective $f(x)=\frac{1}{2}x^{T}Qx+a^{T}x+b$ so that

satisfies the discrete Lyapunov equation:

Next by iterating equation (81) over $k$, we immediately obtain

where we used the estimate $\|(A_{Q}^{*})^{k}\|\leq C_{k}^{*}(\rho_{AG}^{*})^{k}$ from the proof of Theorem 5.

Finally, since $\nabla f$ is $L$-Lipschtiz,

Note that our results in $p$-Wasserstein distances would hold if there exists some $p\geq 1$ so that $p$-th moment of the noise is finite. For instance, the $p<2$ case can arise in applications where the noise has heavy tail (see e.g. [SSG19]).

C.2 Proofs of Results in Section 3.2

where $\alpha>0$ is the step size and $\beta$ is the momentum parameter. In the case when $f$ is quadratic and $f(x)=\frac{1}{2}x^{T}Qx+a^{T}x+b$, we can compute that

and we aim to provide an upper bound to the 2-norm of the matrix, that is:

Let us assume that $Q$ has the decomposition

where $D$ is diagonal consisting of eigenvalues $\lambda_{i}$, $1\leq i\leq d$ in increasing order:

which has the same eigenvalues as the matrix:

are $2\times 2$ matrices with eigenvalues:

Next, we upper bound $\|T_{i}^{k}\|$. We consider three cases (1) $\mu<\lambda_{i}<L$; (2) $\lambda_{i}=\mu$; (3) $\lambda_{i}=L$.

(1) Consider the case $\mu<\lambda_{i}<L$. With the choice of $\alpha$ and $\beta$ in (16), we can compute that for those $\mu<\lambda_{i}<L$, we have

where $1\leq i\leq d$. It is known that the $k$-th power of a $2\times 2$ matrix $A$ with distinct eigenvalues $\mu_{\pm}$ is given by

where $I$ is the $2\times 2$ identity matrix [Wil92]. In our context, $A=T_{i}$ and $\mu_{\pm}=\mu_{i,\pm}$, we get

Hence, it follows from (83), (84), (85), (86) and (87) that

(2) Consider the case $\lambda_{i}=\mu$. With the choice of $\alpha$ and $\beta$ in (16), we can compute that for those $\lambda_{i}=\mu$, we have

so we have double eigenvalues and indeed $1+\beta-\alpha\lambda_{i}=2\sqrt{\beta}$, and

and by a direct computation (e.g. induction on $k$), we get:

Finally, we note that the matrix $T_{i}^{k}/(\sqrt{\beta}^{k}k)$ as $k$ goes to infinity converges to the $2\times 2$ matrix

Therefore, the linear dependency of our bound in (90) with respect to $k$ is tight. This behavior is expected due to the fact that $T_{i}^{k}$ has double roots.

(3) Consider the case $\lambda_{i}=L$. With the choice of $\alpha$ and $\beta$ in (16), we can compute that for those $\lambda_{i}=L$, we have

so we have double eigenvalues and indeed $1+\beta-\alpha\lambda_{i}=-2\sqrt{\beta}$, and

and by a direct computation (e.g. induction on $k$), we get:

Finally, combining the three cases (1) $\mu<\lambda_{i}<L$; (2) $\lambda_{i}=\mu$; (3) $\lambda_{i}=L$, we get

and the proof is complete by applying (94). ∎

Before we state the proof of Theorem 11, let us state the following result, which is built on Theorem 9.

Let us consider two couplings $(x_{k}^{(1)})_{k\geq 0}$ and $(x_{k}^{(2)})_{k\geq 0}$ with the common noise $(\varepsilon_{k+1})_{k\geq 0}$ that starts from $x_{0}^{(1)}$ and $x_{0}^{(2)}$:

where $f$ is quadratic and $f(x)=\frac{1}{2}x^{T}Qx+a^{T}x+b$. Then, we have

where $\rho_{HB}$ and $C_{k}$ are defined by (17) and (29) respectively.

It follows from the estimate (94) in the proof of Theorem 9 and the definitions of $\rho_{HB}$ and $C_{k}$ in (17) and (29) that we have

We recall from Lemma 25 that for any coupling $x^{(1)}$ and $x^{(2)}$

Following from the proof of Theorem 4, we can show by constructing a Cauchy sequence that there exists a unique stationary distribution $\pi_{\alpha,\beta}$. Finally, we assume that $(x_{0}^{(1)},x_{-1}^{(1)})$ starts from the given $(x_{0},x_{-1})$ distributed as $\nu_{0,\alpha,\beta}$ and $(x_{0}^{(2)},x_{-1}^{(2)})$ starts from the stationary distribution $\pi_{\alpha,\beta}$ so that their $L_{p}$ distance is exactly the $\mathcal{W}_{p}$ distance. Then we get

and the proof is complete by taking the power $1/p$ in the above equation. ∎

Before we state the proof of Theorem 12, let us spell out $X$ and $V_{HB}(\xi_{0})$ in the statement of Theorem 12 explicitly here. We will show that Theorem 12 holds with $V_{HB}(\xi_{0})$ given by

In the special case $\Sigma=c^{2}I_{d}$ for some constant $c\geq 0$, we obtain

where $\{\lambda_{i}\}_{i=1}^{d}$ are the eigenvalues of $Q$.

where we consider the quadratic objective $f(x)=\frac{1}{2}x^{T}Qx+a^{T}x+b$ so that

satisfies the discrete Lyapunov equation:

Next by iterating equation (103) over $k$, we immediately obtain

where we used the estimate $\|A_{Q}^{k}\|\leq C_{k}\rho_{HB}^{k}$ from the proof of Theorem 9.

Finally, since $\nabla f$ is $L$-Lipschtiz,

The proof of (31) is complete. To show (102), we can adapt the proof technique of [AFGO18, Proposition 3.2] for gradient descent to HB. Without loss of generality, due to the scaling of the Lyapunov equation, we can assume $c=1$. Consider the eigenvalue decomposition $A_{Q}=V\Lambda V^{T}$ where $Q$ is orthogonal and $\Lambda$ is diagonal with $\Lambda(i,i)=\lambda_{i}$. We can write

If we define $Y:=UXU^{-1}$ for the orthogonal matrix $U=P\bar{V}^{T}$, it solves

where the latter matrix $S$ is a $2d\times 2d$ diagonal matrix with entries $S(i,i)=\alpha^{2}$ if $i$ is odd, and zero if $i$ is even. Due to the special structure of $S$ and $A_{M}$, the solution $Y$ has the structure

where $Y_{i}$ solves the $2\times 2$ Lyapunov equation

with scalars $x_{i}$, $y_{i}$ and $w_{i}$, this equation is equivalent to the linear system

Appendix D Proofs of Results in Section 4

Before we proceed to prove the main results in Section 4, let us first show that the weighted total variation distance $d_{\psi}$ upper bounds the standard $1$-Wasserstein distance.

where $\mathcal{W}_{1}$ is the standard $1$-Wasserstein distance and

where $\hat{c}_{0}$ is the smallest positive eigenvalue of

By applying the Kantorovich-Rubinstein duality for the Wasserstein metric (see e.g. [Vil09]), we get

where we used $1+\psi V_{P}(\xi)\geq c_{0}\|\xi-\xi_{\ast}\|$ from Lemma 27. ∎

Let $\xi^{T}=(x^{T},y^{T})$. If $\|\xi-\xi_{\ast}\|\leq 1$, then $c_{0}=1$ works. Otherwise,

For constrained optimization on a compact set $\mathcal{C}$, we have the following result.

For any $\mu_{1},\mu_{2}$ on the product space $\mathcal{C}^{2}:=\mathcal{C}\times\mathcal{C}$,

where $\mathcal{D}_{\mathcal{C}^{2}}$ is the diameter of $\mathcal{C}^{2}$.

The second inequality in Proposition 28 follows from $d_{\psi}(\mu_{1},\mu_{2})\geq 2\|\mu_{1}-\mu_{2}\|_{TV}$. So it suffices to prove the first inequality. We can compute that

Throughout Section 4, the noise $\varepsilon_{k}$ are assumed to satisfy Assumption 2. Our proof of Theorem 13 relies on the geometric ergodicity and convergence theory of Markov chains. Geometric ergodicity and convergence of Markov chains has been well studied in the literature. Harris’ ergodic theorem of Markov chains essentially states that a Markov chain is ergodic if it admits a small set that is visited infinitely often [Har56]. Such a result often relies on finding an appropriate Lyapunov function [MT93]. The transition probabilities converge exponentially fast towards the unique invariant measure, and the prefactor is controlled by the Lyapunov function [MT93]. Computable bounds for geometric convergence rates of Markov chains has been obtained in e.g. [MT94, HM11]. In the following, we state the results from [HM11]. Before we proceed, let us introduce some definitions and notations.

There exists some constant $\eta\in(0,1)$ and a probability measure $\nu$ so that

Let us recall the definition of the weighted total variation distance:

It is noted in [HM11] that $d_{\psi}$ has the following alternative expression. Define the weighted supremum norm for any $\psi>0$:

and its associated dual metric $d_{\psi}$ on probability measures:

It is also noted in [HM11] that $d_{\psi}$ can also be expressed as:

If the drift condition (Assumption 29) and minorization condition (Assumption 30) hold, then there exists $\bar{\eta}\in(0,1)$ and $\psi>0$ so that

If the drift condition (Assumption 29) and minorization condition (Assumption 30) hold, then $\mathcal{P}$ admits a unique invariant measure $\mu_{\ast}$, i.e. $\mathcal{P}\mu_{\ast}=\mu_{\ast}$.

The drift condition has indeed been obtained in [AFGO18]. The AG method follows the dynamics

Next, let us prove that the drift condition holds. The proof is mainly built on Corollary 4.2. and Lemma 4.5. in [AFGO18].

By Corollary 4.2. and its proof in [AFGO18] (In [AFGO18], the noise are assumed to be independent. But a closer look at the proof of Corollary 4.2. reveals that our Assumption 2 suffices), we have

A closer look at the proof of Corollary 4.2. in [AFGO18] reveals that the following equality also holds:

When $f\in\mathcal{S}_{\mu,L}$ is strongly convex, Lemma 4.5. in [AFGO18] states that for any $\rho\in(0,1)$,

where $X:=\rho X_{1}+(1-\rho)X_{2}$, where

Taking expectation w.r.t. the noise $\varepsilon_{k+1}$ only in (127), we get

With the definition of $\rho_{\alpha,\beta}$, $P_{\alpha,\beta}$ by Lemma 21, we get

Then, combining (115) and (136), applying (137) and the definition of $V_{P_{\alpha,\beta}}$, we get

In the special case $(\alpha,\beta)=(\alpha_{AG},\beta_{AG})$, we obtain the following result.

Given $(\alpha,\beta)=(\alpha_{AG},\beta_{AG})$.

By letting $(\alpha,\beta)=(\alpha_{AG},\beta_{AG})$ in Lemma 33, we get

for some probability measure $\nu_{2}$. Let us define:

We define $\nu_{2}$ such that $\nu_{2}(A\times B)=0$ for any $B$ that does not contain $B_{R}$, and $\nu_{2}(A\times B)=\nu_{1}(A)$ for some probability measure $\nu_{1}$ and for any $B$ that contains $B_{R}$.Then, it suffices for us to show that

where $\nu(x)$ is the probability density function for some probability measure $\nu_{1}(\cdot)$.

For any $\eta\in(0,1)$, there exists some $R>0$ such that

Note that for sufficiently large $M$, $\int_{\|x-x_{\ast}\|\leq M}p(\xi_{\ast},x)dx$ can get arbitrarily close to $1$. Fix $M$, by the continuity of $p(\xi,x)$ in both $\xi$ and $x$, we can find $\eta^{\prime}\in(0,1)$ such that uniformly in $\|x-x_{\ast}\|\leq M$,

which can be arbitrarily close to $1$ if we take $R>0$ to be sufficiently small. In particular, if we fix $\eta\in(0,1)$, then we can take $M>0$ such that

and similarly with fixed $\eta$ and $M$, we take $R>0$ such that uniformly in $\|x-x_{\ast}\|\leq M$,

Finally, we are ready to state the proof of Theorem 13 and Proposition 14.

According to the proof of Lemma 35, for any fixed $\eta>0$, we can define:

where $\bar{\eta}=(1-(\eta-\eta_{0}))\vee(2+R\psi\gamma_{0})/(2+R\psi)$ and $\psi=\eta_{0}/K_{\alpha,\beta}$, where $\eta_{0}\in(0,\eta)$ and $\gamma_{0}\in(\gamma_{\alpha,\beta}+2K_{\alpha,\beta}/R,1)$. In particular, we can choose

where $\psi:=\frac{\eta}{2K_{\alpha,\beta}}$ so that

Let us recall that $\gamma=\rho=1-\frac{1}{\sqrt{\kappa}}$ and $K=\frac{\sigma^{2}}{L}$. Recall that $\gamma_{0}$ satisfies $\gamma_{0}\in(\gamma+2K/R,1)$ and let us assume that $K$ is sufficiently small so that $K\leq\frac{R}{4\sqrt{\kappa}}$, then we can take

We also recall that $\psi=\eta_{0}/K$ and

We have discussed before that we can take $\eta$ to be arbitrarily close to $1$ by taking $M$ sufficiently large, and for fixed $M$ take $R$ sufficiently small. Let us take

If we take $K<R\eta_{0}=\frac{R}{2\sqrt{\kappa}}$, then

Hence, we can take $K\leq\frac{R}{4\sqrt{\kappa}}$, that is,

Finally, we want to take $R>0$ and $M>0$ such that

It is easy to see that we can take $M$ so that

and take $R$ such that for any $\|x-x_{\ast}\|\leq M$,

Recall that $\nu_{k,\alpha,\beta}$ denotes the law of the iterates $\xi_{k}$. By Lemma 32, the Markov chain $\xi_{k}$ admits a unique invariant distribution $\pi_{\alpha,\beta}$. By letting $\mu_{1}=\nu_{0,\alpha,\beta}$ and $\mu_{2}=\pi_{\alpha,\beta}$, we conclude that

Finally, let us prove (33). Given $(\alpha,\beta)=(\alpha_{AG},\beta_{AG})$, we have $\rho_{\alpha,\beta}=1-\frac{1}{\sqrt{\kappa}}$, $\alpha=\frac{1}{L}$. It follows from Lemma 34 and its proof that

By induction on $k$, we can show that for every $k$,

By the definition of $V_{P}$, it follows that

In Proposition 14, the amount of noise that can be tolerated is limited. Nevertheless, in applications where the gradient is estimated from noisy measurements, such results would be applicable if the noise level is mild [BWBZ13].

If the noise $\varepsilon_{k}$ are i.i.d. Gaussian $\mathcal{N}(0,\Sigma)$, then conditional on $x_{k}=x_{k-1}=x_{\ast}$ in the AG method, with stepsize $\alpha=1/L$, $x_{k+1}$ is distributed as $\mathcal{N}(x_{\ast},L^{-2}\Sigma)$ with $\Sigma\preceq L^{2}I_{d}$. Therefore, for $\gamma>0$ sufficiently small,

By Chebychev’s inequality, letting $\gamma=1/2$, for any $m\geq 0$, we get

Conditional on $(x_{k}^{T},x_{k-1}^{T})^{T}=\xi=(\xi_{(1)}^{T},\xi_{(2)}^{T})^{T}$, where $V_{P}(\xi)\leq r$ for some $r>0$, then, $x_{k+1}$ is Gaussian distributed:

Thus, uniformly in $\|x-x_{\ast}\|\leq M$,

Note that $V_{P_{AG}}(\xi)\leq r$ implies that

Thus, uniformly in $\|x-x_{\ast}\|\leq M$,

For the remaining of the proof, without loss of generality assume that $\mu=\Theta(1)$ and $L=\Theta(\kappa)$.Given two scalar-valued functions $f$ and $g$, we say $f=\Theta(g)$, if the ratio $f(x)/g(x)$ lies in an interval $[c_{1},c_{2}]$ for every $x$ and some$c_{1},c_{2}>0$. It is straightforward to see from the Taylor expansion of $M$ that $M=O(\kappa^{-1/8})$ and

D.2 Proofs of Results in Section A

Consider the constrained optimization problem

where $\varepsilon_{k}$ is the random gradient error satisfying Assumption 2, $\alpha,\beta>0$ are the stepsize and momentum parameter and the projection onto the convex compact set $C$ with diameter $\mathcal{D}_{\mathcal{C}}$ can be written as

Due to the non-expansiveness property of the projection operator, we have (see e.g. [CW05, Lemma 2.4])

Following a similar approach to [HL17, FRMP17], we reformulate the projected AG iterations as a linear dynamical system as

where $G_{M}:=\max_{x\in\mathcal{C}}\|\nabla f(x)\|$.

In particular, with $\rho=1-\frac{1}{\sqrt{\kappa}}$, $\beta=\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$, $\alpha=\frac{1}{L}$ where $\kappa=\frac{L}{\mu}$. Then (148) holds with the matrix

We follow the proof technique of [FRMP17] for deterministic proximal AG which is based on [Nes04, Lemma 2.4] and adapt this proof technique to accelerated stochastic projected gradient. Defining the error at step $k$

where in the first inequality we used the fact that the gradient of $f$ is $L$-smooth which implies that

(see e.g. [Bub14]) and second inequality follows from Jensen’s inequality.Finally, the last step is a consequence of (143) and Assumption 2 on the noise. It follows from a similar computation that

We note that the matrices $X_{1}$ and $X_{2}$ can be written as

where $A,B,C,E$ are defined by (147). Using [FRMP17, eqn. (36)–(37)] and Lemma 38, we have

Plugging these into (155) and (156), we obtain

Taking conditional expectations and inserting (161)–(162),

Using the notations as in the proof of Lemma 37, we have the following two inequalities:

Recall that f satisfies following inequalities,

This proves (159). Finally, (160) can also be obtained if we take $x=x_{*}$ and follow similar steps. ∎

Given $\alpha=\frac{1}{L}$, $\beta=\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$, where $\kappa=L/\mu$, we have

and with $\alpha=\frac{1}{L}$, $\beta=\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$, we have

The proof is similar to the proof of Theorem 13 and the proof of (33). We obtain

As in the proof of Proposition 14, we can take

Finally, the proof of (39) is similar as the proof of (37). We obtain

Appendix E Numerical Illustrations

In this section, we illustrate some of our theoretical results over some simple functions with numerical experiments. On the left panel of Figure 1, we compare ASG for the quadratic objective $f(x)=x^{2}/2$ in dimension one with additive i.i.d. Gaussian noise on the gradients for different noise levels $\sigma\in\{0.01,0.1,1,2\}$. The plots show performance with respect to expected suboptimality using $10^{4}$ sample paths. As expected, the performance deteriorates when $\sigma$ increases. The fact that the performance stabilizes after a certain number of iterations supports the claim that a stationary distribution exists, a claim that was proved in Theorem 4. In the middle panel, we repeat the experiment in dimension $d=10$ over the quadratic objective $f(x)=\frac{1}{2}x^{T}Qx$, where $Q$ is a diagonal matrix with diagonal entries $Q_{ii}=1/i$. We observe similar patterns.

Finally, on the right panel of Figure 1, we estimate the distribution of $f(x_{k})$ for $k\in\{5,25,125,625\}$. For this purpose, we plot the histograms of $f(x_{k})$ over $10^{4}$ sample paths for every fixed $k$. We observe that the histograms for $k=125$ and $625$ are similar, illustrating the fact that ASG admits a stationary distribution.