Abstract
Most stochastic optimization methods use gradients once before discardingthem. While variance reduction methods have shown that reusing past gradientscan be beneficial when there is a finite number of datapoints, they do noteasily extend to the online setting. One issue is the staleness due to usingpast gradients. We propose to correct this staleness using the idea of implicitgradient transport (IGT) which transforms gradients computed at previousiterates into gradients evaluated at the current iterate without using theHessian explicitly. In addition to reducing the variance and bias of ourupdates over time, IGT can be used as a dropin replacement for the gradientestimate in a number of wellunderstood methods such as heavy ball or Adam. Weshow experimentally that it achieves stateoftheart results on a wide rangeof architectures and benchmarks. Additionally, the IGT gradient estimatoryields the optimal asymptotic convergence rate for online stochasticoptimization in the restricted setting where the Hessians of all componentfunctions are equal.
Quick Read (beta)
Reducing the variance in online optimization by transporting past gradients
Abstract
Most stochastic optimization methods use gradients once before discarding them. While variance reduction methods have shown that reusing past gradients can be beneficial when there is a finite number of datapoints, they do not easily extend to the online setting. One issue is the staleness due to using past gradients. We propose to correct this staleness using the idea of implicit gradient transport (IGT) which transforms gradients computed at previous iterates into gradients evaluated at the current iterate without using the Hessian explicitly. In addition to reducing the variance and bias of our updates over time, IGT can be used as a dropin replacement for the gradient estimate in a number of wellunderstood methods such as heavy ball or Adam. We show experimentally that it achieves stateoftheart results on a wide range of architectures and benchmarks. Additionally, the IGT gradient estimator yields the optimal asymptotic convergence rate for online stochastic optimization in the restricted setting where the Hessians of all component functions are equal.^{1}^{1} 1 Opensource implementation available at: https://github.com/seba1511/igt.pth
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 [email protected] PierreAntoine Manzagol Google Brain Montréal, QC [email protected] Reza Babanezhad University of British Columbia Vancouver, BC [email protected] Ioannis Mitliagkas Mila, Université de Montréal Montréal, QC [email protected] Nicolas Le Roux Mila, Google Brain Montréal, QC [email protected]
noticebox[b]Preprint. Under review.\[email protected]
1 Introduction
We wish to solve the following minimization problem:
${\theta}^{\ast}$  $=\mathrm{arg}\underset{\theta}{\mathrm{min}}{E}_{x\sim p}[f(\theta ,x)],$  (1) 
where we only have access to samples $x$ and to a firstorder 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 .^{2}^{2} 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 ; ShalevShwartz2013sdca ; 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:16595 . While these methods converge for a constant stepsize in the stochastic case^{3}^{3} 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 plugin 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 updates^{4}^{4} 4 This is slightly different from the standard formulation but equivalent for constant ${\gamma}_{t}$.:
${v}_{t}$  $={\gamma}_{t}{v}_{t1}+(1{\gamma}_{t})g({\theta}_{t},{x}_{t}),{v}_{0}=g({\theta}_{0},{x}_{0})$  
${\theta}_{t+1}$  $={\theta}_{t}{\alpha}_{t}{v}_{t}.$ 
When ${\gamma}_{t}=\gamma $, this leads to ${\theta}_{t+1}={\theta}_{t}{\alpha}_{t}\left({\gamma}^{t}g({\theta}_{0},{x}_{0})+(1\gamma ){\displaystyle \sum _{i=1}^{t}}{\gamma}^{ti}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 ${\mathrm{\Delta}}_{t}={\theta}_{t}{\theta}^{\ast}$ is a random variable. Denoting ${\u03f5}_{i}$ the noise at timestep $i$, i.e. $g({\theta}_{i},{x}_{i})=g({\theta}_{i})+{\u03f5}_{i}$, and writing ${\mathrm{\Delta}}_{t}E[{\mathrm{\Delta}}_{t}]=\alpha {\sum}_{i=0}^{t}{N}_{i,t}{\u03f5}_{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 ${\u03f5}_{i}$ on the iterates. Figure 1 shows the impact of ${\u03f5}_{i}$ on ${\mathrm{\Delta}}_{t}E[{\mathrm{\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 ${\u03f5}_{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 $ti$ 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 ${\u03f5}_{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\mathrm{min}\{\frac{1}{16\widehat{\kappa}},\frac{1}{8N}\}\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. ^{5}^{5} 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 nonconverging 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},\mathrm{\dots},{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}}_{t1}({\theta}_{t1})$ of ${g}_{t1}({\theta}_{t1})$; b) We can compute the gradient $g(\theta ,{x}_{t})$ at any location $\theta $. We shall seek a $\theta $ such that
$\frac{t}{t+1}}{\widehat{g}}_{t1}({\theta}_{t1})+{\displaystyle \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\mathit{}\mathrm{(}\mathrm{\cdot}\mathrm{,}x\mathrm{)}$ 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)
${g}_{t}({\theta}_{t})$  $={\displaystyle \frac{t}{t+1}}{g}_{t1}({\theta}_{t})+{\displaystyle \frac{1}{t+1}}g({\theta}_{t},{x}_{t})$  
$\approx {\displaystyle \frac{t}{t+1}}{\widehat{g}}_{t1}({\theta}_{t1})+{\displaystyle \frac{1}{t+1}}g({\theta}_{t}+t({\theta}_{t}{\theta}_{t1}),{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}_{t1})$. This extrapolation step is reminiscent of Nesterov’s acceleration with the difference that the factor in front of ${\theta}_{t}{\theta}_{t1}$, $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:
${\theta}_{t+1}$  $={\displaystyle \frac{2t+1}{t+1}}{\theta}_{t}{\displaystyle \frac{t}{t+1}}{\theta}_{t1}{\displaystyle \frac{\alpha}{t+1}}g({\theta}_{t}+t({\theta}_{t}{\theta}_{t1}),{x}_{t}).$  (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\mathit{}\mathrm{(}\theta \mathrm{,}x\mathrm{)}\mathrm{=}g\mathit{}\mathrm{(}\theta \mathrm{)}\mathrm{+}\u03f5$ with $\u03f5$ a random i.i.d. noise with covariance bounded by $B\mathit{}I$, then Eq. IGT with stepsize $\alpha \mathrm{=}\mathrm{1}\mathrm{/}L$ leads to iterates ${\theta}_{t}$ satisfying
$E[{\parallel {\theta}_{t}{\theta}^{\ast}\parallel}^{2}]$  $\le {\left(1{\displaystyle \frac{1}{\kappa}}\right)}^{2t}{\parallel {\theta}_{0}{\theta}^{\ast}\parallel}^{2}+{\displaystyle \frac{d{\alpha}^{2}B{\overline{\nu}}_{0}^{2}}{t}},$ 
with $\nu \mathrm{=}\mathrm{(}\mathrm{2}\mathrm{+}\mathrm{2}\mathit{}\mathrm{log}\mathit{}\kappa \mathrm{)}\mathit{}\kappa $ for every $t\mathrm{>}\mathrm{2}\mathit{}\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 plugin 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 dropin 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 (Nonstochastic).
In the nonstochastic case, where $B\mathrm{=}\mathrm{0}$, variance is equal to $\mathrm{0}$ and HeavyballIGT achieves the accelerated linear rate $O\mathit{}\mathrm{\left(}{\mathrm{\left(}\frac{\sqrt{\kappa}\mathrm{}\mathrm{1}}{\sqrt{\kappa}\mathrm{+}\mathrm{1}}\mathrm{\right)}}^{t}\mathrm{\right)}$ using the known, optimal heavy ball tuning, $\mu \mathrm{=}{\mathrm{\left(}\frac{\sqrt{\kappa}\mathrm{}\mathrm{1}}{\sqrt{\kappa}\mathrm{+}\mathrm{1}}\mathrm{\right)}}^{\mathrm{2}}$, $\alpha \mathrm{=}{\mathrm{(}\mathrm{1}\mathrm{+}\sqrt{\mu}\mathrm{)}}^{\mathrm{2}}\mathrm{/}L$.
Proposition 3.3 (Online, stochastic).
When $B\mathrm{>}\mathrm{0}$, there exist constant hyperparameters $\alpha \mathrm{>}\mathrm{0}$, $\mu \mathrm{>}\mathrm{0}$ such that ${\mathrm{\parallel}E\mathit{}\mathrm{[}{\theta}_{t}\mathrm{}{\theta}^{\mathrm{\ast}}\mathrm{]}\mathrm{\parallel}}^{\mathrm{2}}$ converges to zero linearly, and the variance is $\stackrel{\mathrm{~}}{O}\mathit{}\mathrm{(}\mathrm{1}\mathrm{/}t\mathrm{)}$.
The pseudocode can be found in Algorithm 3.3.
{algorithmic}[1] \ProcedureHeavyballIGTStepsize $\alpha $, Momentum $\mu $, Initial parameters ${\theta}_{0}$ \State${v}_{0}\leftarrow g({\theta}_{0},{x}_{0})\mathit{\hspace{1em}},{w}_{0}\leftarrow \alpha {v}_{0}\mathit{\hspace{1em}},{\theta}_{1}\leftarrow {\theta}_{0}+{w}_{0}$ \For$t=1,\mathrm{\dots},T1$ \State${\gamma}_{t}\leftarrow \frac{t}{t+1}$ \State${v}_{t}\leftarrow {\gamma}_{t}{v}_{t1}+(1{\gamma}_{t})g({\theta}_{t}+\frac{{\gamma}_{t}}{1{\gamma}_{t}}({\theta}_{t}{\theta}_{t1}),{x}_{t})$ \State${w}_{t}\leftarrow \mu {w}_{t1}\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:16595 .
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(t1)}{1+c(t1)}\left(1\frac{1}{c}\sqrt{\frac{1c}{t(t1)}}\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 ${\u03f5}_{t}$, where ${\u03f5}_{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 ResNet56 model DBLP:journals/corr/HeZRS15 on the CIFAR10 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 (HBITA). 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.
ImageNet image classification
We also consider the task of training a ResNet50 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).
IMDb sentiment analysis
We train a bidirectional LSTM on the IMDb Large Movie Review Dataset for 200 epochs. maas2011learning We observe that while the training convergence is comparable to HB, HBITA 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 AdamITA, which performs similarly to Adam.
6.2 Reinforcement learning
Linearquadratic regulator
We cast the classical linearquadratic 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 nonconvexity 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 trustregion methods in largescale applications. schulman2015trust ; openai2018learning
6.3 Metalearning
Modelagnostic metalearning
We now investigate the use of IGT in the modelagnostic metalearning (MAML) setting. finn2017model We replicate the 5 ways classification setup with 5 adaptation steps on tasks from the MiniImagenet dataset ravi2016optimization . This setting is interesting because of the many sources contributing to noise in the gradient estimates: the stochastic metagradient depends on the product of 5 stochastic Hessians computed over only 10 data samples, and is averaged over only 4 tasks. We substitute the metaoptimizer 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 metagradient counterpart, both in terms of convergence rate and final accuracy. Those results are also reflected in the final test accuracies where AdamITA ($65.16\%$) performs best, followed by HBITA ($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
 [1] The TensorFlow Authors. Tensorflow official resnet model. 2018.
 [2] 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.
 [3] Francis Bach and Eric Moulines. Nonstronglyconvex smooth stochastic approximation with convergence rate o (1/n). In Advances in Neural Information Processing Systems, pages 773–781, 2013.
 [4] Sébastien Bubeck. Convex optimization: Algorithms and complexity. Foundations and Trends® in Machine Learning, 8(34):231–357, 2015.
 [5] Aaron Defazio, Francis Bach, and Simon LacosteJulien. Saga: A fast incremental gradient method with support for nonstrongly convex composite objectives. In Advances in Neural Information Processing Systems, pages 1646–1654, 2014.
 [6] 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.
 [7] Aymeric Dieuleveut, Nicolas Flammarion, and Francis Bach. Harder, better, faster, stronger convergence rates for leastsquares regression. The Journal of Machine Learning Research, 18(1):3520–3570, 2017.
 [8] 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.
 [9] Chelsea Finn, Pieter Abbeel, and Sergey Levine. Modelagnostic metalearning for fast adaptation of deep networks. In Proceedings of the 34th International Conference on Machine LearningVolume 70, pages 1126–1135. JMLR. org, 2017.
 [10] Nicolas Flammarion and Francis Bach. From averaging to acceleration, there is only a stepsize. In Conference on Learning Theory, pages 658–695, 2015.
 [11] 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.
 [12] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. CoRR, abs/1512.03385, 2015.
 [13] 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.
 [14] Thomas Hofmann, Aurelien Lucchi, Simon LacosteJulien, and Brian McWilliams. Variance reduced stochastic gradient descent with neighbors. In Advances in Neural Information Processing Systems, pages 2305–2313, 2015.
 [15] 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.
 [16] Prateek Jain, Sham M. Kakade, Rahul Kidambi, Praneeth Netrapalli, and Aaron Sidford. Parallelizing stochastic gradient descent for least squares regression: Minibatching, averaging, and model misspecification. Journal of Machine Learning Research, 18(223):1–42, 2018.
 [17] Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In Advances in Neural Information Processing Systems, pages 315–323, 2013.
 [18] Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
 [19] Alex Krizhevsky. Learning multiple layers of features from tiny images. Technical report, 2009.
 [20] Alex Krizhevsky. Learning multiple layers of features from tiny images. Technical report, Citeseer, 2009.
 [21] Huibert Kwakernaak. Linear optimal control systems, volume 1. 1972.
 [22] Simon LacosteJulien, 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.
 [23] Nicolas Le Roux. Anytime tail averaging. arXiv preprint arXiv:1902.05083, 2019.
 [24] 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.
 [25] Yann LeCun and Corinna Cortes. Mnist handwritten digit database. AT&T Labs [Online]. Available: http://yann. lecun. com/exdb/mnist, 2010.
 [26] 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.
 [27] 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.
 [28] Julien Mairal. Optimization with firstorder surrogate functions. In Proceedings of The 30th International Conference on Machine Learning, pages 783–791, 2013.
 [29] 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 inhand manipulation. CoRR, abs/1808.00177, 2018.
 [30] Brendan O’Donoghue and Emmanuel Candes. Adaptive restart for accelerated gradient schemes. Foundations of computational mathematics, 15(3):715–732, 2015.
 [31] 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.
 [32] Boris T Polyak. Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics, 4(5):1–17, 1964.
 [33] Boris T Polyak and Anatoli B Juditsky. Acceleration of stochastic approximation by averaging. SIAM Journal on Control and Optimization, 30(4):838–855, 1992.
 [34] Sachin Ravi and Hugo Larochelle. Optimization as a model for fewshot learning. 2016.
 [35] Herbert Robbins and Sutton Monro. A stochastic approximation method. Annals of Mathematical Statistics, 22(3):400–407, 1951.
 [36] 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 FeiFei. ImageNet Large Scale Visual Recognition Challenge. International Journal of Computer Vision (IJCV), 115(3):211–252, 2015.
 [37] Mark Sandler, Andrew Howard, Menglong Zhu, Andrey Zhmoginov, and LiangChieh Chen. Mobilenetv2: Inverted residuals and linear bottlenecks. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 4510–4520, 2018.
 [38] Mark Schmidt. Convergence rate of stochastic gradient with constant step size. 2014.
 [39] Mark Schmidt, Nicolas Le Roux, and Francis Bach. Convergence rates of inexact proximalgradient methods for convex optimization. In Advances in Neural Information Processing Systems 24, 2011.
 [40] Mark Schmidt, Nicolas Le Roux, and Francis Bach. Minimizing finite sums with the stochastic average gradient. Mathematical Programming, 162(12):83–112, 2017.
 [41] John Schulman, Sergey Levine, Philipp Moritz, Michael I. Jordan, and Pieter Abbeel. Trust region policy optimization. CoRR, abs/1502.05477, 2015.
 [42] Shai ShalevShwartz and Tong Zhang. Stochastic dual coordinate ascent methods for regularized loss. Journal of Machine Learning Research, 14(1):567–599, February 2013.
 [43] Christopher J Shallue, Jaehoon Lee, Joe Antognini, Jascha SohlDickstein, Roy Frostig, and George E Dahl. Measuring the effects of data parallelism on neural network training. arXiv preprint arXiv:1811.03600, 2018.
 [44] 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.
 [45] Sharan Vaswani, Francis Bach, and Mark Schmidt. Fast and faster convergence of sgd for overparameterized models and an accelerated perceptron. In Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics, 2019.
 [46] HoiTo Wai, Wei Shi, Angelia Nedic, and Anna Scaglione. Curvatureaided incremental aggregated gradient method. arXiv preprint arXiv:1710.08936, 2017.
 [47] Ronald J Williams. Simple statistical gradientfollowing algorithms for connectionist reinforcement learning. Machine learning, 8(34):229–256, 1992.
 [48] 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.
 [49] Jian Zhang and Ioannis Mitliagkas. Yellowfin and the art of momentum tuning. In SysML, 2019.
 [50] 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 preprocessing steps,

•
training / validation / testing splits,

•
a description of the hyperparameter 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 [20] consists 50k training and 10k testing images, partitioned over 10 classes. We download and preprocess the images using the TensorFlow models package, available at the following URL: https://github.com/tensorflow/models
Model
Hyperparameters
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 [43]. 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[48] 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, 3e1, 1e1, 3e2, 1e2) for SGD, HB and the IGT variants and (1e2, 3e3, 1e3, 3e4, 1e4) 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 [43]. Note that selecting the hyperparameter combination reaching the lowest training loss yields qualitatively identical curves.
The resulting hyperparameters 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

•
HBITA stepsize 0.03, decay 0.1, tail fraction 18
Infrastructure and Runs
The experiments were run using P100 GPUs (single GPU).
Additional Results
We provide all learning curves for the methods comparison presented in the main manuscript in figure 6.
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 hyperparameters are:

•
ATA stepsize 0.3, decay 0.01, tail fraction 18

•
HBATA stepsize 0.03, decay 0.01, tail fraction 18
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 minibatch 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.
A.2 ImageNet image classification
Dataset
We use the 2015 edition of the ImageNet LargeScale Visual Recognition Challenge (ImageNet) [36] dataset. This dataset consists of 1.2M images partitioned into 1’000 classes. We use the preprocessing 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 [13].
Hyperparameters
We used the same setup and approach as for the CIFAR10 experiments. The setup trains for 90 epochs using minibatches 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 hyperparameters 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

•
HBITA stepsize 0.03, decay 0.01, tail fraction 90, decay step 0.9
Infrastructure and Runs
We ran these experiments using Google TPUv2.
Additional Results
We include error and accuracy curves for training and testing in Figure 9.
A.3 IMDb sentiment analysis
Dataset
The Internet Movie Database (IMDb) [27] 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 preprocessed with torchtext package, available at the following URL: https://github.com/pytorch/text More specifically, we tokenize the text at the wordlevel using the spaCy package, and embed the tokens using the 100dimensional GloVe 6B [31] distributed representations.
Model
The model consists of an embedding lookuptable, followed by a bidirectional LSTM with dropout, and then by a fullyconnected 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 lookuptable initilized with the GloVe vectors. The model is trained to minimize the BCEWithLogitsLoss with a minibatch size of 64.
Hyperparameters
For each method, we used a gridsearch 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 hyperparameters are displayed in the following table for each method.
HB  Adam  ASGD  HBIGT  HBITA  

$\alpha $  0.032  0.0005  0.064  0.128  0.064 
$\mu $  0.95  0.95  n/a  0.9  0.9 
$\xi $  n/a  n/a  100  n/a  n/a 
$\kappa $  n/a  n/a  ${10}^{5}$  n/a  n/a 
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.
Additional Results
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 1520 epochs.
A.4 Linearquadratic regulator
Setup
Our linearquadratic regulator [21] implements the following equations. The cost functional is evaluated at every timestep $h$ and is given by
$C({s}_{h},{a}_{h})={s}_{h}^{\top}Q{s}_{h}+{a}_{h}^{\top}R{a}_{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
${s}_{h+1}=A{s}_{h}+B{a}_{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}=K{s}_{h}+{\u03f5}_{h}^{a}$, where ${\u03f5}_{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 [47] given by
${\nabla}_{K}\mathbb{E}{\displaystyle \sum _{h}^{10}}C({s}_{h},{a}_{h})=\mathbb{E}\left({\displaystyle \sum _{h}^{10}}{\nabla}_{K}\mathrm{log}{\pi}_{K}({a}_{h}{s}_{h})\right)\left({\displaystyle \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.
Hyperparameters
Due to the simplicity of the considered methods, the only hyperparameter is the stepsize. For each method, we choose the stepsize from a logarithmicallyspaced 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 i75820K 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 fullgradient method (GD) about 10h, and computing the solution of the Riccati equation takes less than 5 seconds.
Additional Results
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.
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}_{t1})$ 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 Modelagnostic metalearning
Dataset
We use the MiniImagenet dataset [34] in our modelagnostic metalearning (MAML) [9] 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 metagradient 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 [9]. 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.
Hyperparameters
We only tune the metastepsize for the MAML experiment. We set the momentum constant to 0.9, the adaptationstepsize to 0.01, and average the metagradient of 4 tasks per iterations. Due to the reduced variance in the gradients, we found it necessary to increase the $\u03f5$ of AdamITA to 0.01.
For each method, we search over stepsize values on a logarithmicallyspaced grid and select those values that maximize validation accuracy after 10k metaiterations. These values are reported in Table 2.
HB  Adam  HBITA  AdamITA  
$\alpha $  0.008  0.001  0.008  0.0005 
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.
Additional Results
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: AdamITA reaches $65.16\%$, HBITA $64.57\%$, Adam $63.70\%$, and HB $63.08\%$.
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 HBIGT 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.
We investigate three different scenarios: (a) linearmnist: a logistic regression model on MNIST, (b) mnist: a modified version of LeNet5 [25] on MNIST and (c) cifarsmall: the MobileNetv2 architecture [37] consisting of 19 convolutional layers on CIFAR10. All models are trained with a minibatch size of 64, while the remaining hyperparameters are available in Tables 4, 5, and 6.
For each of the above tasks, models are trained for 200 epochs. We compare the following methods:

•
HB: the heavy ball method [33],

•
Adam [18],

•
ASGD [15],

•
SVRG [17],

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

•
HBIGT: the heavy ball using the IGT as a plugin estimator, and

•
HBITA: same as HBIGT 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, HBIGT performs on par with HBITA 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 HBIGT is not competitive with stateoftheart methods such as Adam or ASGD. HBITA, on the other hand, with its smaller reliance on the assumption, once again performs much better than all other methods. In fact, HBITA 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 TailAveraging mechanism is again apparent: without it, HeavyballIGT 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 HeavyballITA are competitive with the ones discovered by other optimization algorithms.
LinearMNIST  MNIST  CIFAR10  IMDb  

Heavyball  92.52 $\pm $ 0.04  99.08 $\pm $ 0.07  91.55 $\pm $ 0.25  86.90 $\pm $ 0.67 
Adam  92.57 $\pm $ 0.10  98.99 $\pm $ 0.05  89.36 $\pm $ 0.75  85.62 $\pm $ 0.63 
ASGD  92.47 $\pm $ 0.08  99.15 $\pm $ 0.07  91.45 $\pm $ 0.20  87.31 $\pm $ 0.21 
SVRG  92.51 $\pm $ 0.04  99.06 $\pm $ 0.08  86.84 $\pm $ 0.17  n/a 
SGDdec  92.52 $\pm $ 0.06  99.11 $\pm $ 0.06  87.53 $\pm $ 0.23  86.73 $\pm $ 0.34 
HeavyballIGT  92.47 $\pm $ 0.10  99.00 $\pm $ 0.05  12.05 $\pm $ 0.21  86.61 $\pm $ 0.23 
HeavyballITA  92.50 $\pm $ 0.10  99.19 $\pm $ 0.02  90.37 $\pm $ 0.31  87.26 $\pm $ 0.24 
HB  Adam  ASGD  HBIGT  HBITA  

$\alpha $  0.0128  0.0002  $0.0032$  $0.0032$  $0.0016$ 
$\mu $  0.1  0.95  n/a  0.9  0.1 
$\xi $  n/a  n/a  $10$  n/a  n/a 
$\kappa $  n/a  n/a  ${10}^{4}$  n/a  n/a 
HB  Adam  ASGD  HBIGT  HBITA  

$\alpha $  0.0064  0.0016  0.0128  0.0032  0.0032 
$\mu $  0.9  0.95  n/a  0.95  0.95 
$\xi $  n/a  n/a  $10$  n/a  n/a 
$\kappa $  n/a  n/a  ${10}^{4}$  n/a  n/a 
HB  Adam  ASGD  HBIGT  HBITA  

$\alpha $  0.0512  0.0512  0.1024  0.0128  0.0512 
$\mu $  0.95  0.9  n/a  0.9  0.1 
$\xi $  n/a  n/a  100  n/a  n/a 
$\kappa $  n/a  n/a  ${10}^{5}$  n/a  n/a 
Appendix C Proofs
C.1 Transport formula
${g}_{t}({\theta}_{t})$  $={\displaystyle \frac{t}{t+1}}{g}_{t1}({\theta}_{t})+{\displaystyle \frac{1}{t+1}}g({\theta}_{t},{x}_{t})$  
$={\displaystyle \frac{t}{t+1}}\left({g}_{t1}({\theta}_{t1})+H({\theta}_{t}{\theta}_{t1})\right)+{\displaystyle \frac{1}{t+1}}g({\theta}_{t},{x}_{t})$  (Quadratic $f$)  
$={\displaystyle \frac{t}{t+1}}{g}_{t1}({\theta}_{t1})+{\displaystyle \frac{1}{t+1}}\left(g({\theta}_{t},{x}_{t})+tH({\theta}_{t}{\theta}_{t1})\right)$  
$={\displaystyle \frac{t}{t+1}}{g}_{t1}({\theta}_{t1})+{\displaystyle \frac{1}{t+1}}g({\theta}_{t}+t({\theta}_{t}{\theta}_{t1}),{x}_{t})$  (Identical Hessians)  
$\approx {\displaystyle \frac{t}{t+1}}{\widehat{g}}_{t1}({\theta}_{t1})+{\displaystyle \frac{1}{t+1}}g({\theta}_{t}+t({\theta}_{t}{\theta}_{t1}),{x}_{t}).$  (${\widehat{g}}_{t1}$ is an approximation) 
Appendix D Proof of Prop. 3.1
In this proof, we assume that $g$ is a stronglyconvex quadratic function with hessian $H$.
At timestep $t$, we have access to a stochastic gradient $g(\theta ,{x}_{t})=g({\theta}_{t})+{\u03f5}_{t}$ where the ${\u03f5}_{t}$ are i.i.d. with variance $C\u2aaf{\sigma}^{2}H$.
We first prove a simple lemma:
Lemma D.1.
If ${v}_{\mathrm{0}}\mathrm{=}g\mathit{}\mathrm{(}{\theta}_{\mathrm{0}}\mathrm{)}\mathrm{+}{\u03f5}_{\mathrm{0}}$ and, for $t\mathrm{>}\mathrm{0}$, we have
${v}_{t}$  $={\displaystyle \frac{t}{t+1}}{v}_{t1}+{\displaystyle \frac{1}{t+1}}g\left({\theta}_{t}+t({\theta}_{t}{\theta}_{t1})\right)+{\displaystyle \frac{1}{t+1}}{\u03f5}_{t},$ 
then
${v}_{t}$  $=g({\theta}_{t})+{\displaystyle \frac{1}{t+1}}{\displaystyle \sum _{i=0}^{t}}{\u03f5}_{i}.$ 
Proof.
Per our assumption, this is true for $t=0$. Now let us prove the result by induction. Assume this is true for $t1$. Then we have:
${v}_{t}$  $={\displaystyle \frac{t}{t+1}}{v}_{t1}+{\displaystyle \frac{1}{t+1}}g({\theta}_{t}+t({\theta}_{t}{\theta}_{t1}))+{\displaystyle \frac{1}{t+1}}{\u03f5}_{t}$  
$={\displaystyle \frac{t}{t+1}}g({\theta}_{t1})+{\displaystyle \frac{1}{t+1}}{\displaystyle \sum _{i=0}^{t1}}{\u03f5}_{i}$  
$\mathrm{\hspace{1em}\hspace{1em}}+{\displaystyle \frac{1}{t+1}}g({\theta}_{t}+t({\theta}_{t}{\theta}_{t1}))+{\displaystyle \frac{1}{t+1}}{\u03f5}_{t}$  (recurrence assumption)  
$={\displaystyle \frac{t}{t+1}}g({\theta}_{t1})+{\displaystyle \frac{1}{t+1}}{\displaystyle \sum _{i=0}^{t1}}{\u03f5}_{i}$  
$\mathrm{\hspace{1em}\hspace{1em}}+g({\theta}_{t}){\displaystyle \frac{t}{t+1}}g({\theta}_{t1})+{\displaystyle \frac{1}{t+1}}{\u03f5}_{t}$  ($g$ is quadratic)  
$=g({\theta}_{t})+{\displaystyle \frac{1}{t+1}}{\displaystyle \sum _{i=0}^{t}}{\u03f5}_{i}.$ 
This concludes the proof. ∎
Lemma D.2.
Let us assume we perform the following iterative updates:
${v}_{t}$  $={\displaystyle \frac{t}{t+1}}{v}_{t1}+{\displaystyle \frac{1}{t+1}}g({\theta}_{t}+t({\theta}_{t}{\theta}_{t1}))+{\displaystyle \frac{1}{t+1}}{\u03f5}_{t}$  
${\theta}_{t+1}$  $={\theta}_{t}\alpha {v}_{t},$ 
starting from ${v}_{\mathrm{0}}\mathrm{=}g\mathit{}\mathrm{(}{\theta}_{\mathrm{0}}\mathrm{)}\mathrm{+}{\u03f5}_{\mathrm{0}}$. Then, denoting ${\mathrm{\Delta}}_{t}\mathrm{=}{\theta}_{t}\mathrm{}{\theta}^{\mathrm{\ast}}$, we have
${\mathrm{\Delta}}_{t}$  $={(I\alpha H)}^{t}{\mathrm{\Delta}}_{0}\alpha {\displaystyle \sum _{i=0}^{t1}}{N}_{i,t}{\u03f5}_{i}$ 
with
${N}_{i,0}$  $=0$  
${N}_{i,t}$  $$ 
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 $t1$. Using Lemma D.1, we have
${v}_{t1}$  $=g({\theta}_{t1})+{\displaystyle \frac{1}{t}}{\displaystyle \sum _{i=0}^{t1}}{\u03f5}_{i}$ 
and thus, using $g({\theta}_{t1})=H{\mathrm{\Delta}}_{t1}$,
${\mathrm{\Delta}}_{t}$  $={\mathrm{\Delta}}_{t1}\alpha {v}_{t1}$  
$={\mathrm{\Delta}}_{t1}\alpha H{\mathrm{\Delta}}_{t1}{\displaystyle \frac{\alpha}{t}}{\displaystyle \sum _{i=0}^{t1}}{\u03f5}_{i}$  
$=(I\alpha H){\mathrm{\Delta}}_{t1}{\displaystyle \frac{\alpha}{t}}{\displaystyle \sum _{i=0}^{t1}}{\u03f5}_{i}$  
$={(I\alpha H)}^{t}{\mathrm{\Delta}}_{0}\alpha {\displaystyle \sum _{i=0}^{t2}}(I\alpha H){N}_{i,t1}{\u03f5}_{i}{\displaystyle \frac{\alpha}{t}}{\displaystyle \sum _{i=0}^{t1}}{\u03f5}_{i}$  (recurrence assumption)  
$={(I\alpha H)}^{t}{\mathrm{\Delta}}_{0}\alpha {\displaystyle \sum _{i=0}^{t1}}{N}_{i,t}{\u03f5}_{i}$ 
with
${N}_{i,t}$  $$ 
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 [30].
Lemma D.3.
Denote ${r}_{h}\mathrm{=}\mathrm{1}\mathrm{}\alpha \mathit{}h$. We assume $\alpha \mathrm{\le}\frac{\mathrm{1}}{L}$. Then, for any $i$ and any $t$, we have
${N}_{i,t}$  $\ge 0$  (Positivity)  
${N}_{i,t}$  $=0\mathit{\hspace{1em}}\mathit{\text{if}}t\le i$  (Zerostart)  
${N}_{i,t}$  $$  (Constant bound)  
${N}_{i,t}$  $\le {\displaystyle \frac{\mathrm{max}\{1+{r}_{h},2\mathrm{log}\left(\frac{2}{i(1{r}_{h})}\right)\}}{t(1{r}_{h})}}\mathit{\hspace{1em}}\mathit{\text{if}}{\displaystyle \frac{2}{1{r}_{h}}}\le t.$  (Decreasing bound) 
Proof.
The Zerostart case $i\ge 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$,
${N}_{i,t}$  $={r}_{h}{N}_{i,t1}+{\displaystyle \frac{1}{t}}$  
$\le {N}_{i,t1}+{\displaystyle \frac{1}{t}}.$ 
Thus, ${N}_{i,t}{N}_{i,t1}\le \frac{1}{t}$. Summing these inequalities, we get a telescopic sum and, finally:
${N}_{i,t}$  $\le {\displaystyle \sum _{j=i+1}^{t}}{\displaystyle \frac{1}{j}}$  
$\le {\displaystyle {\int}_{x=i}^{t}}{\displaystyle \frac{dx}{x}}$  
$=\mathrm{log}\left({\displaystyle \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
${N}_{0,t}$  $\le 1+\mathrm{log}t.$ 
In the remainder, we shall keep the $\mathrm{log}\left(\frac{t}{i}\right)$ bound for simplicity. The upper bound on the righthand 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
${N}_{i,\frac{2}{1{r}_{h}}}$  $\le \mathrm{log}\left({\displaystyle \frac{\frac{2}{1{r}_{h}}}{i}}\right)$  
$=\mathrm{log}\left({\displaystyle \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
${N}_{i,\frac{2}{1{r}_{h}}}$  $\le \mathrm{log}\left({\displaystyle \frac{2}{i(1{r}_{h})}}\right){\displaystyle \frac{\frac{2}{1{r}_{h}}}{\frac{2}{1{r}_{h}}}}$  
$={\displaystyle \frac{\mathrm{log}\left(\frac{2}{i(1{r}_{h})}\right)}{1{r}_{h}}}{\displaystyle \frac{2}{\frac{2}{1{r}_{h}}}}$  
$\le {\displaystyle \frac{\mathrm{max}\{1+{r}_{h},2\mathrm{log}\left(\frac{2}{i(1{r}_{h})}\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
${N}_{i,t}$  $\le {\displaystyle \frac{\mathrm{max}\{1+{r}_{h},2\mathrm{log}\left(\frac{1}{i(1{r}_{h})}\right)\}}{t(1{r}_{h})}}$  
$\le {\displaystyle \frac{{\nu}_{i}}{t}},$ 
with ${\nu}_{i}=\frac{\mathrm{max}\{1+{r}_{h},2\mathrm{log}\left(\frac{1}{i(1{r}_{h})}\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,t1}\le \frac{{\nu}_{i}}{t1}$, then ${N}_{i,t}\le \frac{{\nu}_{i}}{t}$. Assume that there is such at $t$. Then
${N}_{i,t}$  $={r}_{h}{N}_{i,t1}+{\displaystyle \frac{1}{t}}$  
$\le {\displaystyle \frac{{r}_{h}{\nu}_{i}}{t1}}+{\displaystyle \frac{1}{t}}$  
$={\displaystyle \frac{{r}_{h}t{\nu}_{i}+t1}{t(t1)}}$  
$={\displaystyle \frac{(t1){\nu}_{i}+({r}_{h}1)t{\nu}_{i}+{\nu}_{i}+t1}{t(t1)}}$  
$={\displaystyle \frac{{\nu}_{i}}{t}}+{\displaystyle \frac{({r}_{h}1)t{\nu}_{i}+{\nu}_{i}+t1}{t(t1)}}.$ 
We now shall prove that $({r}_{h}1)t{\nu}_{i}+{\nu}_{i}+t1=[({r}_{h}1){\nu}_{i}+1]t+{\nu}_{i}1$ is negative. First, we have that
$({r}_{h}1){\nu}_{i}+1$  $=1\mathrm{max}\{1+{r}_{h},2\mathrm{log}\left({\displaystyle \frac{1}{i(1{r}_{h})}}\right)\}$  
$\le 0.$ 
Then,
$[({r}_{h}1){\nu}_{i}+1]t+{\nu}_{i}1\le 0$  $\u27fat\ge {\displaystyle \frac{{\nu}_{i}1}{(1{r}_{h}){\nu}_{i}1}}$ 
since $({r}_{h}1){\nu}_{i}+1\le 0$. Thus, the property is true for every $t\ge \frac{{\nu}_{i}1}{(1{r}_{h}){\nu}_{i}1}$. In addition, we have
${\nu}_{i}$  $\ge {\displaystyle \frac{1+{r}_{h}}{1{r}_{h}}}$  
${\nu}_{i}(1{r}_{h})$  $\ge 1+{r}_{h}$  
$2{\nu}_{i}(1{r}_{h})2$  $\ge {\nu}_{i}(1{r}_{h})1+{r}_{h}$  
$\frac{2}{1{r}_{h}}$  $\ge {\displaystyle \frac{{\nu}_{i}1}{{\nu}_{i}(1{r}_{h})1}},$ 
and the property is also true for every $t\ge \frac{2}{1{r}_{h}}$. This concludes the proof. ∎
Finally, we can prove the Proposition 3.1:
Proof.
The expectation of ${\mathrm{\Delta}}_{t}$ is immediate using Lemma D.2 and the fact that the ${\u03f5}_{i}$ are independent, zeromean noises. The variance is equal to $V[{\mathrm{\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
${\overline{\nu}}_{i}$  $={\displaystyle \frac{\mathrm{max}\{2\alpha \mu ,2\mathrm{log}\left(\frac{1}{i\alpha \mu}\right)\}}{\alpha \mu}}\mathit{\hspace{1em}}\text{for}i0$  
${\overline{\nu}}_{0}$  $={\displaystyle \frac{2+2\mathrm{log}\left(\frac{1}{\alpha \mu}\right)}{\alpha \mu}},$ 
which is, for every $i$, the maximum ${\nu}_{i}$ across all dimensions. We get
$V[{\mathrm{\Delta}}_{t}]$  $\le d{\alpha}^{2}B{\displaystyle \sum _{i=0}^{t}}{\displaystyle \frac{{\overline{\nu}}_{i}^{2}}{{t}^{2}}}$  
$\le d{\alpha}^{2}B{\displaystyle \sum _{i=0}^{t}}{\displaystyle \frac{{\overline{\nu}}_{0}^{2}}{{t}^{2}}}\mathit{\hspace{1em}}\text{since}{\nu}_{i}\ge {\nu}_{i+1}\forall i$  
$\le {\displaystyle \frac{d{\alpha}^{2}B{\overline{\nu}}_{0}^{2}}{t}}.$ 
Since we have
$E[{\theta}_{t}{\theta}^{\ast}]$  $={(I\alpha H)}^{t}({\theta}_{0}{\theta}^{\ast}),$ 
we get
$E[{\parallel {\theta}_{t}{\theta}^{\ast}\parallel}^{2}]$  $={\parallel E[{\theta}_{t}{\theta}^{\ast}]\parallel}^{2}+V[{\mathrm{\Delta}}_{t}]$  
$\le {({\theta}_{0}{\theta}^{\ast})}^{\top}{(I\alpha H)}^{2t}({\theta}_{0}{\theta}^{\ast})+{\displaystyle \frac{d{\alpha}^{2}B{\overline{\nu}}_{0}^{2}}{t}}$  
$\le {\left(1{\displaystyle \frac{1}{\kappa}}\right)}^{2t}{\parallel {\theta}_{0}{\theta}^{\ast}\parallel}^{2}+{\displaystyle \frac{d{\alpha}^{2}B{\overline{\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\mathit{}\mathrm{(}\theta \mathrm{,}x\mathrm{)}\mathrm{=}g\mathit{}\mathrm{(}\theta \mathrm{)}\mathrm{+}\u03f5$ with $\u03f5$ a random uncorrelated noise with covariance bounded by $B\mathit{}I$.
Then, Algorithm 3.3 leads to iterates ${\theta}_{t}$ satisfying
$E[{\theta}_{t}{\theta}^{\ast}]$  $=\left(\begin{array}{c}\hfill I\hfill \\ \hfill 0\hfill \end{array}\right){A}^{t}\left(\begin{array}{c}\hfill E[{\theta}_{1}{\theta}^{\ast}]\hfill \\ \hfill E[{\theta}_{0}{\theta}^{\ast}]\hfill \end{array}\right)$  (5) 
where
$A$  $=\left(\begin{array}{cc}\hfill I\alpha H+\mu I\hfill & \hfill \mu I\hfill \\ \hfill I\hfill & \hfill 0\hfill \end{array}\right)$  (6) 
governs the dynamics of this bias. In particular, when its spectral radius, $\rho \mathit{}\mathrm{(}A\mathrm{)}$ is less than $\mathrm{1}$, the iterates converge linearly to ${\theta}^{\mathrm{\ast}}$.
In a similar fashion, the variance dynamics of HeavyballIGT are governed by the matrix
${D}_{i}$  $=\left(\begin{array}{ccc}\hfill {(1\alpha {h}_{i}+\mu )}^{2}+2{\alpha}^{2}{h}_{i}^{2}\hfill & \hfill {\mu}^{2}\hfill & \hfill 2\mu {(1\alpha {h}_{i}+\mu )}^{2}\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\alpha {h}_{i}+\mu \hfill & \hfill 0\hfill & \hfill \mu \hfill \end{array}\right)$ 
If the spectral radius of ${D}_{i}$, $\rho \mathit{}\mathrm{(}{D}_{i}\mathrm{)}$, is strictly less than 1 or all $i$, then there exist constants ${t}_{\mathrm{0}}\mathrm{>}\mathrm{0}$ and $C\mathrm{>}\mathrm{0}$ for which
$\mathrm{Var}({\theta}_{t})$  $\le 2{\alpha}^{2}dBC{\displaystyle \frac{\mathrm{log}(t)}{t}},\text{\mathit{f}\mathit{o}\mathit{r}}t>{t}_{0}$ 
where $B$ is a bound on the variance of noise variables ${\u03f5}_{i}$.
Lemma E.2 (IGT estimator as true gradient plus noise average).
If ${v}_{\mathrm{0}}\mathrm{=}g\mathit{}\mathrm{(}{\theta}_{\mathrm{0}}\mathrm{)}\mathrm{+}{\u03f5}_{\mathrm{0}}$ and for $t\mathrm{>}\mathrm{0}$ we have
$${v}_{t}=\frac{t}{t+1}{v}_{t1}+\frac{1}{t+1}g({\theta}_{t}+t({\theta}_{t}{\theta}_{t1}))+\frac{1}{t+1}{\u03f5}_{t},$$ 
then
$${v}_{t}=g({\theta}_{t})+\frac{1}{t+1}\sum _{i=0}^{t}{\u03f5}_{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 ${\u03f5}_{\mathrm{0}}\mathrm{,}{\u03f5}_{\mathrm{1}}\mathrm{,}\mathrm{\dots}\mathrm{,}{\u03f5}_{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}+{\u03f5}_{0}\right]=g({\theta}_{0}).$$ 
For the inductive case, we can write
$\mathbb{E}\left[{v}_{t}\right]$  $=\mathbb{E}\left[{\displaystyle \frac{t}{t+1}}{v}_{t1}+{\displaystyle \frac{1}{t+1}}\widehat{g}({\theta}_{t}+t({\theta}_{t}{\theta}_{t1}))\right]$  
$=\mathbb{E}\left[{\displaystyle \frac{t}{t+1}}{v}_{t1}+{\displaystyle \frac{1}{t+1}}{g}_{t}+{\displaystyle \frac{t}{t+1}}{g}_{t}{\displaystyle \frac{t}{t+1}}{g}_{t1}+{\displaystyle \frac{1}{t+1}}{\u03f5}_{t}\right]$  
$={\displaystyle \frac{t}{t+1}}\mathbb{E}\left[{v}_{t1}{g}_{t1}\right]+\mathbb{E}\left[{g}_{t}\right]+{\displaystyle \frac{t}{t+1}}\mathbb{E}\left[{\u03f5}_{t}\right]$  
$=\mathbb{E}\left[{g}_{t}\right]=g(\mathbb{E}\left[{\theta}_{t}\right]).$ 
Where, in the third equality, $\mathbb{E}\left[{v}_{t1}{g}_{t1}\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]\le 2{h}^{2}\mathrm{Var}\left[{\theta}_{t}{\theta}^{\star}\right]+\frac{2B}{t},$$ 
where $B$ is the variance of the homoscedastic noise ${\u03f5}_{t}$.
Proof.
$\mathrm{Var}\left[{v}_{t}\right]$  $=\mathrm{Var}\left[{g}_{t}+{\displaystyle \frac{1}{t+1}}{\displaystyle \sum _{i=0}^{t}}{\u03f5}_{i}\right]$  
$=\mathrm{Var}\left[h{\theta}_{t}\right]+\mathrm{Var}\left[{\displaystyle \frac{1}{t+1}}{\displaystyle \sum _{i=0}^{t}}{\u03f5}_{i}\right]$  
$\mathrm{\hspace{1em}\hspace{1em}}+2\mathrm{C}\mathrm{o}\mathrm{v}[h{\theta}_{t},{\displaystyle \frac{1}{t+1}}{\displaystyle \sum _{i=0}^{t}}{\u03f5}_{i}]$  
$\le 2\mathrm{V}\mathrm{a}\mathrm{r}\left[h{\theta}_{t}\right]+2\mathrm{V}\mathrm{a}\mathrm{r}\left[{\displaystyle \frac{1}{t+1}}{\displaystyle \sum _{i=0}^{t}}{\u03f5}_{i}\right]$  
$=2{h}^{2}\mathrm{Var}\left[{\theta}_{t}{\theta}^{\star}\right]+2{\displaystyle \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 HeavyballIGT. We use the quadratic assumption to decouple the vector dynamics of HeavyballIGT into independent scalar dynamics. If the Hessian, $H$, has eigenvalues $L\ge {h}_{1}\ge {h}_{2}\ge \mathrm{\dots}\ge {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}_{\mathrm{0}}\mathrm{=}g\mathit{}\mathrm{(}{\theta}_{\mathrm{0}}\mathrm{)}\mathrm{+}{\u03f5}_{\mathrm{0}}$ and ${w}_{\mathrm{0}}\mathrm{=}\mathrm{0}$, performing the following iterative updates (HeavyballIGT, Algorithm 3.3):
$${v}_{t}=\frac{t}{t+1}{v}_{t1}+\frac{1}{t+1}g(\theta +t({\theta}_{t}{\theta}_{t1}))+\frac{1}{t+1}{\u03f5}_{t},$$ 
$${w}_{t+1}=\mu {w}_{t}+\alpha {v}_{t},{\theta}_{t+1}={\theta}_{t}{w}_{t+1}$$ 
results in
$${\mathrm{\Delta}}_{t}={A}^{t}{\mathrm{\Delta}}_{0}\alpha \sum _{i=0}^{t1}{N}_{i,t}\left[\begin{array}{c}\hfill {\u03f5}_{i}\hfill \\ \hfill 0\hfill \end{array}\right]$$ 
where $$,
${\mathrm{\Delta}}_{t}=\left[\begin{array}{c}\hfill {\theta}_{t}{\theta}^{\ast}\hfill \\ \hfill {\theta}_{t1}{\theta}^{\ast}\hfill \end{array}\right]$ and $A\mathrm{=}\mathrm{\left(}\begin{array}{cc}\hfill \mathrm{1}\mathrm{}\alpha \mathit{}h\mathrm{+}\mu \hfill & \hfill \mathrm{}\mu \hfill \\ \hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill \end{array}\mathrm{\right)}\mathrm{.}$
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 $t1$, and start by using Lemma E.2:
${\mathrm{\Delta}}_{t}$  $=A{\mathrm{\Delta}}_{t1}{\displaystyle \frac{\alpha}{t}}{\displaystyle \sum _{i=0}^{t1}}\left[\begin{array}{c}\hfill {\u03f5}_{i}\hfill \\ \hfill 0\hfill \end{array}\right]$  
$=A({A}^{t1}{\mathrm{\Delta}}_{0}\alpha {\displaystyle \sum _{i=0}^{t2}}{N}_{i,t}\left[\begin{array}{c}\hfill {\u03f5}_{i}\hfill \\ \hfill 0\hfill \end{array}\right]){\displaystyle \frac{\alpha}{t}}{\displaystyle \sum _{i=0}^{t1}}\left[\begin{array}{c}\hfill {\u03f5}_{i}\hfill \\ \hfill 0\hfill \end{array}\right]$  (Inductive assumption)  
$={A}^{t}{\mathrm{\Delta}}_{0}\alpha ({\displaystyle \sum _{i=0}^{t2}}A{N}_{i,t}\left[\begin{array}{c}\hfill {\u03f5}_{i}\hfill \\ \hfill 0\hfill \end{array}\right]+{\displaystyle \frac{1}{t}}{\displaystyle \sum _{i=0}^{t1}}\left[\begin{array}{c}\hfill {\u03f5}_{i}\hfill \\ \hfill 0\hfill \end{array}\right])$  
$={A}^{t}{\mathrm{\Delta}}_{0}\alpha {\displaystyle \sum _{i=0}^{t1}}{N}_{i,t}\left[\begin{array}{c}\hfill {\u03f5}_{i}\hfill \\ \hfill 0\hfill \end{array}\right]$  (Def. of ${N}_{i,t}$) 
∎
Lemma E.6 (Evolution of variance).
Let ${U}_{t}\mathrm{=}\mathrm{Var}\mathit{}\mathrm{\left[}{\theta}_{t}\mathrm{\right]}$ and ${V}_{t}\mathrm{=}\mathrm{Cov}\mathit{}\mathrm{[}{\theta}_{t}\mathrm{,}{\theta}_{t\mathrm{}\mathrm{1}}\mathrm{]}$, where ${\theta}_{t}$ is the $t$th iterate of HeavyballIGT on a 1dimensional quadratic function with curvature $h$. The following matrix describes the variance dynamics of HeavyballIGT.
$$D=\left(\begin{array}{ccc}\hfill {(1\alpha h+\mu )}^{2}+2{\alpha}^{2}{h}^{2}\hfill & \hfill {\mu}^{2}\hfill & \hfill 2\mu {(1\alpha h+\mu )}^{2}\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\alpha h+\mu \hfill & \hfill 0\hfill & \hfill \mu \hfill \end{array}\right)$$  (7) 
If the spectral radius of $D$, $\rho \mathit{}\mathrm{(}D\mathrm{)}$, is strictly less than 1, then there exist constants ${t}_{\mathrm{0}}\mathrm{>}\mathrm{0}$ and $C\mathrm{>}\mathrm{0}$ for which
$\mathrm{Var}({\theta}_{t})$  $\le 2{\alpha}^{2}BC{\displaystyle \frac{\mathrm{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 [49]. We start by expanding ${U}_{t+1}$ as follows.
${U}_{t+1}$  $=\mathbb{E}\left[{({\theta}_{t+1}{\overline{\theta}}_{t+1})}^{2}\right]$  
$=\mathbb{E}\left[{({\theta}_{t}\alpha {v}_{t}+\mu ({\theta}_{t}{\theta}_{t1}){\overline{\theta}}_{t}+\alpha {g}_{t}\mu ({\overline{\theta}}_{t}{\overline{\theta}}_{t1}))}^{2}\right]$  
$=\mathbb{E}[({\theta}_{t}\alpha {g}_{t}+\mu ({\theta}_{t}{\theta}_{t1}){\overline{\theta}}_{t}+\alpha {g}_{t}$  
${\mathrm{\hspace{1em}\hspace{1em}}\mu ({\overline{\theta}}_{t}{\overline{\theta}}_{t1})+\alpha ({g}_{t}{v}_{t}))}^{2}]$  
$=\mathbb{E}\left[{((1\alpha h+\mu )({\theta}_{t}{\overline{\theta}}_{t})\mu ({\theta}_{t1}{\overline{\theta}}_{t1}))}^{2}\right]$  
$\mathrm{\hspace{1em}\hspace{1em}}+{\alpha}^{2}\mathbb{E}\left[{({g}_{t}{v}_{t})}^{2}\right]$  
$\le \mathbb{E}\left[{((1\alpha h+\mu )({\theta}_{t}{\overline{\theta}}_{t})\mu ({\theta}_{t1}{\overline{\theta}}_{t1}))}^{2}\right]$  
$\mathrm{\hspace{1em}\hspace{1em}}+{\alpha}^{2}\left(2{h}^{2}\mathbb{E}\left[{({\theta}_{t}{\overline{\theta}}_{t})}^{2}\right]+{\displaystyle \frac{2B}{t+1}}\right)$  
$\le \left[{(1\alpha h+\mu )}^{2}+2{\alpha}^{2}{\mu}^{2}\right]\mathbb{E}\left[{({\theta}_{t}{\overline{\theta}}_{t})}^{2}\right]$  
$\mathrm{\hspace{1em}\hspace{1em}}2\mu (1\alpha h+\mu )\mathbb{E}\left[({\theta}_{t}{\overline{\theta}}_{t})({\theta}_{t1}{\overline{\theta}}_{t1})\right]$  
$\mathrm{\hspace{1em}\hspace{1em}}+{\mu}^{2}\mathbb{E}\left[{({\theta}_{t1}{\overline{\theta}}_{t1})}^{2}\right]+{\alpha}^{2}{\displaystyle \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}$.
${V}_{t}$  $=\mathbb{E}\left[({\theta}_{t}{\overline{\theta}}_{t})({\theta}_{t1}{\overline{\theta}}_{t1})\right]$  
$=\mathbb{E}[((1\alpha h+\mu )({\theta}_{t1}{\overline{\theta}}_{t1})\mu ({\theta}_{t2}\overline{{\theta}_{t2}})+\alpha ({g}_{t}{v}_{t}))$  
$\mathrm{\hspace{1em}\hspace{1em}}({\theta}_{t1}{\overline{\theta}}_{t1})]$  
$=(1\alpha h+\mu )\mathbb{E}\left[{({\theta}_{t1}{\overline{\theta}}_{t1})}^{2}\right]$  
$\mathrm{\hspace{1em}\hspace{1em}}\mu \mathbb{E}\left[({\theta}_{t1}{\overline{\theta}}_{t1})({\theta}_{t2}{\overline{\theta}}_{t2})\right]$ 
From the above expressions, we obtain
$\left(\begin{array}{c}\hfill {U}_{t+1}\hfill \\ \hfill {U}_{t}\hfill \\ \hfill {V}_{t+1}\hfill \end{array}\right)$  $\le D\left(\begin{array}{c}\hfill {U}_{t}\hfill \\ \hfill {U}_{t1}\hfill \\ \hfill {V}_{t}\hfill \end{array}\right)+\left(\begin{array}{c}\hfill {\alpha}^{2}\frac{2B}{t+1}\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right)$  
$\le 2{\alpha}^{2}B{\displaystyle \sum _{i=0}^{t}}{D}^{i}\left(\begin{array}{c}\hfill \frac{1}{t+1i}\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right)$  
$\le 2{\alpha}^{2}B\left({\displaystyle \sum _{i=0}^{s1}}{D}^{i}\left(\begin{array}{c}\hfill \frac{1}{t+1i}\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right)+{\displaystyle \sum _{i=s}^{t}}{D}^{i}\left(\begin{array}{c}\hfill \frac{1}{t+1i}\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right)\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
${\left(\begin{array}{c}\hfill 1\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right)}^{T}{\displaystyle \sum _{i=0}^{s1}}{D}^{i}\left(\begin{array}{c}\hfill \frac{1}{t+1i}\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right)$  $\le {C}^{\prime}{\displaystyle \sum _{i=0}^{s1}}{\displaystyle \frac{1}{t+1i}}$  
$\le {\displaystyle \frac{{C}^{\prime}s}{t+2s}}$ 
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, $$ and
${\left(\begin{array}{c}\hfill 1\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right)}^{T}{\displaystyle \sum _{i=s}^{t}}{D}^{i}\left(\begin{array}{c}\hfill \frac{1}{t+1i}\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right)\le $  ${\left(\begin{array}{c}\hfill 1\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right)}^{T}{\displaystyle \sum _{i=s}^{t}}{D}^{i}\left(\begin{array}{c}\hfill 1\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right)$  
$\le $  ${C}^{\prime \prime}{\displaystyle \sum _{i=s}^{t}}{(\rho (D)+\zeta )}^{s}$  
$=$  ${C}^{\prime \prime}(ts+1){(\rho (D)+\zeta )}^{s}$ 
Let ${\rho}^{\prime}=\rho (D)+\zeta $ and $s=\lceil 2{\mathrm{log}}_{1/{\rho}^{\prime}}t\rceil $. Then ${(\rho (D)+\zeta )}^{s}=1/{t}^{2}$, and putting the above two bounds together,
${U}_{t+1}$  $\le 2{\alpha}^{2}B\left({\displaystyle \frac{2{C}^{\prime}{\mathrm{log}}_{1/{\rho}^{\prime}}t}{t+22{\mathrm{log}}_{1/{\rho}^{\prime}}t}}+{C}^{\prime \prime}{\displaystyle \frac{t2{\mathrm{log}}_{1/{\rho}^{\prime}}t+1}{{t}^{2}}}\right)$  
$\le 2{\alpha}^{2}BC{\displaystyle \frac{\mathrm{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
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}=\left(\begin{array}{cc}\hfill 1\alpha {h}_{i}+\mu \hfill & \hfill \mu \hfill \\ \hfill 1\hfill & \hfill 0\hfill \end{array}\right).$$ 
To prove that $$ it suffices to prove that $$ 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 $$ and ${\mu}_{0}=0$. Using hyperparameters $(\alpha ,{\mu}_{0})$ one gets the results for gradient descent without momentum. In particular $$, and the spectral radius of $D$ is $$.
We will argue that there exists $\mu >0$, such that $$, and the spectral radius of $D$ is $$. Then the previous lemma implies that bias converges linearly, and variance is $O(\mathrm{log}(t)/t)$.
To argue the existence of $\mu >0$, we will perform eigenvalue perturbation analysis using the BauerFike theorem. Note that $A(\alpha ,\mu ,h)=A(\alpha ,{\mu}_{0},h)+\mu {\mathrm{\Delta}}_{A}$ where
$${\mathrm{\Delta}}_{A}=\left(\begin{array}{cc}\hfill 1\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right).$$ 
Similarly, $D(\alpha ,\mu ,h)\approx D(\alpha ,{\mu}_{0},h)+\mu {\mathrm{\Delta}}_{D}$ where
$${\mathrm{\Delta}}_{D}=\left(\begin{array}{ccc}\hfill 2(1\alpha h)\hfill & \hfill 0\hfill & \hfill 2(1\alpha h)\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right).$$ 
This last approximate inequality is a firstorder 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 BauerFike theorem to bound the eigenvalues of $D(\alpha ,\mu ,h)$. Consider the eigendecomposition $D(\alpha ,{\mu}_{0},h)=V\mathrm{\Lambda}{V}^{1}$. We can compute
$$V=\left(\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill \frac{12\alpha h+3{\alpha}^{2}{h}^{2}}{1\alpha h}\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill \frac{1}{1\alpha h}\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right)$$ 
and
$${V}^{1}=\left(\begin{array}{ccc}\hfill \frac{1\alpha h}{12\alpha h+3{\alpha}^{2}{h}^{2}}\hfill & \hfill \frac{1}{12\alpha h+3{\alpha}^{2}{h}^{2}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right).$$ 
Note that because we assume $$ we get $1\alpha h>0$. Also, $12\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 ${\mathrm{\Delta}}_{D}$ is also finite. The BauerFike 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 \le {\parallel V\parallel}_{p}{\parallel {V}^{1}\parallel}_{p}{\parallel \mu {\mathrm{\Delta}}_{D}\parallel}_{p},$$ 
for any $p$norm. Since by construction $$, the above means that there exists a sufficiently small, but strictly positive value of $\mu $, such that $$. 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 $$.
∎