Accelerated Doubly Stochastic Gradient Algorithm for Large-scale Empirical Risk Minimization

1 Introduction

In this paper, we consider the following problem:

where F(x) = 1n ∑ni=1 fi(x) is the average of n smooth convex component functions fi's and P(x) is a blockseparable, possibly non-smooth, convex regularization function. In machine learning applications, many tasks can be naturally phrased as the above problem, e.g. the Empirical Risk Minimization (ERM) problem [Zhang and Xiao, 2015; Friedman et al., 2001]. However, the latest explosive growth of data introduces unprecedented computational challenges of scalability and storage bottleneck. In consequence, algorithms with fast convergence, small memory footprints, and low per-iteration complexity have been ardently pursued in the recent years. Most successful practices adopt stochastic strategies that incorporate randomness into solving procedures. They followed two parallel tracks. That is, in each iteration, gradient is estimated by either a mini batch of samples or a small block of variable coordinates, commonly designated by a random procedure.

∗Corresponding author.

Accessing only a random mini batch of samples in each iteration, classical Stochastic Gradient Descent (SGD) and its Variance Reduction (VR) variants, such as SVRG [Johnson and Zhang, 2013], SAGA [Defazio et al., 2014], and SAG [Schmidt et al., 2013], have gained increasing attention in the last quinquennium. To further improve the performance, a significant amount of efforts have been made towards reincarnating Nesterov's optimal accelerated convergence rate in SGD type methods [Zhang and Xiao, 2015; Frostig et al., 2015; Lin et al., 2015a; Shalev-Shwartz and Zhang, 2014; Hu et al., 2009; Lan, 2012; Nitanda, 2014]. Recently, a direct accelerated version of SVRG called Katyusha is proposed to obtain optimal convergence results without compromising the low per-iteration sample access [Allen-Zhu, 2016; Woodworth and Srebro, 2016]. However, when datasets are high dimensional, SGD type methods may still suffer from the large memory footprint due to its full vector operation in each iteration, which lefts plenty of scope to push further.

Updating variables only on a randomly selected small block of variable coordinates in each iteration is another important strategy that can be adopted solely to reduce the memory reference. The most noteworthy endeavors include the Randomized Block Coordinate Descent (RBCD) methods [Nesterov, 2012; Richt´arik and Tak´aˇc, ; Lee and Sidford, 2013; Wright, 2015]. Now accelerated versions of RBCD type methods, such as APPROX [Fercoq and Richt´arik, 2015] and APCG [Lin et al., 2015b], have also made their debuts. A drawback of RBCD type methods is that all samples have to be accessed in each iteration. When the number of samples is huge, they can be still quite inefficient.

Doubly stochastic algorithms, simultaneously utilizing the idea of randomness from sample choosing and coordinate selection perspective, have emerged more recently. Zhao et al. carefully combine ideas from VR and RBCD and propose a method called Mini-batch Randomized Block Coordinate Descent with variance reduction (MRBCD), which achieves linear convergence in strongly convex case [Zhao et al., 2014]. Zhang and Gu propose the Accelerated Stochastic Block Coordinate Descent (ASBCD) method, which incorporates RBCD into SAGA and uses a non-uniform probability when sampling data points [Zhang and Gu, 2016]. These methods avoid both full dataset assess and full vector operation in each iteration, and therefore are amenable to solving (1) when n and d are large at the same time. However, none of themmeets the optimal convergence rate [Woodworth and Srebro, 2016], and hence can still be accelerated.

To bridge the gap, we introduce a multi-momentum technique in this paper and devise a method called Accelerated Doubly Stochastic Gradient algorithm (ADSG). Our method enjoys an accelerated convergence rate, superior to existing doubly stochastic methods, without compromising the low sample access and small periteration complexity. Specifically, our contributions are listed as follows.

1. With two novel coupling steps, we incorporate three momenta into ADSG. These two steps enable the acceleration of doubly stochastic optimization procedure. Additionally, we devise an efficient implementation of ADSG for ERM problem.
2. We prove that ADSG has an accelerated convergence rate and the overall computational complexity to obtain an ϵ-accurate solution is O(log(1/ϵ)/√ϵ) in the general convex case and O((n+ √nκ) log(1/ϵ)) in the strongly convex case, where κ stands for the condition number.

Further, to show the efficiency of ADSG in practice, we conduct learning tasks on huge datasets with more than 10M samples and 1M features. The results demonstrate superior computational efficiency of our approach compared to the state-of-the-art.

2 Preliminary

2.1 Notation & Assumptions

We assume that the variable x ∈ Rd can be equally partitioned into B blocks for simplicity, and we let Ω = d/B be the block size. We denote the coordinates in the lth block of x by [x]l ∈ RΩ and the rest by [x]\l. The regularization P(x) is assumed to be block separable with respect to the partition of x, i.e. P(x) = ∑Bl=1 Pl([x]l). Many important functions qualify such separability assumption [Tibshirani, 1996; Simon et al., 2013]. The proximal operator of a convex function g is defined as proxg(y) = argminx g(x) + 12 ‖x - y‖2 where we use ‖ · ‖ to denote the Euclidean norm. We say x is an ϵ-accurate solution of to Problem (1) if FP(x)-FP(x∗) ≤ ϵ, where x∗ is the optimal. Further, we define ϵ0 to be FP(x0), i.e. the objective value at the initial point x0. We now give some assumptions.

Assumption I Each function fi is L-smooth, i.e. ∀x, y ∈ Rd , f (y) ≤ f (x) + 〈∇f (x), y - x〉 + L2 ‖x - y‖2 . Assumption II The average function F is LB-smooth with respect to each block, i.e. for any l and ∀x, h ∈ Rd such that [h]\l ≡ 0,f(x+h) ≤ f(x)+〈∇f(x), h〉+ LB 2 ‖h‖2 . Assumption III Function P is µ-strongly convex, i.e. ∀x, y ∈ Rd , and g ∈ ∂P(x), P(y) ≥ P(x) + 〈g, y - x〉 + µ2 ‖x - y‖2 .

We define κ = L+LB µ as the condition number.

2.2 Doubly Stochastic Algorithms for Convex Composite Optimization

A doubly stochastic method that incorporates both VR and CD technique usually consists of two nested loops. More specifically, at the beginning of each outer loop, the exact full gradient at some snapshot point ˜x is calculated. Then, each iteration of the following inner loop estimates the partial gradient based on a minibatch component functions and modifies it by the obtained exact gradient to perform block coordinate descent. Take MRBCD for example. As noted in Table 1, its per-iteration vector operation complexity and sample access are O(Ω) and O(1) respectively, which is amenable to solving large-sample-high-dimension problems. However, its overall computational complexity to achieve an ϵ-accurate solution scales linearly with respect to the condition number κ in strongly convex case and depends on 1/ϵ (up to a log factor) in general convex case. Such complexity leaves much room to improve in doubly stochastic methods, especially when a highly accurate solution to an ill-conditioned problem is sought (small ϵ, large κ).

2.3 Momentum Accelerating Technique

Accelerating the first order algorithm with momentum is by no mean new but has always been considered difficult [Polyak, 1964; Allen-Zhu and Orecchia, 2014]. Nesterov is the first to prove the accelerated convergence rate for deterministic smooth convex optimization [Nesterov, 1983; 1998]. When randomness is involved in data point sampling, noise in the stochastic gradient is further accumulated by the momentum, making the acceleration much harder. Efforts have been made to overcome such difficulty in various ways: with an outer-inner loop manner [Lin et al., 2015a], from a dual perspective [Shalev-Shwartz and Zhang, 2014], or under a primal-dual framework [Zhang and Xiao, 2015]. However, it was not until recently that an optimal method called Katyusha is proposed [Allen-Zhu, 2016; Woodworth and Srebro, 2016]. In a parallel line, at the first attempt to accelerate the RBCD, the momentum step forces full vector operation in each iteration and hence compromises the advantage of RBCD [Nesterov, 2012]. With continuing efforts, Lin et al. proposed the APCG method for both strongly convex and general convex problems, avoiding full vector operation completely [Lin et al., 2015b]. While these previous mentioned analyses deal with the randomness only in data point sampling or coordinate choosing, our analysis, due to doubly stochastic nature of ADSG, has to consider how randomness affects the momenta from both sample and feature perspective, and hence is more difficult.

2.4 Empirical Risk Minimization

We focus on Empirical Risk Minimization (ERM) with linear predictor, an important class of smooth convex problems. Specifically, each fi in Problem (1) is of the form fi(x) = φi(a⊤i x), where ai is the feature vector of the ith sample and φi(·) : R→R is some smooth convex function. Let A = [a1 . . . an]⊤ be the data matrix. In real applications, A is usually very sparse and we define its sparsity to be ρ = nnz(A) nd . Note that our convergence analysis applies to any convex smooth fi.

3 Methodology

We present the proposed ADSG in algorithm 1, and discuss about some crucial details in this section.

Input: The input of ADSG varies for strongly convex and general convex problems. In the former case, µ should be the strongly convex parameter, and we set α2,s = 1 2B min{1, √ nκ } and α3,s = 1 2B . In the latter case, µ is set to 0, and we set α2,s = 2 s+4B and α3,s = 1 2B . When µ = 0, line 14 adopts uniform probability for sampling. Main Body: Our algorithm is divided into epochs. Four variables xk, yk, zk, and ˜xs are maintained throughout. At the beginning of each epoch, the full gradient at the snapshot point ˜xs is computed. Updating steps are taken in the follow-up m inner loops, where we randomly select a mini-batch I of size b and a feature block l to construct a mixed stochastic gradient at

Algorithm 1 ADSG I

Input: n, x0, S, b,B, µ, {α2,s, α3,s}Ss=0 1: z0 = ˜x0 ←x0,m←Bn, α1,s ← 1-α2,s -α3,s; 2: for s ← 0 to S do 3: ¯Ls ← L Bα3,s +LB; 4: ηs ← 1 ¯Lsα2,sB , θs ← 1 + µ ¯LsB2α2,s+(B-1)µ ; 5: ˜∇s ←∇F(˜xs ); 6: for j ← 1 to m do 7: k ← sm + j; 8: yk ← α1,sxk-1 + α2,szk-1 + α3,s ˜xs ; 9: sample mini batch I of size b and feature block l; 10: [vk]l ← [ ˜∇s ]l + 1b ∑i∈I ([∇fi(yk)]l - [∇fi(˜xs )]l ); 11: [zk]l←proxηsPl([zk-1]l -ηs[vk]l), [zk]\l ← [zk-1]\l; 12: xk ← yk + α2,sB(zk - zk-1); 13: end for 14: sample σ from {1, . . . ,m} with probability θσ-1 s ∑m-1 i=0 θi-1 s ; 15: ˜xs+1 ←xsm+σ; 16: end for

point yk and perform proximal coordinate descent on the auxiliary variable zk.

Momenta: In sharp contrast to existing doubly stochastic algorithms, two coupling steps are added in ADSG to accelerate the convergence. In line 8, yk is constructed as the convex combination of xk-1, zk-1, and the snapshot point ˜xs . Here, zk-1 acts as a historical momentum that adds weight to the previous stochastic gradient {vt}t≤k (note that zk is simply the linear combination of all {vt}t≤k when P ≡ 0). This historical momentum is used in many deterministic accelerated methods, e.g. [Beck and Teboulle, 2009]. ˜xs serves as a negative momentum that ensures the "gradient" variable yk not drifting away from ˜xs and prevents the variance introduced by the randomness from surging. Such negative momentum is recently proposed in [Allen-Zhu, 2016], but with a different weight α3,s = 12 . Note that such weight is crucial to the convergence, see Theorem 1 and 2. In line 12, we have Elxk = α1xk-1 + α2zk-1 + α3 ˜xs + α2(˜zk - zk-1), where ˜zk = proxηsP(zk-1 - ηsvk). Hence, α2,sB(zk - zk-1), in expectation, adds extra weight on the most recent progress ˜zk - zk-1 and is called momentum in expectation for this reason. These three momenta are the key to the accelerated convergence of ADSG.

4 Convergence Analysis

In this section, we present the convergence rate of ADSG in both strongly convex case and general convex case.

4.1 Strongly Convex Case

If κ = O(n), we can see from Table 1 that existing doubly stochastic algorithm like MRBCD requires only O(log 1/ϵ) passes over the whole dataset to achieve an ϵaccurate solution. However, facing ill-conditioned problems where κ > n2 , the performance of such method decays faster than the deterministic APG method due to its linear dependence on the condition number κ. The following theorem shows that ADSG enjoys an accelerated convergence rate and depends only on √κ.

Theorem 1. Suppose Assumption I-III are satisfied. Set α2,s = 1 2B min{1, √ nκ }, α3,s = 1 2B and set the mini batch size to be 1. If κ > 8B, for s ≥ 0, we have

and therefore ADSG takes O((1 + √ κn ) log ϵ0 ϵ ) outer loops to achieve an ϵ-accurate solution.

4.2 General Convex Case

We present both direct and indirect methods to solve Problem (1) in the general convex case. We make the assumption that the optimal x∗ is bounded ‖x∗‖2 ≤ D, which is common in the literature.

Indirect Method: One option is to obtain an ϵ/2accurate solution ¯x to the following strongly convex, but perturbed stand-in of Problem (1):

Such strategy is reasonable in that ϵ2 ≥ FPϵ (¯x)-FPϵ (x∗) = FP(¯x)-FP(x∗)+ ϵ 2D ‖¯x‖2 - ϵ 2D ‖x∗‖2 ⇒ FP(¯x)-FP(x∗) ≤ ϵ, and hence ¯x is an ϵ-accurate solution to Problem (1). Furthermore, since the condition number of Problem (2) is DL ϵ , Theorem 1 implies the following corollary.

Corollary 1. Suppose the optimal point x∗ is bounded ‖x∗‖2 ≤D and suppose ¯x is an ϵ/2-accurate solution to Problem (2). Then ¯x is an ϵ-accurate solution to Problem (1). Further, if Assumption I and II are satisfied, then it takes O((1+√ LD nϵ ) log ϵ0 ϵ ) outer loops to compute ¯x in ADSG.

Direct Method: With different parameter setting, ADSG solves Problem (1) explicitly in general convex case, yielding a direct method with accelerated convergence rate. In this way, we avoid the extra log ϵ0 ϵ factor in the indirect method.

Theorem 2. Suppose Assumption I and II are satisfied. Set α2,s = 2 s+4B , α3,s = 1 2B , for s ≥ 0, we have

and therefore ADSG takes O(B√ ϵ0 ϵ + LD nϵ ) outer loops to achieve an ϵ-accurate solution.

The extra B in Theorem 2 leaves room for improvement. The following theorem shows the dependence on B can be alleviated when P ≡ 0.

Theorem 3. Suppose Assumption I and II are satisfied and P ≡ 0. Denote L′ = L + BLB. Setting α2,s = 2 s+4 , α3,s = 12 , and α1,s = 1-α2,s -α3,s, then for s ≥ 0,

Therefore, ADSG takes O(√ ϵ0 ϵ + L′D nϵ ) outer loops to achieve an ϵ-accurate solution.

4.3 Proof for Strongly Convex Case

In this section, we give two key lemmas for our analysis. The full proofs of Theorem (1-3) are deferred to the long version of this paper.

The idea of the first lemma is to express xk as the convex combination of {˜xi }si=0 and {zl}kl=0 .

Lemma 1. In Algorithm I, by setting α2,0 = α3,0 = 1/2B, for k = sm+ j ≥ 1, we have

where γ 0 0 = 1, γ 0 1 = 1 2 - 1 2B , γ 1 1 = 1 2 , β0 0 = 0, βs 0 = α3,s, λs(s+1)m = βsm, λik+1 = α1,sλik ,

and

Additionally, we have ∑s-1 i=0 λik +βs j +∑kl=0 γ lk = 1. If all α1,s ≥ B-1 B , then each entry in this sum is non-negative for all k ≥ 1, i.e. xk is a convex combination of {˜xi }si=0 and {zl}kl=0 .

The second lemma analyzes ADSG in one iteration. Before proceeding, we first define a few terms: (1) ˆP(xk) def = ∑s-1 i=0 λikP(˜xi )+βs j P(˜xs )+∑kl=0 γ lkP(zl) ≥ P(xk), where the inequality uses the convexity of P; (2) d(xk) def = FP(xk) - FP(x∗); (3) ˜ ds def = FP(˜xs ) - FP(x∗); and (4) dk def = (F(xk) + ˆP(xk)) - FP(x∗). Additionally, we have 0 ≤ d(xk) ≤ dk and d(x0) = d0.

Lemma 2. In the sth epoch of Algorithm 1, we have

with ρs = α2 2,sB2 ¯Ls + (B - 1)µα2,s, θs = 1 + µα2,s ρs , Aj = ρs‖x∗ -zsm+j‖2 /2, and ˆ dj = dsm+j.

5 Efficient Implementation

While ADSG has an accelerated convergence rate, naively implementing Algorithm 1 requires O(d) computation in each inner loop due to the two coupling steps, which compromises the low per-iteration complexity enjoyed by RBCD type methods. Such quandary strikes all existing accelerated RBCD algorithm [Lin et al., 2015b; Nesterov, 2012; Fercoq and Richt´ arik, 2015; Lee and Sidford, 2013]. To bypass this dilemma, we cast Algorithm 1 in an equivalent but more practical form, ADSG II, with the inner loop complexity reduced to O(Ω). ADSG II uses three auxiliary functions {¯xk, ¯yk, ¯zk} which are defined as

with k = sm+ j. The following proposition shows the equivalence between Algorithm 1 and 2.

Algorithm 2 ADSG II

Input: n, x0, S, b,B, µ, {α2,s, α3,s}Ss=0 1: ˆz00←⃗0, ξ= ˙x0 ←x0, k←0,m←Bn, α1,s←1-α2,s-α3,s; 2: for s ← 0 to S do 3: γs ← α2,s α2,s+α3,s , βs-1 ← 1, βs 0 ←α1,s, ¯Ls ← L Bα3,s +LB; 4: ηs ← 1 ¯Lsα2,sB , θs ← 1 + µ ¯LsB2α2,s+(B-1)µ ; 5: ˙∇s ← ∇f ( ˙xs ); 6: us0 ← ξ -γsˆzs0 - ˙xs ; 7: for j ← 1 to m do 8: k ← sm + j; 9: sample mini batch I of size b and feature block l; 10: [ ˙vk]l ← [ ˙∇s ]l + 1b ∑i∈I ([∇fi(¯yk)]l - [∇fi(˜xs )]l ); 11: [ˆzsj]l←proxηsPl([¯zk]l-ηs[ ˙vk]l)-[ ˙xs ]l, [ˆzsj]\l←[ˆzsj-1]\l; 12: usj ← usj-1 + α2,sB-γs βs j-1 (ˆzsj - ˆzsj-1); βs j ← α1,sβs j-1; 13: end for 14: sample σ from {1, . . . ,m} with probability θσ-1 s ∑m-1 i=0 θi-1 s ; 15: ˙xs+1 ← ¯xsm+σ; 16: ˆzs+1 0 ← ¯zk - ˙xs+1 , ξ = ¯xk; 17: end for

Proposition 1. If Algorithm 1 has the same input as Algorithm 2, its iterates xk, yk, and zk equal to ¯xk, ¯yk, ¯zk respectively for all k.

Proof. First, we prove that if at the beginning of the sth epoch, ¯zsm = zsm, ¯xsm = xsm, and ˙xs = ˜xs stand, then the proposition stand in the following iterations in that epoch. We prove with induction. Assume that the equivalence holds till the (k - 1)th iteration. In the kth iteration, for ¯yk we have

since α1,sβs j-2 = βs j-1 and α1,sγs+α2,s = γs. For ¯zk, by induction we have [¯zk]\l = [¯zk-1]\l = [zk-1]\l = [zk]\l. Additionally, since ˙∇s = ˜∇s and yk = ¯yk, we have [vk]l = [ ˙vk]l and thus [¯zk]l = [ˆzsj + ˙xs ]l =proxηPl ([¯zk-1η ˙vk]l)=proxηPl ([zk-1-ηvk]l) = [zk]l. For ¯xk, we have

by the updating rule of uk in line 12 in ADSG II.

We then show that at the beginning of each epoch ¯zsm = zsm, ¯xsm = xsm, and ˙xs = ˜xs stand. For z0, we have ¯z0 = z0 from the initialization. For zsm, s ≥ 1, we have zsm = ¯zsm = ˆzs0 + ˙xs = ¯zsm, where the first equation is from the induction in previous epoch, and the second equation is from the definition of ˆzs0 in line 16 in ADSG II. For x0, we clearly have x0 = ¯x0 by the initialization. For xsm, s ≥ 1, we have xsm = ξ = βs-1us0 +γsˆzs0 + ˙xs = ¯xsm, where the first equation is from the induction in previous epoch, and the second equation is from line 6 in ADSG II. ˙xs = ˜xs because xk = ¯xk for all k in that epoch. Thus we have the result.

5.1 Overall Computational Complexity

We discuss the detailed implementation of ADSG II when solving ERM problems. For simplicity, the mini batch size b is set to 1. In line 10, [∇fi(¯yk)]l = ∇φi(a⊤i ¯yk)[ai]l, where we compute each term in a⊤i ¯yk = βj-1a⊤i usj-1 +γsa⊤i ˆzsj-1 +a⊤i ˙xs separately. This can be done with O(ρd). Line 12 takes O(Ω) because only the lth block is updated. Consequently, we have O(ρd +Ω) from every inner loop and the per-epoch complexity of Algorithm 2 is O(BnΩ + Bnρd). When ρ is small, Ω dominates ρd and it becomes O(dn).

Combining the above per-epoch complexity analysis and the convergence rate, the overall computational complexity of ADSG is O(d(n + √n/κ) log 1/ϵ) in strongly convex case, and O(d(n+ √nLD/ϵ) log 1/ϵ) in general convex case.

6 Experiments

In this section, we conduct several experiments to show the efficiency of ADSG on huge-scale real problems. We use FP(xt )-FP(x∗ ) FP(x0)-FP(x∗ ) to measure the suboptimality.

Problems We conduct experiments of three ERM problems, namely l1-logistic regression, l1l2-logistic regression, and ridge regression. Four large scale datasets from LibSVM [Chang and Lin, 2011] are used: kdd2010raw, avazu-app, new20.binary, and url-combined. Their statistics are given in Table 2.

Algorithms Katyusha [Allen-Zhu, 2016], MRBCD [Zhao et al., 2014] (ASBCD has similar performance), and SVRG [Johnson and Zhang, 2013] are included for comparison. MRBCD and ADSG adopts the same block parameter B. All methods use the same mini batch size b. We use the default inner loop count described in the original paper for SVRG, MRBCD, and Katyusha. We tune the step size to give the best performance. For SVRG and MRBCD, it is usually 1/L.

6.1 l1-Logistic Regression

In this problem, we set fi(x) = log(1 + exp(-yia⊤i x)) and set P(x) = λ‖x‖1, where {ai, yi}ni=1 are data points. We present suboptimality vs. Evaluated Partial Gradients (EPG) and suboptimality vs. time in the first two rows of Figure 1, along with the parameter λ used in experiments. The result shows that (i) Katyusha and ADSG have superior convergence rate over nonaccelerated SVRG and MRBCD, (ii) ADSG and MRBCD, as doubly stochastic methods, enjoy a better time efficiency than Katyusha and SVRG, and (iii) ADSG has the best performance among all competitors.

6.2 l1l2-Logistic Regression

fi(x) is log(1+exp(-yia⊤i x)) and P(x) is µ2 ‖x‖2 +λ‖x‖1 in l1l2-logistic regression. The log suboptimality vs. EPG and log suboptimality vs. time plots are given in the third and fourth rows of Figure 1. We observe similar phenomenon here as in the previous experiment.

6.3 Ridge Regression

We set fi(x) = 1 2 ‖a⊤i x-yi‖2 and P(x) = µ2 ‖x‖2 . In the last two rows of Figure 1, we give the log suboptimality vs. EPG and log suboptimality vs. time plots. ADSG shows exceptional computational efficiency.

7 Conclusion

An accelerated doubly stochastic algorithm called ADSG is proposed in this paper. We give its convergence analyses, and compare our algorithm to the state-of-the-art in large scale ERM problems. The result is promising.

Acknowledgments

This work is partially supported by National Natural Science Foundation of China (Grant No: 61472347, 61672376, and 61672449)