Title: QUASI-HYPERBOLIC MOMENTUM AND ADAM FOR DEEP LEARNING

Abstract: Momentum-based acceleration of stochastic gradient descent (SGD) is widely used in deep learning. We propose the quasi-hyperbolic momentum algorithm (QHM) as an extremely simple alteration of momentum SGD, averaging a plain SGD step with a momentum step. We describe numerous connections to and identities with other algorithms, and we characterize the set of two-state optimization algorithms that QHM can recover. Finally, we propose a QH variant of Adam called QHAdam, and we empirically demonstrate that our algorithms lead to significantly improved training in a variety of settings, including a new state-of-theart result on WMT16 EN-DE. We hope that these empirical results, combined with the conceptual and practical simplicity of QHM and QHAdam, will spur interest from both practitioners and researchers. Code is immediately available.

1. INTRODUCTION:

Stochastic gradient descent (SGD) serves as the optimizer of choice for many recent advances in deep learning across domains (Krizhevsky et al., 2012; He et al., 2016a; Gehring et al., 2017) . SGD for deep learning is typically augmented with either the "heavy ball" momentum technique of Polyak (1964) or the accelerated gradient of Nesterov (1983) . In the deterministic setting, these methods provably yield faster convergence in fairly general settings. In the stochastic setting, these methods lose many theoretical advantages. However, due to its implicit gradient averaging, momentum can confer the benefit of variance reduction, applying less noisy parameter updates than plain SGD. Recent work has explicitly shown the use of momentum as a variance reducer (Roux et al., 2018) .

Algorithms Starting with gradient variance reduction as an informal and speculative motivation, we introduce the quasi-hyperbolic momentum (QHM) optimization algorithm in Section 3. Put as simply as possible, QHM's update rule is a weighted average of momentum's and plain SGD's update rule. We later propose a similar variant of Adam (QHAdam) in Section 5.

Connecting the dots QHM is simple yet expressive. In Section 4, we connect QHM with plain SGD, momentum, Nesterov's accelerated gradient, PID control algorithms (Recht, 2018; An et al., 2018) , synthesized Nesterov variants (Lessard et al., 2016) , noise-robust momentum (Cyrus et al., 2018 ), Triple Momentum (Scoy et al., 2018) , and least-squares acceleration of SGD (Kidambi et al., 2018) . Such connections yield reciprocal benefits -these algorithms aid in analyzing QHM, and conversely QHM recovers many of these algorithms in a more efficient and conceptually simpler manner. We then characterize the set of optimization algorithms that QHM recovers.

Practical validation and considerations:

In Section 6, we empirically demonstrate that QHM and QHAdam provide superior optimization in a variety of deep learning settings. We provide both comprehensive parameter sweep analyses on smaller models and case studies on large real-world models. We demonstrate improvements on strong (sometimes state-of-the-art) models simply by swapping out the vanilla algorithms with the QH counterpart. Notably, taking the WMT16 EN-DE translation model of Ott et al. (2018) , we achieve a 40% improvement in stability, along with a new state-of-the-art result of 29.45 BLEU. We then offer some practical tips for QHM and QHAdam.

Miscellany We provide errata for Kingma & Ba (2015) , Recht (2018) , and Kidambi et al. (2018) . We also offer evidence that momentum often yields negligible improvement over plain SGD.

We emphasize QHM and QHAdam's efficiency and conceptual simplicity. QHM has no extra overhead vs. Nesterov's accelerated gradient, and QHAdam has very little overhead vs. Adam. Also, both algorithms are easily understood as an interpolation between two other well-known algorithms, so they are accessible to practitioners and can be tuned starting with existing practical intuitions. We believe that this contributes strongly to the algorithms' practical promise.

2. PRELIMINARIES:

We begin with notation and a brief review of stochastic gradient descent (SGD) and momentum.

Primitives In this paper, θ ∈ R p denotes a vector of model parameters. L(θ) : R p → R denotes a loss function to be minimized via θ. L(θ) : R p → R denotes an approximator of the loss function (e.g. over a minibatch). ∇L denotes the gradient of function L. Unless otherwise specified, all vector operations are element-wise. We use g, a, s, v, w ∈ R p as auxiliary buffers, and g is typically the "momentum buffer". θ, L(•), and all buffers are subscriptable by t, the optimization step.

Optimization algorithms We consider optimization algorithms that perform a sequence of steps (indexed by t), updating θ at each step towards minimizing L(θ). For brevity, we write algorithms as "update rules", which describe the algorithm's behavior during a single step t, rather than as full pseudocode. Update rules take this basic form (optionally with one or more auxiliary steps):

θ t+1 ← θ t -[. . .]

Plain SGD The SGD algorithm, parameterized by learning rate α ∈ R, uses the update rule:

θ t+1 ← θ t -α • ∇ Lt (θ t )

Momentum The momentum algorithm, parameterized by α ∈ R and β ∈ R, uses the update rule:

g t+1 ← β • g t + (1 -β) • ∇ Lt (θ t ) (1) θ t+1 ← θ t -α • g t+1

(2) where g is commonly called the "momentum buffer". Note that β = 0 recovers plain SGD.

The exponential discount factor β controls how slowly the momentum buffer is updated. In the stochastic setting, β also controls the variance of a normalized momentum buffer. A common rule of thumb for momentum is β = 0.9 (Ruder, 2016). 2 In contrast to common formulations of momentum (Polyak, 1964; Sutskever et al., 2013) , we normalize, or "dampen", the momentum buffer g by (1 -β) in (1). This serves both to remove dependence of the update step magnitude on β, and to allow the interpretation of g as a weighted average of past gradients (and thus a gradient estimator). Of course, this also shrinks the updates by a factor of 1 -β vs. common formulations; this is easily reversible with a corresponding increase to α.

3. ALGORITHM: QUASI-HYPERBOLIC MOMENTUM (QHM):

In this section, we propose and discuss the quasi-hyperbolic momentum (QHM) algorithm.

QHM update rule:

QHM, parameterized by α ∈ R, β ∈ R, and ν ∈ R, uses the update rule:

g t+1 ← β • g t + (1 -β) • ∇ Lt (θ t )

(3)

EQUATION

Section 7.1 provides a recommended rule of thumb (ν = 0.7 and β = 0.999).

Alternative PID setting:

In Appendix E, we briefly discuss another PID setting by An et al. (2018) and relate the resulting optimization algorithm to QHM. In short, the setting is degenerate as the P, I, and D terms are linearly dependent. Thus, QHM can recover the resulting PID control optimizer.

4.3. SYNTHESIZED NESTEROV VARIANTS (SNV):

Section 6 of Lessard et al. (2016) describes a "synthesized Nesterov variant" algorithm, which we call "SNV" for convenience. This algorithm is used to analyze and improve optimizer robustness under "relative deterministic noise" (i.e. multiplicative noise of the gradient).

Update rule SNV, parameterized by γ, β 1 , β 2 ∈ R, uses the update rule:5 

ξ t+1 ← ξ t -γ • ∇ Lt (θ t ) + β 1 (ξ t -ξ t-1 ) θ t+1 ← ξ t+1 + β 2 (ξ t+1 -ξ t )

Connection with QHM We fully relate QHM and SNV in Appendix C.4. To summarize, QHM and SNV recover each other. By extension, QHM recovers the Robust Momentum method, which is a specific parameterization of SNV (Cyrus et al., 2018) . Moreover, since Robust Momentum recovers the Triple Momentum of Scoy et al. (2018) , QHM also recovers Triple Momentum.

Efficiency SNV is costlier than QHM, requiring 2 auxiliary buffers and 5 scaled vector additions. Jain et al. (2017) and Kidambi et al. (2018) point out various failures of momentum and NAG in the setting of stochastic least squares optimization. This motivates their proposal of the AccSGD algorithm, which yields faster convergence over momentum and NAG in certain least-squares regression settings. Here, we discuss the formulation of Kidambi et al. (2018) .

4.4. ACCSGD:

Update rule AccSGD, parameterized by δ > 0, κ > 1, ξ ≤ √ κ, and < 1, uses the update rule:

wt+1 ← 2 ξ κ • wt + 1 - 2 ξ κ w t - κδ • ∇ Lt (θ t ) θ t+1 = w t+1 ← ξ κ + ξ w t -δ • ∇ Lt (θ t ) + κ κ + ξ • wt+1

Connection with QHM We fully relate QHM and AccSGD in Appendix C.5. To summarize, QHM recovers AccSGD. In the reverse direction, AccSGD does not recover QHM; specifically, we disprove the claim in Kidambi et al. (2018) that AccSGD recovers NAG. Since QHM recovers NAG, AccSGD cannot fully recover QHM.

Efficiency AccSGD, like QHM, requires 1 auxiliary buffer. Computationally, AccSGD is costlier, requiring 2 in-place scalar-vector multiplications and 4 scaled vector additions per update step.

4.5. DISCUSSION:



Theoretical convergence results:

We note that various convergence results follow simply via these connections. In the deterministic (full-batch) case, since QHM recovers Triple Momentum, QHM also recovers the global linear convergence rate of 1 -1/ √ κ for strongly convex, smooth loss functions. 6 For first-order methods, this is the fastest known global convergence rate for such functions. In the stochastic (minibatch) case, QHM's recovery of AccSGD gives QHM the same convergence results as in Kidambi et al. (2018) 's least-squares regression setting, of O( √ κ • log κ • log 1 ) iterations for -approximation of the minimal loss.

Unifying two-state optimization algorithms These connections demonstrate that many two-state optimization algorithms are functionally similar or equivalent to each other. However, they are often implemented inefficiently and their parameterizations can be inaccessible to practitioners. QHM yields a highly accessible and efficient version of these algorithms. (Polyak, 1964) subfamily better recovered by QHM with ν = 1 NAG (Nesterov, 1983) subfamily same recovered by QHM with ν = β PID (Recht, 2018) parent worse QHM's β restricts PID's k P /k D PID (An et al., 2018) bijective worse degenerate; either "PI" or "PD" SNV (Lessard et al., 2016) bijective worse used in handling multiplicative noise Robust M. (Cyrus et al., 2018) subfamily worse SNV w/ convergence guarantees Triple M. (Scoy et al., 2018) subfamily worse "fastest" for str. convex, smooth L(•) AccSGD (Kidambi et al., 2018) subfamily worse acceleration for least-squares SGD

In Appendix D, we characterize the set of two-state optimization algorithms recoverable by QHM.

Our hope here is to provide future work with a routine conversion to QHM so that they may leverage the accessibility and efficiency benefits, as well as the many connections to other algorithms.

Many-state optimization algorithms Going beyond a single momentum buffer, it is possible to recover many-state algorithms by linearly combining many momentum buffers (with different discount factors) in the update rule. However, we found in preliminary experiments that using multiple momentum buffers yields negligible value over using a single slow-decaying momentum buffer and setting an appropriate immediate discount -that is, using QHM with high β and appropriate ν.

We note that the Aggregated Momentum (AggMo) algorithm (Lucas et al., 2018) precisely performs this linear combination of multiple momentum buffers. While AggMo takes a simple average of the buffers, an extended variant of AggMo allows for other linear combinations. This extended AggMo can be viewed as a many-state generalization of two-state algorithms (including QHM), recovering them when two buffers are used. Appendix H provides a supplemental discussion and empirical comparison of QHM and AggMo, corroborating our preliminary experiments' findings.

5. ALGORITHM: QHADAM:

The Adam optimizer (Kingma & Ba, 2015) has enabled many compelling results in deep learning (Xu et al., 2015; Vaswani et al., 2017; Yu et al., 2018) . We propose to replace both of Adam's moment estimators with quasi-hyperbolic terms, and we name the resulting algorithm QHAdam.

QHAdam update rule:

QHAdam, parameterized by α, ≥ 0, β 1 , β 2 ∈ [0, 1), and ν 1 , ν 2 ∈ R, uses the update rule:

g t+1 ← β 1 • g t + (1 -β 1 ) • ∇ Lt (θ t ) g t+1 ← 1 -β t+1 1 -1 • g t+1 s t+1 ← β 2 • v t + (1 -β 2 )(∇ Lt (θ t )) 2 s t+1 ← 1 -β t+1 2 -1 • v t+1 θ t+1 ← θ t -α   (1 -ν 1 ) • ∇ Lt (θ t ) + ν 1 • g t+1 (1 -ν 2 )(∇ Lt (θ t )) 2 + ν 2 • s t+1 +  

Note that only the last expression differs from vanilla Adam. In fact, QHAdam recovers Adam when ν 1 = ν 2 = 1. Moreover, modulo bias correction, QHAdam recovers RMSProp (Hinton et al., 2012) when ν 1 = 0 and ν 2 = 1, and NAdam (Dozat, 2016) when ν 1 = β 1 and ν 2 = 1. We note that Adam has inspired many variants such as AMSGrad (Reddi et al., 2018) and AdamW (Loshchilov & Hutter, 2017) , which can be analogously modified. 

6. EXPERIMENTS:

We perform two categories of experiments: parameter sweeps and case studies. For brevity, all experimental settings are summarized in Table 2 and comprehensively detailed in Appendix I.

6.1. PARAMETER SWEEPS:

With parameter sweeps, we aim to comprehensively study the various parameterizations of the QH algorithms using relatively small models. We train for 90 epochs with size-64 minibatches. For QHM, we initialize α = 1 and decay it 10-fold every 30 epochs. The sweep grid for QHM (encapsulating various parameterizations of plain SGD, momentum, and NAG) is:

ν ∈ {0, 0.25, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.98, 0.99, 0.995, 0.998, 0.999, 0.9995, 1} β ∈ {0, 0.25, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.98, 0.99, 0.995, 0.998, 0.999, 0.9995}

For QHAdam, we fix α = 10 -3 , = 10 -8 , ν 2 = 1, and β 2 = 0.999, and sweep over ν 1 and β 1 .

"Default" ν and β values Motivated by the popular momentum/NAG "default" of β = 0.9, we select a QH "default" of ν = 0.7 and β = 0.999 based on preliminary experimentation on the MNIST dataset (LeCun, 1998) along with the intuitions from Appendix A. In the following figures, we show these defaults along with the globally optimal parameterizations.

Results Fig. 1 presents selected results of these sweep experiments (full results in Appendix J). Perhaps the most immediate observation is that the QH algorithms improve both training and validation metrics. Even the hardcoded default ν = 0.7 and β = 0.999 handily outperforms the optimal parameterization of NAG or Adam in all settings. In some settings, there remains a large gap between the QH and vanilla algorithms at the end of training. In other settings, the gap shrinks to smaller levels. However, even for these latter settings, the QH algorithm converges much faster, suggesting that a more aggressive learning rate schedule can significantly reduce training time.

What about plain SGD? We note that in most of these experiments, there is little difference between the performance of plain SGD and NAG (particularly when compared to QHM). Although not shown in the figures, there is also little difference between plain SGD and momentum. This indicates that the benefit of momentum and NAG (in the common, unnormalized formulations) comes in large part from the increase in effective step size. We thus suspect that much of the folk wisdom about momentum's benefits for SGD should instead be folk wisdom about using sensible learning rates. In contrast, QHM provides significant benefits without changing the effective step size.

6.2. CASE STUDIES:

With case studies, we apply the QH algorithms to diverse settings, with (currently or recently) stateof-the-art models. Our case studies cover image recognition, language modeling, reinforcement learning, and neural machine translation. Each case study features a baseline setting and a QH setting, which are identical modulo the optimizer used. Results are presented in Fig. 2 and Table 3 .

Image recognition (RN152-ImageNet-QHM) We train a ResNet152 model (He et al., 2016a) on the ILSVRC2012 dataset (Russakovsky et al., 2015) . The baseline setting is nearly identical to the size-256 minibatch baseline in Goyal et al. (2017) , using NAG with β = 0.9 and a decaying learning rate schedule. The QH setting swaps out NAG for QHM, with ν = 0.7 and β = 0.999.7 Running 3 seeds, QHM plainly trains much faster than NAG, and QHM converges to a marginally superior validation error as well.8 

Language modeling (FConvLM-WikiText103-QHM) Deep learning for NLP often features "spiky" gradient distributions (e.g. encountering rare words). We train a FConv language model (Dauphin et al., 2016) on the WikiText-103 dataset (Merity et al., 2016) . The baseline setting precisely follows the original paper, using NAG with β = 0.99. The QH setting swaps out NAG for QHM, with ν = 0.98 and β = 0.998. 7 We suspect that high β improves stability in the presense of spiky gradients, and QHM's ν allows the use of high β.

Running 10 seeds, QHM outperforms the NAG baseline on validation perplexity by half a point.

Reinforcement learning (TD3-MuJoCo-QHAdam) Reinforcement learning presents a challenging task for gradient-based optimization, since the objective L is not stationary. QH algorithms provide a natural way of upweighting the most recent gradient. Here, we apply the TD3 algorithm (Fujimoto et al., 2018) to various MuJoCo environments (Todorov et al., 2012) . The baseline precisely follows Fujimoto et al. (2018) 's setup, which uses Adam with β 1 = 0.9 and β 2 = 0.999.

The QH setting swaps out Adam for QHAdam, with ν 1 = 0.9 and other parameters identical.9 

Running 10 seeds, QHAdam yields improvements in average reward on four environments out of seven tested, and virtually ties on another.

Neural machine translation (TF-WMT16ENDE-QHAdam) Many state-of-the-art neural machine translation (NMT) models are fragile to train. As in language modeling, the gradient distribution is often "spiky"; thus, Adam training often fails to converge due to a very small number of large parameter updates. 10 Here, we empirically demonstrate that QHAdam improves both performance and robustness by using ν 2 to control the maximum per-step update. We train a large transformer model (Vaswani et al., 2017) on the WMT16 English-German dataset. The baseline setting precisely follows the state-of-the-art setup of Ott et al. (2018) , using β 1 = 0.9 and β 2 = 0.98 for Adam. The QH setting swaps out Adam for QHAdam, with ν 1 = 0.8, β 1 = 0.95, ν 2 = 0.7, and β 2 = 0.98.11 

Running 10 seeds, the Adam baseline explodes on 4 seeds. QHAdam is more robust, converging for all seeds. Ultimately, QHAdam yields a new state-of-the-art-result of 29.45 BLEU. Thus, we improve both the stability and performance of the state-of-the-art with a simple optimizer swap.



7.1. PRACTICAL SUGGESTIONS:

We offer some practical suggestions for deep learning practitioners, particularly those who default to momentum, NAG, or Adam with β = 0.9 as a rule of thumb:

• Consider using QHM or QHAdam, instead of momentum, NAG, or Adam.

• While QHM parameters should be tuned when feasible, a decent rule of thumb is to set ν = 0.7 and β = 0.999. QHAdam parameter selection is somewhat more situational, although as discussed in Section 5, ν 2 = 1 and β 2 unchanged is usually reasonable when replacing a stable Adam optimizer with QHAdam. • Be mindful of learning rate differences between (unnormalized) momentum/NAG and QHM. Convert learning rates from the former to the latter via multiplication by (1 -β) -1 . For example, momentum/NAG with α = 0.1 and β = 0.9 should be replaced by QHM with α = 1. This conversion is unnecessary for Adam, as it already normalizes all buffers.

7.2. FUTURE WORK:

This paper has only scratched the surface when it comes to empirical evaluation of QHM and QHAdam. Future work could apply the algorithms to other well-studied tasks and architectures, both to assess the extent of their performance gains in diverse domains, and to further develop insights into hyperparameter choice.

Effective hyperparameter autotuning methods can improve the practicality of any optimization algorithm. Thus, a useful direction for future work is to create an effective ν, β adapter, possibly based on techniques such as YellowFin (Zhang et al., 2017) or via continuous-time optimal control analysis, as in Li et al. (2017) . Moreover, learning rate adaptation techniques such as Hypergradient Descent (Baydin et al., 2018) can be applied to both QHM and QHAdam.

Future work could develop convergence results for QHAdam. Convergence results for QHM in a reasonably general stochastic setting would also be appealing, although we are not aware of compelling analogous results for momentum or NAG.

Finally, momentum has been studied in the distributed, asynchronous setting, with some noting that the delays in asynchronous SGD are, in some sense, akin to adding momentum (Mitliagkas et al., 2016) . As a result, the optimal momentum constant β shrinks as more asynchronous workers are added to optimization. It would be interesting to extend these results to QHM, especially to disentagle the implicit effects of asynchrony on ν and β.

7.3. CONCLUSION:

QHM and QHAdam are computationally cheap, intuitive to interpret, and simple to implement. They can serve as excellent replacements for momentum/NAG and Adam in a variety of settings.

In particular, they enable the use of high exponential discount factors (i.e. β) through the use of immediate discounting (i.e. ν). QHM recovers numerous other algorithms in an efficient and accessible manner. Parameter sweep experiments and case studies demonstrate that the QH algorithms can handily outpace their vanilla counterparts. We hope that practitioners and researchers will find these algorithms both practically useful and interesting as a subject of further study.

A DISCOUNTED SUM ESTIMATORS (ORIGINAL VARIANCE REDUCTION MOTIVATION FOR QHM):

We now provide an interpretation of the momentum buffer as a discounted sum estimator, seeking to motivate the QHM algorithm from a variance reduction perspective.

A.1 DISCOUNTED SUMS For a discount function δ : N ≥0 → R and a sequence of vectors x 0...t ∈ R p , we define a discounted sum DS δ (x 0...t ) as:

DS δ (x 0...t ) = t i=0 δ(i) • x t-i When t i=0 δ(i)

= 1 for all t ≥ 0, we call this a discounted sum average. When ∞ i=0 δ(i) = 1, we call this a discounted sum average (modulo initialization bias).

A.2 EXPONENTIAL DISCOUNTING AND EWMA:

For β ∈ (-1, 1), we define the exponential discount function δ EXP,β as:

δ EXP,β (i) = (1 -β)β i

and the exponentially weighted moving average EWMA β (x 0...t ) as:

EWMA β (x 0...t ) = DS δ EXP,β (x 0...t ) = (1 -β) • t i=0 β i • x t-i

The EWMA is a discounted sum average (modulo initialization bias), so it can be viewed as an estimator of the expectation of a random variable x if x 0...t ∼ x. Note that the momentum buffer g t from (1) is precisely an EWMA -specifically, g t = EWMA β (∇ L0...t (θ 0...t )).

It is well known that the exponential discount function is the only time-consistent (commonly "memoryless"), discount function -i.e. for any i, τ ≥ 0, the ratio d(i + τ )/d(i) depends only on τ . This is precisely why the EWMA can be tracked with no auxiliary memory -for example, as in momentum's update rule.

A.3 EWMA AS VARIANCE REDUCTION:

We now provide the following fact about the covariance of the EWMA when x 0...t are random variables. Fact A.1 (Limit covariance of EWMA). Assume that x 0...t are independent random vectors, each with the covariance matrix Σ. Then:

lim t→∞ Cov[EWMA β (x 0...t )] = (1 -β) 2 1 -β 2 • Σ Proof. Corollary of Fact A.2.

This means that arbitrary variance reduction of the EWMA is possible by increasing β. For example, β = 0.9 implies that the covariance is reduced to 1 19 • Σ, and β = 0.99 implies that the covariance is reduced to 1 199 • Σ. This provides an intuitive explanation of momentum as a variance reduction technique. Assuming that the momentum buffer is normalized (and thus interpretable as an estimator of the gradient), applying momentum will reduce the variance of the update steps, with higher β leading to more variance reduction.

However, the flip side is that higher β induces more bias (informally, "staleness") in the momentum buffer with respect to the true gradient, as the momentum buffer becomes extremely slow to update. Thus, the question arises: can we achieve variance reduction while guaranteeing that recent gradients contribute significantly to the update step? For this, we must introduce time-inconsistency.

A.4 (PURE) HYPERBOLIC DISCOUNTING AND HWMA Hyperbolic discounting, first proposed by Chung & Hernstein (1961) , is the classical timeinconsistent discount function in consumer choice. It is commonly used to model individual behaviors such as impatience. We consider its use in the setting of stochastic optimization (in place of the EWMA buffer of momentum).

For constants c, k > 0, we define the hyperbolic discount function as:12 δ H,c,k (i) = c 1 + ki and the hyperbolic weighted moving average HWMA c,k (x 0...t ) as:

HWMA c,k (x 0...t ) = DS δ H,c,k (x 0...t ) = c • t i=0 1 1 + ki • x t-i

Note that the hyperbolic discount function is time-inconsistent, since:

δ H,c,k (i + τ ) δ H,c,k (i) = 1 + ki 1 + k(i + τ )

depends on both i and τ .

Unlike the EWMA, the HWMA is not a discounted sum average -in fact, ∞ i=0 δ H,c,k (i) = ∞ holds regardless of choice of c or k. Thus, to use an HWMA of gradients in an optimization algorithm, c (or the learning rate α) must be decayed at a logarithmic rate. More concerning, however, is the computational inefficiency of the HWMA; specifically, the sum must be recomputed from scratch at each iteration from all past gradients. This is unacceptable for use in most practical applications.

However, in preliminary stochastic optimization experiments, we did observe a marked benefit of HWMA over EWMA (i.e. momentum), limiting the number of past gradients used for tractability. This indicates that time-inconsistency might be a useful property to have in a stochastic optimizer.

A.5 QUASI-HYPERBOLIC DISCOUNTING AND QHWMA Quasi-hyperbolic discounting, proposed by Phelps & Pollak (1968) and popularized in consumer choice by Laibson (1997) , seeks to qualitatively approximate the time-inconsistency of hyperbolic discounting by applying a discontinuous "upweighting" of the current step. Its tractability has resulted in much wider adoption in consumer choice vs. pure hyperbolic discounting, and we find that it is also more suited for use in practical optimization.

For constants ν ∈ R and β ∈ (-1, 1), we define the quasi-hyperbolic discount function as:13 

δ QH,ν,β (i) = 1 -νβ i = 0 ν(1 -β)β i i > 0

and the quasi-hyperbolic weighted moving average QHWMA ν,β (x 0...t ) as:

EQUATION

The QHWMA, like the EWMA, is a discounted sum average (modulo initialization bias), so it can also be viewed as an estimator under the same assumptions.

When ν = 1, the QHWMA is precisely the EWMA (with identical β), and the quasi-hyperbolic discount function is precisely the exponential discount function (and thus time-consistent). When ν = 1, the quasi-hyperbolic discount function, like the hyperbolic discount function, is timeinconsistent since:

δ QH,ν,β (i + τ ) δ QH,ν,β (i) = ν(1-β) 1-νβ β τ i = 0 β τ i > 0

depends on both i and τ ; specifically, i = 0 yields a different ratio than i > 0.

Note from (5) that the QHWMA is a ν-weighted average of the EWMA (with identical β) and x 0 . This means that the QHWMA can be easily computed online by simply keeping track of the EWMA, thus requiring no additional memory.

A.6 VARIANCE OF QHWMA

We now characterize the variance of a QHWMA using this fact:

Fact A.2 (Limit covariance of QHWMA). Assume that x 0...t are independent random vectors, each with the covariance matrix Σ. Then:

lim t→∞ Cov[QHWMA ν,β (x 0...t )] = ρ • Σ

where ρ is defined as:

ρ = (1 -νβ) 2 + [νβ(1 -β)] 2 1 -β 2 Proof. Provided in Appendix G.

ρ is essentially a scaling factor for the covariance of the QHWMA. It can be verified that ρ decreases (thus inducing variance reduction) with both increasing β and increasing ν.

A.7 MOTIVATING QHM:

This leads to our motivation for QHM, which simply replaces the EWMA momentum buffer with a QHWMA. Starting with any momentum parameterization (ν = 1 and β ∈ (0, 1)), β can be increased towards variance reduction (i.e. lowering ρ). Then, ν can be decreased to make the QH-WMA less biased as a gradient estimator, thus mitigating the aforementioned "staleness" problem. Note, however, that since decreasing ν will also increase ρ, we cannot simply decrease ν to zero. Specifically, any ν < 1 imposes a tight lower bound of (1 -ν) 2 on ρ, regardless of choice of β.

A.8 MOMENTUM AND QHM UPDATE RULES:

For completeness, we explicitly write the update rules for the momentum and QHM algorithms.

Momentum The momentum update rule is:

EQUATION

which can be efficiently written using an auxiliary buffer g t as:

g t+1 ← β • g t + (1 -β) • ∇ Lt (θ t ) (1, revisited from Section 2) θ t+1 ← θ t -α • g t+1

(2, revisited from Section 2)

QHM The QHM update rule is:

EQUATION

which can be efficiently written using an auxiliary buffer g t as:

g t+1 ← β • g t + (1 -β) • ∇ Lt (θ t ) (3, revisited from Section 3) θ t+1 ← θ t -α (1 -ν) • ∇ Lt (θ t ) + ν • g t+1 (4

, revisited from Section 3) QHM vs. momentum Comparing ( 2) and ( 4), QHM may seem at first glance identical to momentum with discount factor νβ. However, replacing the β in ( 6) with νβ yields:

θ t+1 ← θ t -α (1 -νβ) • t i=0 (νβ) i • ∇ Lt-i (θ t-i ) ← θ t -α   (1 -νβ) • ∇ Lt (θ t ) + (1 -νβ) differs from (7) • t i=1 (νβ) i differs from (7) •∇ Lt-i (θ t-i )   

which plainly differs from (7) -most notably, in the exponential discount factor (νβ) for past gradients. Thus, momentum with discount factor νβ does not recover QHM.

A.9 RELATED WORK IN VARIANCE REDUCTION:

Section 4 presents numerous connections to other optimization algorithms that shed light on both deterministic and stochastic convergence properties of QHM. However, we do not formally analyze the convergence properties of QHM from a variance reduction standpoint; this remains future work.

Here, we briefly discuss other work in variance reduction.

Finite sums Recently, much effort has been devoted towards reducing the variance of the stochastic gradients used in optimization algorithms. Perhaps the most widely-studied setting is the "finite sum", or offline, stochastic optimization setting. Methods analyzed in the finite-sum setting include SAG (Schmidt et al., 2013) , SAGA (Defazio et al., 2014) , SVRG (Johnson & Zhang, 2013) , FSVRG (Shang et al., 2017) , Katyusha (Allen-Zhu, 2017), and others. We do not comment in detail on the finite sum setting due to its limited practical applicability to large-scale deep learning; for a fuller discussion of such methods, see Kidambi et al. (2018) .

Momentum as variance reduction Some work in variance reduction has drawn an explicit connection to momentum. For example, Roux et al. (2018) propose a method involving Bayesian updates of gradient estimates, which induces adaptive gradient averaging. The authors note that this method boils down to momentum with an adaptive β.

B PID CONTROLLERS AND OPTIMIZATION:

We follow Recht (2018) in describing the connection between PID control and gradient-based optimization.

Continuous PID We slightly adapt the setting from Aaström & Hägglund (1995) . t denotes time. There is a setpoint (i.e. target state), r(t), and a process variable (i.e. current state), y(t). The error of the system is defined as e(t) def = r(t) -y(t). A "controller" outputs a control signal u(t), usually towards the goal of making the error zero. The controller's choice of u(t) affects y(t) in some unspecified manner.

A PID controller, parameterized by k P , k I , and k D , uses the control function:

u(t) = k P [e(t)] + k I t 0 e(t )dt + k D de(t) dt (8)

Here, the terms in the square brackets are typically referred to as the P, I, and D terms, respectively.

Discrete approximation:

In discrete time, the setpoint, process variable, and error are trivially discretized as r t , y t , and e t def = r t -y t , respectively. The I term, which we label w t , is discretized as:

w t ← w t-1 + e t

The D term, which we label v t , could be discretized as v t = e t -e t-1 (first differences). However, a low-pass filter is often applied to mitigate noise, thus resulting in:

v t ← βv t-1 + (1 -β)(e t -e t-1 )

We simplify exposition by considering e -1 , w -1 , and v -1 to be 0.

Finally, the PID control function ( 8) is trivially discretized as:

u t = k P • e t + k I • w t + k D • v t (9)

Optimization Recht (2018) relates optimization to PID control as follows:

y t = ∇ Lt (θ t ) ; r t = 0 ; e t def = r t -y t = -∇ Lt (θ t ) ; u t = θ t+1 -θ 0

That is, the process variable is the stochastic gradient, the controller's goal is to make this gradient zero, and the controller achieves this by choosing the next step's model parameters according to the update rule θ t+1 ← u t + θ 0 . The update rule for a PID control optimizer is thus:

e t ← -∇ Lt (θ t ) v t ← βv t-1 + (1 -β)(e t -e t-1 ) w t ← w t-1 + e t θ t+1 ← θ 0 + k P • e t + k I • w t + k D • v t

Recht demonstrates that PID in this setting encapsulates gradient descent, momentum, and NAG; for example, gradient descent is recovered when k P = k D = 0 and k I = α.

Intuition Finally, to provide some additional intuition, we can state the following fact about the D term (v t ): Fact B.1 (D term is gradient and momentum). v t can be written as:

v t = 1 -β β -∇ Lt (θ t ) + (1 -β) • t i=0 β i • ∇ Lt-i (θ t-i ) Proof. Provided in Appendix G.

Thus, the D term is simply a weighted sum of an EWMA of gradients (i.e. momentum buffer) and the current gradient, and a PID control optimizer's output is simply a weighted sum of the momentum buffer, the current gradient, and the sum of all past gradients.

C CONNECTIONS TO OTHER ALGORITHMS (IN-DEPTH):

This appendix presents a deeper theoretical treatment of Section 4.2 through Section 4.4, deriving and discussing connections between QHM and various other optimization algorithms.

C.1 LINEAR OPERATOR ANALYSIS:

Along the lines of Lessard et al. (2016) , we consider optimizers as linear operators, interrupted by a nonlinear step (the gradient evaluation). In this setting, optimizers have b internal state buffers, which we write as a stacked vector S t ∈ R b•p . Optimizers accept the current optimizer state (S t ) and gradient (∇ Lt (θ t )), and they produce the new optimizer state (S t+1 ) and parameters (θ t ) using a square matrix T ∈ R (b+2)p×(b+2)p .14 

Update rule For convenience, we impose the restriction that the output θ t can only depend on the state S t . Then, for analytical purposes, the optimizer can be written as the following update rule:

[S t+1 ∅ ∅] ← T S t ∇ Lt (θ t ) 0 p [∅ ∅ θ t+1 ] ← T [S t+1 0 p 0 p ]

where 0 p denotes the size-p zero vector, ∅ denotes throwaway values, and [v 1 v 2 . . .] denotes vector stacking.

Coordinate-wise decomposition Since we only consider optimizers that act coordinate-wise (except for the gradient evaluation), we can write T as the Kronecker product of a coordinate-wise transition matrix A ∈ R (b+2)×(b+2) and the identity matrix

I p . That is, T = A ⊗ I p .

Then, for t > 0, we can write θ t in terms of the initial state S 0,{1...b} and all past gradients, using the last row of various matrix powers of A:

EQUATION

C.2 QHM

The internal state of QHM includes two buffers: g t (momentum buffer) and θ t (model parameters).

The transition matrix T QHM , mapping from g t θ t ∇ Lt (θ t ) 0 p to [g t+1 θ t+1 0 p θ t ] , is:

A QHM =    β 0 1 -β 0 -ανβ 1 -α(1 -νβ) 0 0 0 0 0 0 1 0 0    T QHM = A QHM ⊗ I p

For n > 0, routine computation yields the last row of the (n + 1)-th matrix power:

A n+1 QHM 4 = -ανβ(1-β n ) 1-β 1 -α(1 -νβ n ) 0

Applying (10), the optimizer state θ t can be written as:

θ t = (A t+1 QHM ) 4,1 • g 0 + (A t+1 QHM ) 4,2 • θ 0 + t i=1 (A i+1 QHM ) 4,3 • ∇ Lt-i (θ t-i ) = - ανβ(1 -β t ) 1 -β g 0 + θ 0 -α • t i=1 (1 -νβ i ) • ∇ Lt-i (θ t-i )

In the typical case of g 0 = 0, we have:

EQUATION

C.3 PID Recht (2018) draws a strong connection between gradient-based optimization and PID control. We regurgitate the excellent exposition (with minor modifications) in Appendix B.

Update rule A PID control optimizer, parameterized by k P , k I , k D ∈ R, uses the update rule:

e t ← -∇ Lt (θ t ) v t ← β • v t-1 + (1 -β)(e t -e t-1 ) w t ← w t-1 + e t θ t+1 ← θ 0 + k P • e t + k I • w t + k D • v t

Coordinate-wise decomposition The internal state of a PID control optimizer includes four buffers: e t-1 (P term), w t-1 (I term), v t-1 (D term), and θ 0 (initial parameters). 15The transition matrix T PID , mapping from e t-1 w t-1 v t-1 θ 0 ∇ Lt (θ t ) 0 p to [e t w t v t θ 0 0 p θ t ] , is:

A PID =        0 0 0 0 -1 0 0 1 0 0 -1 0 -(1 -β) 0 β 0 -(1 -β) 0 0 0 0 1 0 0 0 0 0 0 0 0 k P k I k D 0 0 0        T PID = A PID ⊗ I p

For n > 0, routine computation yields the last row of the (n + 1)-th matrix power:

EQUATION

where:

A n+1 PID 6,5 = -[k P + k I + (1 -β)k D ] n = 1 -k I -(1 -β) 2 β n-2 k D n > 1

Applying (10), the optimizer state θ t can be written as:

θ t = A t+1 PID 6,1 • e -1 + A t+1 PID 6,2 • w -1 A t+1 PID 6,3 • v -1 + A t+1 PID 6,4 • θ 0 + t i=1 A i+1 PID 6,5 • ∇ Lt-i (θ t-i ) = - k D β t-1 (1 -β) β • e -1 + k I • w -1 + k D β n • v -1 + θ 0 + t i=1 A i+1 PID 6,5 • ∇ Lt-i (θ t-i )

Relationship with QHM In the typical case of e -1 = w -1 = v -1 = 0 p , we have:

θ t = θ 0 + t i=1 A i+1 PID 6,5 • ∇ Lt-i (θ t-i )

Then, equating with (11), we have that QHM is PID with:16 

k P = - ανβ 1 -β k I = α k D = ανβ 2 (1 -β) 2

or that PID is QHM with:

α = k I ν = k 2 P k D • k I β = k D k D -k P

Viewing β as a constant, the following restriction holds on the PID coefficients that QHM can recover:

k D k P = - β 1 -β

This restriction is looser than those for plain SGD (which has the additional restriction k P = k D = 0), momentum (which has the additional restriction k P /k I = k D /k P ), and NAG (which has the additional restriction k P /k I = βk D /k P ).

Viewing β as a hyperparameter, QHM can recover all PID coefficients except when k I = 0 (i.e. P, D, or PD controller), or k P = 0 = k D (i.e. PI controller).

To summarize, PID is a superfamily of QHM. Viewing β as a constant, QHM imposes a restriction on the ratio between k P and k D . Viewing β as a free variable, however, QHM can recover nearly all PID coefficients. C.4 SNV Section 6 of Lessard et al. (2016) describes a "synthesized Nesterov variant" algorithm, which we call "SNV" for convenience. This algorithm is used to analyze and improve optimizer robustness under "relative deterministic noise" (i.e. multiplicative noise of the gradient).17 

Update rule SNV, parameterized by γ, β 1 , β 2 ∈ R, uses the update rule:18 

ξ t+1 ← ξ t -γ • ∇ Lt (θ t ) + β 1 (ξ t -ξ t-1 ) θ t+1 ← ξ t+1 + β 2 (ξ t+1 -ξ t )

Coordinate-wise decomposition The internal state of a SNV optimizer includes two buffers: ξ t and ξ t-1 .

The transition matrix T SNV , mapping from ξ t ξ t-1 ∇ Lt (θ t ) 0 p to [ξ t+1 ξ t 0 p θ t ] , is:

A SNV =    1 + β 1 -β 1 -γ 0 1 0 0 0 0 0 0 0 1 + β 2 -β 2 0 0    T SNV = A SNV ⊗ I p

For n > 0, routine computation gives us the last row of the (n + 1)-th matrix power:

A n+1 SNV 4 = 1 1-β1 (1 + χ n ) -1 1-β1 (β 1 + χ n ) - γ β1(1-β1) (β 1 + χ n ) 0

where:

χ n = β n 1 (β 2 (1 -β 1 ) -β 1 )

Applying ( 10), the optimizer state θ t can be written as:

θ t = T t+1 SNV 4,1 • w0 + T t+1 SNV 4,2 • w 0 + t i=1 T i+1 SNV 4,3 • ∇ Lt-i (θ t-i ) = 1 1 -β 1 (1 + χ t ) • ξ 0 - 1 1 -β 1 (β 1 + χ t ) • ξ -1 - t i=1 γ β 1 (1 -β 1 ) (β 1 + χ i ) • ∇ Lt-i (θ t-i )

Relationship with QHM Initialize ξ 0 = ξ -1 = θ 0 . The optimizer state θ t is:

θ t = θ 0 - t i=1 γ β 1 (1 -β 1 ) (β 1 + χ i ) • ∇ Lt-i (θ t-i )

Then, equating with (11), we have that QHM is SNV with:

γ = α(1 -β) β 1 = β β 2 = (1 -ν) β 1 -β or that SNV is QHM with: α = γ 1 -β 1 ν = 1 - 1 -β 1 β 1 β 2 β = β 1

To summarize, QHM and SNV recover each other. By extension, QHM recovers the Robust Momentum method, which is a specific parameterization of SNV (Cyrus et al., 2018) . Moreover, since Robust Momentum recovers the Triple Momentum of Scoy et al. (2018) , QHM also recovers Triple Momentum.

C.5 ACCSGD Jain et al. (2017) and Kidambi et al. (2018) point out various failures of momentum and NAG in the setting of stochastic least squares optimization. This motivates their proposal of the AccSGD algorithm, which yields faster convergence over momentum and NAG in certain least-squares regression settings. Here, we discuss the formulation of Kidambi et al. (2018) .

Update rule AccSGD, parameterized by δ > 0, κ > 1, ξ ≤ √ κ, and < 1, uses the update rule:

EQUATION

Coordinate-wise decomposition The internal state of an AccSGD optimizer includes two buffers: wt (a buffer) and w t (the iterate, identical to θ t ).

The transition matrix T AccSGD , mapping from wt w t ∇ Lt (θ t ) 0 p to [ wt+1 w t+1 0 p θ t ] , is:

A AccSGD =      1 - 2 ξ κ 2 ξ κ -δ ξ 0 ξ κ+ ξ 1 - 2 ξ κ κ+ 3 ξ 2 /κ κ+ ξ -δ κ+ 2 ξ 2 κ+ ξ 0 0 0 0 0 0 1 0 0      T AccSGD = A AccSGD ⊗ I p

For n > 0, routine computation gives us the last row of the (n + 1)-th matrix power:

EQUATION

where:

χ = κ -2 ξ κ + ξ

Applying (10), the optimizer state θ t can be written as:

θ t = T t+1 AccSGD 4,1 • w0 + T t+1 AccSGD 4,2 • w 0 + t i=1 T i+1 AccSGD 4,3 • ∇ Lt-i (θ t-i ) = κ -2 ξ -κ -2 ξ χ t κ(1 + ) • w0 + (κ + ξ) + κ -2 ξ χ t κ(1 + ) • w 0 - t i=1 δ (1 + ξ) -( ξ -1)χ i 1 + • ∇ Lt-i (θ t-i )

Relationship with QHM Fix ∈ (0, 1), and initialize w0 = w 0 = θ 0 . The optimizer state θ t is:

θ t = θ 0 - t i=1 δ (1 + ξ) -( ξ -1)χ i 1 + • ∇ Lt-i (θ t-i )

Then, equating with (11), we have that QHM is AccSGD with:

δ = α(1 -ν) κ = (β + )( ν + 1) (1 -ν)(1 -β) ξ = ν + 1 (1 -ν)

or that AccSGD is QHM with:

α = δ (1 + ξ) 1 + ν = ξ -1 (1 + ξ) β = κ -2 ξ κ + ξ

AccSGD cannot recover NAG Based on the above analysis, NAG (i.e. ν = β) is recovered 19 This disproves the claim in Kidambi et al. (2018) that AccSGD recovers NAG when ξ = √ κ. In fact, we demonstrate that AccSGD cannot recover NAG at all. For ∈ (0, 1) and the aforementioned value of ξ, we have that ξ > √ κ:

when ξ = 1 2 1 -+ 4κ + (1 -) 2 .

ξ = 1 2 1 -+ 4κ + (1 -) 2 > 1 2 4κ + (1 -) 2 > 2 √ κ 2 > √ κ

Since AccSGD requires that ξ ≤ √ κ and that ∈ (0, 1)20 , AccSGD cannot recover NAG.

To summarize, QHM recovers AccSGD. In the reverse direction, AccSGD does not recover QHM; specifically, we disprove the claim in Kidambi et al. (2018) that AccSGD recovers NAG. Since QHM recovers NAG, AccSGD cannot fully recover QHM.

D GENERAL TWO-STATE OPTIMIZER:

This appendix describes a generic two-state optimizer ("TSO") where one of the states is the iterate (θ t ) and the other is an auxiliary buffer (a t ). The optimizer is parameterized by h, k, l, m, q, z ∈ R, and the update rule is:

a t+1 ← h • a t + k • θ t + l • ∇ Lt (θ t ) θ t+1 ← m • a t + q • θ t + z • ∇ Lt (θ t )

We can write this as a transition matrix T TSO ∈ R 3×3 :

T TSO = h k l m q z 0 0 0

To simplify further derivations we diagonalize T TSO as:

φ = (h -q) 2 + 4km ψ = km -hq Q TSO =   lq-kz ψ h-q-φ 2m h-q+φ 2m hz-lm ψ 1 1 1 0 0   Λ TSO =   0 0 0 0 1 2 (h + q -φ) 0 0 0 1 2 (h + q + φ)   T TSO = Q TSO Λ TSO Q -1

TSO:

Relationship with QHM If ψ = 0, φ = 0, φ ≤ 1, 1 2 (h + q + φ) = 1, h -q + φ = 0, and 1 -l = 1 2 (h + q -φ), then QHM implements the TSO optimizer with:

g 0 = - (2 -(h + q -φ))ψ (h + q -φ)(lq -kz) • a 0 α = 1 2ψφ [(h -q -φ)(lm -hz) + 2m(lq -kz)] ν = 2m(lq -kz) (h -q -φ)(lm -hz) + 2m(lq -kz) β = 1 2 (h + q -φ)

Proof. We can write down the unrolled TSO update rule for θ t , as follows:

θ t ← (T t TSO ) 2,1 • a 0 + (T t TSO ) 2,2 • θ 0 + t-1 i=0 (T t-i TSO ) 2,3 • ∇ Li (θ i )

Similarly, for QHM we can define a transition matrix T QHM ∈ R 3×3 that advances state g t θ t ∇ Lt (θ t ) as:

T QHM = β 0 (1 -β) -ανβ 1 -α(1 -νβ) 0 0 0

Thus, the unrolled update rule for QHM takes the following form:

θ t ← (T t QHM ) 2,1 • g 0 + (T t QHM ) 2,2 • θ 0 + t-1 i=0 (T t-i QHM ) 2,3 • ∇ Li (θ i )

Now we match the corresponding coefficients in both of the update rules to establish dependencies:

(T t-i QHM ) 2,3 • ∇ Li (θ i ) = (T t-i TSO ) 2,3 • ∇ Li (θ i ) ∀i ∈ [0, t -1] (T t QHM ) 2,1 • g 0 + (T t QHM ) 2,2 • θ 0 = (T t TSO ) 2,1 • a 0 + (T t TSO ) 2,2

• θ 0 By solving the first equation we can establish values for α, β, and ν:

(T t-i QHM ) 2,3 = (T t-i TSO ) 2,3 ∀i ∈ [0, t -1] (T t-i QHM ) 2,3 = (Q TSO Λ t-i TSO Q -1 TSO ) 2,3 ∀i ∈ [0, t -1] -α(1 -νβ t-i ) = 1 2ψφ (Λ t-i TSO ) 2,2 [(h -q + φ)(lm -hz) + 2m(lq -kz)] - 1 2ψφ (Λ t-i TSO ) 3,3 [(h -q -φ)(lm -hz) + 2m(lq -kz)] ∀i ∈ [0, t -1]

Per our assumption (Λ TSO ) 3,3 = 1 2 (h + q + φ) = 1 and h -q + φ = 0, we can recover the following relationships:

β = (Λ TSO ) 2,2 = 1 2 (h + q -φ) α = 1 2ψφ [(h -q -φ)(lm -hz) + 2m(lq -kz)] ν = 2m(lq -kz) (h -q -φ)(lm -hz) + 2m(lq -kz)

We can solve the second equation to find g 0 :

(T t QHM ) 2,1 • g 0 + (T t QHM ) 2,2 • θ 0 = (T t TSO ) 2,1 • a 0 + (T t TSO ) 2,2 • θ 0 -αβν 1 -β t 1 -β • g 0 + θ 0 = m φ (Λ t TSO ) 3,3 -(Λ t TSO ) 2,2 • a 0 + 1 2φ (Λ t TSO ) 2,2 (h -q + φ) -(Λ t TSO ) 3,3 (h -q -φ) • θ 0

Given that (Λ TSO ) 3,3 = 1 and h -q + φ = 0 =⇒ h -q -φ = -2φ, we can simplify:

-αβν 

1 -β t 1 -β • g 0 + θ 0 = m φ (1 -β t ) • a 0 + θ 0 g 0 = - (1 -β)m αβνφ • a 0 g 0 = - (2 -(h + q -φ))ψ (h + q -φ)(lq -kz) • a 0 as desired.

e t ← -∇ Lt (θ t ) v t ← β • v t-1 -(1 -β)(e t -e t-1 ) w t ← β • w t-1 + r • e t u t ← w t + k D • v t θ t+1 ← θ t + u t

Discussion This setting departs somewhat from typical PID control, in that the signal u t controls the derivative of the controller's output (i.e. θ t+1 -θ t ) rather than the output itself (i.e. θ t+1 -θ 0 ). To avoid parameter blowup, this formulation necessitates the addition of exponential decay to the I term, with discount factor β. 22The I term thus becomes the momentum buffer. However, recall from Fact B.1 that the D term is a weighted sum of the momentum buffer and the P term. It follows that the D term is a weighted sum of the P and I terms, and that this setting is degenerate (either "PI" or "PD").

As a consequence, the proposed PID algorithm of An et al. (2018) is less expressive than that of Recht (2018) . Specifically, applying Fact B.1 demonstrates a mapping into QHM:

α = r 1 -β ν = 1 - k D (1 -β) 2 rβ β = unchanged

Efficiency This PID control optimizer is costlier than QHM. It requires 2 auxiliary buffers of memory. Computationally, it requires 2 in-place scalar-vector multiplications and 5 scaled vector additions per update step.

F (QH)ADAM'S UPDATE BOUND:

This appendix elaborates on Adam and QHAdam's stability properties through the lens of a step size upper bound.

It is well known that the training process for deep learning models can often "explode" due to a very small number of large parameter updates. With Adam, these large updates can occur if there exist parameters whose stochastic gradients are almost always near zero but incur rare "spikes". 23 . This is because the square root of the second moment estimate, used in normalizing the gradient for the update step, will be far below the magnitude of these spikes.

There are three main ways to address this instability:

• Firstly, one can simply decrease the learning rate α. However, this may be undesirable due to slower training.

• Secondly, one can increase the hyperparameter. However, the appropriate setting of depends on the exact magnitudes of these gradient spikes, which is often unknown. Setting too high effectively turns Adam into SGD. Thus, setting often reduces to guesswork.

• Thirdly, one can clip gradients. However, the appropriate magnitude of the gradient clipping also depends on the stochastic gradient distribution. Thus, this solution also involves a fair amount of guesswork.

However, Adam does provide a useful guarantee -unlike SGD, Adam has an upper bound on the per-step update (Kingma & Ba, 2015) . This upper bound is independent of the gradient distribution (or even temporal correlation), depending only on the hyperparameters α, β 1 , and β 2 . Thus, no matter the gradient distribution, Adam will restrict the magnitude of the per-step updates to some known constant. Kingma & Ba (2015) intuitively describe this bound as "establishing a trust region around the current parameter value".

We show that the step size upper bound claimed in Section 2.1 of Kingma & Ba (2015) is incorrect, by providing the correct tight bound for both Adam and QHAdam. We then demonstrate that with QHAdam, one can lower the maximum per-step update (and thus improve stability) simply by lowering ν 2 to be below 1.

F.1 SETTING:

We make two simplifications. Firstly, we fix = 0.24 Secondly, we remove the bias correction of the moment estimators (i.e. we use g t+1 ← g t+1 and s t+1 ← s t+1 ).

In this setting, QHAdam applies the following update rule:

EQUATION

where:

EQUATION

F.2 IMPLICIT UPDATE BOUND:

We now bound QHAdam's update (before scaling by α) by a constant dependent only on β 1 , β 2 , ν 1 , and ν 2 :

Fact F.1 (QHAdam tight upper bound on update). Assume that st+1 is nonzero at each coordinate and that 0 < β 1 < √ β 2 < 1. Then, the following per-coordinate tight upper bound holds:

gt+1 √ st+1 ∞ ≤ (1 -ν 1 β 1 ) 2 1 -ν 2 β 2 + [ν 1 β 1 (1 -β 1 )] 2 1 - β 2 1 β2 t ν 2 (1 -β 2 )(β 2 -β 2 1 )

Proof. Firstly and without loss of generality, we can treat the gradients as single coordinates x i ∈ R.

That is,

x i = ∇ Li (θ i ) ∈ R for i ∈ [0, t].

We perform the following simplification of ( 12) and ( 13):

EQUATION

We now wish to find the values of x i that maximize g2 t+1 st+1 . Applying ( 14) and ( 15), these values are characterized by the following first-order conditions:

x i =    st+1 gt+1 ν1(1-β1)β i 1 ν2(1-β2)β i 2 i < t st+1 gt+1 1-ν1β1 1-ν2β2 i = t (16)

Since the quantity g2 t+1 st+1 is invariant to scalar multiplication of all x i , we can simplify (16) to:

x i = ν1(1-β1)β i 1 ν2(1-β2)β i 2 i < t 1-ν1β1 1-ν2β2 i = t (17)

Plugging the values of x i from ( 17) into ( 14) and (15) yields:

max x0...t g2 t+1 st+1 = (1 -ν 1 β 1 ) 2 1 -ν 2 β 2 + [ν 1 β 1 (1 -β 1 )] 2 1 - β 2 1 β2 t ν 2 (1 -β 2 )(β 2 -β 2 1 )

The desired result follows immediately.

Limit case Consider the limit case of t → ∞. Then, the bound in Fact F.1 simplifies to:

EQUATION

For vanilla Adam (i.e. ν 1 = ν 2 = 1), (18) simplifies further to:

EQUATION

Note that since the bound in ( 19) is tight, this result contradicts the claim in Section 2.1 of Kingma & Ba (2015) that Adam's per-coordinate step size is bounded above by α • max{1, (1β 1 )/ √ 1 -β 2 }. 25 In the following discussion, we use the correct bounds from ( 18) and (19).

F.3 DISCUSSION:

The recommended vanilla Adam setting of β 2 = 0.999 in Kingma & Ba (2015) makes the right-hand side of ( 19) to be large, and various work has employed Adam with a significantly lower β 2 ; e.g. 0.98 (Vaswani et al., 2017; Ott et al., 2018) . 26 Decreasing β 2 is undesirable, often slowing down training. 27 Moving from Adam to QHAdam, an alternative solution is to decrease ν 2 to be below 1. This decreases the right-hand side of ( 18), up to a point, and thus imposes a tighter constraint on the magnitudes of updates than the vanilla Adam setting of ν 2 = 1. Fig. 3 shows an example of this phenomenon using a fixed ν 1 , β 1 , and β 2 . 

:

https://github.com/facebookresearch/qhoptim/

Additionally, Kingma & Ba (2015) recommends β1 = 0.9 for Adam, and β1 = 0.9 is the default for Adam in both the PyTorch and TensorFlow frameworks(Paszke et al., 2017;Abadi et al., 2015).

The appendix also sheds some light on the nomenclature, adopted from the hyperbolic discounting work pioneered byChung & Hernstein (1961),Phelps &Pollak (1968), andLaibson (1997) in consumer choice. We caution that "quasi-hyperbolic" does not directly relate to the geometry of hyperbolas.

Appendix C presents a deeper theoretical treatment of Section 4.2 through Section 4.4, which are largely narrative.

The learning rate is α in the original paper; we use γ to avoid confusion with QHM's α.

Here, κ is the ratio between the Lipschitz constant of ∇L(•) and the strong convexity parameter of L(•).

Here, we did not sweep over alternate parameterizations.

We also train the model using plain SGD, again finding that plain SGD performs nearly as well as NAG throughout training. Although not shown, plain SGD in fact performs better than momentum. The validation loss curves for plain SGD, momentum, and NAG are indistinguishable throughout training, suggesting that momentum/NAG is not needed inGoyal et al. (2017).

Here, we tried higher values of β1. Significantly increasing β1 was not fruitful for either algorithm.

Refer to Appendix F for a more detailed theoretical treatment.

Here, we tried two other parameterizations (higher β1) with marginal success.

Slightly adapted from the original formulation.

Significantly adapted from (but still equivalent to) the original formulation.

Typically, the output would be θt+1; however, we lag one step so that the parameters can be specified in terms of the current-step state. This is done purely for notational convenience.

The offset of -1 in the P, I, and D term subscripts is purely for convenience.

This is inconsistent with Recht (2018)'s derivation by a factor of (1 -β) in kP , kD and (1 -β)/β in kI . While (1 -β) is explainable as the difference between a normalized and unnormalized momentum buffer, we suspect that the extra factor of β in kI is a mistake in the original derivation.

This is similar to the vanishing gradient assumption used inLoizou & Richtárik (2017), which removes the need to consider variance reduction of the gradient.

The learning rate is α in the original paper; we use γ to avoid confusion with QHM's α.

Empirical simulations confirm this finding.

The recommended value for is 0.7(Kidambi et al., 2018).

The exponential discount is α in the original paper; we use β to avoid confusion.

One final oddity is that the D term (vt) calculates the negation of the derivative.

Somewhat more concretely, a stochastic gradient is "spiky" if its distribution has extremely large higherorder cumulants relative to its typical magnitudes.

24 The analysis with free follows the same approach but is somewhat messier.

The difference can be rather large -for example, for the recommended Adam parameters of β1 = 0.9 and β2 = 0.999,Kingma & Ba (2015) claim an upper bound of 3.16 • α, while (19) implies a tight bound of 7.27 • α.

We performed experiments on these models indicating that increasing β2 far beyond 0.98 led to training explosion. We suspect that these instability issues are especially prevalent in settings with rare inputs or labels, such as machine translation.

In proposing the AdamNC algorithm,Reddi et al. (2018) suggests that β2 should be high to capture a sufficiently long history of past gradients.

https://github.com/AtheMathmo/AggMo

These learning rates may seem high, but recall that the effective step size is identical to that of "typical", unnormalized momentum/NAG with α ∈ {0.1, 0.01, 0.001} and β = 0.9.

https://github.com/sfujim/TD3

https://github.com/facebookresearch/qhoptim/releases/download/ emptytag/qhoptim_parameter_sweep_data.tar.gz

Figures and Tables:

Figure 1: Selected parameter sweep results (full results in Appendix J). Top row shows train loss. Bottom row shows validation error. "Best" refers to the optimal parameterization within the sweep, with respect to the metric. Shaded bands indicate ±1 standard deviation.

Figure 2: Case study results. Top row: RN152-Imagenet-QHM (left), FConvLM-WikiText103-QHM (center), and TF-WMT16ENDE-QHAdam (right). Bottom row: TD3-MuJoCo-QHAdam Ant (left), HalfCheetah (center), Hopper (right). Shaded bands indicate ±1 standard deviation.

Figure 3: Bound from (18), fixing ν 1 = 0.8, β 1 = 0.95, and β 2 = 0.98, and varying ν 2 .

Summary of connections between QHM and other optimization algorithms

Table content: <table><tr><td>ALGORITHM</td><td colspan="2">RELATION  *  EFF.  †</td><td>BRIEF NOTES</td></tr><tr><td>Plain SGD</td><td>subfamily</td><td colspan="2">better recovered by QHM with ν = 0</td></tr><tr><td>Momentum</td><td/><td/></tr></table>

Summary of experimental settingsPractical notes We leave formal convergence analysis to future work. However, a couple of informal points are worth mentioning. Firstly, ν 1 and β 1 can be reasoned about in a similar manner to QHM's ν and β. Secondly, when replacing Adam with QHAdam, setting ν 2 = 1 and β 2 unchanged is usually reasonable if the original Adam training is stable. Thirdly, when Adam training is not stable, it is possible that setting ν 2 < 1 can improve stability by imposing a tighter step size bound -Appendix F elaborates, disproving the step size bound claimed inKingma & Ba (2015).

Table content: <table><tr><td>SHORT NAME</td><td/><td>MODEL</td><td colspan="2">DATASET/TASK OPTIMIZER</td></tr><tr><td>Logistic-EMNIST-QHM PS</td><td colspan="3">logistic regression EMNIST digits</td><td>QHM</td></tr><tr><td>Logistic-EMNIST-QHAdam PS</td><td colspan="3">logistic regression EMNIST digits</td><td>QHAdam</td></tr><tr><td>MLP-EMNIST-QHM PS</td><td colspan="3">3-layer tanh MLP EMNIST digits</td><td>QHM</td></tr><tr><td>MLP-EMNIST-QHAdam PS</td><td colspan="3">3-layer tanh MLP EMNIST digits</td><td>QHAdam</td></tr><tr><td>RN18-CIFAR10-QHM PS</td><td colspan="2">PreActResNet18</td><td>CIFAR10</td><td>QHM</td></tr><tr><td>RN50-ImageNet-QHM PS</td><td colspan="2">ResNet50</td><td>ILSVRC2012</td><td>QHM</td></tr><tr><td>RN152-ImageNet-QHM CS</td><td colspan="2">ResNet152</td><td>ILSVRC2012</td><td>QHM</td></tr><tr><td colspan="3">FConvLM-WikiText103-QHM CS FConvLM</td><td>WikiText-103</td><td>QHM</td></tr><tr><td>TD3-MuJoCo-QHAdam CS</td><td>TD3</td><td/><td>MuJoCo</td><td>QHAdam</td></tr><tr><td>TF-WMT16ENDE-QHAdam CS</td><td colspan="2">Transformer</td><td colspan="2">WMT16 EN-DE QHAdam</td></tr><tr><td colspan="2">PS Parameter sweep experiment (Section 6.1).</td><td colspan="2">CS Case study (Section 6.2).</td></tr><tr><td colspan="4">Efficiency QHAdam incurs four extra scaled vector additions over Adam.</td></tr></table>

Case study results

Table content: <table><tr><td>SHORT NAME</td><td>QH</td><td>BASELINE</td></tr><tr><td>RN152-ImageNet-QHM ERR</td><td>0.2128±0.0005</td><td>0.2137±0.0011</td></tr><tr><td>FConvLM-WikiText103-QHM PPL</td><td>34.45±0.17</td><td>34.92±0.33</td></tr><tr><td>TD3-MuJoCo-QHAdam HalfCheetah AR</td><td colspan="2">10001.66±1106.50 9829.95±872.29</td></tr><tr><td>TD3-MuJoCo-QHAdam Hopper AR</td><td>2948.14±1078.70</td><td>2401.39±1566.00</td></tr><tr><td>TD3-MuJoCo-QHAdam Walker2d AR</td><td>4282.15±707.78</td><td>4356.95±860.73</td></tr><tr><td>TD3-MuJoCo-QHAdam Ant AR</td><td>4657.12±881.03</td><td>4303.32±1279.86</td></tr><tr><td>TD3-MuJoCo-QHAdam InvPendulum AR</td><td>986.30±59.70</td><td>950.27±186.63</td></tr><tr><td>TD3-MuJoCo-QHAdam InvDoublePendulum AR</td><td>8890.59±1941.56</td><td>9333.49±20.02</td></tr><tr><td>TD3-MuJoCo-QHAdam Reacher AR</td><td>-3.99±0.55</td><td>-3.98±0.56</td></tr><tr><td>TF-WMT16ENDE-QHAdam BLEU</td><td>29.45±0.06  *</td><td>29.17±0.07</td></tr></table>

E AN ALTERNATIVE PID SETTINGWe very briefly comment onAn et al. (2018)'s PID control setting.Update ruleAn et al. (2018)'s PID control optimizer, parameterized by r, k D , β > 0, uses the following update rule:21

Table content: <table/>

Peer review (iclr2019):
 The authors introduce a class of quasi-hyperbolic algorithms that mix SGD with SGDM (or similar with Adam) and show improved empirical results. They also prove theoretical convergence of the methods and motivate the design well. The paper is well-written and contained the necessary references. Although I did feel that the authors could have better compared their method against the recent AggMom (Aggregated Momentum: Stability Through Passive Damping by Lucas et al.). Seems like there are a few similarities there. I enjoyed reading this paper and endorse it for acceptance. The theoretical results presented and easy to follow and state the assumptions clearly. I appreciated the fact that the authors aimed to keep the paper self-contained in its theory. The numerical experiments are thorough and fair. The authors test the algorithms on an extremely wide set of problems ranging from image recognition (including CIFAR and ImageNet), natural language processing (including the state-of-the-art machine translation model), and reinforcement learning (including MuJoCo). I have not seen such a wide comparison in any paper proposing training algorithms before. Further, the numerical experiments are well-designed and also fair. The hyperparameters are chosen carefully, and both training and validation errors are presented. I also appreciate that the authors made the code available during the reviewing phase. Out of curiosity, I ran the code on some of my workflows and found that there was some improvement in performance as well.