Accelerated Incremental Gradient Descent using Momentum Acceleration with Scaling Factor

1 Introduction

Recently, stochastic/incremental first-order methods have received extensive attention due to their low per-iteration cost and the ability to handle large-scale problems including unconstrained/constrained composite convex minimization [Allen-Zhu, 2018; Liu et al., 2017; Shang et al., 2018a]. In particular, the research on stochastic variance reduced gradient descent methods (e.g., SAG [Roux et al., 2012], SVRG [Johnson and Zhang, 2013], SDCA [Shalev-Shwartz and Zhang, 2013], SAGA [Defazio et al., 2014a]) and their proximal variants (e.g., Prox-SVRG [Xiao and Zhang, 2014] and VR-SGD [Shang et al., 2018a]) has made exciting progress, e.g., linear convergence for strongly convex (SC) problems. These methods use past gradients to progressively reduce the variance of stochastic gradient estimators, which leads to a revolution in the area of first-order optimization [Shang et al., 2019]. In this paper, we consider the following composite convex minimization problem in many problems of machine learning, statistics, and operations research, such as regularized empirical risk minimization (ERM).

∗Corresponding author

where f(x)= 1n ∑ni=1 fi(x) is the finite average of n convex component functions fi(x) :Rd →R, and h(x) is a simple but possibly non-smooth convex function.

For solving SC problems of the formulation (1), the oracle complexity (total number of gradient evaluations to find an ϵ-suboptimal solution) of the stochastic variance reduction methods mentioned above (including SVRG and SAGA) is O((n+κ) log(1/ϵ)), where κ is the condition number of Problem (1), and n is the number of samples, while the oracle complexity of accelerated deterministic methods, e.g., FISTA [Beck and Teboulle, 2009], is O(n√κ log(1/ϵ)). Obviously, the oracle complexities show that those stochastic variance reduction methods always converge faster than deterministic methods (including their accelerated variants) as long as κ ≤O(n2 ). However, there is still a gap between the complexity of those stochastic variance reduction methods and the upper bound provided in [Woodworth and Srebro, 2016].

More recently, there is a surge of interests in accelerating stochastic variance reduction gradient optimization. The acceleration techniques in accelerated methods mainly include the Nesterov's acceleration technique [Nitanda, 2014; Frostig et al., 2015; Lin et al., 2015; Murata and Suzuki, 2017], the choice of growing epoch length [Mahdavi et al., 2013; Allen-Zhu and Yuan, 2016; Shang et al., 2017], and the momentum acceleration tricks [Shang et al., 2018b; Allen-Zhu, 2018; Zhou et al., 2018; Hien et al., 2019]. [Lin et al., 2015] proposed a Catalyst framework for accelerating some stochastic variance reduction algorithms (including SVRG and SAGA) and proved that their accelerated variants achieve an oracle complexity of O((n+ √nκ) log(κ) log(1/ϵ)) for SC problems. In particular, as the accelerated variants of SVRG, Katyusha [Allen-Zhu, 2018] and MiG [Zhou et al., 2018] achieve the best-known oracle complexity of O((n+ √nκ) log(1/ϵ)) for SC problems, which is identical to the upper bound in [Woodworth and Srebro, 2016].

As an accelerated variant of SAGA, point-SAGA [Defazio, 2016] requires a proximal operator oracle of each fi and can attain the same oracle complexity as Katyusha. However, the proximal operator oracle for point-SAGA may not be efficiently calculated in practice. As we all know that a large amount of work has been done for accelerating SVRG, while the notable incremental gradient method, SAGA, does not have a direct accelerated variant until recently, except for SSNM [Zhou et al., 2019], which can also obtain the best-known oracle complexity. However, the memory complexity of SSNM is always O(nd). Therefore, we will propose a simple and direct accelerated variant for SAGA.

Contributions: Existing accelerated methods mentioned above can attain the theoretical upper complexity bounds provided in [Woodworth and Srebro, 2016]. We ask the following question in this paper: Using only gradient information, can we further improve the convergence rates of those methods such as Katyusha [Allen-Zhu, 2018], point-SAGA [Defazio, 2016] and SSNM [Zhou et al., 2019]? In this paper, we propose a novel accelerated incremental gradient descent (AIGD) algorithm to push towards the convergence rates.

• We design a general momentum acceleration scheme for the direct acceleration of SAGA, in which we introduce a novel momentum to replace the Nesterov's momentum and Katyusha momentum used in [Allen-Zhu, 2018].
• We prove that AIGD achieves a linear convergence rate and the oracle complexity of O((n+√nκ) log(1/ϵ)) for strongly convex problems, which is identical to the convergence results of existing accelerated algorithms such as Katyusha [Allen-Zhu, 2018], point-SAGA [Defazio, 2016], and SSNM [Zhou et al., 2019].
• In particular, our convergence results also show that AIGD can slightly improve the convergence rates of Katyusha, point-SAGA and SSNM for the case of n≥ L/µ, as shown in Table 1. It means that this study can partly answer the above-mentioned question.
• We also discuss some subtle differences between AIGD and existing accelerated incremental algorithms such as point-SAGA and SSNM, which imply that AIGD is more suitable for solving various large-scale machine learning problems.
• AIGD can also be extended to the non-convex setting. Our experimental results further verify that AIGD is usually faster than the state-of-the-art accelerated stochastic/incremental methods.

2 Preliminaries and Notations

Throughout this paper, the norm ‖·‖ is the standard Euclidean norm. We denote by ∇f(x) the full gradient of f(x) if it is differentiable, or ∂f(x) a sub-gradient of f(x) if f(x) is only Lipschitz continuous. We mostly focus on the case of Problem (1) when each fi(x) is L-smooth.

Assumption 1. Each fi(·) is L-smooth, i.e., there exists a constant L>0 such that for all x, y∈Rd ,

Assumption 2. The function h(·) is µ-strongly convex, i.e., there exists µ>0 such that for all x, y∈Rd ,

where ∂h(x) is the subdifferential of h(·) at x.

Recently, there emerges a stream of studies on stochastic variance reduced methods, such as [Zhang et al., 2013; Johnson and Zhang, 2013; Xiao and Zhang, 2014; Defazio et

Algorithm 1 AIGD for Strongly Convex Objectives

• 2: Pick ik uniformly at random from {1, 2, . . . , n}, and update yk by (2);
• 3: ˜∇fik(yk)=∇fik(yk)-∇fik(φk-1 ik )+ 1n ∑nj=1∇fj (φk-1 j );
• 4: Update zk by solving (4);
• 5: Update xk by (3);
• 6: Take φkik = xk, calculate and store ∇fik(φkik ) in the table;
• 7: end for

Output: x = xK.

al., 2014a; Shang et al., 2019]. The two most popular choices for stochastic gradient estimators are the SVRG estimator in [Zhang et al., 2013; Johnson and Zhang, 2013] and the SAGA estimator in [Defazio et al., 2014a]. The main update steps of SAGA [Defazio et al., 2014a] are

where η>0 is a learning rate, and proxhγ (·) is a proximal operator. More details can be found in [Defazio et al., 2014a]. We start with an initial vector x0 ∈Rd and the known derivatives ∇fj(φ0j ) with φ0j =x0 for each j ∈{1, 2, . . . , n}. These derivatives are stored in a table data structure of length n. At the k-th iteration, we choose randomly sample i, and take φki =xk and store the gradient ∇fi(φki ) in the table, and all other entries in the table remain unchanged.

3 Accelerated Incremental Gradient Descent

In section, we propose a simple and direct accelerated incremental gradient descent (AIGD) algorithm.

3.1 ANovel Momentumwith Scaling Factor

In this subsection, we design a novel momentum acceleration scheme for incremental gradient descent optimization, as shown in Algorithm 1. In particular, a scaling factor is introduced into our momentum acceleration scheme to improve the convergence rate. For solving Problem (1), our main update rules with momentum acceleration are

where θ is the momentum parameter, β is a scaling factor using which we can improve the convergence rate, and zk is the solution to the following problem:

where ˜∇fi(yk) :=∇fi(yk)-∇fi(φk-1 i )+ 1n ∑nj=1∇fj (φk-1 j ). Note that in our algorithm, we also use h(z/β) instead of h(z) used in SNNM. The parameters θ, η and β are set in Algorithm 1.

Remark 1. There is an important difference between our method and SSNM. That is, we use yk =θzk-1+(1-βθ)xk-1 instead of yk =θzk-1+(1-θ)φk-1 i used in SSNM. Our setting has the following two advantages: We only need to store the gradient table, while SSNM requires to store both the "point" and gradient tables. The second advantage is that our method removes the independent sample assumption used in SSNM. This means that our algorithm has a much weaker convergence condition than SSNM.

3.2 Comparison with SAGA, Point-SAGA and SSNM

• This paper uses F(x)-F(x⋆ ) as a convergence critical rule instead of the Lyapunov function used in SAGA, point-SAGA and SSNM, making our algorithm easier to extend for solving structure optimization problems, such as graph-guided fuzed Lasso [Kim et al., 2009] and generalized Lasso [Tibshirani and Taylor, 2011].
• In particular, we design a novel momentum scheme for accelerating SAGA instead of that in SSNM, such that our algorithm only needs to store the gradient table, while SSNM requires to store both the "point" and gradient tables. Therefore, the memory complexity of SSNM is always O(nd), as analyzed in [Zhou et al., 2019]. This is a disadvantage of SSNM when the objective is a linear model, e.g., linear logistic regression and ridge regression. In contrast, our algorithm only requires an O(d) memory complexity to simply store a scalar to represent the gradient of each component function.
• Our AIGD method is a direct accelerated incremental gradient method, while point-SAGA [Defazio, 2016] requires the proximal operator oracle of each component function. However, the proximal operator may not be efficiently computed in practice, which makes pointSAGA not suitable for many real-world problems.
• Moreover, the update rules of our algorithm are more elegant than Katyusha and MiG, both of which require a tricky weighted averaged scheme at the end of each inner loop. In particular, our algorithm can further improve the convergence rate of Katyusha, MiG, point-SAGA and SSNM for the case of n≥3κ/4, where κ=L/µ.

3.3 Extensions and Complexity Analysis

AIGD can be extended to non-smooth and Lipschitz continuous settings. By using adaptive regularization and smoothing techniques as in [Allen-Zhu and Hazan, 2016; Allen-Zhu, 2018], one can get a new and smooth function, which approximates the original function. That is, our AIGD method can attain at least the same complexity bounds as Katyusha, i.e., O(n log(1/ϵ)+L√n/(µϵ)) for strongly convex and Lipschitz continuous problems and O(n log(1/ϵ)+ √nL/ϵ) for non-strongly convex and Lipstize continuous problems.

Each iteration of AIGD computes the stochastic gradients ∇fik(yk) and ∇fik(φk-1 ik ), which is the same as existing stochastic methods such as SVRG [Johnson and Zhang, 2013] and Katyusha [Allen-Zhu, 2018]. In ERM problems, the loss function fi(·) takes the form fi(aTi x) for ai, where ai is a constant vector. With such a structure, we apply the widely used scheme as in [Roux et al., 2012; Defazio et al., 2014b] and only need to store the scalar aTikφk-1 ik for ∇fik(φk-1 ik ) rather than the vector in each iteration. Therefore, the storage cost of AIGD can be reduced from O(nd) to O(n).

4 Convergence Analysis

In this section, we analyze the convergence property of our AIGD algorithm. Before giving our main convergence result, we first give and prove the following lemmas.

Lemma 1 (Variance Upper Bound). Let x⋆ be the optimal solution of Problem (1). Suppose each fi(·) is convex and L-smooth, and let ˜∇k := ˜∇fik(yk), then we have

Proof. By the definition of ˜∇k and taking expectation over the random choice of ik, we have

where the first inequality follows from the facts that Eik[∇fik(φk-1 ik )] = 1n ∑nj=1∇fj (φk-1 j ), Eik[∇fik(yk)] = ∇f(yk), E[fik(φk-1 ik )] = 1n ∑nj=1 fj (φk-1 j ) and Eik[‖Eik [x]- x‖2 ] = Eik[‖x‖2 ] - ‖Eik [x]‖2 ≤ Eik[‖x‖2 ]; the third inequality holds due to Eik[∇fik(yk)] = ∇f(yk), Eik[∇fik(x⋆ )] = ∇f(x⋆ ), and ‖∇fj(x) - ∇fj(y)‖2 ≤ 2L[fj(y)-fj(x)+〈∇fj(x), x-y〉] for all x, y∈Rd .

Lemma 2. Suppose each fi(·) is convex and L-smooth, then we have

where E denotes the expectation with respect to all randomness, and Fj(·)=fj(·)+h(·).

Proof.

where ∂h(x⋆ ) is a subgradient of h(·) and ∂h(x⋆ )+∇f(x⋆ )= 0. Let F ′ j (x⋆ ) :=∇fj(x⋆ )+∂h(x⋆ ), and using the result in SAGA [Defazio et al., 2014a], we have

Let F ′ (x⋆ ) :=∇f(x⋆ )+∂h(x⋆ ), and we have F ′ (x⋆ ) = 0. Thus, the following result holds:

Therefore, we have

This completes the proof.

Theorem 1. Suppose h(·) is µ-strongly convex, and let {(xk, yk, zk)} be the sequence generated by Algorithm 1. Then Algorithm 1 has the following geometric convergence in expectation:

where

That is, Algorithm 1 achieves an ϵ-suboptimal solution using at most

Remark 2. Theorem 1 shows that AIGD achieves the oracle complexity of O((n+√nκ) log(1/ϵ)) for strongly convex problems, which is the same as the best-known result in [Defazio, 2016; Allen-Zhu, 2018], and also matches the upper complexity bound provided in [Woodworth and Srebro, 2016]. In particular, in the case of n ≥ 3L/4µ, the convergence rate of both SAGA in [Defazio et al., 2014a] and Finito [Defazio et al., 2014b] is (1- 1 2n )n ≈ exp(-1/2) = 0.606, and Katyusha [Allen-Zhu, 2018] can obtain the rate of 1/1.5 = 0.667, while the convergence rate of AIGD is (1- 3 4n )n ≈0.517. That is, our AIGD algorithm improves the best-known convergence rates.

Before giving the proof of Theorem 1, we first give the following property.

Property 1. Given any w1, w2, w3, w4 ∈ Rd , then we have

Proof of Theorem 1:

Proof. From the optimality condition of (4) in Algorithm 1 with respect to zk with η = 1 Lαθ , we have

θ〈 ˜∇k+β∂˜h(zk)+Lαθ(zk-zk-1), z-zk〉≥0, for any z ∈Rd , where ˜h(z) :=h(z/β), and ∂˜h(z) =∂h(z/β)/β. Using the above inequality with z=βx⋆ , then we have

where the second inequality holds due to the strong convexity of h(·); the last inequality holds due to the update rule xk = θzk+(1-βθ)xk-1, Property 1 and the convexity of h(·).

Since f(·) is L-smooth, and by the update rule yk = θzk-1 + (1-βθ)xk-1, xk = θzk + (1-βθ)xk-1, and taking expectation over the random choice of ik, we have

where the inequality a≤ follows from the Young's inequality 〈∇f(yk)- ˜∇k, xk -yk〉 ≤ 1 2La ‖∇f (yk) - ˜∇k‖2 + La 2 ‖xk - yk‖2 ; the inequality b≤ follows from the inequality (5), and Ak-1 := 1n ∑nj=1(fj (φk-1 j )-f(x⋆ )-〈∇fj(x⋆ ), φk-1 j -x⋆ 〉); the inequality c≤ holds due to the fact that

where the first equation follows from Lemma 1 and Eik [˜∇k] = ∇f(yk), and the last inequality holds due to the convexity of the function f(·) and βθ - 2/(α-1) > 0. By using Lemma 2 and the above analysis, we have

Using the result in SAGA [Defazio et al., 2014a], we have

Thus, 1n E[Gk-1] = Pk-(1- 1n )Pk-1, whereGk :=F(xk)- F(x⋆ ), Pk := 1n ∑nj=1Fj (φkj )-F (x∗ ) and E[Pk] = Σk-1 j=0 (1/n)k-j-1 1/nPj. Let Qk :=‖βx⋆ -zk‖2 , and using (6) and multiplying each side of the above inequality by γK-k , 0 < γ ≤ 1, and summing it over k=1, · · · ,K - 1, we have

Next, we consider the following two cases.

Thus, (7) is rewritten as follows:

where τ =3/2-1/(4n)+1/8, since P0=G0 and ‖x⋆-x0‖≤ 1µ G0, we have (1-βθ)G0+( 12- 1 2n+ 2 α-1 )P0+Lαθ2 2 Q0=τG0. Using (6) with k=K and α=5, we have

Since γ > (1-βθ) and 5Lθ2 2 < γθ(Lαθ+µ/β) 2 , and by using the above two inequalities, we have

That is,

Case 2: n ≥ 3L 4µ . We set γ = 1 1+3/(4n) = 1-3/(4n+3). With η= 3 4µn , α=8n/β and θ= 1 Lαη , and using the similar derivation, we have the following result

This completes the proof.

5 Experimental Results

In this section, we evaluate the performance of our algorithm for justifying our theoretical results.

We conducted many experiments of the strongly convex logistic regression problem on the two real-world data sets: Covtype (581,012 examples and 54 features) and a9a (32,562 examples and 123 features). Note that each feature vector of the two data sets was scaled down by the average Euclidean norm of the whole data set as in [Allen-Zhu and Yuan, 2016]. For SVRG and Katyusha, we set the epoch size m = 2n, as suggested in [Johnson and Zhang, 2013; Allen-Zhu, 2018]. Figs. 1 and 2 show how the objective gap (i.e., F(xK)-F(x⋆ )) of all these algorithms decreases for logistic regression with different regularization parameters λ=10-4 , 10-7 , 10-8 . All the results show that the accelerated methods (i.e., Katyusha, SSNM and our algorithm) usually perform much better than the non-accelerated methods, SVRG and SAGA, especially with relatively small regularization parameters, e.g., λ=10-7 . SAGA and Katyusha achieve similar performance for the well-conditioned case, while Katyusha is significantly faster than SAGA for the illconditioned case (i.e., the case of small regularization parameters). Our algorithm outperforms the other methods (including the accelerated algorithms, Katyusha and SSNM) in terms of running time. This further justifies the effectiveness of our momentum acceleration technique for accelerating SAGA.

6 Conclusions

In this paper, we proposed a novel accelerated incremental gradient descent algorithm with the proposed momentum acceleration technique. Unlike the existing accelerated algorithms such as point-SAGA and SSNM, our algorithm is a direct accelerated method, and requires significantly less memory than SSNM for various large-scale linear models. Moreover, we provided the convergence property of our algorithm for solving strongly convex problems, which shows that our algorithm attains the best-known oracle complexity, and an improved convergence rate for the case of n≥L/µ. Our algorithm can be extended to non-convex setting and constrained composite convex minimization setting, especially, the problems of the general ADMM form as in [Liu et al., 2017].

Acknowledgments

This work was supported by Project supported Foundation for the Innovative Research Groups of the National Natural Science Foundation of China (No. 61621005), the Major Research Plan of the National Natural Science Foundation of China (Nos. 91438201 and 91438103), the National Natural Science Foundation of China (Nos. 61876220, 61876221, 61836009, U1701267, 61871310, 61573267, 61502369 and 61473215), the Science Foundation of Xidian University (Nos. 10251180018 and 10251180019), the Program for Cheung Kong Scholars and Innovative Research Team in University (No. IRT 15R53), the Fund for Foreign Scholars in University Research and Teaching Programs (No. B07048), and the Key Special Project of China High Resolution Earth Observation System-Young Scholar Innovation Fund.