### Reducing the variance in online optimization by transporting past gradients

• 2019-06-08 22:02:28
• Sébastien M. R. Arnold, Pierre-Antoine Manzagol, Reza Babanezhad, Ioannis Mitliagkas, Nicolas Le Roux
• 24

# Reducing the variance in online optimization by transporting past gradients

Sébastien M. R. Arnold
University of Southern California
Los Angeles, CA
[email protected]
&Pierre-Antoine Manzagol
Montréal, QC
[email protected]
University of British Columbia
Vancouver, BC
[email protected]
&Ioannis Mitliagkas
Mila, Université de Montréal
Montréal, QC
[email protected]
&Nicolas Le Roux
Montréal, QC
[email protected]
Work done while at Mila.
###### Abstract

Reducing the variance in online optimization by transporting past gradients

Sébastien M. R. Arnold thanks: Work done while at Mila. University of Southern California Los Angeles, CA Pierre-Antoine Manzagol Google Brain Montréal, QC Reza Babanezhad University of British Columbia Vancouver, BC Ioannis Mitliagkas Mila, Université de Montréal Montréal, QC Nicolas Le Roux Mila, Google Brain Montréal, QC

\@float

noticebox[b]Preprint. Under review.

## 1 Introduction

We wish to solve the following minimization problem:

 $\displaystyle\theta^{\ast}$ $\displaystyle=\arg\min_{\theta}E_{x\sim p}[f(\theta,x)]\;,$ (1)

where we only have access to samples $x$ and to a first-order oracle that gives us, for a given $\theta$ and a given $x$, the derivative of $f(\theta,x)$ with respect to $\theta$, i.e. $\frac{\partial f(\theta,x)}{\partial\theta}=g(\theta,x)$. It is known robbins1951stochastic that, when $f$ is smooth and strongly convex, there is a converging algorithm for Problem 1 that takes the form $\theta_{t+1}=\theta_{t}-\alpha_{t}g(\theta_{t},x_{t})$, where $x_{t}$ is a sample from $p$. This algorithm, dubbed stochastic gradient (SG), has a convergence rate of $O(1/t)$ (see for instance bubeck2015convex ), within a constant factor of the minimax rate for this problem. When one has access to the true gradient $g(\theta)=E_{x\sim p}[g(\theta,x)]$ rather than just a sample, this rate dramatically improves to $O(e^{-\nu t})$ for some $\nu>0$.

In addition to hurting the convergence speed, noise in the gradient makes optimization algorithms harder to tune. Indeed, while full gradient descent is convergent for constant stepsize $\alpha$, and also amenable to line searches to find a good value for that stepsize, the stochastic gradient method from robbins1951stochastic with a constant stepsize only converges to a ball around the optimum schmidt2014convergence .22 2 Under some conditions, it does converge linearly to the optimum (e.g., Vaswani19, ) Thus, to achieve convergence, one needs to use a decreasing stepsize. While this seems like a simple modification, the precise decrease schedule can have a dramatic impact on the convergence speed. While theory prescribes $\alpha_{t}=O(t^{-\alpha})$ with $\alpha\in(1/2,1]$ in the smooth case, practictioners often use larger stepsizes like $\alpha_{t}=O(t^{-1/2})$ or even constant stepsizes.

When the distribution $p$ has finite support, Eq. 1 becomes a finite sum and, in that setting, it is possible to achieve efficient variance reduction and drive the noise to zero, allowing stochastic methods to achieve linear convergence rates leroux2012stochastic ; johnson2013accelerating ; zhang2013linear ; mairal2013optimization ; Shalev-Shwartz2013sdca ; defazio2014saga . Unfortunately, the finite support assumption is critical to these algorithms which, while valid in many contexts, does not have the broad applicability of the standard SG algorithm. Several works have extended these approaches to the online setting by applying these algorithms while increasing $N$ babanezhad2015wasting ; hofmann2015variance but they need to revisit past examples mutiple times and are not truly online.

Another line of work reduces variance by averaging iterates polyak1992acceleration ; lacoste2012simpler ; bach2013non ; flammarion2015averaging ; dieuleveut2017harder ; dieuleveut2017bridging ; JMLR:v18:16-595 . While these methods converge for a constant stepsize in the stochastic case33 3 Under some conditions on $f$., their practical speed is heavily dependent on the fraction of iterates kept in the averaging, a hyperparameter that is thus hard to tune, and they are rarely used in deep learning.

Our work combines two existing ideas and adds a third: a) At every step, it updates the parameters using a weighted average of past gradients, like in SAG (leroux2012stochastic, ; schmidt2017minimizing, ), albeit with a different weighting scheme; b) It reduces the bias and variance induced by the use of these old gradients by transporting them to “equivalent” gradients computed at the current point, similar to gower2017tracking ; c) It does so implicitly by computing the gradient at a parameter value different from the current one. The resulting gradient estimator can then be used as a plug-in replacement of the stochastic gradient within any optimization scheme. Experimentally, both SG using our estimator and its momentum variant outperform the most commonly used optimizers in deep learning.

## 2 Momentum and other approaches to dealing with variance

Stochastic variance reduction methods use an average of past gradients to reduce the variance of the gradient estimate. At first glance, it seems like their updates are similar to that of momentum polyak1964some , also known as the heavy ball method, which performs the following updates44 4 This is slightly different from the standard formulation but equivalent for constant $\gamma_{t}$.:

 $\displaystyle v_{t}$ $\displaystyle=\gamma_{t}v_{t-1}+(1-\gamma_{t})g(\theta_{t},x_{t}),\qquad v_{0}% =g(\theta_{0},x_{0})$ $\displaystyle\theta_{t+1}$ $\displaystyle=\theta_{t}-\alpha_{t}v_{t}\;.$

When $\gamma_{t}=\gamma$, this leads to $\displaystyle\theta_{t+1}=\theta_{t}-\alpha_{t}\left(\gamma^{t}g(\theta_{0},x_% {0})+(1-\gamma)\sum_{i=1}^{t}\gamma^{t-i}g(\theta_{i},x_{i})\right)$. Hence, the heavy ball method updates the parameters of the model using an average of past gradients, bearing similarity with SAG leroux2012stochastic , albeit with exponential instead of uniform weights.

Interestingly, while momentum is a popular method for training deep networks, its theoretical analysis in the stochastic setting is limited (sutskever2013importance, ), except in the particular setting when the noise converges to 0 at the optimum (loizou2017momentum, ). Also surprising is that, despite the apparent similarity with stochastic variance reduction methods, current convergence rates are slower when using $\gamma>0$ in the presence of noise Schmidt11_inexact , although this might be a limitation of the analysis.

### 2.1 Momentum and variance

We propose here an analysis of how, on quadratics, using past gradients as done in momentum does not lead to a decrease in variance. If gradients are stochastic, then $\Delta_{t}=\theta_{t}-\theta^{\ast}$ is a random variable. Denoting $\epsilon_{i}$ the noise at timestep $i$, i.e. $\displaystyle g(\theta_{i},x_{i})=g(\theta_{i})+\epsilon_{i}$, and writing $\Delta_{t}-E[\Delta_{t}]=\alpha\sum_{i=0}^{t}N_{i,t}\epsilon_{i}$, with $N_{i,t}$ the impact of the noise of the $i$-th datapoint on the $t$-th iterate, we may now analyze the total impact of each $\epsilon_{i}$ on the iterates. Figure 1 shows the impact of $\epsilon_{i}$ on $\Delta_{t}-E[\Delta_{t}]$ as measured by $N_{i,t}^{2}$ for three datapoints ($i=1$, $i=25$ and $i=50$) as a function of $t$ for stochastic gradient ($\gamma=0$, left) and momentum ($\gamma=0.9$, right). As we can see, when using momentum, the variance due to a given datapoint first increases as the noise influences both the next iterate (through the parameter update) and the subsequent updates (through the velocity). Due to the weight $1-\gamma$ when a point is first sampled, a larger value of $\gamma$ leads to a lower immediate impact of the noise of a given point on the iterates. However, a larger $\gamma$ also means that the noise of a given gradient is kept longer, leading to little or no decrease of the total variance (dashed blue curve). Even in the case of stochastic gradient, the noise at a given timestep carries over to subsequent timesteps, even if the old gradients are not used for the update, as the iterate itself depends on the noise.

At every timestep, the contribution to the noise of the 1st, the 25th and the 50th points in Fig. 1 is unequal. If we assume that the $\epsilon_{i}$ are i.i.d., then the total variance would be minimal if the contribution from each point was equal. Further, one can notice that the impact of datapoint $i$ is only a function of $t-i$ and not of $t$. This guarantees that the total noise will not decrease over time.

To address these two points, one can increase the momentum parameter over time. In doing so, the noise of new datapoints will have a decreasing impact on the total variance as their gradient is multiplied by $1-\gamma_{t}$. Figure 0(c) shows the impact $N_{i,t}^{2}$ of each noise $\epsilon_{i}$ for an increasing momentum $\gamma_{t}=1-\frac{1}{t}$. The peak of noise for $i=25$ is indeed lower than that of $i=1$. However, the variance still does not go to 0. This is because, as the momentum parameter increases, the update is an average of many gradients, including stale ones. Since these gradients were computed at iterates already influenced by the noise over previous datapoints, that past noise is amplified, as testified by the higher peak at $i=1$ for the increasing momentum. Ultimately, increasing momentum does not lead to a convergent algorithm in the presence of noise when using a constant stepsize.

### 2.2 SAG and Hessian modelling

The impact of the staleness of the gradients on the convergence is not limited to momentum. In SAG, for instance, the excess error after $k$ updates is proportional to $\left(1-\min\left\{\frac{1}{16\widehat{\kappa}},\frac{1}{8N}\right\}\right)^{k}$, compared to the excess error of the full gradient method which is $\left(1-\frac{1}{\kappa}\right)^{k}$ where $\kappa$ is the condition number of the problem. 55 5 The $\widehat{\kappa}$ in the convergence rate of SAG is generally larger than the $\kappa$ in the full gradient algorithm. The difference between the two rates is larger when the minimum in the SAG rate is the second term. This happens either when $\widehat{\kappa}$ is small, i.e. the problem is well conditioned and a lot of progress is made at each step, or when $N$ is large, i.e. there are many points to the training set. Both cases imply that a large distance has been travalled between two draws of the same datapoint.

Recent works showed that correcting for that staleness by modelling the Hessian wai2017curvature ; gower2017tracking leads to improved convergence. As momentum uses stale gradients, the velocity is an average of current and past gradients and thus can be seen as an estimate of the true gradient at a point which is not the current one but rather a convex combination of past iterates. As past iterates depend on the noise of previous gradients, this bias in the gradients amplifies the noise and leads to a non-converging algorithm. We shall thus “transport” the old stochastic gradients $g(\theta_{i},x_{i})$ to make them closer to their corresponding value at the current iterate, $g(\theta_{t},x_{i})$. Past works did so using the Hessian or an explicit approximation thereof, which can be expensive and difficult to compute and maintain. We will resort to using implicit transport, a new method that aims at compensating the staleness of past gradients without making explicit use of the Hessian.

## 3 Converging optimization through implicit gradient transport

Before showing how to combine the advantages of both increasing momentum and gradient transport, we demonstrate how to transport gradients implicitly. This transport is only exact under a strong assumption that will not hold in practice. However, this result will serve to convey the intuition behind implicit gradient transport. We will show in Section 4 how to mitigate the effect of the unsatisfied assumption.

### 3.1 Implicit gradient transport

Let us assume that we received samples $x_{0},\ldots,x_{t}$ in an online fashion. We wish to approach the full gradient $g_{t}(\theta_{t})=\frac{1}{t+1}\sum_{i=0}^{t}g(\theta_{t},x_{i})$ as accurately as possible. We also assume here that a) We have a noisy estimate $\widehat{g}_{t-1}(\theta_{t-1})$ of $g_{t-1}(\theta_{t-1})$; b) We can compute the gradient $g(\theta,x_{t})$ at any location $\theta$. We shall seek a $\theta$ such that

 $\displaystyle\frac{t}{t+1}\widehat{g}_{t-1}(\theta_{t-1})+\frac{1}{t+1}g(% \theta,x_{t})\approx g_{t}(\theta_{t})\;.$

To this end, we shall make the following assumption:

###### Assumption 3.1.

All individual functions $f(\cdot,x)$ are quadratics with the same Hessian $H$.

This is the same assumption as (flammarion2015averaging, , Section 4.1). Although it is unlikely to hold in practice, we shall see that our method still performs well when that assumption is violated.

Under Assumption 3.1, we then have (see details in Appendix)

 $\displaystyle g_{t}(\theta_{t})$ $\displaystyle=\frac{t}{t+1}g_{t-1}(\theta_{t})+\frac{1}{t+1}g(\theta_{t},x_{t})$ $\displaystyle\approx\frac{t}{t+1}\widehat{g}_{t-1}(\theta_{t-1})+\frac{1}{t+1}% g(\theta_{t}+t(\theta_{t}-\theta_{t-1}),x_{t})\;.$

Thus, we can transport our current estimate of the gradient by computing the gradient on the new point at a shifted location $\theta=\theta_{t}+t(\theta_{t}-\theta_{t-1})$. This extrapolation step is reminiscent of Nesterov’s acceleration with the difference that the factor in front of $\theta_{t}-\theta_{t-1}$, $t$, is not bounded.

### 3.2 Combining increasing momentum and implicit gradient transport

We now describe our main algorithm, Implicit Gradient Transport (IGT). IGT uses an increasing momentum $\gamma_{t}=\frac{t}{t+1}$. At each step, when updating the velocity, it computes the gradient of the new point at an extrapolated location so that the velocity $v_{t}$ is a good estimate of the true gradient $g(\theta_{t})$.

We can rewrite the updates to eliminate the velocity $v_{t}$, leading to the update:

 $\displaystyle\theta_{t+1}$ $\displaystyle=\frac{2t+1}{t+1}\theta_{t}-\frac{t}{t+1}\theta_{t-1}-\frac{% \alpha}{t+1}g\left(\theta_{t}+t(\theta_{t}-\theta_{t-1}),x_{t}\right)\;.$ (IGT)

We see in Fig. 0(d) that IGT allows a reduction in the total variance, thus leading to convergence with a constant stepsize. This is captured by the following proposition:

###### Proposition 3.1.

If $f$ is a quadratic function with positive definite Hessian $H$ with largest eigenvalue $L$ and condition number $\kappa$ and if the stochastic gradients satisfy: $g(\theta,x)=g(\theta)+\epsilon$ with $\epsilon$ a random i.i.d. noise with covariance bounded by $BI$, then Eq. IGT with stepsize $\alpha=1/L$ leads to iterates $\theta_{t}$ satisfying

 $\displaystyle E[\|\theta_{t}-\theta^{\ast}\|^{2}]$ $\displaystyle\leq\left(1-\frac{1}{\kappa}\right)^{2t}\|\theta_{0}-\theta^{\ast% }\|^{2}+\frac{d\alpha^{2}B\bar{\nu}_{0}^{2}}{t}\;,$

with $\nu=(2+2\log\kappa)\kappa$ for every $t>2\kappa$.

The proof of Prop. 3.1 is provided in the appendix.

Despite this theoretical result, two limitations remain: First, Prop. 3.1 shows that IGT does not improve the dependency on the conditioning of the problem; Second, the assumption of equal Hessians is unlikely to be true in practice, leading to an underestimation of the bias. We address the conditioning issue in the next section and the assumption on the Hessians in Section 4.

### 3.3 IGT as a plug-in gradient estimator

We demonstrated that the IGT estimator has lower variance than the stochastic gradient estimator for quadratic objectives. IGT can also be used as a drop-in replacement for the stochastic gradient in an existing, popular first order method: the heavy ball (HB). This is captured by the following two propositions:

###### Proposition 3.2 (Non-stochastic).

In the non-stochastic case, where $B=0$, variance is equal to $0$ and Heavyball-IGT achieves the accelerated linear rate $O\big{(}\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{t}\big{)}$ using the known, optimal heavy ball tuning, $\mu=\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{2}$, $\alpha=(1+\sqrt{\mu})^{2}/L$.

###### Proposition 3.3 (Online, stochastic).

When $B>0$, there exist constant hyperparameters $\alpha>0$, $\mu>0$ such that $\|E[\theta_{t}-\theta^{\ast}]\|^{2}$ converges to zero linearly, and the variance is $\tilde{O}(1/t)$.

The pseudo-code can be found in Algorithm 3.3.

{algorithm}

Heavyball-IGT {algorithmic} \ProcedureHeavyball-IGTStepsize $\alpha$, Momentum $\mu$, Initial parameters $\theta_{0}$ \State$v_{0}\leftarrow g(\theta_{0},x_{0})\quad,\quad w_{0}\leftarrow-\alpha v_{0}% \quad,\quad\theta_{1}\leftarrow\theta_{0}+w_{0}$ \For$t=1,\ldots,T-1$ \State$\gamma_{t}\leftarrow\frac{t}{t+1}$ \State$v_{t}\leftarrow\gamma_{t}v_{t-1}+(1-\gamma_{t})g\left(\theta_{t}+\frac{\gamma_% {t}}{1-\gamma_{t}}(\theta_{t}-\theta_{t-1}),x_{t}\right)$ \State$w_{t}\leftarrow\mu w_{t-1}-\alpha v_{t}$ \State$\theta_{t+1}\leftarrow\theta_{t}+w_{t}$ \EndFor \Statereturn $\theta_{T}$ \EndProcedure

## 4 IGT and Anytime Tail Averaging

So far, IGT weighs all gradients equally. This is because, with equal Hessians, one can perfectly transport these gradients irrespective of the distance travelled since they were computed. In practice, the individual Hessians are not equal and might change over time. In that setting, the transport induces an error which grows with the distance travelled. We wish to average a linearly increasing number of gradients, to maintain the $O(1/t)$ rate on the variance, while forgetting about the oldest gradients to decrease the bias. To this end, we shall use anytime tail averaging leroux2019anytime , named in reference to the tail averaging technique used in optimization JMLR:v18:16-595 .

Tail averaging is an online averaging technique where only the last points, usually a constant fraction $c$ of the total number of points seen, is kept. Maintaining the exact average at every timestep is memory inefficient and anytime tail averaging performs an approximate averaging using $\gamma_{t}=\frac{c(t-1)}{1+c(t-1)}\left(1-\frac{1}{c}\sqrt{\frac{1-c}{t(t-1)}}\right)$. We refer the reader to leroux2019anytime for additional details.

## 5 Impact of IGT on bias and variance in the ideal case

To understand the behaviour of IGT when Assumption 3.1 is verified, we minimize a strongly convex quadratic function with Hessian $Q\in\mathbb{R}^{100\times 100}$ with condition number $1000$, and we have access to the gradient corrupted by noise $\epsilon_{t}$, where $\epsilon_{t}\sim N(0,0.3\cdot I_{100})$. In that scenario where all Hessians are equal and implicit gradient transport is exact, Fig. 1(a) confirms the $O(1/t)$ rate of IGT with constant stepsize while SGD and HB only converge to a ball around the optimum.

To further understand the impact of IGT, we study the quality of the gradient estimate. Standard stochastic methods control the variance of the parameter update by scaling it with a decreasing stepsize, which slows the optimization down. With IGT, we hope to have a low variance while maintaining a norm of the update comparable to that obtained with gradient descent. To validate the quality of our estimator, we optimized a quadratic function using IGT, collecting iterates $\theta_{t}$. For each iterate, we computed the squared error between the true gradient and either the stochastic or the IGT gradient. In this case where both estimators are unbiased, this is the trace of the noise covariance of our estimators. The results in Figure 1(b) show that, as expected, this noise decreases linearly for IGT and is constant for SGD.

We also analyse the direction and magnitude of the gradient of IGT on the same quadratic setup. Figure 1(c) displays the cosine similarity between the true gradient and either the stochastic or the IGT gradient, as a function of the distance to the optimum. We see that, for the same distance, the IGT gradient is much more aligned with the true gradient than the stochastic gradient is, confirming that variance reduction happens without the need for scaling the estimate.

## 6 Experiments

While Section 5 confirms the performance of IGT in the ideal case, the assumption of identical Hessians almost never holds in practice. In this section, we present results on more realistic and larger scale machine learning settings. All experiments are extensively described in the Appendix A and additional baselines compared in Appendix B.

### 6.1 Supervised learning

##### CIFAR10 image classification

We first consider the task of training a ResNet-56 model DBLP:journals/corr/HeZRS15 on the CIFAR-10 image classification dataset cifar10 . We use TF official models code and setup TfResnet , varying only the optimizer: SGD, HB, Adam and our algorithm with anytime tail averaging both on its own (ITA) and combined with Heavy Ball (HB-ITA). We tuned the step size for each algorithm by running experiments using a logarithmic grid. To factor in ease of tuning wilson2017marginal , we used Adam’s default parameter values and a value of 0.9 for HB’s parameter. We used a linearly decreasing stepsize as it was shown to be simple and perform well shallue2018measuring . For each optimizer we selected the hyperparameter combination that is fastest to reach a consistently attainable target train loss shallue2018measuring . Selecting the hyperparameter combination reaching the lowest training loss yields qualitatively identical curves. Figure 3 presents the results, showing that IGT with the exponential anytime tail average performs favourably, both on its own and combined with Heavy Ball: the learning curves show faster improvement and are much less noisy. Figure 3: Resnet-56 on CIFAR10. Left: Train loss. Center: Train accuracy. Right: Test accuracy.
##### ImageNet image classification

We also consider the task of training a ResNet-50 modelDBLP:journals/corr/HeZRS15 on the larger ImageNet dataset ILSVRC15 . The setup is similar to the one used for CIFAR10 with the difference that we trained using larger minibatches (1024 instead of 128). In Figure 4, one can see that IGT is as fast as Adam for the train loss, faster for the train accuracy and reaches the same final performance, which Adam does not. We do not see the noise reduction we observed with CIFAR10, which could be explained by the larger batch size (see Appendix A.1). Figure 4: ResNet-50 on ImageNet. Left: Train loss. Center: Train accuracy. Right: Test accuracy.
##### IMDb sentiment analysis

We train a bi-directional LSTM on the IMDb Large Movie Review Dataset for 200 epochs. maas2011learning We observe that while the training convergence is comparable to HB, HB-ITA performs better in terms of validation and test accuracy. In addition to the baseline and IGT methods, we also train a variant of Adam using the ITA gradients, dubbed Adam-ITA, which performs similarly to Adam. Figure 5: Validation curves for different large-scale machine learning settings. Shading indicates one standard deviation computed over three random seeds. Left: Reinforcement learning via policy gradient on a LQR system. Right: Meta-learning using MAML on Mini-Imagenet.

### 6.2 Reinforcement learning

We cast the classical linear-quadratic regulator (LQR) kwakernaak1972linear as a policy learning problem to be optimized via gradient descent. This setting is extensively described in Appendix A. Note that despite their simple linear dynamics and a quadratic cost functional, LQR systems are notoriously difficult to optimize due to the non-convexity of the loss landscape. fazel2018global

The left chart in Figure 5 displays the evaluation cost computed along training and averaged over three random seeds. The first method (Optimal) indicates the cost attained when solving the algebraic Riccati equation of the LQR – this is the optimal solution of the problem. SGD minimizes the costs using the REINFORCE williams1992simple gradient estimator, averaged over 600 trajectories. ITA is similar to SGD but uses the ITA gradient computed from the REINFORCE estimates. Finally, GD uses the analytical gradient by taking the expectation over the policy.

We make two observations from the above chart. First, ITA initially suffers from the stochastic gradient estimate but rapidly matches the performance of GD. Notably, both of them converge to a solution significantly better than SGD, demonstrating the effectiveness of the variance reduction mechanism. Second, the convergence curve is smoother for ITA than for SGD, indicating that the ITA iterates are more likely to induce similar policies from one iteration to the next. This property is particularly desirable in reinforcement learning as demonstrated by the popularity of trust-region methods in large-scale applications. schulman2015trust ; openai2018learning

### 6.3 Meta-learning

##### Model-agnostic meta-learning

We now investigate the use of IGT in the model-agnostic meta-learning (MAML) setting. finn2017model We replicate the 5 ways classification setup with 5 adaptation steps on tasks from the Mini-Imagenet dataset ravi2016optimization . This setting is interesting because of the many sources contributing to noise in the gradient estimates: the stochastic meta-gradient depends on the product of 5 stochastic Hessians computed over only 10 data samples, and is averaged over only 4 tasks. We substitute the meta-optimizer with each method, select the stepsize that maximizes the validation accuracy after 10K iterations, and use it to train the model for 100K iterations.

The right graph of Figure 5 compares validation accuracies for three random seeds. We observe that methods from the IGT family significantly outperform their stochastic meta-gradient counter-part, both in terms of convergence rate and final accuracy. Those results are also reflected in the final test accuracies where Adam-ITA ($65.16\%$) performs best, followed by HB-ITA ($64.57\%$), then Adam ($63.70\%$), and finally HB ($63.08\%$).

## 7 Conclusion and open questions

We proposed a simple optimizer which, by reusing past gradients and transporting them, offers excellent performance on a variety of problems. While it adds an additional parameter, the ratio of examples to be kept in the tail averaging, it remains competitive across a wide range of such values. Further, by providing a higher quality gradient estimate that can be plugged in any existing optimizer, we expect it to be applicable to a wide range of problems. As the IGT is similar to momentum, this further raises the question on the links between variance reduction and curvature adaptation. Whether there is a way to combine the two without using momentum on top of IGT remains to be seen.

#### Acknowledgments

The authors would like to thank Liyu Chen for his help with the LQR experiments and Fabian Pedregosa for insightful discussions.

## References

•  The TensorFlow Authors. Tensorflow official resnet model. 2018.
•  Reza Babanezhad, Mohamed Osama Ahmed, Alim Virani, Mark Schmidt, Jakub Konec̆ný, and Scott Sallinen. Stop wasting my gradients: Practical SVRG. In Advances in Neural Information Processing Systems, 2015.
•  Francis Bach and Eric Moulines. Non-strongly-convex smooth stochastic approximation with convergence rate o (1/n). In Advances in Neural Information Processing Systems, pages 773–781, 2013.
•  Sébastien Bubeck. Convex optimization: Algorithms and complexity. Foundations and Trends® in Machine Learning, 8(3-4):231–357, 2015.
•  Aaron Defazio, Francis Bach, and Simon Lacoste-Julien. Saga: A fast incremental gradient method with support for non-strongly convex composite objectives. In Advances in Neural Information Processing Systems, pages 1646–1654, 2014.
•  Aymeric Dieuleveut, Alain Durmus, and Francis Bach. Bridging the gap between constant step size stochastic gradient descent and markov chains. arXiv preprint arXiv:1707.06386, 2017.
•  Aymeric Dieuleveut, Nicolas Flammarion, and Francis Bach. Harder, better, faster, stronger convergence rates for least-squares regression. The Journal of Machine Learning Research, 18(1):3520–3570, 2017.
•  Maryam Fazel, Rong Ge, Sham M. Kakade, and Mehran Mesbahi. Global convergence of policy gradient methods for linearized control problems. CoRR, abs/1801.05039, 2018.
•  Chelsea Finn, Pieter Abbeel, and Sergey Levine. Model-agnostic meta-learning for fast adaptation of deep networks. In Proceedings of the 34th International Conference on Machine Learning-Volume 70, pages 1126–1135. JMLR. org, 2017.
•  Nicolas Flammarion and Francis Bach. From averaging to acceleration, there is only a step-size. In Conference on Learning Theory, pages 658–695, 2015.
•  Robert M Gower, Nicolas Le Roux, and Francis Bach. Tracking the gradients using the hessian: A new look at variance reducing stochastic methods. arXiv preprint arXiv:1710.07462, 2017.
•  Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. CoRR, abs/1512.03385, 2015.
•  Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Identity mappings in deep residual networks. In European conference on computer vision, pages 630–645. Springer, 2016.
•  Thomas Hofmann, Aurelien Lucchi, Simon Lacoste-Julien, and Brian McWilliams. Variance reduced stochastic gradient descent with neighbors. In Advances in Neural Information Processing Systems, pages 2305–2313, 2015.
•  Prateek Jain, Sham M Kakade, Rahul Kidambi, Praneeth Netrapalli, and Aaron Sidford. Accelerating stochastic gradient descent for least squares regression. In Conference On Learning Theory, pages 545–604, 2018.
•  Prateek Jain, Sham M. Kakade, Rahul Kidambi, Praneeth Netrapalli, and Aaron Sidford. Parallelizing stochastic gradient descent for least squares regression: Mini-batching, averaging, and model misspecification. Journal of Machine Learning Research, 18(223):1–42, 2018.
•  Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In Advances in Neural Information Processing Systems, pages 315–323, 2013.
•  Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
•  Alex Krizhevsky. Learning multiple layers of features from tiny images. Technical report, 2009.
•  Alex Krizhevsky. Learning multiple layers of features from tiny images. Technical report, Citeseer, 2009.
•  Huibert Kwakernaak. Linear optimal control systems, volume 1. 1972.
•  Simon Lacoste-Julien, Mark Schmidt, and Francis Bach. A simpler approach to obtaining an o (1/t) convergence rate for the projected stochastic subgradient method. arXiv preprint arXiv:1212.2002, 2012.
•  Nicolas Le Roux. Anytime tail averaging. arXiv preprint arXiv:1902.05083, 2019.
•  Nicolas Le Roux, Mark Schmidt, and Francis Bach. A stochastic gradient method with an exponential convergence rate for finite training sets. In Advances in Neural Information Processing Systems, pages 2663–2671, 2012.
•  Yann LeCun and Corinna Cortes. Mnist handwritten digit database. AT&T Labs [Online]. Available: http://yann. lecun. com/exdb/mnist, 2010.
•  Nicolas Loizou and Peter Richtárik. Momentum and stochastic momentum for stochastic gradient, newton, proximal point and subspace descent methods. arXiv preprint arXiv:1712.09677, 2017.
•  Andrew L. Maas, Raymond E. Daly, Peter T. Pham, Dan Huang, Andrew Y. Ng, and Christopher Potts. Learning word vectors for sentiment analysis. In Proceedings of the 49th Annual Meeting of the Association for Computational Linguistics: Human Language Technologies, pages 142–150, Portland, Oregon, USA, June 2011. Association for Computational Linguistics.
•  Julien Mairal. Optimization with first-order surrogate functions. In Proceedings of The 30th International Conference on Machine Learning, pages 783–791, 2013.
•  OpenAI, Marcin Andrychowicz, Bowen Baker, Maciek Chociej, Rafal Józefowicz, Bob McGrew, Jakub W. Pachocki, Jakub Pachocki, Arthur Petron, Matthias Plappert, Glenn Powell, Alex Ray, Jonas Schneider, Szymon Sidor, Josh Tobin, Peter Welinder, Lilian Weng, and Wojciech Zaremba. Learning dexterous in-hand manipulation. CoRR, abs/1808.00177, 2018.
•  Brendan O’Donoghue and Emmanuel Candes. Adaptive restart for accelerated gradient schemes. Foundations of computational mathematics, 15(3):715–732, 2015.
•  Jeffrey Pennington, Richard Socher, and Christopher Manning. Glove: Global vectors for word representation. In Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP), pages 1532–1543, 2014.
•  Boris T Polyak. Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics, 4(5):1–17, 1964.
•  Boris T Polyak and Anatoli B Juditsky. Acceleration of stochastic approximation by averaging. SIAM Journal on Control and Optimization, 30(4):838–855, 1992.
•  Sachin Ravi and Hugo Larochelle. Optimization as a model for few-shot learning. 2016.
•  Herbert Robbins and Sutton Monro. A stochastic approximation method. Annals of Mathematical Statistics, 22(3):400–407, 1951.
•  Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, Alexander C. Berg, and Li Fei-Fei. ImageNet Large Scale Visual Recognition Challenge. International Journal of Computer Vision (IJCV), 115(3):211–252, 2015.
•  Mark Sandler, Andrew Howard, Menglong Zhu, Andrey Zhmoginov, and Liang-Chieh Chen. Mobilenetv2: Inverted residuals and linear bottlenecks. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 4510–4520, 2018.
•  Mark Schmidt. Convergence rate of stochastic gradient with constant step size. 2014.
•  Mark Schmidt, Nicolas Le Roux, and Francis Bach. Convergence rates of inexact proximal-gradient methods for convex optimization. In Advances in Neural Information Processing Systems 24, 2011.
•  Mark Schmidt, Nicolas Le Roux, and Francis Bach. Minimizing finite sums with the stochastic average gradient. Mathematical Programming, 162(1-2):83–112, 2017.
•  John Schulman, Sergey Levine, Philipp Moritz, Michael I. Jordan, and Pieter Abbeel. Trust region policy optimization. CoRR, abs/1502.05477, 2015.
•  Shai Shalev-Shwartz and Tong Zhang. Stochastic dual coordinate ascent methods for regularized loss. Journal of Machine Learning Research, 14(1):567–599, February 2013.
•  Christopher J Shallue, Jaehoon Lee, Joe Antognini, Jascha Sohl-Dickstein, Roy Frostig, and George E Dahl. Measuring the effects of data parallelism on neural network training. arXiv preprint arXiv:1811.03600, 2018.
•  Ilya Sutskever, James Martens, George Dahl, and Geoffrey Hinton. On the importance of initialization and momentum in deep learning. In International conference on machine learning, pages 1139–1147, 2013.
•  Sharan Vaswani, Francis Bach, and Mark Schmidt. Fast and faster convergence of sgd for over-parameterized models and an accelerated perceptron. In Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics, 2019.
•  Hoi-To Wai, Wei Shi, Angelia Nedic, and Anna Scaglione. Curvature-aided incremental aggregated gradient method. arXiv preprint arXiv:1710.08936, 2017.
•  Ronald J Williams. Simple statistical gradient-following algorithms for connectionist reinforcement learning. Machine learning, 8(3-4):229–256, 1992.
•  Ashia C Wilson, Rebecca Roelofs, Mitchell Stern, Nathan Srebro, and Benjamin Recht. The marginal value of adaptive gradient methods in machine learning. arXiv preprint arXiv:1705.08292, 2017.
•  Jian Zhang and Ioannis Mitliagkas. Yellowfin and the art of momentum tuning. In SysML, 2019.
•  Lijun Zhang, Mehrdad Mahdavi, and Rong Jin. Linear convergence with condition number independent access of full gradients. In Advances in Neural Information Processing Systems, pages 980–988, 2013.

## Appendix A Experimental Details

This section provides additional information regarding the experiments included in the main text.

For each experimental setting we strive to follow the reproducibility checklist, and provide:

• a description and citation of the dataset,

• a description of pre-processing steps,

• training / validation / testing splits,

• a description of the hyper-parameter search process and chosen values for each method,

• the exact number of evaluation runs,

• a clear definition of statistics reported, and

• a description of the computing infrastructure.

### A.1 CIFAR10 image classification

##### Dataset

The CIFAR10 dataset  consists 50k training and 10k testing images, partitioned over 10 classes. We download and pre-process the images using the TensorFlow models package, available at the following URL: https://github.com/tensorflow/models

##### Model

We use a residual convolutional network  with 56 layers as defined in the models package. Specifically, we use the second version whose blocks are built as a batch normalization, then a ReLU activation, and then a convolutional layer. 

##### Hyper-parameters

We use the exact setup from https://github.com/tensorflow/models/officials/resnet. As such, training is carried out with minibatches of 128 examples for 182 epochs and the training data is augmented with random crops and horizontal flips. Also note this setup multiplies the step size by the size of the minibatch. One deviation from the setup is our use of a linearly decaying learning rate instead of an explicit schedule. The linearly decaying learning rate schedule is simple and was shown to perform well . This schedule is specified using two parameters: the decay rate, a multiplier specifying the final step size (0.1 or 0.01), and the decay step, specifying the step at which the fully decayed rate is reached (always set to 90% of the training steps). To factor in ease of tuning we used Adam’s default parameter values and a value of 0.9 for HB’s parameter. We used IGT with the exponential Anytime Tail Averaging approach. For the tail fraction, we tried two values: the number of epochs and a tenth of that number (180 and 18). We ran using the following learning rate: (1e0, 3e-1, 1e-1, 3e-2, 1e-2) for SGD, HB and the IGT variants and (1e-2, 3e-3, 1e-3, 3e-4, 1e-4) for Adam. We ran a grid search over the base learning rate and its decay rate with a single run per combination. For each optimizer we selected the hyperparameter combination that is fastest to reach a consistently attainable target train loss of 0.2 . Note that selecting the hyperparameter combination reaching the lowest training loss yields qualitatively identical curves.

The resulting hyper-parameters are:

• SGD stepsize 0.3, decay 0.01

• HB stepsize 0.03, decay 0.01

• Adam stepsize 0.001, decay 0.01

• ITA stepsize 0.3, decay 0.01, tail fraction 18

• HB-ITA stepsize 0.03, decay 0.1, tail fraction 18

##### Infrastructure and Runs

The experiments were run using P100 GPUs (single GPU).

We provide all learning curves for the methods comparison presented in the main manuscript in figure 6. Figure 6: Convergence and accuracy curves along training for the CIFAR10 experiments comparing baseline methods to ours. Left: Training. Right: Testing.
##### Ablation study: importance of IGT correction

We confirm the importance of the implicit gradient transport correction by running an experiment in which an increasing momentum is used without transport. The results appear in figure 7.

The resulting hyper-parameters are:

• ATA stepsize 0.3, decay 0.01, tail fraction 18

• HB-ATA stepsize 0.03, decay 0.01, tail fraction 18 Figure 7: Convergence and accuracy curves along training for the CIFAR10 experiments comparing the use of ATA combined with our proposed implicit transport mechanism. Left: Training. Right: Testing.
##### Effect of the batch size

We look into the effect of the batch size. To do so, we plot the number of steps required to reach a reliably attainable training loss of 0.4 as a function of the batch size. We ran using the following mini-batch sizes: 32, 128, 512 and 2048. Running with larger minibatches led to out of memory errors on our single GPU setup. The results presented in figure 8 show the benefit of IGT lowers as the batch size increases. Note that Adam’s ability to keep benefiting from larger batch sizes is consistent with previous observations. Figure 8: Effect of mini-batch size on the number of steps to reach a target training loss.

### A.2 ImageNet image classification

##### Dataset

We use the 2015 edition of the ImageNet Large-Scale Visual Recognition Challenge (ImageNet)  dataset. This dataset consists of 1.2M images partitioned into 1’000 classes. We use the pre-processing and loading utilities of the TensorFlow models package, available at the following URL: https://github.com/tensorflow/models

##### Model

Our model is again a large residual network, consisting of 50 layers. Similar to our CIFAR10 experiments above, we use the implementation described in .

##### Hyper-parameters

We used the same setup and approach as for the CIFAR-10 experiments. The setup trains for 90 epochs using mini-batches of 1024 examples. We used a larger grid for the learning rate schedule: decay (0.1, 0.01, 0.001) and decay step fraction (0.7, 0.8, 0.9).

The resulting hyper-parameters are:

• SGD stepsize 0.3, decay 0.01, decay step 0.8

• HB stepsize 0.03, decay 0.001, decay step 0.9

• Adam stepsize 0.0001, decay 0.01, decay step 0.9

• ITA stepsize 0.3, decay 0.01, tail fraction 90, decay step 0.9

• HB-ITA stepsize 0.03, decay 0.01, tail fraction 90, decay step 0.9

##### Infrastructure and Runs

We ran these experiments using Google TPUv2.

We include error and accuracy curves for training and testing in Figure 9. Figure 9: Convergence and accuracy curves along training for the ImageNet experiments. Left: Training. Right: Testing.

### A.3 IMDb sentiment analysis

##### Dataset

The Internet Movie Database (IMDb)  consists of 25’000 training and 25’000 test movie reviews. The objective is binary sentiment classification based on the review’s text. We randomly split the training set in two folds of 17’536 and 7’552 reviews, the former being used for training and the latter for testing. The data is downloaded, splitted, and pre-processed with torchtext package, available at the following URL: https://github.com/pytorch/text More specifically, we tokenize the text at the word-level using the spaCy package, and embed the tokens using the 100-dimensional GloVe 6B  distributed representations.

##### Model

The model consists of an embedding lookup-table, followed by a bi-directional LSTM with dropout, and then by a fully-connected layer. The LSTM uses 256 hidden units and the dropout rate is set to 0.5. The whole model consists of 3.2M trainable parameters, with the embedding lookup-table initilized with the GloVe vectors. The model is trained to minimize the BCEWithLogitsLoss with a mini-batch size of 64.

##### Hyper-parameters

For each method, we used a grid-search to find the stepsize minimizing validation error after 15 epochs. The grid starts at 0.00025 and doubles until reaching 0.1024, so as to ensure that no chosen value lies on its boundaries. When applicable, the momentum factor is jointly optimized over values 0.1 to 0.95. The final hyper-parameters are displayed in the following table for each method.

##### Infrastructure and Runs

All IMDb experiments use a single NVIDIA GTX 1080, with PyTorch v0.3.1.post2, CUDA 8.0, and cuDNN v7.0.5. We run each final configurations with 5 different random seeds and always report the mean tendency $\pm$ one standard deviation. Each run lasts approximately three hours and thirty minutes.

In addition to the results reported in the main text, we include training, validation, and testing curves for each method in Figure 10. Shading indicates the one standard deviation interval. Note that our focus is explicitly on optimization: in the specific case of IMDb, training for 200 epochs is completely unnecessary from a generalization standpoint as performance degrades rapidly after 15-20 epochs. Figure 10: Convergence and accuracy curves along training for the IMDb experiments. Left: Convergence. Right: Accuracy.

### A.4 Linear-quadratic regulator

##### Setup

Our linear-quadratic regulator  implements the following equations. The cost functional is evaluated at every timestep $h$ and is given by

 $\displaystyle C(s_{h},a_{h})=s_{h}^{\top}Qs_{h}+a_{h}^{\top}Ra_{h},$ (2)

for random symmetric positive definite matrices $Q\in\mathbb{R}^{20\times 20}$ and $R\in\mathbb{R}^{12\times 12}$ each with condition number 3. The initial state $s_{0}\sim\mathcal{N}(0,3\cdot I_{20})$ is sampled around the origin, and the subsequent states evolve according to

 $\displaystyle s_{h+1}=As_{h}+Ba_{h},$ (3)

where entries of $A\in\mathbb{R}^{20\times 20},B\in\mathbb{R}^{20\times 12}$ are independently sampled from a Normal distribution and then scaled such that the matrix has unit Frobenius norm. The actions are given by the linear stochastic policy $a_{h}=Ks_{h}+\epsilon_{h}^{a}$, where $\epsilon_{h}^{a}\sim\mathcal{N}(0,I)$ and $K$ are the parameters to be optimized.

Gradient methods in this manuscript optimize the sum of costs using the REINFORCE estimate  given by

 $\displaystyle\nabla_{K}\mathbb{E}\sum_{h}^{10}C(s_{h},a_{h})=\mathbb{E}\left(% \sum_{h}^{10}\nabla_{K}\log\pi_{K}(a_{h}|s_{h})\right)\left(\sum_{h}^{10}C(s_{% h},a_{h})\right).$ (4)

In our experiments, the above expectation is approximated by the average of 600 trajectory rollouts. Due to the noisy dynamics of the system, it is possible for the gradient norm to explode leading to numerical unstabilities – especially when using larger stepsizes. To remedy this issue, we simply discard such problematic trajectories from the gradient estimator.

For each training iteration, we first gather 600 trajectories used for learning and then 600 more used to report evaluation metrics.

##### Hyper-parameters

Due to the simplicity of the considered methods, the only hyper-parameter is the stepsize. For each method, we choose the stepsize from a logarithmically-spaced grid so as to minimize the evaluation cost after 600 iterations on a single seed. Incidentally, the optimal stepsize for GD, SGD, and ITA is 0.0002.

##### Infrastructure and Runs

We use an Intel Core i7-5820K CPU to run the LQR experiments. All methods are implemented using numpy v1.15.4. We present results averaged over 3 random seeds, and also report the standard deviation. For stochastic gradient methods (SGD, ITA) training for 20K iterations takes about 3h, for full-gradient method (GD) about 10h, and computing the solution of the Riccati equation takes less than 5 seconds.

In addition to the evaluation cost reported in the main text, we also include the cost witnessed during training (and used for optimization) in Figure 11. Figure 11: LQR costs along training iterations. Left: Costs used for learning. Right: Costs used for evaluation.

We notice that the training cost curve of ITA is not as smooth as the evaluation one. Similarly, the observed learning costs never reach as good a minima as the evaluation ones. This phenomena is easily clarified: during learning, ITA esimates the gradient using the shifted parameters $K_{t}+\frac{\gamma_{t}}{1-\gamma_{t}}(K_{t}-K_{t-1})$ as opposed to the true parameters $K_{t}$. Those shifted parameters are not subject to a reduced variance, hence explaining the observed noise in the cost as well as deteriorated convergence.

### A.5 Model-agnostic meta-learning

##### Dataset

We use the Mini-Imagenet dataset  in our model-agnostic meta-learning (MAML)  experiments. This dataset comprises 64 training, 12 validation, and 24 test classes. For each of train, validation, and test sets, tasks are constructed by sampling 5 classes from their respective split, and further sampling 5 images per class. Images are downsampled to 84x84x3 tensors of RGB values. For more details, please refer to the official code repository – which we carefully replicated – at the following URL: https://github.com/cbfinn/maml

Our implementation departs in two ways from the reference. First, we train our models for 100k iterations as opposed to 60k and only use 5 image samples to compute a meta-gradient whereas the reference implementation uses 15. Second, we only use 5 adaptation steps at evaluation time, when the reference uses 10.

##### Model

The model closely replicates the convolutional neural network of MAML . It consists of 4 layers, each with 32 3x3 kernels, followed by batch normalization and ReLU activations. For specific implementation details, we refer the reader to the above reference implementation.

##### Hyper-parameters

We only tune the meta-stepsize for the MAML experiment. We set the momentum constant to 0.9, the adaptation-stepsize to 0.01, and average the meta-gradient of 4 tasks per iterations. Due to the reduced variance in the gradients, we found it necessary to increase the $\epsilon$ of Adam-ITA to 0.01.

For each method, we search over stepsize values on a logarithmically-spaced grid and select those values that maximize validation accuracy after 10k meta-iterations. These values are reported in Table 2.

##### Infrastructure and Runs

Each MAML experiment is run on a single NVIDIA GTX TITAN X, with PyTorch v1.1.0, CUDA 8.0, and cuDNN v7.0.5. We run each configuration with 3 different random seeds and report the mean tendency $\pm$ one standard deviation. Each run takes approximately 36 hours, and we evaluate the validation and testing accuracy every 100 iteration.

We complete the MAML validation curves from the main manuscript with training and testing accuracy curves in Figure 12. Moreover, we recall the final test accuracies for each method: Adam-ITA reaches $65.16\%$, HB-ITA $64.57\%$, Adam $63.70\%$, and HB $63.08\%$. Figure 12: Training, validation, and testing accuracies for the MAML experiments along training. Shading indicates the 1 standard deviation interval. Left: Training. Center: Validation. Right: Testing.

## Appendix B Additional Experiments

This section presents additional experiments to the ones included in the main text.

### B.1 Baselines comparisons

While experiments in Section 5 highlighted properties of IGT and HB-IGT when the assumption of identical, constant Hessians was verified, we now turn to more realistic scenarios where individual functions are neither quadratic nor have the same Hessian to compare our proposed methods against popular baselines for the online stochastic optimization setting. We target optimization benchmarks and focus on training loss minimization. Figure 13: Training loss curves for different optimization algorithms on several popular benchmarks. For each method, the hyper-parameters are tuned to minimize the training error after 15 epochs. Algorithms using the IGT gradient estimates tend to outperform their stochastic gradient counter-parts. Left: Logistic regression on MNIST. Center: LeNet5 on MNIST. Right: MobileNetv2 on CIFAR10.

We investigate three different scenarios: (a) linear-mnist: a logistic regression model on MNIST, (b) mnist: a modified version of LeNet5  on MNIST and (c) cifar-small: the MobileNetv2 architecture  consisting of 19 convolutional layers on CIFAR10. All models are trained with a mini-batch size of 64, while the remaining hyper-parameters are available in Tables 45, and 6.

For each of the above tasks, models are trained for 200 epochs. We compare the following methods:

• HB: the heavy ball method ,

• ASGD ,

• SVRG ,

• SGD-dec: stochastic gradient method with an exponential learning rate schedule and exponential constant $0.999$,

• HB-IGT: the heavy ball using the IGT as a plug-in estimator, and

• HB-ITA: same as HB-IGT but using the anytime tail averaging to forget the oldest gradients.

The hyperparameters of each method, and in particular the stepsize, are tuned independently according to a logarithmic grid so as to minimize the mean training error after epoch 15 on one seed. We then use those parameters on 5 random seeds and report the mean and standard deviation of the performance.

Figure 13 shows the training curves for the five algorithms in the three settings.

First, we observe that, for the logistic regression, HB-IGT performs on par with HB-ITA and far better than all the other methods, even though the assumption on the Hessians is violated. When using a ConvNet, however, we see that HB-IGT is not competitive with state-of-the-art methods such as Adam or ASGD. HB-ITA, on the other hand, with its smaller reliance on the assumption, once again performs much better than all other methods. In fact, HB-ITA not only converges to a lower final train error but also has a faster initial rate.

While our focus is on optimization, we also report generalization metrics in Table 3. For each algorithm, we computed the best mean accuracy after each epoch on the test set and report this value together with its standard deviation. The importance of the Anytime Tail-Averaging mechanism is again apparent: without it, Heavyball-IGT is unable to improve for more than a few epochs on the CIFAR10 validation set, regardless of the stepsize choice. On the other hand, it is evident from those results that the solutions found by Heavyball-ITA are competitive with the ones discovered by other optimization algorithms.

## Appendix C Proofs

### C.1 Transport formula

 $\displaystyle g_{t}(\theta_{t})$ $\displaystyle=\frac{t}{t+1}g_{t-1}(\theta_{t})+\frac{1}{t+1}g(\theta_{t},x_{t})$ $\displaystyle=\frac{t}{t+1}\left(g_{t-1}(\theta_{t-1})+H(\theta_{t}-\theta_{t-% 1})\right)+\frac{1}{t+1}g(\theta_{t},x_{t})$ (Quadratic $f$) $\displaystyle=\frac{t}{t+1}g_{t-1}(\theta_{t-1})+\frac{1}{t+1}\left(g(\theta_{% t},x_{t})+tH(\theta_{t}-\theta_{t-1})\right)$ $\displaystyle=\frac{t}{t+1}g_{t-1}(\theta_{t-1})+\frac{1}{t+1}g(\theta_{t}+t(% \theta_{t}-\theta_{t-1}),x_{t})$ (Identical Hessians) $\displaystyle\approx\frac{t}{t+1}\widehat{g}_{t-1}(\theta_{t-1})+\frac{1}{t+1}% g(\theta_{t}+t(\theta_{t}-\theta_{t-1}),x_{t})\;.$ ($\widehat{g}_{t-1}$ is an approximation)

## Appendix D Proof of Prop. 3.1

In this proof, we assume that $g$ is a strongly-convex quadratic function with hessian $H$.

At timestep $t$, we have access to a stochastic gradient $g(\theta,x_{t})=g(\theta_{t})+\epsilon_{t}$ where the $\epsilon_{t}$ are i.i.d. with variance $C\preceq\sigma^{2}H$.

We first prove a simple lemma:

###### Lemma D.1.

If $v_{0}=g(\theta_{0})+\epsilon_{0}$ and, for $t>0$, we have

 $\displaystyle v_{t}$ $\displaystyle=\frac{t}{t+1}v_{t-1}+\frac{1}{t+1}g\left(\theta_{t}+t(\theta_{t}% -\theta_{t-1})\right)+\frac{1}{t+1}\epsilon_{t}\;,$

then

 $\displaystyle v_{t}$ $\displaystyle=g(\theta_{t})+\frac{1}{t+1}\sum_{i=0}^{t}\epsilon_{i}\;.$
###### Proof.

Per our assumption, this is true for $t=0$. Now let us prove the result by induction. Assume this is true for $t-1$. Then we have:

 $\displaystyle v_{t}$ $\displaystyle=\frac{t}{t+1}v_{t-1}+\frac{1}{t+1}g(\theta_{t}+t(\theta_{t}-% \theta_{t-1}))+\frac{1}{t+1}\epsilon_{t}$ $\displaystyle=\frac{t}{t+1}g(\theta_{t-1})+\frac{1}{t+1}\sum_{i=0}^{t-1}% \epsilon_{i}$ $\displaystyle\qquad+\frac{1}{t+1}g(\theta_{t}+t(\theta_{t}-\theta_{t-1}))+% \frac{1}{t+1}\epsilon_{t}$ (recurrence assumption) $\displaystyle=\frac{t}{t+1}g(\theta_{t-1})+\frac{1}{t+1}\sum_{i=0}^{t-1}% \epsilon_{i}$ $\displaystyle\qquad+g(\theta_{t})-\frac{t}{t+1}g(\theta_{t-1})+\frac{1}{t+1}% \epsilon_{t}$ ($g$ is quadratic) $\displaystyle=g(\theta_{t})+\frac{1}{t+1}\sum_{i=0}^{t}\epsilon_{i}\;.$

This concludes the proof. ∎

###### Lemma D.2.

Let us assume we perform the following iterative updates:

 $\displaystyle v_{t}$ $\displaystyle=\frac{t}{t+1}v_{t-1}+\frac{1}{t+1}g(\theta_{t}+t(\theta_{t}-% \theta_{t-1}))+\frac{1}{t+1}\epsilon_{t}$ $\displaystyle\theta_{t+1}$ $\displaystyle=\theta_{t}-\alpha v_{t}\;,$

starting from $v_{0}=g(\theta_{0})+\epsilon_{0}$. Then, denoting $\Delta_{t}=\theta_{t}-\theta^{\ast}$, we have

 $\displaystyle\Delta_{t}$ $\displaystyle=(I-\alpha H)^{t}\Delta_{0}-\alpha\sum_{i=0}^{t-1}N_{i,t}\epsilon% _{i}$

with

 $\displaystyle N_{i,0}$ $\displaystyle=0$ $\displaystyle N_{i,t}$ $\displaystyle=(I-\alpha H)N_{i,t-1}+1_{i
###### Proof.

The result is true for $t=0$. We now prove the result for all $t$ by induction. Let us assume this is true for $t-1$. Using Lemma D.1, we have

 $\displaystyle v_{t-1}$ $\displaystyle=g(\theta_{t-1})+\frac{1}{t}\sum_{i=0}^{t-1}\epsilon_{i}$

and thus, using $g(\theta_{t-1})=H\Delta_{t-1}$,

 $\displaystyle\Delta_{t}$ $\displaystyle=\Delta_{t-1}-\alpha v_{t-1}$ $\displaystyle=\Delta_{t-1}-\alpha H\Delta_{t-1}-\frac{\alpha}{t}\sum_{i=0}^{t-% 1}\epsilon_{i}$ $\displaystyle=(I-\alpha H)\Delta_{t-1}-\frac{\alpha}{t}\sum_{i=0}^{t-1}% \epsilon_{i}$ $\displaystyle=(I-\alpha H)^{t}\Delta_{0}-\alpha\sum_{i=0}^{t-2}(I-\alpha H)N_{% i,t-1}\epsilon_{i}-\frac{\alpha}{t}\sum_{i=0}^{t-1}\epsilon_{i}$ (recurrence assumption) $\displaystyle=(I-\alpha H)^{t}\Delta_{0}-\alpha\sum_{i=0}^{t-1}N_{i,t}\epsilon% _{i}$

with

 $\displaystyle N_{i,t}$ $\displaystyle=(I-\alpha H)N_{i,t-1}+1_{i

This concludes the proof. ∎

For the following lemma, we will assume that the Hessian is diagonal and will focus on one dimension with eigenvalue $h$. Indeed, we know that there are no interactions between the eigenspaces and that we can analyze each of them independently .

###### Lemma D.3.

Denote $r_{h}=1-\alpha h$. We assume $\alpha\leq\frac{1}{L}$. Then, for any $i$ and any $t$, we have

 $\displaystyle N_{i,t}$ $\displaystyle\geq 0$ (Positivity) $\displaystyle N_{i,t}$ $\displaystyle=0\quad\textrm{ if }t\leq i$ (Zero-start) $\displaystyle N_{i,t}$ $\displaystyle\leq\log\left(\frac{2}{i(1-r_{h})}\right)\quad\textrm{ if }i (Constant bound) $\displaystyle N_{i,t}$ $\displaystyle\leq\frac{\max\left\{1+r_{h},2\log\left(\frac{2}{i(1-r_{h})}% \right)\right\}}{t(1-r_{h})}\quad\textrm{if }\frac{2}{1-r_{h}}\leq t\;.$ (Decreasing bound)
###### Proof.

The Zero-start case $i\geq t$ is immediate from the recursion of Lemma D.2. The Positivity property of $N_{i,t}$ is also immediate from the recursion since the stepsize $\alpha$ is such that $r_{h}=1-\alpha h$ is positive.

We now turn to the Constant bound property. We have, for $t>i$,

 $\displaystyle N_{i,t}$ $\displaystyle=r_{h}N_{i,t-1}+\frac{1}{t}$ $\displaystyle\leq N_{i,t-1}+\frac{1}{t}\;.$

Thus, $N_{i,t}-N_{i,t-1}\leq\frac{1}{t}$. Summing these inequalities, we get a telescopic sum and, finally:

 $\displaystyle N_{i,t}$ $\displaystyle\leq\sum_{j=i+1}^{t}\frac{1}{j}$ $\displaystyle\leq\int_{x=i}^{t}\frac{dx}{x}$ $\displaystyle=\log\left(\frac{t}{i}\right)\;.$

This bound is trivial in the case $i=0$. In that case, we keep the first term in the sum separate and get

 $\displaystyle N_{0,t}$ $\displaystyle\leq 1+\log t\;.$

In the remainder, we shall keep the $\log\left(\frac{t}{i}\right)$ bound for simplicity. The upper bound on the right-hand size is increasing with $t$ and its value for $t=\frac{2}{1-r_{h}}$ is thus an upper bound for all smaller values of $t$. Replacing $t$ with $\frac{2}{1-r_{h}}$ leads to

 $\displaystyle N_{i,\frac{2}{1-r_{h}}}$ $\displaystyle\leq\log\left(\frac{\frac{2}{1-r_{h}}}{i}\right)$ $\displaystyle=\log\left(\frac{2}{i(1-r_{h})}\right)\;.$

This proves the third inequality.

We shall now prove the Decreasing bound by induction. This bound states that, for $t$ large enough, each $N_{i,t}$ decreases as $O(1/t)$. Using the second and third inequalities, we have

 $\displaystyle N_{i,\frac{2}{1-r_{h}}}$ $\displaystyle\leq\log\left(\frac{2}{i(1-r_{h})}\right)\frac{\frac{2}{1-r_{h}}}% {\frac{2}{1-r_{h}}}$ $\displaystyle=\frac{\log\left(\frac{2}{i(1-r_{h})}\right)}{1-r_{h}}\frac{2}{% \frac{2}{1-r_{h}}}$ $\displaystyle\leq\frac{\max\left\{1+r_{h},2\log\left(\frac{2}{i(1-r_{h})}% \right)\right\}}{\frac{2}{1-r_{h}}(1-r_{h})}\ .$

The maximum will help us prove the last property. Thus, for $t=\frac{2}{1-r_{h}}$, we have

 $\displaystyle N_{i,t}$ $\displaystyle\leq\frac{\max\left\{1+r_{h},2\log\left(\frac{1}{i(1-r_{h})}% \right)\right\}}{t(1-r_{h})}$ $\displaystyle\leq\frac{\nu_{i}}{t}\;,$

with $\nu_{i}=\frac{\max\left\{1+r_{h},2\log\left(\frac{1}{i(1-r_{h})}\right)\right% \}}{(1-r_{h})}$. The Decreasing bound is verified for $t=\frac{2}{1-r_{h}}$.

We now show that if, for any $t>\frac{2}{1-r_{h}}$, we have $N_{i,t-1}\leq\frac{\nu_{i}}{t-1}$, then $N_{i,t}\leq\frac{\nu_{i}}{t}$. Assume that there is such at $t$. Then

 $\displaystyle N_{i,t}$ $\displaystyle=r_{h}N_{i,t-1}+\frac{1}{t}$ $\displaystyle\leq\frac{r_{h}\nu_{i}}{t-1}+\frac{1}{t}$ $\displaystyle=\frac{r_{h}t\nu_{i}+t-1}{t(t-1)}$ $\displaystyle=\frac{(t-1)\nu_{i}+(r_{h}-1)t\nu_{i}+\nu_{i}+t-1}{t(t-1)}$ $\displaystyle=\frac{\nu_{i}}{t}+\frac{(r_{h}-1)t\nu_{i}+\nu_{i}+t-1}{t(t-1)}\;.$

We now shall prove that $(r_{h}-1)t\nu_{i}+\nu_{i}+t-1=[(r_{h}-1)\nu_{i}+1]t+\nu_{i}-1$ is negative. First, we have that

 $\displaystyle(r_{h}-1)\nu_{i}+1$ $\displaystyle=1-\max\left\{1+r_{h},2\log\left(\frac{1}{i(1-r_{h})}\right)\right\}$ $\displaystyle\leq 0\;.$

Then,

 $\displaystyle[(r_{h}-1)\nu_{i}+1]t+\nu_{i}-1\leq 0$ $\displaystyle\Longleftrightarrow t\geq\frac{\nu_{i}-1}{(1-r_{h})\nu_{i}-1}$

since $(r_{h}-1)\nu_{i}+1\leq 0$. Thus, the property is true for every $t\geq\frac{\nu_{i}-1}{(1-r_{h})\nu_{i}-1}$. In addition, we have

 $\displaystyle\nu_{i}$ $\displaystyle\geq\frac{1+r_{h}}{1-r_{h}}$ $\displaystyle\nu_{i}(1-r_{h})$ $\displaystyle\geq 1+r_{h}$ $\displaystyle 2\nu_{i}(1-r_{h})-2$ $\displaystyle\geq\nu_{i}(1-r_{h})-1+r_{h}$ $\displaystyle\frac{2}{1-r_{h}}$ $\displaystyle\geq\frac{\nu_{i}-1}{\nu_{i}(1-r_{h})-1}\;,$

and the property is also true for every $t\geq\frac{2}{1-r_{h}}$. This concludes the proof. ∎

Finally, we can prove the Proposition 3.1:

###### Proof.

The expectation of $\Delta_{t}$ is immediate using Lemma D.2 and the fact that the $\epsilon_{i}$ are independent, zero-mean noises. The variance is equal to $V[\Delta_{t}]=\alpha^{2}B\sum_{i=0}^{t}N_{i,t}^{2}$. While our analysis was only along one eigenspace of the Hessian with associated eigenvalue $h$, we must now sum over all dimensions. We will thus define

 $\displaystyle\bar{\nu}_{i}$ $\displaystyle=\frac{\max\left\{2-\alpha\mu,2\log\left(\frac{1}{i\alpha\mu}% \right)\right\}}{\alpha\mu}\quad\textrm{ for }i>0$ $\displaystyle\bar{\nu}_{0}$ $\displaystyle=\frac{2+2\log\left(\frac{1}{\alpha\mu}\right)}{\alpha\mu}\;,$

which is, for every $i$, the maximum $\nu_{i}$ across all dimensions. We get

 $\displaystyle V[\Delta_{t}]$ $\displaystyle\leq d\alpha^{2}B\sum_{i=0}^{t}\frac{\bar{\nu}_{i}^{2}}{t^{2}}$ $\displaystyle\leq d\alpha^{2}B\sum_{i=0}^{t}\frac{\bar{\nu}_{0}^{2}}{t^{2}}% \quad\textrm{since }\nu_{i}\geq\nu_{i+1}\;\forall i$ $\displaystyle\leq\frac{d\alpha^{2}B\bar{\nu}_{0}^{2}}{t}\;.$

Since we have

 $\displaystyle E[\theta_{t}-\theta^{\ast}]$ $\displaystyle=(I-\alpha H)^{t}(\theta_{0}-\theta^{\ast})\;,$

we get

 $\displaystyle E[\|\theta_{t}-\theta^{\ast}\|^{2}]$ $\displaystyle=\|E[\theta_{t}-\theta^{\ast}]\|^{2}+V[\Delta_{t}]$ $\displaystyle\leq(\theta_{0}-\theta^{\ast})^{\top}(I-\alpha H)^{2t}(\theta_{0}% -\theta^{\ast})+\frac{d\alpha^{2}B\bar{\nu}_{0}^{2}}{t}$ $\displaystyle\leq\left(1-\frac{1}{\kappa}\right)^{2t}\|\theta_{0}-\theta^{\ast% }\|^{2}+\frac{d\alpha^{2}B\bar{\nu}_{0}^{2}}{t}\;.$

This concludes the proof. ∎

## Appendix E Proof of Proposition 3.2 and Proposition 3.3

In this section we list and prove all lemmas used in the proofs of Proposition 3.2 and Proposition 3.3; all lemmas are stated in the same conditions as the proposition.

We start the following proposition:

###### Proposition E.1.

Let $f$ be a quadratic function with positive definite Hessian $H$ with largest eigenvalue $L$ and condition number $\kappa$ and if the stochastic gradients satisfy $g(\theta,x)=g(\theta)+\epsilon$ with $\epsilon$ a random uncorrelated noise with covariance bounded by $BI$.

Then, Algorithm 3.3 leads to iterates $\theta_{t}$ satisfying

 $\displaystyle E[\theta_{t}-\theta^{\ast}]$ $\displaystyle=\begin{pmatrix}I\\ 0\end{pmatrix}A^{t}\begin{pmatrix}E[\theta_{1}-\theta^{\ast}]\\ E[\theta_{0}-\theta^{\ast}]\end{pmatrix}$ (5)

where

 $\displaystyle A$ $\displaystyle=\begin{pmatrix}I-\alpha H+\mu I&-\mu I\\ I&0\end{pmatrix}$ (6)

governs the dynamics of this bias. In particular, when its spectral radius, $\rho(A)$ is less than $1$, the iterates converge linearly to $\theta^{\ast}$.

In a similar fashion, the variance dynamics of Heavyball-IGT are governed by the matrix

 $\displaystyle D_{i}$ $\displaystyle=\begin{pmatrix}(1-\alpha h_{i}+\mu)^{2}+2\alpha^{2}h_{i}^{2}&\mu% ^{2}&-2\mu(1-\alpha h_{i}+\mu)^{2}\\ 1&0&0\\ 1-\alpha h_{i}+\mu&0&-\mu\end{pmatrix}$

If the spectral radius of $D_{i}$, $\rho(D_{i})$, is strictly less than 1 or all $i$, then there exist constants $t_{0}>0$ and $C>0$ for which

 $\displaystyle\mathrm{Var}{(\theta_{t})}$ $\displaystyle\leq 2\alpha^{2}dBC\frac{\log(t)}{t},\quad\textrm{for}\ t>t_{0}$

where $B$ is a bound on the variance of noise variables $\epsilon_{i}$.

###### Lemma E.2 (IGT estimator as true gradient plus noise average).

If $v_{0}=g(\theta_{0})+\epsilon_{0}$ and for $t>0$ we have

 $v_{t}=\frac{t}{t+1}v_{t-1}+\frac{1}{t+1}g(\theta_{t}+t(\theta_{t}-\theta_{t-1}% ))+\frac{1}{t+1}\epsilon_{t},$

then

 $v_{t}=g(\theta_{t})+\frac{1}{t+1}\sum_{i=0}^{t}\epsilon_{i}.$

This lemma is already proved in the previous section for the IGT estimator (Lemma D.1) and is just repeated here for completeness. We will use this result in the next few lemmas.

###### Lemma E.3 (The IGT gradient estimator is unbiased on quadratics).

For the IGT gradient estimator, $v_{t}$, corresponding to parameters $\theta_{t}$ we have

 $\mathbb{E}\left[v_{t}\right]=g(\mathbb{E}\theta_{t}),$

where the expectation is over all gradient noise vectors $\epsilon_{0},\epsilon_{1},\ldots,\epsilon_{t}$.

###### Proof.

The proof proceeds by induction. The base case holds as we have

 $\mathbb{E}\left[v_{0}\right]=\mathbb{E}\left[g_{0}+\epsilon_{0}\right]=g(% \theta_{0}).$

For the inductive case, we can write

 $\displaystyle\mathbb{E}\left[v_{t}\right]$ $\displaystyle=\mathbb{E}\left[\frac{t}{t+1}v_{t-1}+\frac{1}{t+1}\hat{g}(\theta% _{t}+t(\theta_{t}-\theta_{t-1}))\right]$ $\displaystyle=\mathbb{E}\left[\frac{t}{t+1}v_{t-1}+\frac{1}{t+1}g_{t}+\frac{t}% {t+1}g_{t}-\frac{t}{t+1}g_{t-1}+\frac{1}{t+1}\epsilon_{t}\right]$ $\displaystyle=\frac{t}{t+1}\mathbb{E}\left[v_{t-1}-g_{t-1}\right]+\mathbb{E}% \left[g_{t}\right]+\frac{t}{t+1}\mathbb{E}\left[\epsilon_{t}\right]$ $\displaystyle=\mathbb{E}\left[g_{t}\right]=g(\mathbb{E}\left[\theta_{t}\right]).$

Where, in the third equality, $\mathbb{E}\left[v_{t-1}-g_{t-1}\right]=0$ by the inductive assumption, and the last equality because the gradient of a quadratic function is linear. ∎

###### Lemma E.4 (Bounding the IGT gradient variance).

Let $v_{t}$ be the IGT gradient estimator. Then

 $\mathrm{Var}\left[v_{t}\right]\leq 2h^{2}\mathrm{Var}\left[\theta_{t}-\theta^{% \star}\right]+\frac{2B}{t},$

where $B$ is the variance of the homoscedastic noise $\epsilon_{t}$.

###### Proof.
 $\displaystyle\mathrm{Var}\left[v_{t}\right]$ $\displaystyle=\mathrm{Var}\left[g_{t}+\frac{1}{t+1}\sum_{i=0}^{t}\epsilon_{i}\right]$ $\displaystyle=\mathrm{Var}\left[h\theta_{t}\right]+\mathrm{Var}\left[\frac{1}{% t+1}\sum_{i=0}^{t}\epsilon_{i}\right]$ $\displaystyle\qquad+2\mathrm{Cov}\left[h\theta_{t},\frac{1}{t+1}\sum_{i=0}^{t}% \epsilon_{i}\right]$ $\displaystyle\leq 2\mathrm{Var}\left[h\theta_{t}\right]+2\mathrm{Var}\left[% \frac{1}{t+1}\sum_{i=0}^{t}\epsilon_{i}\right]$ $\displaystyle=2h^{2}\mathrm{Var}\left[\theta_{t}-\theta^{\star}\right]+2\frac{% B}{t}$

Now that we have these basic results on the IGT estimator, we can analyze the evolution of the bias and variance of Heavyball-IGT. We use the quadratic assumption to decouple the vector dynamics of Heavyball-IGT into independent scalar dynamics. If the Hessian, $H$, has eigenvalues $L\geq h_{1}\geq h_{2}\geq\ldots\geq h_{n}=L/\kappa$, then we can assume without loss of generality that $H$ is diagonal with $H_{ii}=h_{i}$.

###### Lemma E.5 (Evolution of bias for scalar quadratic).

Assume that the Hessian, second derivative, is h.

Starting with $v_{0}=g(\theta_{0})+\epsilon_{0}$ and $w_{0}=0$, performing the following iterative updates (Heavyball-IGT, Algorithm 3.3):

 $v_{t}=\frac{t}{t+1}v_{t-1}+\frac{1}{t+1}g(\theta+t(\theta_{t}-\theta_{t-1}))+% \frac{1}{t+1}\epsilon_{t},$
 $w_{t+1}=\mu w_{t}+\alpha v_{t},\qquad\theta_{t+1}=\theta_{t}-w_{t+1}$

results in

 $\Delta_{t}=A^{t}\Delta_{0}-\alpha\sum_{i=0}^{t-1}N_{i,t}\begin{bmatrix}% \epsilon_{i}\\ 0\end{bmatrix}$

where $N_{j,0}=0_{2\times 2},\qquad N_{i,t}=AN_{i,t-1}+1_{i,

$\Delta_{t}=\begin{bmatrix}\theta_{t}-\theta^{\ast}\\ \theta_{t-1}-\theta^{\ast}\end{bmatrix}$ and $A=\begin{pmatrix}1-\alpha h+\mu&-\mu\\ 1&0\end{pmatrix}.$

###### Proof.

The proof proceeds by induction. First notice that for $t=0$ the equality naturally holds. We make the inductive assumption that it holds for $t-1$, and start by using Lemma E.2:

 $\displaystyle\Delta_{t}$ $\displaystyle=A\Delta_{t-1}-\frac{\alpha}{t}\sum_{i=0}^{t-1}\begin{bmatrix}% \epsilon_{i}\\ 0\end{bmatrix}$ $\displaystyle=A(A^{t-1}\Delta_{0}-\alpha\sum_{i=0}^{t-2}N_{i,t}\begin{bmatrix}% \epsilon_{i}\\ 0\end{bmatrix})-\frac{\alpha}{t}\sum_{i=0}^{t-1}\begin{bmatrix}\epsilon_{i}\\ 0\end{bmatrix}$ (Inductive assumption) $\displaystyle=A^{t}\Delta_{0}-\alpha(\sum_{i=0}^{t-2}AN_{i,t}\begin{bmatrix}% \epsilon_{i}\\ 0\end{bmatrix}+\frac{1}{t}\sum_{i=0}^{t-1}\begin{bmatrix}\epsilon_{i}\\ 0\end{bmatrix})$ $\displaystyle=A^{t}\Delta_{0}-\alpha\sum_{i=0}^{t-1}N_{i,t}\begin{bmatrix}% \epsilon_{i}\\ 0\end{bmatrix}$ (Def. of $N_{i,t}$)

###### Lemma E.6 (Evolution of variance).

Let $U_{t}=\mathrm{Var}\left[\theta_{t}\right]$ and $V_{t}=\mathrm{Cov}\left[\theta_{t},\theta_{t-1}\right]$, where $\theta_{t}$ is the $t$-th iterate of Heavyball-IGT on a 1-dimensional quadratic function with curvature $h$. The following matrix describes the variance dynamics of Heavyball-IGT.

 $D=\begin{pmatrix}(1-\alpha h+\mu)^{2}+2\alpha^{2}h^{2}&\mu^{2}&-2\mu(1-\alpha h% +\mu)^{2}\\ 1&0&0\\ 1-\alpha h+\mu&0&-\mu\end{pmatrix}$ (7)

If the spectral radius of $D$, $\rho(D)$, is strictly less than 1, then there exist constants $t_{0}>0$ and $C>0$ for which

 $\displaystyle\mathrm{Var}{(\theta_{t})}$ $\displaystyle\leq 2\alpha^{2}BC\frac{\log(t)}{t}$

, where $B$ is a bound on the variance of the noise.

###### Proof.

The proof (and lemma) is similar to the proof of Lemma 9 in . We start by expanding $U_{t+1}$ as follows.

 $\displaystyle U_{t+1}$ $\displaystyle=\mathbb{E}\left[(\theta_{t+1}-\bar{\theta}_{t+1})^{2}\right]$ $\displaystyle=\mathbb{E}\left[(\theta_{t}-\alpha v_{t}+\mu(\theta_{t}-\theta_{% t-1})-\bar{\theta}_{t}+\alpha g_{t}-\mu(\bar{\theta}_{t}-\bar{\theta}_{t-1}))^% {2}\right]$ $\displaystyle=\mathbb{E}[(\theta_{t}-\alpha g_{t}+\mu(\theta_{t}-\theta_{t-1})% -\bar{\theta}_{t}+\alpha g_{t}$ $\displaystyle\qquad-\mu(\bar{\theta}_{t}-\bar{\theta}_{t-1})+\alpha(g_{t}-v_{t% }))^{2}]$ $\displaystyle=\mathbb{E}\left[((1-\alpha h+\mu)(\theta_{t}-\bar{\theta}_{t})-% \mu(\theta_{t-1}-\bar{\theta}_{t-1}))^{2}\right]$ $\displaystyle\qquad+\alpha^{2}\mathbb{E}\left[(g_{t}-v_{t})^{2}\right]$ $\displaystyle\leq\mathbb{E}\left[((1-\alpha h+\mu)(\theta_{t}-\bar{\theta}_{t}% )-\mu(\theta_{t-1}-\bar{\theta}_{t-1}))^{2}\right]$ $\displaystyle\qquad+\alpha^{2}\left(2h^{2}\mathbb{E}\left[(\theta_{t}-\bar{% \theta}_{t})^{2}\right]+\frac{2B}{t+1}\right)$ $\displaystyle\leq\left[(1-\alpha h+\mu)^{2}+2\alpha^{2}\mu^{2}\right]\mathbb{E% }\left[(\theta_{t}-\bar{\theta}_{t})^{2}\right]$ $\displaystyle\qquad-2\mu(1-\alpha h+\mu)\mathbb{E}\left[(\theta_{t}-\bar{% \theta}_{t})(\theta_{t-1}-\bar{\theta}_{t-1})\right]$ $\displaystyle\qquad+\mu^{2}\mathbb{E}\left[(\theta_{t-1}-\bar{\theta}_{t-1})^{% 2}\right]+\alpha^{2}\frac{2B}{t+1}\;.$

Where the fourth equality is obtained since we know that the IGT gradient estimator is unbiased, i.e. $\mathbb{E}\left[g_{t}-v_{t}\right]=0$. The first inequality stems from Lemma E.4. We similarly expand $V_{t}$.

 $\displaystyle V_{t}$ $\displaystyle=\mathbb{E}\left[(\theta_{t}-\bar{\theta}_{t})(\theta_{t-1}-\bar{% \theta}_{t-1})\right]$ $\displaystyle=\mathbb{E}\left[\left((1-\alpha h+\mu)(\theta_{t-1}-\bar{\theta}% _{t-1})-\mu(\theta_{t-2}-\bar{\theta_{t-2}})+\alpha(g_{t}-v_{t})\right)\right.$ $\displaystyle\qquad\left.(\theta_{t-1}-\bar{\theta}_{t-1})\right]$ $\displaystyle=(1-\alpha h+\mu)\mathbb{E}\left[(\theta_{t-1}-\bar{\theta}_{t-1}% )^{2}\right]$ $\displaystyle\qquad-\mu\mathbb{E}\left[(\theta_{t-1}-\bar{\theta}_{t-1})(% \theta_{t-2}-\bar{\theta}_{t-2})\right]$

From the above expressions, we obtain

 $\displaystyle\begin{pmatrix}U_{t+1}\\ U_{t}\\ V_{t+1}\end{pmatrix}$ $\displaystyle\leq D\begin{pmatrix}U_{t}\\ U_{t-1}\\ V_{t}\end{pmatrix}+\begin{pmatrix}\alpha^{2}\frac{2B}{t+1}\\ 0\\ 0\end{pmatrix}$ $\displaystyle\leq 2\alpha^{2}B\sum_{i=0}^{t}D^{i}\begin{pmatrix}\frac{1}{t+1-i% }\\ 0\\ 0\end{pmatrix}$ $\displaystyle\leq 2\alpha^{2}B\left(\sum_{i=0}^{s-1}D^{i}\begin{pmatrix}\frac{% 1}{t+1-i}\\ 0\\ 0\end{pmatrix}+\sum_{i=s}^{t}D^{i}\begin{pmatrix}\frac{1}{t+1-i}\\ 0\\ 0\end{pmatrix}\right)$

where an inequality of vectors implies the corresponding elementwise inequalities.

If the spectral radius of $D$, $\rho(D)$ is strictly less than 1, then there exists constant $C^{\prime}>0$ such that

 $\displaystyle\begin{pmatrix}1\\ 0\\ 0\end{pmatrix}^{T}\sum_{i=0}^{s-1}D^{i}\begin{pmatrix}\frac{1}{t+1-i}\\ 0\\ 0\end{pmatrix}$ $\displaystyle\leq C^{\prime}\sum_{i=0}^{s-1}\frac{1}{t+1-i}$ $\displaystyle\leq\frac{C^{\prime}s}{t+2-s}$

If the spectral radius of $D$, $\rho(D)$, is strictly less than 1, then there exists constant $\zeta>0$ and constant $C^{\prime\prime}(\zeta)>0$ such that, $\rho(D)+\zeta<1$ and

 $\displaystyle\begin{pmatrix}1\\ 0\\ 0\end{pmatrix}^{T}\sum_{i=s}^{t}D^{i}\begin{pmatrix}\frac{1}{t+1-i}\\ 0\\ 0\end{pmatrix}\leq$ $\displaystyle\begin{pmatrix}1\\ 0\\ 0\end{pmatrix}^{T}\sum_{i=s}^{t}D^{i}\begin{pmatrix}1\\ 0\\ 0\end{pmatrix}$ $\displaystyle\leq$ $\displaystyle C^{\prime\prime}\sum_{i=s}^{t}(\rho(D)+\zeta)^{s}$ $\displaystyle=$ $\displaystyle C^{\prime\prime}(t-s+1)(\rho(D)+\zeta)^{s}$

Let $\rho^{\prime}=\rho(D)+\zeta$ and $s=\lceil 2\log_{1/\rho^{\prime}}t\rceil$. Then $(\rho(D)+\zeta)^{s}=1/t^{2}$, and putting the above two bounds together,

 $\displaystyle U_{t+1}$ $\displaystyle\leq 2\alpha^{2}B\left(\frac{2C^{\prime}\log_{1/\rho^{\prime}}t}{% t+2-2\log_{1/\rho^{\prime}}t}+C^{\prime\prime}\frac{t-2\log_{1/\rho^{\prime}}t% +1}{t^{2}}\right)$ $\displaystyle\leq 2\alpha^{2}BC\frac{\log(t+1)}{t+1}$

where the last inequality holds for $t>t_{0}$ for some $t_{0}$ and some constant $C>0$.

We can now prove Proposition E.1.

###### Proof of Proposition E.1.

The bias statement of the proposition follows directly from taking an expectation on the bound of Lemma B.4, and the variance statement from summing up the $d$ different variance terms given for each scalar component by Lemma B.5. ∎

### E.1 Proof of Proposition 3.2

This Proposition follows from the observation that, in the noiseless case, $\epsilon_{t}=0$ in our model. In that case, Lemma E.3 shows that Heavyball-IGT reduces to the heavy ball, and the rest follows from the optimal tuning of the heavy ball .

### E.2 Proof of Proposition 3.3

###### Proof.

Like we did in previous proofs, we can assume without loss of generality that the Hessian, $H$, is diagonal with elements $h_{i}$. For a diagonal $H$, matrix $A$ can be permuted to be block diagonal with blocks

 $A_{i}=\begin{pmatrix}1-\alpha h_{i}+\mu&\mu\\ 1&0\\ \end{pmatrix}.$

To prove that $\rho(A)<1$ it suffices to prove that $\rho(A_{i})<1$ for all $i$. For the rest of the proof we will focus on the dynamics along a single eigendirection with curvature $h_{i}$. The rest of this proof used $D$ to denote $D_{i}$, $A$ to denote $A_{i}$ and $h$ to denote $h_{i}$.

To make explicit the dependence of matrices $A$ and $D$ on hyperparameters and curvature, we write $A(\alpha,\mu,h)$ and $D(\alpha,\mu,h)$. Let $0<\alpha<2/(3h)$ and $\mu_{0}=0$. Using hyperparameters $(\alpha,\mu_{0})$ one gets the results for gradient descent without momentum. In particular $\rho(A(\alpha,\mu_{0},h))=|1-\alpha h|<1$, and the spectral radius of $D$ is $\rho(D(\alpha,\mu_{0},h))=|(1-\alpha h)^{2}+2\alpha^{2}h^{2}|<1$.

We will argue that there exists $\mu>0$, such that $\rho(A(\alpha,\mu,h))<1$, and the spectral radius of $D$ is $\rho(D(\alpha,\mu,h))<1$. Then the previous lemma implies that bias converges linearly, and variance is $O(\log(t)/t)$.

To argue the existence of $\mu>0$, we will perform eigenvalue perturbation analysis using the Bauer-Fike theorem. Note that $A(\alpha,\mu,h)=A(\alpha,\mu_{0},h)+\mu\Delta_{A}$ where

 $\Delta_{A}=\begin{pmatrix}1&-1\\ 0&0\end{pmatrix}.$

Similarly, $D(\alpha,\mu,h)\approx D(\alpha,\mu_{0},h)+\mu\Delta_{D}$ where

 $\Delta_{D}=\begin{pmatrix}2(1-\alpha h)&0&-2(1-\alpha h)\\ 0&0&0\\ 1&0&-1\end{pmatrix}.$

This last approximate inequality is a first-order approximation, in the sense that we are working with arbitrarily small, positive values of $\mu$, and we have kept terms linear in $\mu$ but ignored higher order terms, like $\mu^{2}$.

We will apply the Bauer-Fike theorem to bound the eigenvalues of $D(\alpha,\mu,h)$. Consider the eigendecomposition $D(\alpha,\mu_{0},h)=V\Lambda V^{-1}$. We can compute

 $V=\begin{pmatrix}0&0&\frac{1-2\alpha h+3\alpha^{2}h^{2}}{1-\alpha h}\\ 0&1&\frac{1}{1-\alpha h}\\ 1&0&1\end{pmatrix}$

and

 $V^{-1}=\begin{pmatrix}\frac{1-\alpha h}{1-2\alpha h+3\alpha^{2}h^{2}}&-\frac{1% }{1-2\alpha h+3\alpha^{2}h^{2}}&0\\ 0&1&0\\ 1&0&0\end{pmatrix}.$

Note that because we assume $\alpha<2/(3h)$ we get $1-\alpha h>0$. Also, $1-2\alpha h+3\alpha^{2}h^{2}>0$ regardless of the choice of hyperparameters. This means that matrices $V$ and $V^{-1}$ are singular and of finite norm. The norm of $\Delta_{D}$ is also finite. The Bauer-Fike theorem states that, if $\nu$ is an eigenvalue of $D(\alpha,\mu_{0},h)$, then there exists an eigenvalue $\lambda$ of $D(\alpha,\mu,h)$ such that

 $|\lambda-\nu|\leq\|V\|_{p}\|V^{-1}\|_{p}\|\mu\Delta_{D}\|_{p},$

for any $p$-norm. Since by construction $|\nu|\leq\rho(D(\alpha,\mu_{0},h))<1$, the above means that there exists a sufficiently small, but strictly positive value of $\mu$, such that $\lambda<1$. By repeating this argument for all pairs of eigenvalues, we get the stated result. The same argument can be repeated to prove the existence of a strictly positive $\mu$ such that $\rho(A(\alpha,\mu,h))<1$.