OPTIMAL CONVERGENCE OF ADAPTIVE POLYAK’S HEAVY-BALL METHODS
Year: 2021

Abstract:
The adaptive stochastic gradient descent (SGD) with momentum has been widely adopted in deep learning as well as convex optimization. In practice, the last iterate is commonly used as the final solution. However, the available regret analysis and the setting of constant momentum parameters only guarantee the optimal convergence of the averaged solution. In this paper, we fill this theory-practice gap by investigating the convergence of the last iterate (referred to as individual convergence), which is a more difficult task than convergence analysis of the averaged solution. Specifically, in the constrained convex cases, we prove that the adaptive Polyak’s Heavy-ball (HB) method, in which the step size is only updated using the exponential moving average strategy, attains an individual convergence rate of O( 1 √ t ), as opposed to that of O( log t √ t ) of SGD, where t is the number of iterations. Our new analysis not only shows how the HB momentum and its timevarying weight help us to achieve the acceleration in convex optimization but also gives valuable hints how the momentum parameters should be scheduled in deep learning. Empirical results validate the correctness of our convergence analysis in optimizing convex functions and demonstrate the improved performance of the adaptive HB methods in training deep networks.

text: t ), as opposed to that of O( log t√ t ) of SGD, where t is the number of
iterations. Our new analysis not only shows how the HB momentum and its timevarying weight help us to achieve the acceleration in convex optimization but also gives valuable hints how the momentum parameters should be scheduled in deep learning. Empirical results validate the correctness of our convergence analysis in optimizing convex functions and demonstrate the improved performance of the adaptive HB methods in training deep networks.

heading: 1 INTRODUCTION
text: One of the most popular optimization algorithms in deep learning is the momentum method (Krizhevsky et al., 2012). The first momentum can be traced back to the pioneering work of Polyak’s heavy-ball (HB) method (Polyak, 1964), which helps accelerate stochastic gradient descent (SGD) in the relevant direction and dampens oscillations (Ruder, 2016). Recent studies also find that the HB momentum has the potential to escape from the local minimum and saddle points (Ochs et al., 2014; Sun et al., 2019a). From the perspective of theoretical analysis, HB enjoys a smaller convergence factor than SGD when the objective function is twice continuously differentiable and strongly convex (Ghadimi et al., 2015). In nonsmooth convex cases, with suitably chosen step size, HB attains an optimal convergence rate of O( 1√
t ) in terms of the averaged output (Yang et al., 2016),
where t is the number of iterations.
∗Equal contribution †Corresponding author
To overcome the data-independent limitation of predetermined step size rules, some adaptive gradient methods have been proposed to exploit the geometry of historical data. The first algorithm in this line is AdaGrad (Duchi et al., 2011). The intuition behind AdaGrad is that the seldom-updated weights should be updated with a larger step size than the frequently-updated weights. Typically, AdaGrad rescales each coordinate and estimates the predetermined step size by a sum of squared past gradient values. As a result, AdaGrad has the same convergence rate as vanilla SGD but enjoys a smaller factor especially in sparse learning problems. The detailed analysis of AdaGrad (Mukkamala & Hein, 2017) implies that one can derive similar convergence rates for the adaptive variants of the predetermined step size methods without additional difficulties.
Unfortunately, experimental results illustrate that AdaGrad under-performed when applied to training deep neural newtworks (Wilson et al., 2017). Practical experience has led to the development of adaptive methods that is able to emphasize the more recent gradients. Specifically, an exponential moving average (EMA) strategy was proposed in RMSProp to replace the cumulative sum operation (Tieleman & Hinton, 2012). Adam (Kingma & Ba, 2014), which remains one of the most popular optimization algorithms in deep learning till today, built upon RMSProp together with updating the search directions via the HB momentum. Generally speaking, the gradient-based momentum algorithms that simultaneously update the search directions and learning rates using the past gradients are referred to as the Adam-type methods (Chen et al., 2019). These kinds of methods have achieved several state-of-the-art results on various learning tasks (Sutskever et al., 2013).
Compared with HB and AdaGrad, the main novelty of Adam lies in applying EMA to gradient estimate (first-order) and to element-wise square-of-gradients (second-order), with the momentum parameter β1t and step size parameter β2t (see (6)) (Alacaoglu et al., 2020). However, the use of EMA causes a lot of complexities to the convergence analysis. For example, in the online setting, (Kingma & Ba, 2014) offered a proof that Adam would converge to the optimum. Despite its remarkable practicality, Adam suffers from the non-convergence issue. To overcome its advantages, several variants such as AMSGrad and AdamNC were proposed (Reddi et al., 2018). Unfortunately, the regret bound of AMSGrad in (Reddi et al., 2018) is O( √ log t √ t) for nonsmooth convex problems, as opposed to that of O( √ t) of SGD. On the other hand, EMA uses the current step size in exponential moving averaging while the original HB can use the previous information (Zou et al., 2018). This will lead the update to stagnate when β1t is very close to 1. Fortunately, such a dilemma will not appear in Polyak’s HB method and a simple proof on the convergence of this kind of Adams in smooth cases has been provided (Défossez et al., 2020).
In this paper, we will focus on the adaptive Polyak’s HB method, in which the step size is only updated using EMA. Despite various reported practical performance for the Adam-type methods, there still exist some gaps between theoretical guarantees and empirical success.
• First of all, some important regret bounds have been established to guarantee the performance of online Adam-type algorithms. Nevertheless, the online-to-batch conversion can inevitably lead the solution of the induced stochastic algorithm to take the form of averaging of all the past iterates. In practice, the last iterate is popularly used as the final solution, which has the advantage of readily enforcing the learning structure (Chen et al., 2012). For SGD, the convergence of the last iterate, which is referred to as individual convergence in (Tao et al., 2020b), was posed as an open problem (Shamir, 2012). Only recently, its optimal individual convergence rate is proved to be O( log t√
t ) and O( log tt ) for general and
strongly convex problems respectively (Harvey et al., 2019; Jain et al., 2019). Despite enjoying the optimal averaging convergence (Yang et al., 2016), as far as we know, the individual convergence about the adaptive HB has not been discussed.
• Secondly, the momentum technique is often claimed as an accelerated strategy in machine learning community. However, almost all the theoretical analysis is only limited to the Nesterov’s accelerated gradient (NAG) (Nesterov, 1983) method especially in smooth cases (Hu et al., 2009; Liu & Belkin, 2020), which accelerates the rate of SGD from O( 1t ) to O( 1t2 ). While the individual convergence of HB is also concerned in some papers (Sebbouh et al., 2020; Sun et al., 2019b), the considered problem is limited to smooth and the derived rate is not optimal in convex cases. It is discovered that NAG is capable of accelerating the rate of individual convergence of SGD from O( log t√
t ) to O( 1√ t ) (Tao et al.,
2020a) in nonsmooth convex cases. Nevertheless, there is still a lack of the report about the acceleration of the adaptive HB. • Finally, in practice, almost all the momentum and Adam-type algorithms are often used
with a constant momentum parameter β1t (typically between 0.9 and 0.99). In theory, regret guarantees in the online Adam require a rapidly decaying β1t → 0 schedule, which is also considered in (Sutskever et al., 2013; Orvieto et al., 2019). This gap is recently bridged by getting the same regret bounds as that in (Reddi et al., 2018) with a constant β1t (Alacaoglu et al., 2020). In each state-of-the-art deep learning library (e.g. TensorFlow, PyTorch, and Keras), HB is named as SGD with momentum and β1t is empirically set to 0.9 (Ruder, 2016). Despite its intuition in controlling the number of forgotten past gradients and guarantee in optimal averaging convergence (Yang et al., 2016), how β1t affects individual convergence has not been discussed (Gitman et al., 2019).
The goal of this paper is to close a theory-practice gap when using HB to train the deep neural networks as well as optimize the convex objective functions. Specifically,
• By setting β1t = tt+2 , we prove that the adaptive HB attains an individual convergence rate of O( 1√
t ) (Theorem 5), as opposed to that of O( log t√ t ) of SGD. Our proof is different
from all the existing analysis of averaging convergence. It not only provides a theoretical guarantee for the acceleration of HB but also clarifies how the momentum and its parameter β1t help us to achieve the optimal individual convergence. • If 0 ≤ β1t ≡ β < 1, we prove that the adaptive HB attains optimal averaging convergence
(Theorem 6). To guarantee the optimal individual convergence, Theorem 5 suggests that time-varying β1t can be adopted. Note β1t = tt+2 → 1, thus our new convergence analysis not only offers an interesting explanation why we usually restrict β1t → 1 but also gives valuable hints how the momentum parameters should be scheduled in deep learning.
We mainly focus on the proof of individual convergence of HB (Theorem 3, Appendix A.1). The analysis of averaging convergence (Theorem 4) is simpler. Their extensions to adaptive cases are slightly more complex (Theorem 5 and 6), but it is similar to the proof of AdaGrad (Mukkamala & Hein, 2017) and the details can be found in the supplementary material.

heading: 2 PROBLEM STATEMENT AND RELATED WORK
text: Consider the following optimization problem,
min f(w), s.t. w ∈ Q. (1) where Q ⊆ Rd is a closed convex set and f(w) is a convex function. Denote that w∗ is an optimal solution and P is the projection operator on Q. Generally, the averaging convergence is defined as
f(w̄t)− f(w∗) ≤ (t), (2) where w̄t = 1t ∑t i=1 wi and (t) is the convergence bound about t. By contrast, the individual convergence is described as f(wt)− f(w∗) ≤ (t). (3)
Throughout this paper, we use g(wt) to denote the subgradient of f at wt. Projected subgradient descent (PSG) is one of the most fundamental algorithms for solving problem (1) (Dimitri P. et al., 2003), and the iteration of which is
wt+1 = P [wt − αtg(wt)], where αt > 0 is the step size. To analyze the convergence, we need the following assumption.
Assumption 1. Assume that there exists a number M > 0 such that
‖g(w)‖ ≤M, ∀w ∈ Q. It is known that the optimal bound for the nonsmooth convex problem (1) is O( 1√
t ) (Nemirovsky &
Yudin, 1983). PSG can attain this optimal convergence rate in terms of the averaged output while its optimal individual rate is only O( log t√
t ) (Harvey et al., 2019; Jain et al., 2019).
When Q = RN , the regular HB for solving the unconstrained problem (1) is wt+1 = wt − αtg(wt) + βt(wt −wt−1). (4)
If 0 ≤ βt ≡ β < 1, the key property of HB is that it can be reformulated as (Ghadimi et al., 2015)
wt+1 + pt+1 = wt + pt − αt
1− β g(wt),where pt =
β
1− β (wt −wt−1). (5)
Thus its convergence analysis makes almost no difference to that of PSG. Especially, if αt ≡ α√T , its averaging convergence rate is O( 1√
T ) (Yang et al., 2016), where T is the total number of iterations.
Simply speaking, the regular Adam (Kingma & Ba, 2014) takes the form of
wt+1 = wt − α√ t V − 12 t ĝt,
where ĝ(wt) is a unbiased estimation of g(wt) and ĝt = β1tĝt−1 + (1− β1t)ĝ(wt), Vt = β2tVt−1 + (1− β2t)diag [ ĝ(wt)ĝ(wt) >] . (6)

heading: 3 INDIVIDUAL CONVERGENCE OF HB
text: To solve the constrained problem (1), HB can be naturally reformulated as
wt+1 = PQ[wt − αtg(wt) + βt(wt −wt−1)]. (7) We first prove a key lemma, which extends (5) to the constrained and time-varying cases. Lemma 1. (Dimitri P. et al., 2003) For w ∈ Rd,w0 ∈ Q,
〈w −w0,u−w0〉 ≤ 0, for all u ∈ Q if and only if w0 = P (w). Lemma 2. Let {wt}∞t=1 be generated by HB (7). Let
pt = t(wt −wt−1), βt = t
t+ 2 , αt =
α
(t+ 2) √ t .

heading: Then HB (7) can be reformulated as
text: wt+1 + pt+1 = PQ[wt + pt − α√ t g(wt)]. (8)
Proof. The projection operation can be rewritten as an optimization problem (Duchi, 2018), i.e., wt+1 = PQ[wt − αtg(wt) + βt(wt −wt−1)] is equivalent to
wt+1 = arg min w∈Q {αt〈g(wt),w〉+
1 2 ‖w −wt − βt(wt −wt−1)‖2}. (9)
Then, ∀w ∈ Q, we have 〈wt+1 −wt − βt(wt −wt−1) + αtg(wt),wt+1 −w〉 ≤ 0.
This is 〈wt+1 + pt+1 − (wt + pt) +
α√ t g(wt),wt+1 −w〉 ≤ 0. (10)
Specifically, 〈wt+1 + pt+1 − (wt + pt) +
α√ t g(wt),wt+1 −wt〉 ≤ 0. (11)
From (10) and (11),
〈wt+1 + pt+1 − (wt + pt) + α√ t g(wt),wt+1 −w + (t+ 1)(wt+1 −wt)〉 ≤ 0.
i.e., 〈wt+1 + pt+1 − (wt + pt) +
α√ t g(wt),wt+1 + pt+1 −w〉 ≤ 0.
Using Lemma 1, Lemma 2 is proved.
Due to the non-expansive property of PQ (Dimitri P. et al., 2003), Lemma 2 implies that the convergence analysis for unconstrained problems can be applied to analyze the constrained problems.
Theorem 3. Assume that Q is bounded. Let {wt}∞t=1 be generated by HB (7). Set
βt = t
t+ 2 and αt =
α
(t+ 2) √ t .
Then f(wt)− f(w∗) ≤ O(
1√ t ).
Proof. According to Lemma 2,
‖w∗ − (wt+1 + pt+1)‖2 ≤ ‖w∗ − (wt + pt) + α√ t g(wt)‖2.
‖w∗ − (wt + pt) + α√ t g(wt)‖2
= ‖w∗ − (wt + pt)‖2 + ‖ α√ t g(wt)‖2 + 2〈 α√ t g(wt),w ∗ −wt〉+ 2〈 αt√ t g(wt),wt−1 −wt〉
Note
〈g(wt),w∗ −wt〉 ≤ f(w∗)− f(wt), 〈g(wt),wt−1 −wt〉 ≤ f(wt−1)− f(wt).
Then
(t+ 1)[f(wt)− f(w∗)]
≤ t[f(wt−1)− f(w∗)] + √ t
2α ‖w∗ − (wt + pt)‖2 −
√ t 2α ‖w∗ − (wt+1 + pt+1)‖2 + α 2 √ t ‖g(wt)‖2.
Summing this inequality from k = 1 to t, we obtain
(t+ 1)[f(wt)− f(w∗)]
≤ f(w0)− f(w∗) + t∑
k=1
α
2 √ k ‖g(wk)‖2 + t∑ k=1 [√k 2α (‖w∗ − (wk + pk)‖2 − ‖w∗ − (wk+1 + pk+1)‖2) ] .
Note t∑
k=1
1
2 √ k ‖g(wk)‖2 ≤
√ tM2.
and t∑
k=1
[√k 2 (‖w∗ − (wk + pk)‖2 − ‖w∗ − (wk+1 + pk+1)‖2) ] .
≤ 1 2 ‖w∗ − (w1 + p1)‖2 −
√ t
2 ‖w∗ − (wt+1 + pt+1)‖2 + t∑ k=2 ( √ k 2 − √ k − 1 2 )‖w∗ − (wk + pk)‖2.
Since Q is a bounded set, there exists a positive number M0 > 0 such that
‖w∗ − (wt+1 + pt+1)‖2 ≤M0,∀t ≥ 0.
Therefore
(t+ 1)[f(wt)− f(w∗)] ≤ f(w0)− f(w∗) + α √ tM2 +
√ t 2α M0.
This completes the proof of Theorem 3.
It is necessary to give some remarks about Theorem 3.
• In nonsmooth convex cases, Theorem 3 shows that the individual convergence rate of SGD can be accelerated from O( log t√
t ) to O( 1√ t ) via the HB momentum. The proof here clarifies
how the HB-type momentum wt −wt−1 and its time-varying weight βt help us to derive the optimal individual convergence.
• The convergence analysis in Theorem 3 is obviously different from the regret analysis in all the available papers, this is because the connection between f(wt) − f(w∗) and f(wt−1) − f(w∗) should be established here. It can be seen that seeking an optimal individual convergence is more difficult than the analysis of averaging convergence in many papers such as (Zinkevich, 2003) and (Yang et al., 2016). • We can get a stochastic HB by replacing the subgradient g(wt) in (7) with its unbiased
estimation ĝ(wt). Such substitution will not influence our convergence analysis. This means that we can get E[f(wt)− f(w∗)] ≤ O( 1√t ) under the same assumptions.
If βt remains a constant, we can get the averaging convergence rate, in which the proof of the first part is similar to Lemma 2 and that of the second part is similar to online PSG (Zinkevich, 2003). Theorem 4. Assume that Q is bounded and 0 ≤ βt ≡ β < 1. Let {wt}∞t=1 be generated by HB (7). Set
pt = β
1− β (wt −wt−1) and αt = α√ t .
Then we have
wt+1 + pt+1 = PQ[wt + pt − αt
1− β g(wt)], f(
1
t t∑ k=1 wk)− f(w∗) ≤ O( 1√ t ).
If Q is not bounded, the boundness of sequence ‖w∗ − (wt+1 + pt+1)‖ can not be ensured, which may lead to the failure of Theorem 4. Fortunately, like that in (Yang et al., 2016), E[f( 1T ∑T k=1 wk)− f(w∗)] ≤ O( 1√ T ) still holds, but we need to set αt ≡ α√T .

heading: 4 EXTENSION TO ADAPTIVE CASES
text: It is easy to find that HB (8) is in fact a gradient-based algorithm with predetermined step size α√ t . Thus its adaptive variant with EMA can be naturally formulated as
wt+1 = PQ[wt − αβ1t t √ t V − 12 t ĝ(wt) + β1t(wt −wt−1)]. (12)
where β1t = t
t+ 2 , Vt = β2tVt−1 + (1− β2t)diag
[ ĝ(wt)ĝ(wt) >] . The detailed steps of the adaptive HB are shown in Algorithm 1.
Algorithm 1 Adaptive HB Input: momentum parameters β1t, β2t, constant δ > 0, the total number of iterations T
1: Initialize w0 = 0, V0 = 0d×d 2: repeat 3: ĝt(wt) = ∇ft(wt), 4: Vt = β2tVt−1 + (1− β2t)diag(ĝt(wt)ĝt(wt)>), 5: V̂t = V 1 2 t + δ√ t Id, 6: wt+1 = PQ[wt − αβ1t
t √ t V̂t −1 ĝ(wt) + β1t(wt −wt−1)],
7: until t = T Output: wT
Theorem 5. Assume that Q is a bounded set. Let {wt}∞t=1 be generated by the adaptive HB (Algorithm 1). Denote pt = t(wt −wt−1). Suppose that β1t = tt+2 and 1− 1 t ≤ β2t ≤ 1− γ t for some 0<γ ≤ 1. Then wt+1 + pt+1 = PQ[wt + pt −
α√ t V̂t −1 ĝ(wt)] (13)
E[f(wt)− f(w∗)] ≤ O( 1√ t ).
The proof of (13) is identical to that of Lemma 2. It is easy to find that (13) is an adaptive variant of (8). This implies that the proof of the second part is similar to that of AdaGrad (Mukkamala & Hein, 2017). When 0 ≤ β1t ≡ β < 1, the adaptive variant of HB (7) is
wt+1 = PQ[wt − α√ t V − 12 t ĝ(wt) + β(wt −wt−1)]. (14)
where Vt = β2tVt−1 + (1− β2t)diag(ĝ(wt)ĝ(wt)>).
Similar to the proof of Theorem 5, we can get the following averaging convergence.
Theorem 6. Assume that Q is bounded and 0 ≤ β1t ≡ β < 1 in Algorithm 1. Let {wt}∞t=1 be generated by the adaptive HB (Algorithm 1). Suppose that 1− 1t ≤ β2t ≤ 1− γ t for some 0<γ ≤ 1. Denote pt = β1−β (wt −wt−1). Then
wt+1 + pt+1 = PQ[wt + pt − α (1− β) √ t V̂t −1 ĝ(wt)], E[f( 1 t t∑ k=1 wk)− f(w∗)] ≤ O( 1√ t ).
It is necessary to give some remarks about Theorem 5 and Theorem 6.
• The adaptive HB is usually used with a constant β1t in deep learning. However, according to Theorem 6, the constant β1t only guarantees the optimal data-dependent averaging convergence. The convergence property of the last iterate still remains unknown.
• In order to assure the optimal individual convergence, according to Theorem 5, β1t has to be time-varying. β1t = tt+2 can explain why we usually restrict β1t → 1 in practice. It also offers a new schedule about the selection of momentum parameters in deep learning.

heading: 5 EXPERIMENTS
text: In this section, we present some empirical results. The first two experiments are to validate the correctness of our convergence analysis and investigate the performance of the suggested parameters schedule. For fair comparison, we independently repeat the experiments five times and report the averaged results. The last experiment (Appendix A.4) is to show the effective acceleration of HB over GD in terms of the individual convergence.

heading: 5.1 EXPERIMENTS ON OPTIMIZING GENERAL CONVEX FUNCTIONS
text: This experiment is to optimize hinge loss with the l1-ball constraints. Let τ denotes the radius of the l1-ball. For implementation of the l1 projection operation, we use SLEP package1.
min f(w), s.t. w ∈ {w : ‖w‖1 ≤ τ}. (15)
Datasets: A9a, W8a, Covtype, Ijcnn1, Rcv1, Realsim (available at LibSVM2 website).
Algorithms: PSG (αt = α√t ), HB (αt = α (t+2) √ t , βt = tt+2 ), NAG (Tao et al., 2020a) and adaptive HB (12) (β1t = tt+2 ).
The relative function value f(wt) − f(w∗) v.s. epoch is illustrated in Figure 1. As expected, the individual convergence of the adaptive HB has almost the same behavior as the averaging output of PSG, and the individual output of HB and NAG. Since the three stochastic methods have the optimal convergence for general convex problems, we conclude that the stochastic adaptive HB attains the optimal individual convergence.
1http://yelabs.net/software/SLEP/ 2http://www.csie.ntu.edu.tw/˜cjlin/libsvmtools/datasets/

heading: 5.2 TRAINING DEEP NEURAL NETWORKS
text: These experiments are conducted on 4-layer CNN and ResNet-18 using a server with 2 NVIDIA 2080Ti GPUs.
Datasets: MNIST (60000 training samples, 10000 test samples), CIFAR10 (50000 training samples, 10000 test samples), and CIFAR100 (50000 training samples, 10000 test samples).
Algorithms: Adam (α, β1t ≡ 0.9, β2t ≡ 0.999, = 10−8) (Kingma & Ba, 2014), SGD (αt ≡ α), SGD-momentum (αt ≡ α, βt ≡ 0.9), AdaGrad (αt ≡ α) (Duchi et al., 2011), RMSprop (αt ≡ α, β2t ≡ 0.9, = 10−8) (Tieleman & Hinton, 2012). For our adaptive HB, γ = 0.1 and δ = 10−8. Different from the existing methods, we set β1t = tt+2 and β2t = 1 − γ t in Algorithm 1. Within each epoch, β1t and β2t remain unchanged.
Note that all methods have only one adjustable parameter α, we choose α from the set of {0.1, 0.01, 0.001, 0.0001} for all experiments. Following (Mukkamala & Hein, 2017) and (Wang et al., 2020), we design a simple 4-layer CNN architecture that consists two convolutional layers (32 filters of size 3 × 3), one max-pooling layer (2 × 2 window and 0.25 dropout) and one fully connected layer (128 hidden units and 0.5 dropout). We also use weight decay and batch normalization to reduce over-fitting. The optimal rate is always chosen for each algorithm separately so that one achieves either best training objective or best test performance after a fixed number of epochs.
The loss function is the cross-entropy. The training loss results are illustrated in Figure 2 and 4, and the test accuracy results are presented in Figure 3 and 5. As can be seen, the adaptive HB achieves the improved training loss. Moreover, this improvement also leads to good performance on test accuracy. The experimental results show that our suggested schedule about the momentum parameters could gain improved practical performance even in deep learning tasks.

heading: 6 CONCLUSION
text: In this paper, we prove that the adaptive HB method attains an optimal data-dependent individual convergence rate in the constrained convex cases, which bridges a theory-practice gap in using momentum methods to train the deep neural networks as well as optimize the convex functions. Our new analysis not only clarifies how the HB momentum and its time-varying weight β1t = tt+2 help us to achieve the acceleration but also gives valuable hints how its momentum parameters should be scheduled in deep learning. Empirical results on optimizing convex functions validate the
correctness of our convergence analysis and several typical deep learning experiments demonstrate the improved performance of the adaptive HB.

heading: 7 ACKNOWLEDGEMENTS
text: This work was supported in part by National Natural Science Foundation of China under Grants (62076252, 61673394, 61976213) and in part by Beijing Advanced Discipline Fund.

heading: A SUPPLEMENTARY MATERIAL

heading: A.1 PROOF FOR THEOREM 4
text: Let {wt}∞t=1 be generated by HB (7). Set
pt = β
1− β (wt −wt−1) and αt = α√ t .
Then, ∀w ∈ Q, according to Lemma 1, we have 〈wt+1 −wt − β(wt −wt−1) + αtg(wt),wt+1 −w〉 ≤ 0.
This is 〈 1 1− β (wt+1 −wt)− pt + αt 1− β g(wt),wt+1 −w〉 ≤ 0.
i.e., 〈wt+1 + pt+1 − (wt + pt) +
αt 1− β g(wt),wt+1 −w〉 ≤ 0 (16)
Specifically,
〈wt+1 + pt+1 − (wt + pt) + αt
1− β g(wt), β(wt+1 −wt) 1− β 〉 ≤ 0 (17)
From (16) and (17),
〈wt+1 + pt+1 − (wt + pt) + αt
1− β g(wt),wt+1 + pt+1 −w〉 ≤ 0.
Using Lemma 1, we have
wt+1 + pt+1 = PQ[wt + pt − αt
1− β g(wt)].
Then
‖w∗ − (wt+1 + pt+1)‖2 ≤‖w∗ − (wt + pt) + αt
1− β g(wt)‖2
=‖w∗ − (wt + pt)‖2 + ‖ αt
1− β g(wt)‖2 + 2〈 αt 1− β g(wt),w ∗ −wt〉
+2〈 αtβ (1− β)2 g(wt),wt−1 −wt〉
Note
〈g(wt),w∗ −wt〉 ≤ f(w∗)− f(wt), 〈g(wt),wt−1 −wt〉 ≤ f(wt−1)− f(wt). Then
‖w∗ − (wt+1 + pt+1)‖2
≤‖w∗ − (wt + pt)‖2 + α2t
(1− β)2 ‖g(wt)‖2
+ 2αt 1− β [f(w∗)− f(wt)] + 2αtβ (1− β)2 [f(wt−1)− f(wt)].
Rearrange the inequality, we have
2αt 1− β [f(wt)− f(w∗)] ≤ 2αtβ (1− β)2 [f(wt−1)− f(wt)] + ‖w∗ − (wt + pt)‖2
−‖w∗ − (wt+1 + pt+1)‖2 + α2t
(1− β)2 ‖g(wt)‖2.
i.e.,
f(wt)− f(w∗) ≤ β
1− β [f(wt−1)− f(wt)] + 1− β 2αt [‖w∗ − (wt + pt)‖2
− ‖w∗ − (wt+1 + pt+1)‖2] + αt
2(1− β) ‖g(wt)‖2.
Summing this inequality from k = 1 to t, we obtain t∑
k=1
[f(wk)− f(w∗)]
≤ β 1− β [f(w0)− f(wt)] + 1− β 2α1 ‖w∗ − (w1 + p1)‖2
−1− β 2αt
‖w∗ − (wt+1 + pt+1)‖2 + t∑
k=1
αk 2(1− β) ‖g(wk)‖2
+ t∑ k=2 ‖w∗ − (wk + pk)‖2( 1− β 2αk − 1− β 2αk−1 ).
i.e.,
t∑ k=1 [f(wk)− f(w∗)]
≤ β 1− β [f(w0)− f(wt)] + 1− β 2α ‖w∗ − (w1 + p1)‖2
+ t∑ k=2 ‖w∗ − (wk + pk)‖2( (1− β) √ k 2α − (1− β) √ k − 1 2α )
+ t∑ k=1
α
2(1− β) √ k ‖g(wk)‖2.
(18)
Note t∑
k=1
1
2 √ k ‖g(wk)‖2 ≤
√ tM2. (19)
and since Q is a bounded set, there exists a positive number M0 > 0 such that ‖w∗ − (wt+1 + pt+1)‖2 ≤M0,∀t ≥ 0. (20) From (18)(19)(20) we have,
t∑ k=1 [f(wk)− f(w∗)] ≤ β 1− β [f(w0)− f(wt)] + (1− β) √ tM0 2α + α √ tM2 1− β .
By convexity of f(w), we obtain
f( 1
t t∑ k=1 wk)− f(w∗) ≤ β (1− β)t [f(w0)− f(wt)] + (1− β)M0 2α √ t + αM2 (1− β) √ t .
This completes the proof of Theorem 4.

heading: A.2 PROOF FOR THEOREM 5
text: Notation. For a positive definite matrix H ∈ Rd×d, the weighted `2-norm is defined by ‖x‖2H = x>Hx. The H-weighted projection PHQ (x) of x onto Q is defined by P H Q (x) = argminy∈Q ‖y− x‖2H . We use g(wk) to denote the subgradient of fk(·) at wk. For the diagonal matrix sequence {Mk}tk=1, we use mk,i to denote the i-th element in the diagonal of Mk. We introduce the notation, g1:k,i = (g1,i, g2,i, .., gk,i) >, where gk,i is the i-th element of g(wk).
Lemma 7. (Mukkamala & Hein, 2017) Suppose that 1− 1t ≤ β2t ≤ 1− γ t for some 0<γ ≤ 1, and t ≥ 1, then d∑ i=1 t∑ k=1 g2k,i√ kvk,i + δ ≤ d∑ i=1 2(2− γ) γ ( √ tvt,i + δ).
Proof for Theorem 5. Without loss of generality, we only prove Theorem 5 in the full gradient setting. It can be extended to stochastic cases using the regular technique in (Rakhlin et al., 2011).
Note that the projection operation can be rewritten as an optimization problem (Duchi, 2018), i.e., wt+1 = PQ[wt − αtV̂ −1t g(wt) + β1t(wt −wt−1)] is equivalent to
wt+1 = arg min w∈Q {αt〈V̂ −1t g(wt),w〉+
1 2 ‖w −wt − β1t(wt −wt−1)‖2}. (21)
Then, ∀u ∈ Q, we have
〈wt+1 −wt − βt(wt −wt−1) + αtV̂ −1t g(wt),wt+1 −w〉 ≤ 0.
This is 〈wt+1 + pt+1 − (wt + pt) +
α√ t V̂ −1t g(wt),wt+1 −w〉 ≤ 0. (22)
Specifically, 〈wt+1 + pt+1 − (wt + pt) +
α√ t V̂ −1t g(wt),wt+1 −wt〉 ≤ 0. (23)
From (22) and (23),
〈wt+1 + pt+1 − (wt + pt) + α√ t V̂ −1t g(wt),wt+1 −wt + (t+ 1)(wt+1 −wt)〉 ≤ 0.
i.e., 〈wt+1 + pt+1 − (wt + pt) +
α√ t V̂ −1t g(wt),wt+1 + pt+1 −wt〉 ≤ 0.
Using Lemma 1, we have
wt+1 + pt+1 = P V̂t Q [wt + pt − α√ t V̂ −1t g(wt)].
Then
‖w∗ − (wt+1 + pt+1)‖2V̂t ≤ ‖w ∗ − (wt + pt) + α√ t V̂ −1t g(wt)‖2V̂t
= ‖w∗ − (wt + pt)‖2V̂t + ‖ α√ t g(wt)‖2V̂t
+ 2〈 α√ t g(wt),w ∗ −wt〉+ 2〈 αt√ t g(wt),wt−1 −wt〉.
Note
〈g(wt),w∗ −wt〉 ≤ f(w∗)− f(wt), 〈g(wt),wt−1 −wt〉 ≤ f(wt−1)− f(wt).
Then
(t+ 1)[f(wt)− f(w∗)] ≤t[f(wt−1)− f(w∗)] + √ t
2α ‖w∗ − (wt + pt)‖2V̂t
− √ t
2α ‖w∗ − (wt+1 + pt+1)‖2V̂t +
α
2 √ t ‖g(wt)‖2V̂ −1t .
Summing this inequality from k = 1 to t, we obtain
(t+ 1)[f(wt)− f(w∗)] ≤ f(w0)− f(w∗) + t∑
k=1
α
2 √ k ‖g(wk)‖2V̂ −1k
+ t∑ k=1 [√k 2α (‖w∗ − (wk + pk)‖2V̂k − ‖w ∗ − (wk+1 + pk+1)‖2V̂k) ] .
Using Lemma 7, we have t∑
k=1
α
2 √ k ‖g(wk)‖2V̂ −1k ≤ d∑ i=1 α(2− γ) γ ( √ tvt,i + δ).
Note t∑
k=1
[√k 2α (‖w∗ − (wk + pk)‖2V̂k − ‖w ∗ − (wk+1 + pk+1)‖2V̂k) ] =
d∑ i=1 v̂1,i 2α (w∗i − (w1,i + p1,i))2 − d∑ i=1 √ tv̂t,i 2α (w∗i − (wt+1,i + pt+1,i))2
+ d∑ i=1 t∑ k=2 1 2α ( √ kv̂k,i − √ k − 1v̂k−1,i)(w∗i − (wk,i + pk,i))2.
(24)
Since Q is a bounded set, there exists a positive number M1 > 0 such that
(w∗i − (wt+1,i + pt+1,i))2 ≤M1,∀t ≥ 0, i = 1, 2, ..., d.
and vk,i = β2kvk−1,i + (1− β2k)g2k,i as well as β2k ≥ 1− 1k which implies kβ2k ≥ k − 1, we get √ kv̂k,i = √ kvk,i + δ
= √ kβ2kvk−1,i + k(1− β2k)g2k,i + δ
≥ √ (k − 1)vk−1,i + δ = √ k − 1v̂k−1,i.
Therefore t∑
k=1
[√k 2α (‖w∗ − (wk + pk)‖2V̂k − ‖w ∗ − (wk+1 + pk+1)‖2V̂k) ] ≤
d∑ i=1 v̂1,i 2α M1 + d∑ i=1 t∑ k=2 1 2α ( √ kv̂k,i − √ k − 1v̂k−1,i)M1
= d∑ i=1 v̂1,iM1 2α + d∑ i=1 √ tv̂t,iM1 2α − d∑ i=1 v̂1,iM1 2α
= M1 2α d∑ i=1 ( √ tvt,i + δ).
(25)
Since √ tvt,i = ‖g1:t,i‖, therefore
(t+ 1)[f(wt)− f(w∗)] ≤f(w0)− f(w∗) + M1 2α d∑ i=1 ( √ tvt,i + δ) + d∑ i=1 α(2− γ) γ ( √ tvt,i + δ)
=f(w0)− f(w∗) + ( M1 2α + α(2− γ) γ ) d∑ i=1 (‖g1:t,i‖+ δ).
This proves
f(wt)− f(w∗) ≤ O( 1√ t ).

heading: A.3 PROOF FOR THEOREM 6
text: Let {wt}∞t=1 be generated by the adaptive HB (Algorithm 1). Set
pt = β
1− β (wt −wt−1) and αt = α√ t .
Then, ∀u ∈ Q, according to Lemma 1, we have
〈wt+1 −wt − β(wt −wt−1) + αtV̂ −1t g(wt),wt+1 −w〉 ≤ 0. This is
〈 1 1− β (wt+1 −wt)− pt + αtV̂
−1 t
1− β g(wt),wt+1 −w〉 ≤ 0.
i.e.,
〈wt+1 + pt+1 − (wt + pt) + αtV̂
−1 t
1− β g(wt),wt+1 −w〉 ≤ 0 (26)
Specifically,
〈wt+1 + pt+1 − (wt + pt) + αtV̂
−1 t
1− β g(wt), β(wt+1 −wt) 1− β 〉 ≤ 0 (27)
From (26) and (27),
〈wt+1 + pt+1 − (wt + pt) + αtV̂
−1 t
1− β g(wt),wt+1 + pt+1 −w〉 ≤ 0.
Using Lemma 1, we have
wt+1 + pt+1 = P V̂t Q [wt + pt −
αtV̂ −1 t
1− β g(wt)].
According to Lemma 2,
‖w∗ − (wt+1 + pt+1)‖2V̂t
≤‖w∗ − (wt + pt) + αtV̂
−1 t
1− β g(wt)‖2V̂t
=‖w∗ − (wt + pt)‖2V̂t + ‖ αt 1− β g(wt)‖2V̂t
+2〈 αt 1− β g(wt),w ∗ −wt〉+ 2〈 αtβ (1− β)2 g(wt),wt−1 −wt〉
Note
〈g(wt),w∗ −wt〉 ≤ f(w∗)− f(wt), 〈g(wt),wt−1 −wt〉 ≤ f(wt−1)− f(wt). Then
‖w∗ − (wt+1 + pt+1)‖2V̂t
≤‖w∗ − (wt + pt)‖2V̂t + α2t (1− β)2 ‖g(wt)‖2V̂ −1t
+ 2αt 1− β [f(w∗)− f(wt)] + 2αtβ (1− β)2 [f(wt−1)− f(wt)].
Rearrange the inequality, we have 2αt 1− β [f(wt)− f(w∗)]
≤ 2αtβ (1− β)2 [f(wt−1)− f(wt)] + ‖w∗ − (wt + pt)‖2V̂t
−‖w∗ − (wt+1 + pt+1)‖2V̂t + α2t (1− β)2 ‖g(wt)‖2V̂ −1t .
i.e., f(wt)− f(w∗)
≤ β 1− β [f(wt−1)− f(wt)] + 1− β 2αt [‖w∗ − (wt + pt)‖2V̂t −‖w∗ − (wt+1 + pt+1)‖2V̂t ] + αt 2(1− β) ‖g(wt)‖2V̂ −1t .
Summing this inequality from k = 1 to t, we obtain
t∑ k=1 [f(wk)− f(w∗)]
≤ β 1− β [f(w0)− f(wt)] + 1− β 2α1 ‖w∗ − (w1 + p1)‖2V̂1
−1− β 2αt ‖w∗ − (wt+1 + pt+1)‖2V̂t + t∑
k=1
αk 2(1− β) ‖g(wk)‖2V̂ −1k
+ d∑ i=1 t∑ k=2 (w∗i − (wk,i + pk,i))2( (1− β)v̂k,i 2αk − (1− β)v̂k−1,i 2αk−1 ).
i.e., t∑
k=1
[f(wk)− f(w∗)]
≤ β 1− β [f(w0)− f(wt)] + 1− β 2α ‖w∗ − (w1 + p1)‖2V̂1
+ d∑ i=1 t∑ k=2 (w∗i − (wk,i + pk,i))2 1− β 2α ( √ kv̂k,i − √ k − 1v̂k−1,i)
+ t∑ k=1
α
2(1− β) √ k ‖g(wk)‖2V̂ −1k .
(28)
Using Lemma 7, we have
t∑ k=1
α
2 √ k(1− β) ‖g(wk)‖2V̂ −1k ≤ d∑ i=1 α(2− γ) γ(1− β) ( √ tvt,i + δ) = α(2− γ) γ(1− β) d∑ i=1 (‖g1:t,i‖+ δ). (29)
and since Q is a bounded set, there exists a positive number M0 > 0 such that
‖w∗ − (wt+1 + pt+1)‖2 ≤M0,∀t ≥ 0. (30)
From (28)(29)(30) we have,
t∑ k=1 [f(wk)− f(w∗)]
≤ β 1− β [f(w0)− f(wt)] + d∑ i=1 (1− β)v̂1,iM0 2α + α(2− γ) γ(1− β) d∑ i=1 (‖g1:t,i‖+ δ)
+ d∑ i=1 (1− β)v̂t,i √ tM0 2α − d∑ i=1 (1− β)v̂1,iM0 2α .
i.e., t∑ k=1 [f(wk)− f(w∗)] ≤ β 1− β [f(w0)−f(wt)]+ α(2− γ) γ(1− β) d∑ i=1 (‖g1:t,i‖+δ)+ (1− β)M0 2α d∑ i=1 (‖g1:t,i‖+δ).
By convexity of f(wk), we obtain
f( 1
t t∑ k=1 wk)−f(w∗) ≤ β (1− β)t [f(w0)−f(wt)]+ α(2− γ) γ(1− β)t d∑ i=1 (‖g1:t,i‖+δ)+ (1− β)M0 2αt d∑ i=1 (‖g1:t,i‖+δ).
This completes the proof of Theorem 6.

heading: A.4 EXPERIMENTS ON OPTIMIZING A SYNTHETIC CONVEX FUNCTION
text: A constrained convex optimization problem was constructed in (Harvey et al., 2019) to show that the optimal individual convergence rate of SGD is O( log t√
t ). We will use example to illustrate the
acceleration of HB.
Let Q be unit ball in RT . For i ∈ [T ] and c ≥ 1, define the positive scalar parameters
ai = 1
8c(T − i+ 1) bi =
√ i
2c √ T
Define f : Q→ R and hi ∈ RT for i ∈ [T + 1] by
f(w) = max i∈[T ]
h>i w where hi,j =  aj , 1 ≤ j < i −bj , i = j < T
0, i < j ≤ T
Obviously, the minimum value of f on the unit ball is non-positive because f(0) = 0. It can be proved f(wT ) ≥ log T32c√T . Set c = 2, the function value f(wt) v.s. iteration is illustrated in Figure 6, where the step size of GD is c√
t and the parameters of the constrained HB (7) (α = 8) and AdaHB
(12) (α = 0.08, γ = 0.9, δ = 10−8) are selected according to Theorem 3 and Theorem 5. As expected, the individual convergence of HB is much faster than that of PSG. We thus conclude that HB is an effective acceleration of GD in terms of the individual convergence.

Review summaries:
The authors investigate the convergence of the projected Heavy-ball method (and an adaptive variant) for convex problems with convex constraints. The authors prove 4 results: 2 individual (last iterate) convergence rates and 2 rates using averaging. Notably, in their proofs they require an increasing (from 1/2 to 1) momentum parameter and a decreasing stepsize. Finally, the authors present some experimental results.