From Averaging to Acceleration, There is Only a Step-size

Introduction

Many problems in machine learning are naturally cast as convex optimization problems over a Euclidean space; for supervised learning this includes least-squares regression, logistic regression, and the support vector machine. Faced with large amounts of data, practitioners often favor first-order techniques based on gradient descent, leading to algorithms with many cheap iterations. For smooth problems, two extensions of gradient descent have had important theoretical and practical impacts: acceleration and averaging.

Acceleration techniques date back to Nesterov (1983) and have their roots in momentum techniques and conjugate gradient (Polyak, 1987). For convex problems, with an appropriately weighted momentum term which requires to store two iterates, Nesterov (1983) showed that the traditional convergence rate of $O(1/n)$ for the function values after $n$ iterations of gradient descent goes down to $O(1/n^{2})$ for accelerated gradient descent, such a rate being optimal among first-order techniques that can access only sequences of gradients (Nesterov, 2004). Like conjugate gradient methods for solving linear systems, these methods are however more sensitive to noise in the gradients; that is, to preserve their improved convergence rates, significantly less noise may be tolerated (d’Aspremont, 2008; Schmidt et al., 2011; Devolder et al., 2014).

Averaging techniques which consist in replacing the iterates by the average of all iterates have also been thoroughly considered, either because they sometimes lead to simpler proofs, or because they lead to improved behavior. In the noiseless case where gradients are exactly available, they do not improve the convergence rate in the convex case; worse, for strongly-convex problems, they are not linearly convergent while regular gradient descent is. Their main advantage comes with random unbiased gradients, where it has been shown that they lead to better convergence rates than the unaveraged counterparts, in particular because they allow larger step-sizes (Polyak and Juditsky, 1992; Bach and Moulines, 2011). For example, for least-squares regression with stochastic gradients, they lead to convergence rates of $O(1/n)$, even in the non-strongly convex case (Bach and Moulines, 2013).

In this paper, we show that for quadratic problems, both averaging and acceleration are two instances of the same second-order finite difference equation, with different step-sizes. They may thus be analyzed jointly, together with a non-strongly convex version of the heavy-ball method (Polyak, 1987, Section 3.2). In presence of random zero-mean noise on the gradients, this joint analysis allows to design a novel intermediate algorithm that exhibits the good aspects of both acceleration (quick forgetting of initial conditions) and averaging (robustness to noise).

In this paper, we make the following contributions:

We show in Section 2 that accelerated gradient descent, averaged gradient descent and the heavy-ball method for non-strongly-convex problems may be reformulated as constant parameter second-order difference equation algorithms, where stability of the system is equivalent to convergence at rate $O(1/n^{2})$.

In Section 3, we provide a detailed analysis of the eigenvalues of the corresponding linear dynamical system, showing various oscillatory and non-oscillatory behaviors, together with a sharp stability result with explicit constants.

In Section 4, we consider the situation where noisy gradients are available, where we extend our general convergence result, which suggests an alternative algorithm (i.e., with different step sizes) that exhibits the good aspects of both averaging and acceleration.

In Section 5, we illustrate our results with simulations on synthetic examples.

Second-Order Iterative Algorithms for Quadratic Functions

In this paper we study second-order iterative algorithms of the form:

By letting $\phi_{n}=\theta_{n}-\theta_{*}$ we then have $\phi_{n+1}=A_{n}\phi_{n}+B_{n}\phi_{n-1}$, started from $\phi_{0}=\phi_{1}=\theta_{0}-\theta_{*}$. Thus, we restrict our problem to the study of the convergence of an iterative system to .

In connection with accelerated methods, we are interested in algorithms for which $f(\theta_{n})-f(\theta_{*})=\frac{1}{2}\langle\phi_{n},H\phi_{n}\rangle$ converges to at a speed of $O\left(1/n^{2}\right)$. Within this context we impose that $A_{n}$ and $B_{n}$ have the form :

By letting $\eta_{n}=n\phi_{n}=n(\theta_{n}-\theta_{*})$, we can now study the simple iterative system with constant terms $\eta_{n+1}=A\eta_{n}+B\eta_{n-1}$, started at $\eta_{0}=0$ and $\eta_{1}=\theta_{0}-\theta_{*}$. Showing that the function values remain bounded, we directly have the convergence of $f(\theta_{n})$ to $f(\theta_{*})$ at the speed $O\left(1/n^{2}\right)$. Thus the n-scalability property allows to switch from a convergence problem to a stability problem.

For feasibility concerns the method can only access $H$ through matrix-vector products. Therefore $A$ and $B$ should be polynomials in $H$ and $c$ a polynomial in $H$ times $q$, if possible of low degree. The following theorem clarifies the general form of iterative systems which share these three properties (see proof in Appendix B).

Note that our result prevents $A_{n}$ and $B_{n}$ from being zero, thus requiring the algorithm to strictly be of second order. This illustrates the fact that first-order algorithms as gradient descent do not have the convergence rate in $O(1/n^{2})$.

The recursion in Eq. (3) may be written with gradients of $f$ in multiple ways. In order to preserve the parallel with accelerated techniques, we rewrite it as:

2 Examples

When computing the average online as $\theta_{n+1}=\theta_{n}+\frac{1}{n+1}(\psi_{n+1}-\theta_{n})$ and seeing the average as the main iterate, the algorithm becomes (see proof in Appendix B.2):

This corresponds to Eq. (4) with $\alpha=0$ and $\beta=\gamma$.

For smooth optimization the accelerated literature (Nesterov, 2004; Beck and Teboulle, 2009) uses the step-size $\delta_{n}=1-\frac{3}{n+1}$ and their results are not valid for bigger step-size $\delta_{n}$. However $\delta_{n}=1-\frac{2}{n+1}$ is compatible with the framework of Lan (2012) and is more convenient for our set-up. This corresponds to Eq. (4) with $\alpha=\gamma$ and $\beta=\gamma$. Note that accelerated techniques are more generally applicable, e.g., to composite optimization with smooth functions (Nesterov, 2013; Beck and Teboulle, 2009).

when $\delta_{n}=1-\frac{2}{n+1}$. We note that typically $\delta_{n}$ is constant for strongly-convex problems. This corresponds to Eq. (4) with $\alpha=\gamma$ and $\beta=0$.

Convergence with Noiseless Gradients

We study the convergence of the iterates defined by: $\eta_{n+1}=\left(I-\alpha H\right)\eta_{n}+\left(I-\beta H\right)\left(\eta_{n}-\eta_{n-1}\right)$. This is a second-order iterative system with constant coefficients that it is standard to cast in a linear framework (see, e.g., Ortega and Rheinboldt, 2000). We may rewrite it as:

Thus $\Theta_{n}=F^{n}\Theta_{0}$. Following O’Donoghue and Candes (2013), if we consider an eigenvalue decomposition of $H$, i.e., $H=P\text{Diag}(h)P^{\top}$ with $P$ an orthogonal matrix and $(h_{i})$ the eigenvalues of $H$, sorted in decreasing order: $h_{d}=L\geq h_{d-1}\geq\cdots\geq h_{2}\geq h_{1}=\mu>0$, then Eq. (3) may be rewritten as:

In this section, we consider a fixed $i\in\{1,\dots,d\}$ and study the stability in the corresponding eigenspace. This linear dynamical system may be analyzed by studying the eigenvalues of the $2\times 2$-matrix $F_{i}=\begin{pmatrix}2-(\alpha+\beta)h_{i}&\beta h_{i}-1\\ 1&0\end{pmatrix}$. These eigenvalues are the roots of its characteristic polynomial which is:

To compute the roots of the second-order polynomial, we compute its reduced discriminant:

Depending on the sign of the discriminant $\Delta_{i}$, there will be two real distinct eigenvalues ($\Delta_{i}>0)$, two complex conjugate eigenvalues ($\Delta_{i}<0)$ or a single real eigenvalue ($\Delta_{i}=0)$.

We will now study the sign of $\Delta_{i}$. In each different case, we will determine under what conditions on $\alpha$ and $\beta$ the modulus of the eigenvalues is less than one, which means that the iterates $(\eta^{i}_{n})_{n}$ remain bounded and the iterates $(\theta_{n})_{n}$ converge to $\theta_{*}$. We may then compute function values as $f(\theta_{n})-f(\theta_{\ast})=\frac{1}{2n^{2}}\sum_{i=1}^{d}(\eta_{n}^{i})^{2}h_{i}=\frac{1}{2}\sum_{i=1}^{d}(\phi_{n}^{i})^{2}h_{i}$.

The various regimes are summarized in Figure 1: there is a triangle of values of $(\alpha h_{i},\beta h_{i})$ for which the algorithm remains stable (i.e., the iterates $(\eta_{n})_{n}$ do not diverge), with either complex or real eigenvalues. In the following lemmas (see proof in Appendix C), we provide a detailed analysis that leads to Figure 1.

The discriminant $\Delta_{i}$ is strictly positive and the algorithm is stable if and only if

We then have two real roots $r_{i}^{\pm}=r_{i}\pm\sqrt{\Delta_{i}}$, with $r_{i}=1-(\frac{\alpha+\beta}{2})h_{i}$. Moreover, we have:

Therefore, for real eigenvalues, $(({\phi_{n}^{i}})^{2}h_{i})_{n}$ will converge to at a speed of $O(1/n^{2})$ however the constant $\Delta_{i}$ may be arbitrarily small (and thus the scaling factor arbitrarily large). Furthermore we have linear convergence if the inequalities in the lemmas are strict.

The discriminant $\Delta_{i}$ is stricly negative and the algorithm is stable if and only if

We then have two complex conjugate eigenvalues: $r_{i}^{\pm}=r_{i}\pm\sqrt{-1}\sqrt{-\Delta_{i}}$. Moreover, we have:

with $\rho_{i}=\sqrt{1-\beta h_{i}}$, and $\omega_{i}$ defined through $\sin(\omega_{i})=\sqrt{-\Delta_{i}}/\rho_{i}$ and $\cos(\omega_{i})=r_{i}/\rho_{i}$.

Therefore, for complex eigenvalues, there is a linear convergence if the inequalities in the lemma are strict. Moreover, $(({\phi_{n}^{i}})^{2}h_{i})_{n}$ oscillates to at a speed of $O(1/n^{2})$ even if $h_{i}$ is arbitrarily small.

When the discriminant goes to zero in the explicit formulas of the real and complex cases, both the denominator and numerator of $(({\phi_{n}^{i}})^{2}h_{i})_{n}$ will go to zero. In the limit case, when the discriminant is equal to zero, we will have a double real eigenvalue. This happens for $\beta=2\sqrt{\alpha/h_{i}}-\alpha$. Then the eigenvalue is $r_{i}=1-\sqrt{\alpha h_{i}}$, and the algorithm is stable for $0<\alpha<{4}/{h_{i}}$, we then have $({\phi_{n}^{i}})^{2}h_{i}=h_{i}({\phi_{1}^{i}})^{2}(1-\sqrt{\alpha h_{i}})^{2(n-1)}$. This can be obtained by letting $\Delta_{i}$ goes to in the real and complex cases (see also Appendix C.3).

To conclude the iterate $(\eta_{n}^{i})_{n}=(n(\theta_{n}^{i}-\theta_{\ast}^{i}))_{n}$ will be stable for $\alpha\in[0,4/h_{i}]$ and $\beta\in[0,2/h_{i}-\alpha/2]$. According to the values of $\alpha$ and $\beta$ this iterate will have a different behavior. In the complex case, the roots are complex conjugate with magnitude $\sqrt{1-\beta h_{i}}$. Thus, when $\beta>0$, $(\eta_{n}^{i})_{n}$ will converge to , oscillating, at rate $\sqrt{1-\beta h_{i}}$. In the real case, the two roots are real and distinct. However the product of the two roots is equal to $\sqrt{1-\beta h_{i}}$, thus one will have a higher magnitude and $(\eta_{n}^{i})_{n}$ will converges to at rate higher than in the complex case (as long as $\alpha$ and $\beta$ belong to the interior of the stability region).

Finally, for a given quadratic function $f$, all the $d$ iterates $(\eta_{n}^{i})_{n}$ should be bounded, therefore we must have $\alpha\in[0,4/L]$ and $\beta\in[0,2/L-\alpha/2]$. Then, depending on the value of $h_{i}$, some eigenvalues may be complex or real.

2 Classical examples

For particular choices of $\alpha$ and $\beta$, displayed in Figure 1, the eigenvalues are either all real or all complex, as shown in the table below.

Averaged gradient descent loses linear convergence for strongly-convex problems, because $r_{i}^{+}=1$ for all eigensubspaces. Similarly, the heavy ball method is not adaptive to strong convexity because $\rho_{i}=1$. However, accelerated gradient descent, although designed for non-strongly-convex problems, is adaptive because $\rho_{i}=\sqrt{1-\gamma h_{i}}$ depends on $h_{i}$ while $\alpha$ and $\beta$ do not. These last two algorithms have an oscillatory behavior which can be observed in practice and has been already studied (Su et al., 2014).

Note that all the classical methods choose step-sizes $\alpha$ and $\beta$ either having all the eigenvalues real either complex; whereas we will see in Section 4, that it is significant to combine both behaviors in presence of noise.

3 General bound

Even if the exact formulas in Lemmas 1 and 2 are computable, they are not easily interpretable. In particular when the two roots become close, the denominator will go to zero, which prevents from bounding them easily. When we further restrict the domain of $(\alpha,\beta)$, we can always bound the iterate by the general bound (see proof in Appendix D):

For $\alpha\leq 1/h_{i}$ and $0\leq\beta\leq 2/h_{i}-\alpha$, we have

These bounds are shown by dividing the set of $(\alpha,\beta)$ in three regions where we obtain specific bounds. They do not depend on the regime of the eigenvalues (complex or real); this enables us to get the following general bound on the function values, our main result for the deterministic case.

For $\alpha\leq 1/L$ and $0\leq\beta\leq 2/L-\alpha$:

The first bound $\frac{{\|\theta_{0}-\theta_{*}\|}^{2}}{\alpha n^{2}}$ corresponds to the traditional acceleration result, and is only relevant for $\alpha>0$ (that is, for Nesterov acceleration and the heavy-ball method, but not for averaging). We recover the traditional convergence rate of second-order methods for quadratic functions in the singular case, such as conjugate gradient (Polyak, 1987, Section 6.1).

While the result above focuses on function values, like most results in the non-strongly convex case, the distance to optimum $\|\theta_{n}-\theta_{\ast}\|^{2}$ typically does not go to zero (although it remains bounded in our situation).

When $\alpha=0$ (averaged gradient descent), then the second bound $\frac{4{\|\theta_{0}-\theta_{*}\|}^{2}}{(\alpha+\beta)n}$ provides a convergence rate of $O(1/n)$ if no assumption is made regarding the starting point $\theta_{0}$, while the last bound of Theorem 2 would lead to a bound $\frac{{8\|H^{-1/2}(\theta_{0}-\theta_{*})\|}^{2}}{(\alpha+\beta)^{2}n^{2}}$, that is a rate of $O(1/n^{2})$, only for some starting points.

As shown in Appendix E by exhibiting explicit sequences of quadratic functions, the inverse dependence in $\alpha n^{2}$ and $(\alpha+\beta)n$ in Eq. (10) is not improvable.

Quadratic Optimization with Additive Noise

In many practical situations, the gradient of $f$ is not available for the recursion in Eq. (4), but only a noisy version. In this paper, we only consider additive uncorrelated noise with finite variance.

For quadratic functions, for the reduced variable $\eta_{n}=n\phi_{n}=n(\theta_{n}-\theta_{\ast})$, we get:

Note that algorithms with $\alpha\neq 0$ will have an important level of noise because of the term $n\alpha\varepsilon_{n+1}$. We denote by $\xi_{n+1}=\begin{pmatrix}[n\alpha+\beta]\varepsilon_{n+1}\\ 0\end{pmatrix}$ and we now have the recursion:

2 Convergence result

For a quadratic function $f$ with arbitrarily small eigenvalues and uncorrelated noise with finite covariance, we obtain the following convergence result (see proof in Appendix F); since we will allow the parameters $\alpha$ and $\beta$ to depend on the time we stop the algorithm, we introduce the horizon $N$:

Although we only provide an upper-bound, the proof technique relies on direct moment computations in each eigensubspace with few inequalities, and we conjecture that the scalings with respect to $n$ are tight.

For $\alpha=0$ and $\beta=1/L$ (which corresponds to averaged gradient descent), the second bound leads to $\frac{4{L\|\theta_{0}-\theta_{*}\|}^{2}}{N}+\frac{4\mathop{\rm tr}(C)}{L}$, which is bounded but not converging to zero. We recover a result from Bach and Moulines (2011, Theorem 1).

For $\alpha=\beta=1/L$ (which corresponds to Nesterov’s acceleration), the first bound leads to $\frac{{L\|\theta_{0}-\theta_{*}\|}^{2}}{N^{2}}+\frac{(N+1)\mathop{\rm tr}(C)}{L}$, and our bound suggests that the algorithm diverges, which we have observed in our experiments in Appendix A.

For $\alpha=0$ and $\beta=1/L\sqrt{N}$, the second bound leads to $\frac{{4L\|\theta_{0}-\theta_{*}\|}^{2}}{\sqrt{N}}+\frac{4\mathop{\rm tr}(C)}{L\sqrt{N}}$, and we recover the traditional rate of $1/\sqrt{N}$ for stochastic gradient in the non-strongly-convex case.

When the values of the bias and the variance are known we can choose $\alpha$ and $\beta$ such that the trade-off between the bias and the variance is optimal in our bound, as the following corrollary shows. Note that in the bound below, taking a non zero $\beta$ enables the bias term to be adaptive to hidden strong-convexity.

For $\alpha=\min\left\{\frac{\|\theta_{0}-\theta_{*}\|}{2\sqrt{\mathop{\rm tr}C}N^{3/2}},1/L\right\}$ and $\beta\in[0,\min\{N\alpha,1/L\}]$, we have:

3 Structured noise and least-square regression

When only the noise total variance $\mathop{\rm tr}(C)$ is considered, as shown in Section 4.4, Corollary 2 recover existing (more general) results. Our framework however leads to improved result for structured noise processes frequent in machine learning, in particular in least-squares regression which we now consider but this goes beyond (see, e.g. Bach and Moulines, 2013).

For this particular structured noise we can take advantage of a large $\beta$:

For $\alpha=0$ and $\beta=1/L$ (which corresponds to averaged gradient descent), the second bound leads to $\frac{4{L\|\theta_{0}-\theta_{*}\|}^{2}}{N}+\frac{8\mathop{\rm tr}(CH^{-1})}{N}$. We recover a result from Bach and Moulines (2013, Theorem 1). Note that when $C\preccurlyeq\sigma^{2}H$, $\mathop{\rm tr}(CH^{-1})\leqslant\sigma^{2}d$.

For $\alpha=\beta=1/L$ (which corresponds to Nesterov’s acceleration), the first bound leads to $\frac{{L\|\theta_{0}-\theta_{*}\|}^{2}}{N^{2}}+\mathop{\rm tr}(CH^{-1})$, which is bounded but not converging to zero (as opposed to the the unstructured noise where the algorithm may diverge).

For $\alpha=1/(LN^{a})$ with $0\leq a\leq 1$ and $\beta=1/L$, the first bound leads to $\frac{L\|\theta_{0}-\theta_{*}\|^{2}}{N^{2-a}}+\frac{\mathop{\rm tr}(CH^{-1})}{N^{a}}$. We thus obtain an explicit bias-variance trade-off by changing the value of $a$.

When the values of the bias and the variance are known we can choose $\alpha$ and $\beta$ with an optimized trade-off, as the following corrollary shows:

For $\alpha=\min\left\{\frac{\|\theta_{0}-\theta_{*}\|}{\sqrt{L\mathop{\rm tr}(CH^{-1})}N},1/L\right\}$ and $\beta=\min\left\{N\alpha,1/L\right\}$ we have:

4 Related work

Several authors (Lan, 2012; Hu et al., 2009; Xiao, 2010) have shown that using a step-size proportional to $1/N^{3/2}$ accelerated methods with noisy gradients lead to the same convergence rate of O\big{(}\frac{L\|\theta_{0}-\theta_{*}\|^{2}}{N^{2}}+\frac{\|\theta_{0}-\theta_{*}\|\sqrt{\mathop{\rm tr}(C)}}{\sqrt{N}}\big{)} than in Corollary 2, for smooth functions. Thus, for unstructured noise, our analysis provides insights in the behavior of second-order algorithms, without improving bounds. We get significant improvements for structured noises.

When the noise is structured as in least-square regression and more generally in linear supervised learning, Bach and Moulines (2011) have shown that using averaged stochastic gradient descent with constant step-size leads to the convergence rate of O\big{(}\frac{L\|\theta_{0}-\theta_{0}\|^{2}}{N}+\frac{\sigma^{2}d}{N}\big{)}. It has been highlighted by Défossez and Bach (2014) that the bias term $\frac{L\|\theta_{0}-\theta_{*}\|^{2}}{N}$ may often be the dominant one in practice. Our result in Corollary 3 leads to an improved bias term in $O(1/N^{2})$ with the price of a potentially slightly worse constant in the variance term. However, with optimal constants in Corollary 3, the new algorithm is always an improvement over averaged stochastic gradient descent in all situations. If constants are unknown, we may use $\alpha=1/(LN^{a})$ with $0\leq a\leq 1$ and $\beta=1/L$ and we choose $a$ depending on the emphasis we want to put on bias or variance.

For noisy quadratic problems, the convergence rate nicely decomposes into two terms, a bias term which corresponds to the noiseless problem and the variance term which corresponds to a problem started at $\theta_{\ast}$. For each of these two terms, lower bounds are known. For the bias term, if $N\leq d$, then the lower bound is, up to constants, $L\|\theta_{0}-\theta_{\ast}\|^{2}/N^{2}$ (Nesterov, 2004, Theorem 2.1.7). For the variance term, for the general noisy gradient situation, we show in Appendix H that for $N\leq d$, it is ${(\mathop{\rm tr}C)}/({L\sqrt{N}})$, while for least-squares regression, it is $\sigma^{2}d/N$ (Tsybakov, 2003). Thus, for the two situations, we attain the two lower bounds simultaneously for situations where respectively $L\|\theta_{0}-\theta_{\ast}\|^{2}\leq(\mathop{\rm tr}C)/L$ and $L\|\theta_{0}-\theta_{\ast}\|^{2}\leq d\sigma^{2}$. It remains an open problem to achieve the two minimax terms in all situations.

We also note as shown in Appendix G that in the special case of quadratic functions, the algorithms of Lan (2012); Hu et al. (2009); Xiao (2010) could be unified into our framework (although they have significantly different formulations and justifications in the smooth case).

Experiments

We compare our algorithm to other stochastic accelerated algorithms, that is, AC-SA (Lan, 2012), SAGE (Hu et al., 2009) and Acc-RDA (Xiao, 2010) which are presented in Appendix G. For all these algorithms (and ours) we take the optimal step-sizes defined in these papers. We show results averaged over 10 replications.

We first consider an i.i.d. zero mean noise whose covariance matrix is proportional to $H$. We also consider a variant of our algorithm with an any-time step-size function of $n$ rather than $N$ (for which we currently have no proof of convergence). In Figure 3, we take into account two different set-ups. In the left plot, the variance dominates the bias (with $r=\|\theta_{0}-\theta_{\ast}\|=\sigma$). We see that (a) Acc-GD does not converge to the optimum but does not diverge either, (b) Av-GD and our algorithms achieve the optimal rate of convergence of $O(\sigma^{2}d/n)$, whereas (c) other accelerated algorithms only converge at rate $O(1/\sqrt{n})$. In the right plot, the bias dominates the variance ($r=10$ and $\sigma=0.1$). In this situation our algorithm outperforms all others.

We now see how these algorithms behave for least-squares regressions and the regular (non-homoscedastic) stochastic gradients described in Section 4.3. We consider normally distributed inputs. The covariance matrix $H$ is the same as before. The outputs are generated from a linear function with homoscedatic noise with a signal-to-noise ratio of $\sigma$. We consider $d=20$. We show results averaged over 10 replications. In Figure 4, we consider again a situation where the bias dominates (left) and vice versa (right). We see that our algorithm has the same good behavior than in the homoscedastic noise case and we conjecture that our bounds also hold in this situation.

Conclusion

We have provided a joint analysis of averaging and acceleration for non-strongly-convex quadratic functions in a single framework, both with noiseless and noisy gradients. This allows to define a class of algorithms that can benefit simultaneously of the known improvements of averaging and accelerations: faster forgetting of initial conditions (for acceleration), and better robustness to noise when the noise covariance is proportional to the Hessian (for averaging).

Our current analysis of our class of algorithms in Eq. (4), that considers two different affine combinations of previous iterates (instead of one for traditional acceleration), is limited to quadratic functions; an extension of its analysis to all smooth or self-concordant-like functions would widen its applicability. Similarly, an extension to least-squares regression with natural heteroscedastic stochastic gradient, as suggested by our simulations, would be an interesting development.

This work was partially supported by the MSR-Inria Joint Centre and a grant by the European Research Council (SIERRA project 239993). The authors would like to thank Aymeric Dieuleveut for helpful discussions.

References

Appendix A Additional experimental results

In this appendix, we provide additional experimental results to illustrate our theoretical results.

In Figure 5, we minimize a one-dimensional quadratic function $f(\theta)=\frac{1}{2}\theta^{2}$ for a fixed step-size $\alpha=1/10$ and different step-sizes $\beta$. In the left plot, we compare Acc-GD, HB and Av-GD. We see that HB and Acc-GD both oscillate and that Acc-GD leverages strong convexity to converge faster. In the right plot, we compare the behavior of the algorithm for different values of $\beta$. We see that the optimal rate is achieved for $\beta=\beta_{*}$ defined to be the one for which there is a double coalescent eigenvalue, where the convergence is linear at speed $O(1-\sqrt{\alpha L})^{n}$. When $\beta>\beta_{*}$, we are in the real case and when $\beta<\beta_{*}$ the algorithm oscillates to the solution.

Figure 6 shows interactions between different eigenspaces. In the left plot, we optimize a quadratic function of dimension $d=2$. The first eigenvalue is $L=1$ and the second is $\mu=2^{-8}$. For Av-GD the convergence is of order $O(1/n)$ since the problem is “not” strongly convex (i.e., not appearing as strongly convex since $n\mu$ remains small). The convergence is at the beginning the same for HB and Acc-GD, with oscillation at speed $O(1/n^{2})$, since the small eigenvalue prevents Acc-GD from having a linear convergence. Then for large $n$, the convergence becomes linear for Acc-GD, since $\mu n$ becomes large. In the right plot, we optimize a quadratic function in dimension $d=5$ with eigenvalues from $1$ to $0.1$. We show the function values of the projections of the iterates $\eta_{n}$ on the different eigenspaces. We see that high eigenvalues first dominate, but converge quickly to zero, whereas small ones keep oscillating, and converge more slowly.

In Figure 7, we optimize two $20$-dimensional quadratic functions with different eigenvalues with Av-GD, HB and Acc-GD for a fixed step-size $\gamma=1/10$. In the left plot, the eigenvalues are $1/k^{2}$ and in the right one, they are $1/k^{8}$, for $k=1,\dots,d$. We see that in both cases, Av-GD converges at a rate of $O(1/n)$ and HB at a rate of $O(1/n^{2})$. For Acc-GD the convergence is linear when $\mu$ is large (left plot) and becomes sublinear at a rate of $O(1/n^{2})$ when $\mu$ becomes small (right plot).

A.2 Noisy convergence with unstructured additive noise

We optimize the same quadratic function, but now with noisy gradients. We compare our algorithm to other stochastic accelerated algorithms, that is, AC-SA [Lan, 2012], SAGE [Hu et al., 2009] and Acc-RDA [Xiao, 2010], which are presented in Appendix G. For all these algorithms (and ours) we take the optimal step-sizes defined in these papers. We plot the results averaged over 10 replications.

We consider in Figure 8 an i.i.d. zero mean noise of variance $C=I$. We see that all the accelerated algorithms achieve the same precision whereas Av-GD with constant step-size does not converge and Acc-Gd diverges. However SAGE and AC-SA are anytime algorithms and are faster at the beginning since their step-sizes are decreasing and not a constant (with respect to $n$) function of the horizon $N$.

Appendix B Proofs of Section 2

And in connection with Eq. (16) we can rewrite $P$ and $Q$ as:

Thus for $n=1$, we have $\bar{p}=2$. Then $-\frac{n-1}{n+1}\bar{q}=\frac{2n}{n+1}-1=\frac{n-1}{n+1}$ and $\bar{q}=-1$. Therefore

We let $\bar{A}=-(\bar{P}+\bar{Q})$ and $\bar{B}=\bar{Q}$ so that we have:

B.2 Av-GD as two steps-algorithm

Appendix C Proof of Section 3

The discriminant $\Delta_{i}$ is strictly positive when \big{(}\frac{\alpha+\beta}{2}\big{)}^{2}h_{i}-\alpha>0. This is always true for $\alpha$ strictly negative. For $\alpha$ positive and for $h_{i}\neq 0$, this is true for $|\frac{\alpha+\beta}{2}|>\sqrt{\alpha/h_{i}}$ . Thus the discriminant $\Delta_{i}$ is strictly positive for

Then we determine when the modulus of the eigenvalues is less than one (which corresponds to $-1\leq r_{i}^{-}\leq r_{i}^{+}\leq 1$).

Figure 9 (where we plot all the constraints we have so far) enables to conclude that the discriminant $\Delta_{i}$ is strictly positive and the algorithm is stable when the following three conditions are satisfied:

For any of those $\alpha$ et $\beta$ we will have:

Since $\eta_{0}^{i}=0$, $c_{1}+c_{2}=0$ and for $n=1$, $c_{1}=\eta^{i}_{1}/(r_{i}^{-}-r_{i}^{+})$; we thus have:

C.2 Proof of Lemma 2

The discriminant $\Delta_{i}$ is strictly negative if and only if \big{(}\frac{\alpha+\beta}{2}\big{)}^{2}h_{i}-\alpha<0. This implies $|\frac{\alpha+\beta}{2}|<\sqrt{\alpha/h_{i}}$. The modulus of the eigenvalues is $|r_{i}^{\pm}|^{2}=1-\beta h_{i}$. Thus the discriminant $\Delta_{i}$ is strictly negative and the algorithm is stable for

For any of those $\alpha$ et $\beta$ we have:

with $\rho_{i}=\sqrt{1-\beta h_{i}}$, $\sin(\omega_{i})=\sqrt{-\Delta_{i}}/\rho_{i}$ and $\cos(\omega_{i})=r_{i}/\rho_{i}$. Since $\eta_{0}^{i}=0$, $c_{1}=0$ and we have for $n=1$, $c_{2}=\eta^{i}_{1}/(\sin(\omega_{i})\rho_{i})$. Therefore

C.3 Coalescing eigenvalues

When $\beta=2\sqrt{\alpha/h_{i}}-\alpha$, the discriminant $\Delta_{i}$ is equal to zero and we have a double real eigenvalue:

Thus the algorithm is stable for $\alpha<\frac{4}{h_{i}}$. For any of those $\alpha$ et $\beta$ we have:

This gives with $\eta_{0}^{i}=0$, $c_{1}=0$ and $c_{2}=\eta_{1}^{i}/r$. Therefore

In the presence of coalescing eigenvalues the convergence is linear if $0<\alpha<4/h_{i}$ and $h_{i}>0$, however one might worry about the behavior of $(({\phi_{n}^{i}})^{2}h_{i})_{n}$ when $h_{i}$ becomes small. Using the bound $x^{2}\exp(-x)\leq 1$ for $x\leq 1$, we have for $\alpha<4/h_{i}$:

Therefore we always have the following bound for $\alpha<4/h_{i}$:

Appendix D Proof of Theorem 2

We divide the domain of validity of Theorem 2 in three subdomains as explained in Figure 14. On the domain described in Figure 14 we have a first bound on the iterate $\eta_{n}^{i}$:

For $0\leq\alpha\leq 1/h_{i}$ and $1-\sqrt{1-\alpha h_{i}}<\beta h_{i}<1+\sqrt{1-\alpha h_{i}}$, we have:

And on the domain described Figure 14 we also have:

For $0\leq\alpha\leq 1/h_{i}$ and $\beta\leq\alpha$ we have:

These two lemmas enable us to prove the first bound of Theorem 2 since the domain of this theorem is included in the intersection of the two domains of these lemmas as shown in Figure 14.

Then we have the following bound on domain described in Figure 14:

For $0\leq\alpha\leq 2/h_{i}$ and $0\leq\beta\leq 2/h_{i}-\alpha$, we have:

Since the domain of definition of Theorem 2 is included in the domain of definition of Lemma 5 (as shown in Figure 14), this lemma proves the last two bounds of the theorem.

D.2 Outline of the proofs of the Lemmas

We also prove that $G(\eta^{i}_{n},\eta^{i}_{n-1})$ dominates $c\|\eta^{i}_{n}\|^{2}$ when we want to have a bound on $\|\eta^{i}_{n}\|^{2}$ of the form $\frac{1}{c}G(\eta^{i}_{1},\eta^{i}_{0})=\frac{1}{c}G(\theta_{0}^{i}-\theta_{*}^{i},0)$.

For readability, we remove the index $i$ and take $h_{i}=1$ without loss of generality.

D.3 Proof of Lemma 3

We first consider a quadratic Lyapunov function $\begin{pmatrix}\eta_{n}\\ \eta_{n-1}\end{pmatrix}^{\top}G_{1}\begin{pmatrix}\eta_{n}\\ \eta_{n-1}\end{pmatrix}$ with $G_{1}=\begin{pmatrix}1&\alpha-1\\ \alpha-1&1-\alpha\end{pmatrix}$. We note that $G_{1}$ is symmetric positive semi-definite for $\alpha\leq 1$. We recall $F_{i}=\begin{pmatrix}2-(\alpha+\beta)&\beta-1\\ 1&0\end{pmatrix}$.

For the result to be true we need for $0\leq\alpha\leq 1$ and $1-\sqrt{1-\alpha}<\beta<1+\sqrt{1-\alpha}$ two properties:

This especially shows Eq. (19) for the boundaries of the interval with $x=\pm\sqrt{1-\alpha}$.

This shows that for $0\leq\alpha\leq 1/h_{i}$ and $1-\sqrt{1-\alpha h_{i}}<\beta h_{i}<1+\sqrt{1-\alpha h_{i}}$:

D.4 Proof of Lemma 4

We consider now a second Lyapunov function $G_{2}(\eta_{n},\eta_{n-1})=(\eta_{n}-r\eta_{n-1})^{2}-\Delta(\eta_{n-1})^{2}$. We have:

Where we have used twice $r^{2}-\Delta=(1-\beta)$ and $\eta_{n}=2r\eta_{n-1}-(1-\beta)\eta_{n-2}$. Moreover $G_{2}(\eta_{n},\eta_{n-1})$ can be rewritten as:

Thus for $\alpha+\beta\leq 2$ and $\beta\leq\alpha$ we have:

Therefore for $\alpha+\beta\leq 2/h_{i}$ and $\beta\leq\alpha$, we have:

D.5 Proof of Lemma 5

Moreover for all $u\in$ and $n\geq 1$ we have $1-(1-u)^{n}\leq\sqrt{nu}$, since $1-(1-u)^{n}\leq 1$ and $1-(1-u)^{n}=u\sum(1-u)^{k}\leq nu$. Thus

Therefore for $0\leq\alpha\leq 2/h_{i}$ and $\alpha+\beta\leq 2/h_{i}$ we have:

Appendix E Lower bound

We have the following lower-bound for the bound shown in Corollary 1, which shows that depending on which of the two terms dominates, we may always find a sequence of functions that makes it tight.

Let $L\geq 0$. For all sequences $0\leq\alpha_{n}\leq 1/L$ and $0\leq\beta_{n}\leq 2/L-\alpha_{n}$, such that $\alpha_{n}+\beta_{n}=o(n\alpha_{n})$ there exists a sequence of one-dimensional quadratic functions $(f_{n})_{n}$ with second-derivative less than $L$ such that:

For all sequences $0\leq\alpha_{n}\leq 1/L$ and $0\leq\beta_{n}\leq 2/L-\alpha_{n}$, such that ${n\alpha_{n}}=o({\alpha_{n}+\beta_{n}})$, there exists a sequence of one-dimensional quadratic functions $(g_{n})_{n}$ with second-derivative less than $L$ such that:

For the first lower bound we consider $0\leq\alpha_{n}\leq 1/L$ and $0\leq\beta_{n}\leq 2/L-\alpha$, such that $\alpha_{n}+\beta_{n}=o(n\alpha_{n})$. We define $f_{n}=\pi^{2}/(4\alpha_{n}n^{2})$ and we consider the sequence of quadratic functions $f_{n}(\theta)=\frac{f_{n}\theta^{2}}{2}$. We consider the iterate $(\eta_{n})_{n}$ defined by our algorithm. We will show that

since $\frac{\beta_{n}}{\alpha_{n}n}=o(1)$. Also, $1-\frac{\pi^{2}(\alpha_{n}+\beta_{n})^{2}}{(4\alpha_{n}n)^{2}}=1+o(1)$, since $\alpha_{n}+\beta_{n}=o(n\alpha_{n})$. Moreover

thus $\omega_{n}={\pi}/(2n)+o(1/n)$ and $\sin(n\omega_{n})=1+o(1)$.

We consider now the situation where the second bound is active. Thus we take sequences $(\alpha_{n})$ and $(\beta_{n})$, such that ${n\alpha_{n}}=o({\alpha_{n}+\beta_{n}})$. We define $g_{n}=\frac{2}{n(\alpha_{n}+\beta_{n})}+\frac{4\alpha_{n}}{(\alpha_{n}+\beta_{n})^{2}}$ and consider the sequence of quadratic functions $g_{n}(\theta)=\frac{g_{n}\theta^{2}}{2}$. We will show for the iterate $(\eta_{n})$ defined by our algorithm that:

Thus $(n\Delta_{n})/g_{n}=\left(\frac{\alpha_{n}+\beta_{n}}{2}\right)$ and

Appendix F Proofs of Section 4

We decompose again vectors in an eigenvector basis of $H$ with $\eta_{n}^{i}=p_{i}^{\top}\eta_{n}$ and $\varepsilon_{n}^{i}=p_{i}^{\top}\varepsilon_{n}$:

We denote by $\xi_{n+1}^{i}=\begin{pmatrix}[n\alpha+\beta]\varepsilon_{n+1}^{i}\\ 0\end{pmatrix}$ and we have the reduced equation:

Unfortunately $F_{i}$ is not Hermitian and this formulation will not be convenient for calculus. Without loss of generality, we assume $r_{i}^{-}\neq r_{i}^{+}$ even if it means having $r_{i}^{-}-r_{i}^{+}$ goes to in the final bound. Let $Q_{i}=\begin{pmatrix}r_{i}^{-}&r_{i}^{+}\\ 1&1\end{pmatrix}$ be the transfer matrix of $F_{i}$, i.e., $F_{i}=Q_{i}D_{i}Q_{i}^{-1}$ with $D_{i}=\begin{pmatrix}r_{i}^{-}&0\\ 0&r_{i}^{+}\end{pmatrix}$ and $Q_{i}^{-1}=\frac{1}{r_{i}^{-}-r_{i}^{+}}\begin{pmatrix}1&-r_{i}^{+}\\ -1&r_{i}^{-}\end{pmatrix}$. We can reparametrize the problem in the following way:

with now $D_{i}$ Hermitian (even diagonal).

Thus it is easier to tackle using standard techniques for stochastic approximation [see, e.g., Polyak and Juditsky, 1992, Bach and Moulines, 2011]:

This is a bias-variance decomposition; the left term only depends on the initial condition and the right term only depends on the noise process.

Moreover we have $\Theta_{0}^{i}=\begin{pmatrix}\phi^{i}_{1}/(r_{i}^{-}-r_{i}^{+})\\ -\phi^{i}_{1}/(r_{i}^{-}-r_{i}^{+})\end{pmatrix}$. Thus

This is the bias term we have studied in Section 3.3 which we bound with Theorem 2. The variance term is controlled by the next proposition.

We note that if we restrict $\beta$ to $\beta\leq 3/(2h_{i})-\alpha/2$, then $4-(\alpha+2\beta)h_{i}\geq 1$ and the first bound of Proposition 2 is simplified to $\frac{2(\alpha n+\beta)^{2}}{\alpha\beta n^{2}}\frac{c_{i}}{h_{i}}$. This allows to conclude to prove Theorem 4.

F.2 Proof of Corollary 3

We let $\nu=\frac{\|\theta_{0}-\theta_{*}\|}{\sqrt{L\mathop{\rm tr}(CH^{-1})}}$ and consider three different regimes depending on $\nu$ and $L$.

If $\nu<1/L$, we have $\nu/N<1/L$ and thus $\alpha=\nu/N$ and $\beta=\nu$. Therefore

where we have used $\sqrt{L}\|\theta_{0}-\theta_{*}\|<\sqrt{\mathop{\rm tr}(CH^{-1})}$ since $\nu<1/L$.

If $\nu>1/L$ and $\nu<N/L$, we have $\alpha=\nu/N$ and $\beta=1/L$. Therefore

where we have used $\sqrt{L}\|\theta_{0}-\theta_{*}\|>\sqrt{\mathop{\rm tr}(CH^{-1})}$ since $\nu>1/L$.

If $\nu>N/L$, we have $\alpha=1/L$ and $\beta=1/L$. Therefore

where we have used that the real bound in Proposition 2 is in fact in $(N-1)\alpha+\beta$, (see Lemma 6) and that $\mathop{\rm tr}(CH^{-1})<\frac{L\|\theta_{0}-\theta_{*}\|^{2}}{N^{2}}$ since $\nu>N/L$.

F.3 Proof of Proposition 2

To prove Proposition 2 we will use Lemmas 6, 7 and 8, that are stated and proved in Section F.3.2.

Thus, by bouding $(k-1)\alpha+\beta$ by $(n-1)\alpha+\beta$, we get

Then, we have from Lemma 7 the inequality:

This allows to prove the first part of the bound. The other parts are much simpler and are done in Lemma 8. Thus, adding these bounds gives for $\alpha\leq 1/h_{i}$ and $0\leq\beta\leq 2/h_{i}-\alpha$:

F.3.2 Some technical Lemmas

We first compute an explicit expansion of the noise term as a function of the eigenvalues of the dynamical system.

For all $\alpha\leq 1/h_{i}$ and $0\leq\beta\leq 2/h_{i}-\alpha$ we have

We first turn the Euclidean norm into a trace, using that $\mathop{\rm tr}[AB]=\mathop{\rm tr}[BA]$ for two matrices $A$ and $B$ and that $\mathop{\rm tr}[x]=x$ for a real $x$.

And the first part of Eq. (23) is equal to:

because $D_{i}=\begin{pmatrix}r_{i}^{-}&0\\ 0&r_{i}^{+}\end{pmatrix}$ and $M_{i}=\begin{pmatrix}h_{i}^{1/2}&h_{i}^{1/2}\\ 0&0\end{pmatrix}$. Therefore:

In the following leamma, we bound a certain sum of powers of the roots.

For all $\alpha\leq 1/h_{i}$ and $0\leq\beta\leq 2/h_{i}-\alpha$ we have

We first note that when the two roots become close, the denominator and the numerator will go to zero, which prevents from bounding the numerator easily. We also note that this bound is very tight since the difference between the two terms goes to zero when $n$ goes to infinity.

We first expand the square of the difference of the powers of the roots and compute their sums.

with I_{n}=\bigg{[}\frac{{r_{i}^{+}}^{2n}}{1-{r_{i}^{+}}^{2}}+\frac{{r_{i}^{-}}^{2n}}{1-{r_{i}^{-}}^{2}}-2\frac{(r_{i}^{+}r_{i}^{-})^{n}}{1-(r_{i}^{+}r_{i}^{-})}\bigg{]}.

This sum is therefore equal to the sum of one term we will compute explicitly and one other term which will go to zero. We have for the first term:

with $J_{n}=\frac{I_{n}}{[(r_{i}^{-})-(r_{i}^{+})]^{2}}$.

Then we simplify the first term of this sum using the explicit values of the roots. We recall $r_{i}^{\pm}=r_{i}\pm\sqrt{\Delta_{i}}=1-\frac{\alpha+\beta}{2}h_{i}\pm\sqrt{\left(\frac{\alpha+\beta}{2}\right)^{2}h_{i}^{2}-\alpha h_{i}}$, therefore

Even if $J_{n}$ will be asymptotically small, we want a non-asymptotic bound, thus we will show that $J_{n}$ is always positive.

and using ${r_{i}^{+}}^{2}+{r_{i}^{-}}^{2}\geq 2{r_{i}^{-}}{r_{i}^{-}}$ we have

and using ${r_{i}^{+}}^{2}+{r_{i}^{-}}^{2}\leq 2{r_{i}^{-}}{r_{i}^{-}}$ we have

However we can also bound roughly Eq. (22) using Theorem 2 since we recall we have $\eta_{n}^{i}=\frac{[(r_{i}^{-})^{n-k}-(r_{i}^{+})^{n-k}]^{2}}{(r_{i}^{-}-r_{i}^{+})^{2}}$. This gives us the following lemma which enables to prove the second part of Proposition 2.

For all $\alpha\leq 1/h_{i}$ and $0\leq\beta\leq 2/h_{i}-\alpha$ we have

Appendix G Comparison with additional other algorithms

When the objective function $f$ is quadratic and for correct choices of step-sizes, the AC-SA algorithm of Lan , the SAGE algorithm of Hu et al. and the Accelerated RDA algorithm of Xiao are all equivalent to:

where we use $H_{n}\theta+\varepsilon_{n}$ as an unbiased estimate of the gradient and $\delta_{n}$ as step-size which values will be specified later.

Lan and Hu et al. only consider bounded cases by projecting their iterates on a bounded space. Xiao deals with the unbounded case and prove the following convergence result:

This result is significantly more general than ours since it is valid for composite optimization and general noise on the gradients.

We now present the different algorithms and show they all share the same form.

G.2 AC-SA

AC-SA algorithm with step size $\gamma_{n}$ and $\beta_{n}$ and gradient estimate $H_{n+1}\theta_{n}+\varepsilon_{n+1}$ is equivalent to:

Let the initial points $x_{1}^{ag}=x_{1}$, and the step-sizes $\{\beta_{n}\}_{n\leq 1}$ and $\{\gamma_{n}\}_{n\leq 1}$ be given.

Step 1. Set $x_{n}^{md}=\beta_{n}^{-1}x_{n}+(1-\beta_{n}^{-1})x_{n}^{ag}$,

Step 3. Set $n\rightarrow n+1$ and go to step 1.

When $f$ is quadratic we will have $G(x_{n}^{md},\xi_{n})=H_{n+1}x_{n}^{md}-\varepsilon_{n+1}$, thus $x_{n+1}=x_{n}-\gamma_{n}H_{n+1}x_{n}^{md}+\gamma_{n}\varepsilon_{n+1}$, and:

These give the result for $\theta_{n}=x_{n}^{ag}$. ∎

G.3 SAGE

The algorithm SAGE with step-sizes $L_{n}$ and $\alpha_{n}$ is equivalent to:

Let the initial points $x_{0}=z_{0}=0$, and the step-sizes $\{\beta_{n}\}_{n\leq 1}$ and $\{L_{n}\}_{n\leq 1}$ be given.

Step 1. Set $x_{n}=(1-\alpha_{n})y_{n-1}+\alpha_{n}z_{n-1}$,

Step 3. Set $n\rightarrow n+1$ and go to step 1.

These give the result for $\theta_{n}=y_{n}$. ∎

G.4 Accelerated RDA method

The algorithm AccRDA with step-sizes $\beta$ and $\alpha_{n}$ is equivalent to:

with $\gamma_{n}=\frac{\alpha_{n}\theta_{n}}{L+\beta}$.

We recall the general Accelerated RDA method:

Step 1. Set $A_{n}=A_{n-1}+\alpha_{n}$ and $\theta_{n}=\frac{\alpha_{n}}{A_{n}}$.

Step 2. Compute the query point $u_{n}=(1-\theta_{n})w_{n-1}+\theta_{n}v_{n-1}$

Step 5. Set $w_{n}=(1-\theta_{n})w_{n-1}+\theta_{n}v_{n}$.

Step 6. Set $n\rightarrow n+1$ and go to step 1.

With $\beta_{n}=\beta$ we have $v_{n}=v_{n-1}-\frac{\alpha_{n}}{L+\beta}(H_{n+1}u_{n}+\varepsilon_{n+1})]$ and

Since $v_{n-1}=\theta_{n-1}^{-1}w_{n-1}-\theta_{n-1}^{-1}(1-\theta_{n-1})w_{n-2}$, then

Appendix H Lower bound for stochastic optimization for least-squares

In this section, we show a lower bound for optimization of quadratic functions with noisy access to gradients. We follow very closely the framework of Agarwal et al. and use their notations. The only difference with their Theorem 1 in the different choice of two functions $f_{i}^{+}$ and $f_{i}^{-}$, which we choose to be:

with a non-increasing sequence $(c_{i})$ to be chosen later. The function $g_{\alpha}$ that is optimized is thus:

This function is quadratic and its Hessian has eigenvalues equal to $2c_{i}/d$. Thus, its largest eigenvalue is $2c_{1}/d$, which we choose equal to $L$.

Noisy gradients are obtained by sampling $d$ independent Bernoulli random variables $b_{i}$, $i=1,\dots,d$, with parameters $(\frac{1}{2}+\alpha_{i}\delta)$ and using the gradient of the random function \frac{1}{d}\sum_{i=1}^{d}\big{\{}b_{i}f_{i}^{+}(x)+(1-b_{i})f_{i}^{-}(x)\big{\}}. The variance of the random gradient is equal to

The function $g_{\alpha}$ is minimized for $x=-\alpha\delta r$, and the discrepancy measure between two functions $g_{\alpha}$ and $g_{\beta}$ is greater than

Since the vectors $\alpha,\beta\in\{-1,1\}^{d}$ are so that their Hamming distance $\Delta(\alpha,\beta)\geqslant d/4$ for $\alpha\neq\beta$, we have a discrepancy measure greater than $\frac{3c_{d}r^{2}\delta^{2}}{16}$. Thus, for a an approximate optimality of $\varepsilon=\frac{c_{d}r^{2}\delta^{2}}{38}$, we have, following the proof of Theorem 1 (equation (29)) from Agarwal et al. , for $N$ iterations of any method that accesses a random gradient, we have:

Thus, for $d$ large, we get, up to constants, $\delta^{2}\geqslant 1/N$ and thus $\varepsilon\geqslant\frac{r^{2}c_{d}}{N}$.

For $c_{1}=2Ld$ and $c_{i}=L\sqrt{d}$ for the remaining ones, we get (up to constants):

This leads to the desired result for $N\leqslant d$.