1 Introduction and Main Result
Training machine learning models very often involves minimizing a scalar function of a large number of real parameters. This function is usually called the loss function or cost function in machine learning, and often called the objective function in the optimization literature. Arguably, the most popular class of methods for minimization are firstorder methods based on gradient descent (LeCun et al., 2015; Ruder, 2016; Goodfellow et al., 2016)
. For example, for training (deep) neural networks, these are the preferred class of methods because gradients of the loss function with respect to parameters can be easily calculated via backpropagation
(Rumelhart et al., 1986; LeCun et al., 1998; Nielsen, 2015), whereas calculating second derivatives or the Hessian is expensive.The algorithm for basic gradient descent is the iteration
(1) 
Here is a
dimensional vector, consisting of the parameters of the machine learning model. The superscripts represent the iteration count. The objective function to be minimized,
, is a function known as the loss function or cost function. The positive constant is the learning rate. It is often modified during minimization, i.e., could depend on the iteration count . In this work, we will mainly focus on constant .We consider a commonly used variation of gradient descent:
(2) 
This is generally known as gradient descent with momentum, and also as Polyak’s heavy ball method (Polyak, 1964). The direction of descent at each step is a linear combination of the direction of the gradient and the direction of the previous step. The hyperparameter provides an ‘intertia’. Although some suggestions for varying can be found in the literature (Srinivasan et al., 2018; Chen and Kyrillidis, 2019), the value is so widely used that we will mostly focus on this value unless explicitly stated. Momentum is often used by evaluating the gradient not at the current position in parameter space, but rather at the estimated position at the next iteration:
(3) 
This modification is commonly known as Nesterov’s accelerated momentum (Nesterov, 1983). The original algorithm of Nesterov (1983) looks somewhat different, but the form (3), introduced in Sutskever et al. (2013), is what is more commonly known in current machine learning literature as Nesterov acceleration. The words ‘momentum’ and ‘acceleration’ are both used in a different sense from their meaning in physics/mechanics.
In this paper, we propose and analyze the following modification to the algorithm (3):
(4) 
We have introduced a new hyperparameter, . When , one obtains the heavyball method (2), and when , one obtains Nesterov’s usual acceleration, (3). In this work we show that it can be advantageous to use values of significantly larger than . Instead of using the gradient at an estimated point one step ahead, this is using the gradient at an estimated point several steps ahead. As this is an extension/strengthening of Nesterov’s “acceleration” idea, we refer to as “superacceleration”.
We will first examine the benefit of superacceleration in the ‘toy’ setting for convex optimization: a parabolic loss function in a oneparameter system, . The detailed analysis is in Section 2, but here we provide an summary of the main lessons. The algorithm (4) moves the parameters exponentially toward the minimum (possibly with oscillatory factors). The decay constant in the exponent provides the time scale for minimization. In Figure 1 we show the value of the superaccelearation parameter that minimizes this time scale. The optimal value of depends on the learning rate in the combination . Except for extremely small , the optimal is generally significantly larger than . This provides basic support for our proposal that superaccelerating Nesterov momentum (using , i.e., looking ahead several estimated steps) may be advantageous in training by gradient descent. Of course, the generalization from a quadratic onedimensional problem to a generic machinelearning problem is highly nontrivial. We therefore tested the effect of for twodimensional and dimensional minimization problems (Section 3), and for training neural networks for handwritten digit classification on the MNIST dataset (Section 4). In all cases, we find that extending beyond provides significant speedup in the optimization/training process.
In addition to momentum and acceleration defined above, many other variants and improvements of gradient descent are in use or under discussion in the machine learning community. Listings and discussions of the most common algorithms can be found in (Ruder, 2016; Mehta et al., 2019; Sun et al., 2019; Zhang, 2019). Overviews are also provided in Chapter 8 of Goodfellow et al. (2016), and Chapter 5 of Kochenderfer and Wheeler (2019). Specially prominent are methods that adaptively modify the learning rate , such as AdaGrad (Duchi et al., 2011), RMSProp (Tieleman and Hinton, 2012), Adam (Kingma and Ba, 2014) and Adadelta (Zeiler, 2012). Momentum and acceleration can be readily incorporated into these adaptive schemes (Dozat, 2016; Ruder, 2016; Goodfellow et al., 2016). Hence, superacceleration (arbitrary ) can also be defined. The analysis of superaccelerating one of the modern adaptive algorithms is significantly complicated. However, we provide a preliminary analysis for RMSProp (Section 5) which shows that using should benefit RMSProp as well.
This Article is organized as follows. In Section 2
, we provide a detailed treatment of the oneparameter parabolic problem. The exponential relaxation can be understood by approximating the algorithm by an ordinary differential equation (ODE) which describes the damped harmonic oscillator. Optimizing the algorithm is then seen to be equivalent to tuning parameters to “critical damping”; this allows us to predict the optimal superacceleration, as presented in Figure
1. In Section 3we explore our algorithm in two higherdimensional cases: a synthetic nonquadratic objective function and a 51dimensional quadratic objective function resulting from linear regression. This exploration shows that superacceleration is broadly advantageous, but also reveals possible interesting sideeffects such as nonlinear instabilities. In Section
4the algorithm is applied to the MNIST digit classification problem. Superacceleration is applied to examples using both nonstochastic and stochastic gradient descent. In Section
5 we show how superacceleration is readily incorporated into RMSProp, and the benefits of this combination. Finally, Section 6 provides context and discussion.2 Quadratic loss function for one parameter — damped harmonic oscillators and optimal superacceleration
In this section we analyze the effect of superacceleration on the minmization of a quadratic objective function of a single parameter . In the Introduction (Figure 1) we have provided an ‘executive summary’; here we consider the problem in some detail.
We are interested in gradient descent with superacceleration applied to minimizing the loss function of the form . Here is a single real variable, and is a positive constant. This is the simplest and most natural model potential for a convex function of a single variable. The algorithm for this special case is
(5) 
which can be rewritten as
(6) 
with a negative real number. Written in this form, the iteration equations are reminiscent of difference equaions that are obtained upon discretizing ordinary differential equations (ODEs), as commonly done in numerical solutions of ODEs (Golub and Ortega, 1992; Stoer and Bulirsch, 2002). It turns out that equations (6) can be obtained as approximations to a differential equation representing damped harmonic motion (Qian, 1999; Su et al., 2014; Kovachki and Stuart, 2019). In fact, there are multiple differential equations which, upon discretization, can lead to the scheme of Eq. (6), depending on whether the derivatives are approximated by a forward differnce rule, a backward difference rule, or a midpoint rule.
In Subsection 2.1, we explain how the representation in terms of a damped harmonic oscillator allows us to define a time scale for the optimization process and thus leads to the prediction of optimal values of , as presented in the Introduction. In Subsection 2.2, we will outline several possible ODEs and compare how their predictions match the behavior of the algorithm. Subsection 2.3 discusses the optimal values of the superacceleration predicted by this analysis.
2.1 Time scale for minimization — critical damping
We will argue below that the algorithm (6) minimizes the variable exponentially, with possible oscillations, i.e., the evolution of is welldescribed by the form , where is a continuous counterpart of the iteration count . When the dynamics is of this form, the quantity is the “timescale” for the minimization process, i.e., the inverse of the rate at which the algorithm approaches the minimum . The goal of optimizing the algorithm is thus to choose hyperparameters that will minimize the timescale . Other criteria, such as minimizing the time it takes to come within a small distance from , is often used in the more mathematical literature on optimization. We consider the smallness of the relaxation timescale to be a more practically relevant and physical measure of the performance of an optimization algorithm.
The algorithm (6) is approximated by a differential equation of the form
(7) 
with and . Here the function is a continuous version of . Expressions for and in terms of and are discussed in the next subsection. Here, we consider the consequence of this form of differential equation.
Eq. (7) describes damped harmonic oscillation, wellknown from introductory classical mechanics (Landau and Lifshitz, 1982; Greiner, 2003; Thornton and Marion, 2004). The solution has the form
(8) 
with . For small (‘underdamped’ regime), the solution is oscillatory with an exponentially decreasing envelop. In this regime the relaxation constant can be read off from the exponential factor to be , which decreases with .
For (‘overdamped’ regime), the oscillation frequency becomes imaginary (), so that the dynamics is nonoscillatory. Now, the dynamics consists of two exponentially decaying terms, with rate constants . There are different timescales corresponding to the initial and the longtime relaxation. One may interpret either of as . In either case, increases with in the overdamped regime. (When we plot in Figure 1, we use the initial relaxation timescale, because this seems more relevant to the efficiency of training a machine learning model.)
In Figure 2(a), we show some example trajectories of the algorithm (5) starting from , together with fits to the function , with fit parameters , and . The fits are seen to be excellent, demonstrating that the minimization process of algorithm (5) is welldescribed by exponential relaxation with or without oscillations, as predicted by the ODE approximation (7).
As the relaxation constant decreases in the underdamped regime and increases in the overdamped regime, the point of ‘critical’ damping ( or ) is where is smallest. It turns out that is a monotonically increasing function of (Subsection 2.2). Thus there is an optimal value, , at which is a minimum. This was the result shown and discussed in Figure 1. In Figure 1 we obtained the optimization timescale simply by fitting data from the algorithm to the form . Approximate analytical expressions for can be obtained from the differential equations, by finding the parameter values corresponding to critical damping (Subsection 2.3).
2.2 Damped harmonic oscillators: ODEs approximating Eq. (6)
In the numerical solution of ordinary differential equations describing evolution, say Newton’s equation for the motion of a particle, it is common to distretize time into small steps and represent derivatives as finite differences (Golub and Ortega, 1992; Stoer and Bulirsch, 2002). Here, we reverse the approximation procedure: our starting point, the algorithm (5) or (6
), is a finite difference equation which we relate approximately to an ODE. The procedure is not unique and several slightly different ODEs can serve this purpose. Fortunately, they all have the same form (
7).In Eq. (6), we interpret as the discretized values of a position variable . The ’s are then integer values of . The second line of Eq. (6) then suggests that should be interpreted as the discretized values of velocity, , because the left side of this equation, , is the firstorder forward difference formula for the derivative of , with step size unity. In the first line of Eq. (6), we could identify the difference as the acceleration . We thus have the ODE
(9) 
An inconsistency in this procedure is that both and have been separately identified as . However, the error due to this inconsistency is of the same order as the difference formulas. If one insists on identifying as also in the first line of Eq. (6), one obtains by using the backward difference for the derivative of :
(10) 
Another approach would be to identify as , so that the difference formula used for is now a secondorder approximation. This leads to
(11) 
Since the difference formula for is still firstorder, the relation between this last equation (11) and the original discrete algorithm is still a firstorder approximation.
In Figure 2(b) we compare the oscillatory minimization dynamics of the algorithm (5) or (6) to numerically exact solutions of the three ODEs. It is clear that Eq. (11) performs best as an approximation.
Although it is not of higher order than the other two ODEs, the empirical observation is that Eq. (11) reproduces features of the discrete algorithm (such as the value of ) more faithfully than the other two ODEs, (9) and (10). This is presumably because Eq. (11) incorporates a secondorder discretization approximation.
Of course, there may be additional ODEs which approximate the same finite difference equations, but we confine ourselves to these three. Our purpose is not to explore the space of possible ODE approximations, but rather to show that all reasonable approximations are of the form (7), i.e., describe damped harmonic motion.
2.3 Optimal parameter
The three ODE’s, (9), (10), and (11), each have the form (7) discussed in Subsection 2.1. Each of these equations describe damped harmonic oscillations. In each case, and can be expressed in terms of and . The case of critical damping is obtained from the condition . Solving this condition for gives estimates for the optimal superacceleration:
(12) 
(13) 
and
(14) 
We can obtain the optimal superacceleration parameter directly from the discrete algorithm, by fitting a function of the form to the vs data (identfying with time), and then finding the value of for which the timescale is minimal. In Figure 2(c) we compare the numerically determined with the approximations derived above. The overall shape is described qualitatively by any of the ODE approximations, but the midpoint approximation (11) performs best quantitatively.
2.4 Very small
The value of calculated numerically for the algorithm becomes illdefined when is very small, namely, when (Figure 1). The analytic predictions, (12), (13), (14), give negative answers for very small . This indicates that, in very shallow valleys or when a very small learning parameter is used, the benefit of superacceleration disappears.
At such small values of , any positive is ‘above’ the curve and hence the algorithm is automatically in the overdamped regime, so that the minimization dynamics is nonoscillatory. Moreover, since appears in the algorithm only through the quantity , when the superacceleration parameter does not affect the algorithm at all. We will later see an example of this effect for a highdimensional system (Subsection 3.2).
2.5 Summary of oneparameter analysis
In this Section, we provided details and derivations related to the message presented in the Introduction (Figure 1): that superaccelerating momentumbased gradient descent is beneficial for minimization in the onedimensional parabolic case.
Our analysis shows that either Nesterovaccelerated () or superaccelerated () gradient descent relaxes to the minimum exponentially. We can therefore associate a timescale to the relaxation dynamics. Optimizing our algorithm means minimizing . With continuoustime approximations, these algorithms are welldescribed as damped harmonic oscillators. Since decreases with in the underdamped regime and increases with in the overdamped regime, the corresponding to critical damping provides the minimal relaxation time. Thus there is a welldefined optimal superacceleration . This was presented in the Introduction (Figure 1).
The result of the analysis is that is preferable for all values above a minimum. (Numerically, we find this value to be .) This indicates that superaccelerating the algorithm will give no benefit for minimization along very shallow valleys. We will see examples of this effect later.
3 Higherdimensional examples
We now turn to the performance of superacceleration in the minimization of multivariable loss functions. In our analysis of the onevariable case with quadratic objective function, we found a welldefined optimal value of the superacceleration parameter . It is, of course, impossible to generalize this prediction to an arbitrary multidimensional problem. Even if the dependences of on each direction could be approximated as quadratic (), the constants corresponding to different directions can be wildly different, making it impossible to identify an optimal superacceleration . Nevertheless, we expect that when the learning rate is not forced to be too small, superacceleration would be beneficial for many of the directions of space, so that it is likely to be beneficial for the overall optimization.
In this Section, we test this idea on two synthetic multidimensional objective functions. Subsection 3.1 treats a nonquadratic function of two parameters. In Subsection 3.2 we treat quadratic dimensional problems. In both cases, we show that using provides benefits; at the same time we learn of some possible pitfalls.
3.1 Minimizing 2D function(s)
In Figure 3 we show the effect of superacceleration on gradient descent minimization for the objective function
(15) 
of two variables, . The selected function is nonconvex and its gradient is nonlinear. It is designed to have a valley around which the ‘hills’ have some spatial nonisotropic structure, as seen in the contour plots in Figure 3. The last term ensures that the minimization procedure does not escape to infinity in any direction, and could be thought of as mimicking a regularization term for machine learning. In this Subsection, we consider the effects of adding momentum and (super)acceleration to gradient descent for this potential, starting at the initial point and using a learning rate . We will see that a superacceleration value around provides the best minimization. Of course, the exact performance and the benefits of superacceleration will depend on the initial point and the learning rate, even for this particular objective function. However, we expect that there are broadly applicable lessons to be learned from analyzing this particular case.
In Figure 3 we show the trajectories of gradient descent with(out) momentum and with various values of superacceleration . The first 180 steps are shown for each minimization algorithm. In the top row (ac), trajectories are shown on the  plane. To avoid excessive clutter, each panel shows the trajectories of only two algorithms, one with large empty red circles and one with small blue stars. Successive positions are joined with straight lines for better visualization of the trajectories. In the bottom row, trajectories are characterized using the Euclidean distance to the position of the minimum.
Panel 3(a) shows pure gradient descent without momentum, and also with accelerationless momentum, i.e., the heavy ball algorithm of Eq. (2) corresponding to . Both these algorithms fall short of reaching the minimum in 180 steps. The striking effect of adding momentum is that the trajectory performs wide oscillations and overshoots. It is often stated that adding momentum dampens oscillation (Ruder, 2016; Goodfellow et al., 2016), in apparent contradiction to the observation in panel 3(a). The situation in which momentum has a dampening effect is when the learning rate is so large that pure gradient descent moves from one side of a valley to another in each step, which appears to be the example(s) considered in (Ruder, 2016; Goodfellow et al., 2016). If
is sufficiently small, the situation observed in panel (a) is probably more typical: adding momentum causes strong oscillations, but can speed up the minimization because the step sizes tend to be much larger.
In 3(b) we see that adding acceleration, , has a dampening effect on the trajectory. Also, superacceleration () works significantly better than usual Nesterov acceleration (). This improvement is also clearly seen in panel (d), where the Eucledean distance to the minimum is plotted. The distance falls much faster with superacceleration . (In this particular case, adding momentum or Nesterov acceleration does not improve the minimization compared to pure gradient descent, but this is not generic — typically momentum by itself will improve the algorithm.)
In 3(c) and 3(e) we show the effect of further increasing . For this case, the optimal values appear to be around or . As is further increased, the algorithm tends to get trapped in regions other than the minimum. For example, 3(c) shows that the algorithm shoots commendably fast into the valley, but after that, gets stuck in a region which is not at the minimum. We find that the iteration oscillates between two nearby points. In 3(e), inset, we show the algorithm approaching a similar situation. It appears that the nonlinear map defined by the superaccelearted iteration has attractors composed of points which are not at the minimum of the loss function. The difference equations (map) describing the algorithm is nonlinear due to the nonlinearity of the potential (15). In the algorithm (4), the parameter appears inside the (gradient of the) loss function; hence it is not surprising that nonlinear effects are enhanced when is increased. In the present case, when is increased too far, nonlinear effects are enhanced to the point of developing attractors at locations which are not at the minimum.
Comparing in 3(e) the curves for or , one notices that curve has a dip indicating that the algorithm overshoots the minimum and then returns along the valley (as with smaller values of ), but does not have such an overshooting behavior. This can be understood from our previous oneparameter analysis applied to dynamics within the valley: presumably, corresponds to underdamping, and is closer to criticial damping. This interpretation is meaningful only to the extent that the valley may be approximated as a onedimensional quadratic structure.
To summarize, we have found that the advantages of superacceleration, i.e., increasing beyond , extends beyond the toy oneparameter problem analyzed in Section 2, and holds for nonquadratic potentials such that the gradient is a nonlinear function of parameters. However, in such a potential, for sufficiently large , the nonlinearity may give rise to unwanted trapping effects.
3.2 Linear regression
We now consider applying gradient descent with superacceleration to a manyvariable regression problem.
Consider realvalued observations which are functions of realvalued features, , with , , , …, . Given such sets of feature values and corresponding observations, the regression problem fits the data using a linear model
The parameter space ( space) is dimensional. The loss function to be minimized in linear regression is
Of course, linear regression can be solved by matrix methods which may be more efficient than gradient descent. Nevertheless, this is a standard pedagogical setup for introducing gradient descent. The loss function for linear regression is known to be convex. Here we apply our algorithm, Eq. (4), to this toy problem, to gain intuition about superacceleration in manydimensional situations.
In Figure 4 we present minimization algorithms for , i.e., a dimensional problem. We used data points, i.e., sets of values. The were chosen uniformly in and the values were . All are set to the same value at the start of each iteration.
Of course, the details are modified when other datasets, initialization, or are used. We have explored a number of variants, and observed that the aspects we emphasize below are broadly applicable.
A curious aspect visible in Figure 4 (inset) is that, at late times in the iteration, the iterations with momentum all follow the same path, irrespective of superacceleration . Asymptotically, there appears to be a universal, independent . In contrast, pure gradient descent without momentum trains slower (the loss function decreases slower) than any of the momentumbased procedures. The asymptotic curve is reached differently for different , as seen in the main panel. As is increased, the initial minimization is more effective. In fact, in this particular case, for initial minimization the case is slower even than pure gd, while the case (Nesterov acceleration) is similar to pure gd. The momentumbased methods perform better than pure gd only when superacceleration, , is used.
These dynamical behaviors can be understood by decomposing the loss function into the “eigendirections” of its hessian at the minimum. Since the loss function is quadratic, it can then be regarded as a sum of decoupled oneparameter quadratic objective functions, , where the ’s are linear combinations of the parameters. The qualitative features of Figure 4 can be understood as a sum of the decoupled cases of the oneparameter case treated in Section 2. The steeper directions (with larger ) are minimized faster by the algorithm. For , we expect superacceleration to be beneficial, this explains why larger provides faster initial minimization. At late times, the directions with larger have already been minimized, so that only the smallest value is relevant. The dynamics is then in the regime discussed in Subsection 2.4. Any nonzero puts us in the overdamped region (above the line of Figure 1(b)); hence we have the slow nonoscillatory decay seen at late times. The timescale of this decay is nearly independent of , as appears in the theory only in the combination , which becomes effectively independent when .
Given that the longtime behavior is identical for a range of , one might object that superacceleration provides no benefit in this case. However, for practical machinelearning, the initial rapid minimization is often the most relevant process, and the latetime fine minimization might not even be desirable (e.g., might lead to overfitting). For the initial minimization, Figure 4 shows superacceleration to be very beneficial.
3.3 Lessons learned
From our studies with the two synthetic objective functions in this section, we have seen that superaccelerating Nesterov momentum in general speeds up minimization, but that there could be pitfalls. In Subsection 3.1, we found that at large nonlinearities might cause convergence to spurious fixed points not located at the minimum. In Subsection 3.2 we saw a case where latetime minimization happens in an extremely shallow direction, for which superacceleration is no better (but no worse) than Nesterov acceleration or even pure momentum.
4 MNIST classification with neural networks
Having explored synthetic potentials and datasets in the previous section, we now turn to applying superacceleration to a standard machine learning benchmark problem: handwritten digit recognition using neural networks, using the MNIST dataset (LeCun, 1998; Nielsen, 2015). In Figure 5 we show results obtained with gradient descent without stochasticity, i.e., instead of using minibatches, all of the training data is used for computing the gradient. Thus each step is an epoch.
As usual, the training data consists of 60k images and the test data consists of 10k images. Figure 5 shows the accuracy and the loss function computed on the test data, as training progresses. The input layer has nodes, each input being a pixel value. The output layer has 10 nodes, representing desired classification results for the digits from 0 to 9. The loss function used is the Euclidean distance of the output from the onehot vector corresponding to the digit label, averaged over the images.
To display the robustness of our results, we show the training process for two different network architectures. The panels on the top row, 5(a,b), are for a network with a single hidden layer with 30 nodes, i.e., a –– architecture. The lower row shows data for training a ––– network, with two hidden layers.
We show both the classification accuracy and loss function on the test data. (The curves for the training data are very similar as there is little or no overfitting at these training stages.) In both cases, each training algorithm is started with the same values of parameters. The weights and biases are each drawn randomly from a gaussian distribution, centred at zero and with widths either unity (biases) or normalized by the number of nodes in that layer (weights), as suggested in Chapter 3 of
Nielsen (2015). We run the minimization algorithms with learning rate in each case.Although some details are different for the two networks (top and bottom row of Figure 5), the outstanding feature in both cases is that superacceleration provides clear improvement in the training process. In both networks, provides superior training compared to smaller values of . We have also tested . In the synthetic examples of Section 3, such large values generally led to instabilities or overflow errors. However, in the MNIST classification tasks, superacceleration apears to function well, and even outperforms in the twolayer case. According to our onedimensional analysis in Section 2, would be beyond the optimal value of , but in the MNIST case large values apparently continue to be performant.
Similar to the multidimensional quadratic case studied in Subsection 3.2, we note that, late in the training process, the curves with momentum appear to merge together independent of . In our data up to 400 steps, the curves for have already merged. In analogy to Figure 4, one surmises that the other curves will merge later on during training. This suggests that the loss landscape for MNIST classification has some similarities to that of the regression problem, i.e., to a multidimensional quadratic function with a range of curvatures in different directions.
In Figure 6, we show evaluations of superacceleration using
stochastic gradient descent, which is more common in practical machine learning. The
training set is now stochastically divided into minibatches, and the minibatches are used to
calculate approximations to the gradient. The overall features are similar to the nonstochastic
case of Figure 5:
Larger allows rapid optimization;
hence superacceleration is beneficial.
In the longtime limit, the accuracy or loss tends
to become independent of .
5 Combining with adaptive methods (RMSProp)
The learning rate is a hyperparameter that obviously affects the minimization process strongly. Therefore, it is tempting to modify the effective learning rate based on data acquired during the process. In recent years, many such adaptive algorithms have appeared, and at this time these algorithms have arguably become more common than gradient descent with constant or predetermined . In particular, RMSProp (Tieleman and Hinton, 2012) and Adam (Kingma and Ba, 2014) are often mentioned as default minimization algorithms.
When RMSProp is combined with momentum, it is almost identical to Adam (Ruder, 2016; Goodfellow et al., 2016). Nesterovaccelerating these methods, by evaluating the gradient at the estimated onestepforward position, have also been discussed (Dozat, 2016; Ruder, 2016). Clearly, if the algorithm can be Nesterov accelerated, it can also be superaccelerated, in the sense introduced in the present work. A thorough understanding of the effects and potential benefits of combining adaptive learning rates with superacceleration would require a very extensive study; here we present some basic calculations.
The RMSProp algorithm with momentum and superacceleration is
(16) 
We have taken the RMSProp+momentum algorithm as described in (Ruder, 2016; Goodfellow et al., 2016), and added superacceleration, i.e., the gradient is now evaluated not at the current values of the parameters, but at an estimated forward point . The values and are reasonably standard; we confine ourselves to these values. The ‘division’ in the second line is componentwise division, which is not standard mathematical notation but standard in this field. The intuition behind this algorithm is that directions in space which have seen little change (smaller gradient magnitudes) in recent iterations are favored at the expense of those directions that have enjoyed rapid changes (larger gradient magnitudes). This is supposed to help escape from saddle points and minimize rapidly along shallow directions.
The Adam and Nadam algorithms are based on almost the same ideas (Kingma and Ba, 2014; Dozat, 2016; Ruder, 2016; Goodfellow et al., 2016). The main differences are (1) that the adaptive factor appears in the third equation for rather than in the second equation for ; (2) a bias correction that is important for initial iterations, which we believe does not play much of a role in determining the timescale. We do not currently have a good theoretical understanding of the differences between RMSProp (with momentum) and Adam. It is reasonable to expect that the basic intuitions obtained below for combining superaccelearation with RMSProp should also hold for combining superacceleration with other adaptive algorithms.
We provide below some results for the algorithm (16) for the oneparameter case with a parabolic loss function (Subsection 5.1) and with the twoparameter loss function, Eq. (15), used previously (Subsection 5.2). We reiterate that this should be a regarded as an initial overview and a full analysis exploring various parameter regimes would involve a substantial study in its own right.
5.1 Oneparameter parabolic function
As we did in Section 2 for the nonadaptive superacceleration algorithm, we now examine RMSProp with superacceleration for the simplest model: a loss function of the form , where is a single real variable, and is a positive constant. The algorithm (16) for this special case becomes
(17) 
with the understanding that , and now have singlecomponent real values instead of having many components. (If they are regarded as vectors in , then .) If Eq. (17) is regarded as a map of the variable , then it is clearly a heavily nonlinear map. As in Section 2, we can find ordinary differential equations approximating this discrete algorithm; these ODE’s are nonlinear and therefore not as easy to draw lessons from as was the case with the damped harmonic oscillators. Instead, we present in Figure 7 results from direct analysis of numerical runs of the algorithm (17).
Unlike the nonadaptive case in Section 2, in Eq. (17) the constant , parameterizing the steepness of the parabolic potential, cannot be absorbed into a factor . It is thus an independent parameter, meaning that an exploration of values is necessary for a full understanding of the algorithm. For the purpose of the present preliminary exploration, we restrict to .
Although the damped harmonic oscillator (7) is no longer an excellent approximation, we find that the algorithm has the same overall qualitative features as in Section 2: for small we obtain oscillatory relaxation to the minimum , and the relaxation time scale decreases with increasing , while above a ‘critical’ the relaxation is nonoscillatory and increases with further increase of . Thus there is an optimal value of which minimizes . Despite the qualitative similarity, a function of the form does not provide quantitative fits to the vs data (identfying with time as before). Therefore, we simply fit the function to , in order to extract , effectively averaging over the oscillations, Figure 7(a). This procedure is not theoretically rigorous but suffices to extract a relaxation time scale. As before, plotting against provides us with the optimal superacceleration value . Figure 7(b) shows as a function of for .
From Figure 7(b) we conclude that a superacceleration parameter considerably larger than is beneficial, as long as the learning parameter is not very small. This is very similar to the nonadaptive case.
5.2 Twodimensional function
We have performed some exploratory calculations of RMSProp with superacceleration, algorithm (16), for twodimensional test functions. Overall, the situation is less clear than the nonadaptive case, and the algorithm seems more likely to show classic behaviors associated with nonlinearlities, such as attractions to regions other than the minimum. However, the general trend is still that superacceleration with is beneficial.
In Figure 8, we show the minimization process for the twoparameter function, Eq. 15, explored in Subsection 3.1, for and the same starting point as in Figure 3. As in the nonadaptive case, the minimization becomes faster as is increased beyond , but when is too large, nonlinearities cause the algorithm to converge elsewhere instead of at the minimum. These results suggest that it might be generally advantageous to start with a large for rapid initial minimization, and then at later times reduce so as to avoid nonlinearities.
6 Discussion, Context and Recommendations
Currently, there is no consensus on the ‘best’ variant of gradient descent for use in machine learning, and minimization methods continue to be an active research topic. The most basic improvements of gradient descent are momentum and Nesterov acceleration. There is a large body of current research either analyzing or suggesting modifications to (nonadaptive) momentumbased methods (Wibisono and Wilson, 2015; Wilson et al., 2016; Wibisono et al., 2016; Yuan et al., 2016; Jin et al., 2017; Lucas et al., 2018; Ma and Yarats, 2018; Cyrus et al., 2018; Srinivasan et al., 2018; Kovachki and Stuart, 2019; Chen and Kyrillidis, 2019; Gitman et al., 2019).
The present work can be considered a contribution of this type. We have proposed a simple twist on the idea of accelerated momentum — extending the acceleration idea of ‘looking ahead’ to beyond a single estimated step. We have also presented an analysis of the momentum and Nesterov acceleration algorithms in terms of continuoustime ordinary differential equations (ODEs), and exploited the interpretation of these ODE’s as damped harmonic motion to predict the optimal superacceleration parameter, . The idea that moving toward the critical point is advantageous has appeared before in the literature (Qian, 1999). We have used this idea quantitatively to predict the optimal .
Although adaptive methods (RMSProp, Adam, etc.) are generally considered to be better performant, the relative merits of adaptive versus nonadaptive methods are not yet fully understood (Wilson et al., 2017; Reddi et al., 2019; Choi et al., 2019). Naturally, there is also a significant amount of current research on analyzing and improving adaptive methods (Basu et al., 2018; Zhou et al., 2018; Ma and Yarats, 2018; Barakat and Bianchi, 2018; Luo et al., 2019; Tong et al., 2019; Wu et al., 2019; Zhang et al., 2019). In the present work, we have presented some preliminary results on how our superacceleration ideas can also be incorporated into adaptive methods such as RMSProp or Adam/Nadam.
We have described a variety of cases where superacceleration () speeds up optimization. Our exploration has also revealed two ways in which superacceleration may fail to provide an advantage. First, in a very shallow part of the landscape, the algorithm may be in the regime of ulrasmall , in the language of Section 2. In such a case, the algorithm becomes nearly independent. Another possible pitfall is due to large nonlinearities, which can be amplified by to the extent that the algorithm can get trapped in a nonextremal region. Examples were seen in Subsections 3.1 and 5.2. For the severely nonquadratic synthetic function of Eq. (15), such instabilities or attractors appeared for . Note however that we did not see this type of phenomenon in our experimentation with neural networks (MNIST classification, Section 4), even for much larger values, which suggests that these problems (nonlinear instabilities) might not be serious in practical applications.
In view of these considerations, we can provide the following recommendation for using gradient descent with superacceleration: start the algorithm with a relatively large , say around 5, and then gradually lower to around so as to avoid spurious nonlinear phenomena.
This work opens up some questions which require further investigation. (1) First, a thorough analysis and exploration of superacceleration in combination with the standard adaptive methods — RMSProp, Adam and Adadelta — remains to be performed. In Section 5 we have provided a very preliminary look. Our initial results are promising and suggest that superaccleration might improve these algorithms significantly. In particular, it is possible that in combination with Adadelta, superacceleration might not suffer from the small problem, as the effective learning rate might be increased when the slope (effective ) becomes too small. (2) Second, we have focused on nonstochastic minimization, except for a cursory exploration in Section 4 in the context of MNIST digit classification. A thorough analysis of superacceleration for stochastic gradient descent remains another task for future work.
References
 Barakat and Bianchi (2018) Anas Barakat and Pascal Bianchi. Convergence of the adam algorithm from a dynamical system viewpoint. arXiv preprint arXiv:1810.02263, 2018.
 Basu et al. (2018) Amitabh Basu, Soham De, Anirbit Mukherjee, and Enayat Ullah. Convergence guarantees for rmsprop and adam in nonconvex optimization and their comparison to nesterov acceleration on autoencoders. arXiv preprint arXiv:1807.06766, 2018.
 Chen and Kyrillidis (2019) John Chen and Anastasios Kyrillidis. Decaying momentum helps neural network training. arXiv preprint arXiv:1910.04952, 2019.
 Choi et al. (2019) Dami Choi, Christopher J Shallue, Zachary Nado, Jaehoon Lee, Chris J Maddison, and George E Dahl. On empirical comparisons of optimizers for deep learning. arXiv preprint arXiv:1910.05446, 2019.
 Cyrus et al. (2018) Saman Cyrus, Bin Hu, Bryan Van Scoy, and Laurent Lessard. A robust accelerated optimization algorithm for strongly convex functions. In 2018 Annual American Control Conference (ACC), pages 1376–1381. IEEE, 2018.
 Dozat (2016) Timothy Dozat. Incorporating nesterov momentum into adam. In ICLR Workshop, 2016.
 Duchi et al. (2011) John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of Machine Learning Research, 12(Jul):2121–2159, 2011.
 Gitman et al. (2019) Igor Gitman, Hunter Lang, Pengchuan Zhang, and Lin Xiao. Understanding the role of momentum in stochastic gradient methods. In Advances in Neural Information Processing Systems, pages 9630–9640, 2019.
 Golub and Ortega (1992) G.H. Golub and J.M. Ortega. Scientific Computing and Differential Equations: An Introduction to Numerical Methods. Academic Press, 1992.
 Goodfellow et al. (2016) Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep learning. MIT press, 2016.
 Greiner (2003) W. Greiner. Classical Mechanics: Point Particles and Relativity. Springer New York, 2003.
 Jin et al. (2017) Chi Jin, Praneeth Netrapalli, and Michael I Jordan. Accelerated gradient descent escapes saddle points faster than gradient descent. arXiv preprint arXiv:1711.10456, 2017.
 Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
 Kochenderfer and Wheeler (2019) Mykel J Kochenderfer and Tim A Wheeler. Algorithms for optimization. 2019.
 Kovachki and Stuart (2019) Nikola B Kovachki and Andrew M Stuart. Analysis of momentum methods. arXiv preprint arXiv:1906.04285, 2019.
 Landau and Lifshitz (1982) L.D. Landau and E.M. Lifshitz. Mechanics. Elsevier Science, 1982.
 LeCun et al. (1998) Y. LeCun, L. Bottou, G. Orr, and K. Muller. Efficient backprop. In G. Orr and K. Muller, editors, Neural Networks: Tricks of the trade. Springer, 1998.

LeCun (1998)
Yann LeCun.
The mnist database of handwritten digits.
http://yann. lecun. com/exdb/mnist/, 1998.  LeCun et al. (2015) Yann LeCun, Yoshua Bengio, and Geoffrey Hinton. Deep learning. Nature, 521(7553):436–444, 2015.
 Lucas et al. (2018) James Lucas, Shengyang Sun, Richard Zemel, and Roger Grosse. Aggregated momentum: Stability through passive damping. arXiv preprint arXiv:1804.00325, 2018.
 Luo et al. (2019) Liangchen Luo, Yuanhao Xiong, Yan Liu, and Xu Sun. Adaptive gradient methods with dynamic bound of learning rate. arXiv preprint arXiv:1902.09843, 2019.
 Ma and Yarats (2018) Jerry Ma and Denis Yarats. Quasihyperbolic momentum and adam for deep learning. arXiv preprint arXiv:1810.06801, 2018.

Mehta et al. (2019)
Pankaj Mehta, Marin Bukov, ChingHao Wang, Alexandre GR Day, Clint Richardson,
Charles K Fisher, and David J Schwab.
A highbias, lowvariance introduction to machine learning for physicists.
Physics reports, 2019.  Nesterov (1983) Yurii Nesterov. A method for unconstrained convex minimization problem with the rate of convergence . Soviet Mathematics Doklady, 269:543–547, 1983.
 Nielsen (2015) M.A. Nielsen. Neural Networks and Deep Learning. Determination Press, 2015. URL https://books.google.nl/books?id=STDBswEACAAJ.
 Polyak (1964) Boris T Polyak. Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics, 4:1–17, 1964.
 Qian (1999) Ning Qian. On the momentum term in gradient descent learning algorithms. Neural networks, 12(1):145–151, 1999.
 Reddi et al. (2019) Sashank J Reddi, Satyen Kale, and Sanjiv Kumar. On the convergence of adam and beyond. arXiv preprint arXiv:1904.09237, 2019.
 Ruder (2016) Sebastian Ruder. An overview of gradient descent optimization algorithms. arXiv preprint arXiv:1609.04747, 2016.
 Rumelhart et al. (1986) D Rumelhart, G Hinton, and R Williams. Learning representations by backpropagating error. Nature, 323(9):533–536, 1986.

Srinivasan et al. (2018)
Vishwak Srinivasan, Adepu Ravi Sankar, and Vineeth N Balasubramanian.
Adine: An adaptive momentum method for stochastic gradient descent.
Proceedings of the ACM India Joint International Conference on Data Science and Management of Data  CoDSCOMAD ’18
, 2018.  Stoer and Bulirsch (2002) J. Stoer and R. Bulirsch. Introduction to Numerical Analysis. Springer New York, 2002.
 Su et al. (2014) Weijie Su, Stephen Boyd, and Emmanuel Candes. A differential equation for modeling nesterov’s accelerated gradient method: Theory and insights. In Advances in Neural Information Processing Systems, pages 2510–2518, 2014.
 Sun et al. (2019) Shiliang Sun, Zehui Cao, Han Zhu, and Jing Zhao. A survey of optimization methods from a machine learning perspective. IEEE Transactions on Cybernetics, pages 1–14, 2019.
 Sutskever et al. (2013) 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.
 Thornton and Marion (2004) S.T. Thornton and J.B. Marion. Classical dynamics of particles and systems. Thompson Brooks/Cole, 5th edition, 2004.
 Tieleman and Hinton (2012) Tijmen Tieleman and Geoffrey Hinton. Lecture 6.5rmsprop: Divide the gradient by a running average of its recent magnitude. COURSERA: Neural networks for machine learning, 4(2):26–31, 2012.
 Tong et al. (2019) Qianqian Tong, Guannan Liang, and Jinbo Bi. Calibrating the learning rate for adaptive gradient methods to improve generalization performance. arXiv preprint arXiv:1908.00700, 2019.
 Wibisono and Wilson (2015) Andre Wibisono and Ashia C Wilson. On accelerated methods in optimization. arXiv preprint arXiv:1509.03616, 2015.
 Wibisono et al. (2016) Andre Wibisono, Ashia C Wilson, and Michael I Jordan. A variational perspective on accelerated methods in optimization. proceedings of the National Academy of Sciences, 113(47):E7351–E7358, 2016.
 Wilson et al. (2016) Ashia C Wilson, Benjamin Recht, and Michael I Jordan. A lyapunov analysis of momentum methods in optimization. arXiv preprint arXiv:1611.02635, 2016.
 Wilson et al. (2017) Ashia C Wilson, Rebecca Roelofs, Mitchell Stern, Nati Srebro, and Benjamin Recht. The marginal value of adaptive gradient methods in machine learning. In Advances in Neural Information Processing Systems, pages 4148–4158, 2017.
 Wu et al. (2019) Xiaoxia Wu, Simon S Du, and Rachel Ward. Global convergence of adaptive gradient methods for an overparameterized neural network. arXiv preprint arXiv:1902.07111, 2019.
 Yuan et al. (2016) Kun Yuan, Bicheng Ying, and Ali H. Sayed. On the influence of momentum acceleration on online learning. Journal of Machine Learning Research, 17(192):1–66, 2016.
 Zeiler (2012) Matthew D Zeiler. Adadelta: an adaptive learning rate method. arXiv preprint arXiv:1212.5701, 2012.
 Zhang (2019) Jiawei Zhang. Gradient descent based optimization algorithms for deep learning models training. arXiv preprint arXiv:1903.03614, 2019.
 Zhang et al. (2019) Michael Zhang, James Lucas, Jimmy Ba, and Geoffrey E Hinton. Lookahead optimizer: steps forward, step back. In Advances in Neural Information Processing Systems, pages 9593–9604, 2019.
 Zhou et al. (2018) Dongruo Zhou, Yiqi Tang, Ziyan Yang, Yuan Cao, and Quanquan Gu. On the convergence of adaptive gradient methods for nonconvex optimization. arXiv preprint arXiv:1808.05671, 2018.