On Variance Reduction in Stochastic Gradient Descent and its Asynchronous Variants

Sashank J. Reddi, Ahmed Hefny, Suvrit Sra, Barnabás Póczos, Alex Smola

Introduction

There has been a steep rise in recent work on “variance reduced” stochastic gradient algorithms for convex problems of the finite-sum form:

Although these advances have great value in general, for large-scale problems we still require parallel or distributed processing. And in this setting, asynchronous variants of SGD remain indispensable . Therefore, a key question is how to extend the synchronous finite-sum VR algorithms to asynchronous parallel and distributed settings.

We answer one part of this question by developing new asynchronous parallel stochastic gradient methods that provably converge at a linear rate for smooth strongly convex finite-sum problems. Our methods are inspired by the influential Svrg , S2gd , Sag and Saga family of algorithms. We list our contributions more precisely below.

Contributions. Our paper makes two core contributions: (i) a formal general framework for variance reduced stochastic methods based on discussions in ; and (ii) asynchronous parallel VR algorithms within this framework. Our general framework presents a formal unifying view of several VR methods (e.g., it includes SAGA and SVRG as special cases) while expressing key algorithmic and practical tradeoffs concisely. Thus, it yields a broader understanding of VR methods, which helps us obtain asynchronous parallel variants of VR methods. Under sparse-data settings common to machine learning problems, our parallel algorithms attain speedups that scale near linearly with the number of processors.

As a concrete illustration, we present a specialization to an asynchronous Svrg-like method. We compare this specialization with non-variance reduced asynchronous SGD methods, and observe strong empirical speedups that agree with the theory.

Related work. As already mentioned, our work is closest to (and generalizes) Sag , Saga , Svrg and S2gd , which are primal methods. Also closely related are dual methods such as sdca and Finito , and in its convex incarnation Miso ; a more precise relation between these dual methods and VR stochastic methods is described in Defazio’s thesis . By their algorithmic structure, these VR methods trace back to classical non-stochastic incremental gradient algorithms , but by now it is well-recognized that randomization helps obtain much sharper convergence results (in expectation). Proximal and accelerated VR methods have also been proposed ; we leave a study of such variants of our framework as future work. Finally, there is recent work on lower-bounds for finite-sum problems .

Within asynchronous SGD algorithms, both parallel and distributed variants are known. In this paper, we focus our attention on the parallel setting. A different line of methods is that of (primal) coordinate descent methods, and their parallel and distributed variants . Our asynchronous methods share some structural assumptions with these methods. Finally, the recent work generalizes S2GD to the mini-batch setting, thereby also permitting parallel processing, albeit with more synchronization and allowing only small mini-batches.

A General Framework for VR Stochastic Methods

While our analysis focuses on strongly convex functions, we can extend it to just smooth convex functions along the lines of .

Next, we provide different instantiations of the framework and construct a new algorithm derived from it. In particular, we consider incremental methods Sag , Svrg and Saga , and classic gradient descent GradientDescent for demonstrating our framework.

Figure 1 shows the schedules for the aforementioned algorithms. In case of Svrg, ScheduleUpdate is triggered every mm iterations (here mm denotes precisely the number of inner iterations used in ); so AtA^{t} remains unchanged for the mm iterations and all αit\alpha_{i}^{t} are updated to the current iterate at the mthm^{\text{th}} iteration. For Saga, unlike Svrg, AtA^{t} changes at the ttht^{th} iteration for all t∈[T]t\in[T]. This change is only to a single element of AtA^{t}, and is determined by the index iti_{t} (the function chosen at iteration tt). The update of Sag is similar to Saga insofar that only one of the αi\alpha_{i} is updated at each iteration. However, the update for At+1A^{t+1} is based on it+1i_{t+1} rather than iti_{t}. This results in a biased estimate of the gradient, unlike Svrg and Saga. Finally, the schedule for gradient descent is similar to Sag, except that all the αi\alpha_{i}’s are updated at each iteration. Due to the full update we end up with the exact gradient at each iteration. This discussion highlights how the scheduler determines the resulting gradient method.

To motivate the design of another schedule, let us consider the computational and storage costs of each of these algorithms. For Svrg, since we update AtA^{t} after every mm iterations, it is enough to store a full gradient, and hence, the storage cost is O(d)O(d). However, the running time is O(d)O(d) at each iteration and O(nd)O(nd) at the end of each epoch (for calculating the full gradient at the end of each epoch). In contrast, both Sag and Saga have high storage costs of O(nd)O(nd) and running time of O(d)O(d) per iteration. Finally, GradientDescent has low storage cost since it needs to store the gradient at O(d)O(d) cost, but very high computational costs of O(nd)O(nd) at each iteration.

Svrg has an additional computation overhead at the end of each epoch due to calculation of the whole gradient. This is avoided in Sag and Saga at the cost of additional storage. When mm is very large, the additional computational overhead of Svrg amortized over all the iterations is small. However, as we will later see, this comes at the expense of slower convergence to the optimal solution. The tradeoffs between the epoch size mm, additional storage, frequency of updates, and the convergence to the optimal solution are still not completely resolved.

A straightforward approach to design a new scheduler is to combine the schedules of the above algorithms. This allows us to tradeoff between the various aforementioned parameters of our interest. We call this schedule hybrid stochastic average gradient (Hsag). Here, we use the schedules of Svrg and Saga to develop Hsag. However, in general, schedules of any of these algorithms can be combined to obtain a hybrid algorithm. Consider some S⊆[n]S\subseteq[n], the indices that follow Saga schedule. We assume that the rest of the indices follow an Svrg-like schedule with schedule frequency sis_{i} for all i∈S‾≜[n]∖Si\in\overline{S}\triangleq[n]\setminus S. Figure 2 shows the corresponding update schedule of Hsag. If S=[n]S=[n] then Hsag is equivalent to Saga, while at the other extreme, for S=∅S=\emptyset and si=ms_{i}=m for all i∈[n]i\in[n], it corresponds to Svrg. Hsag exhibits interesting storage, computational and convergence trade-offs that depend on SS. In general, while large cardinality of SS likely incurs high storage costs, the computational cost per iteration is relatively low. On the other hand, when cardinality of SS is small and sis_{i}’s are large, storage costs are low but the convergence typically slows down.

Before concluding our discussion on the general framework, we would like to draw the reader’s attention to the advantages of studying Algorithm 1. First, note that Algorithm 1 provides a unifying framework for many incremental/stochastic gradient methods proposed in the literature. Second, and more importantly, it provides a generic platform for analyzing this class of algorithms. As we will see in Section 3, this helps us develop and analyze asynchronous versions for different finite-sum algorithms under a common umbrella. Finally, it provides a mechanism to derive new algorithms by designing more sophisticated schedules; as noted above, one such construction gives rise to Hsag.

For any positive parameters c,β,κ>1c,\beta,\kappa>1, step size η\eta and epoch size mm, we define the following quantities:

Suppose the probabilities pi∝(1−1κ)m−ip_{i}\propto(1-\frac{1}{\kappa})^{m-i}, and that c,β,κc,\beta,\kappa, step size η\eta and epoch size mm are chosen such that the following conditions are satisfied:

Then, for iterates of Algorithm 1 under the Hsag schedule, we have

As a corollary, we immediately obtain an expected linear rate of convergence for Hsag.

We emphasize that there exist values of the parameters for which the conditions in Theorem 1 and Corollary 1 are easily satisfied. For instance, setting η=1/16(λn+L)\eta=1/16(\lambda n+L), κ=4/λη\kappa=4/\lambda\eta, β=(2λn+L)/L\beta=(2\lambda n+L)/L and c=2/ηnc=2/\eta n, the conditions in Theorem 1 are satisfied for sufficiently large mm. Additionally, in the high condition number regime of L/λ=nL/\lambda=n, we can obtain constant θ<1\theta<1 (say 0.50.5) with m=O(n)m=O(n) epoch size (similar to ). This leads to a computational complexity of O(nlog⁡(1/ϵ))O(n\log(1/\epsilon)) for Hsag to achieve ϵ\epsilon accuracy in the objective function as opposed to O(n2log⁡(1/ϵ))O(n^{2}\log(1/\epsilon)) for batch gradient descent method. Please refer to the appendix for more details on the parameters in Theorem 1.

Asynchronous Stochastic Variance Reduction

We are now ready to present asynchronous versions of the algorithms captured by our general framework. We first describe our setup before delving into the details of these algorithms. Our model of computation is similar to the ones used in Hogwild! and AsySCD . We assume a multicore architecture where each core makes stochastic gradient updates to a centrally stored vector xx in an asynchronous manner. There are four key components in our asynchronous algorithm; these are briefly described below.

Read: Read the iterate xx and compute the gradient ∇fit(x)\nabla f_{i_{t}}(x) for a randomly chosen iti_{t}.

Read schedule iterate: Read the schedule iterate AA and compute the gradients required for update in Algorithm 1.

Update: Update the iterate xx with the computed incremental update in Algorithm 1.

Schedule Update: Run a scheduler update for updating AA.

Each processor repeatedly runs these procedures concurrently, without any synchronization. Hence, xx may change in between Step 1 and Step 3. Similarly, AA may change in between Steps 2 and 4. In fact, the states of iterates xx and AA can correspond to different time-stamps. We maintain a global counter tt to track the number of updates successfully executed. We use D(t)∈[t]D(t)\in[t] and D′(t)∈[t]D^{\prime}(t)\in[t] to denote the particular xx-iterate and AA-iterate used for evaluating the update at the ttht^{\text{th}} iteration. We assume that the delay in between the time of evaluation and updating is bounded by a non-negative integer τ\tau, i.e., t−D(t)≤τt-D(t)\leq\tau and t−D′(t)≤τt-D^{\prime}(t)\leq\tau. The bound on the staleness captures the degree of parallelism in the method: such parameters are typical in asynchronous systems (see e.g., ). Furthermore, we also assume that the system is synchronized after every epoch i.e., D(t)≥kmD(t)\geq km for t≥kmt\geq km. We would like to emphasize that the assumption is not strong since such a synchronization needs to be done only once per epoch.

For the purpose of our analysis, we assume a consistent read model. In particular, our analysis assumes that the vector xx used for evaluation of gradients is a valid iterate that existed at some point in time. Such an assumption typically amounts to using locks in practice. This problem can be avoided by using random coordinate updates as in (see Section 4 of ) but such a procedure is computationally wasteful in practice. We leave the analysis of inconsistent read model as future work. Nonetheless, we report results for both locked and lock-free implementations (see Section 4).

Suppose step size η\eta, epoch size mm are chosen such that the following condition holds:

Then, for the iterates of an asynchronous variant of Algorithm 1 with Svrg schedule and probabilities pi=1/mp_{i}=1/m for all i∈[m]i\in[m], we have

The bound obtained in Theorem 2 is useful when Δ\Delta is small. To see this, as earlier, consider the indicative case where L/λ=nL/\lambda=n. The synchronous version of Svrg obtains a convergence rate of θ=0.5\theta=0.5 for step size η=0.1/L\eta=0.1/L and epoch size m=O(n)m=O(n). For the asynchronous variant of Svrg, by setting η=0.1/2(max⁡{1,Δ1/2τ}L)\eta=0.1/2(\max\{1,\Delta^{1/2}\tau\}L), we obtain a similar rate with m=O(n+Δ1/2τn)m=O(n+\Delta^{1/2}\tau n). To obtain this, set η=ρ/L\eta=\rho/L where ρ=0.1/2(max⁡{1,Δ1/2τ})\rho=0.1/2(\max\{1,\Delta^{1/2}\tau\}) and θs=0.5\theta_{s}=0.5. Then, a simple calculation gives the following:

where c′c^{\prime} is some constant. This follows from the fact that ρ=0.1/2(max⁡{1,Δ1/2τ})\rho=0.1/2(\max\{1,\Delta^{1/2}\tau\}). Suppose τ<1/Δ1/2\tau<1/\Delta^{1/2}. Then we can achieve nearly the same guarantees as the synchronous version, but τ\tau times faster since we are running the algorithm asynchronously. For example, consider the sparse setting where Δ=o(1/n)\Delta=o(1/n); then it is possible to get near linear speedup when τ=o(n1/2)\tau=o(n^{1/2}). On the other hand, when Δ1/2τ>1\Delta^{1/2}\tau>1, we can obtain a theoretical speedup of 1/Δ1/21/\Delta^{1/2}.

We finally provide the convergence result for the asynchronous algorithm in the general case. The proof is complicated by the fact that set AA, unlike in Svrg, changes during the epoch. The key idea is that only a single element of AA changes at each iteration. Furthermore, it can only change to one of the iterates in the epoch. This control provides a handle on the error obtained due to the staleness. Due to space constraints, the proof is relegated to the appendix.

For any positive parameters c,β,κ>1c,\beta,\kappa>1, step size η\eta and epoch size mm, we define the following quantities:

Suppose probabilities pi∝(1−1κ)m−ip_{i}\propto(1-\frac{1}{\kappa})^{m-i}, parameters β,κ\beta,\kappa, step-size η\eta, and epoch size mm are chosen such that the following conditions are satisfied:

Then, for the iterates of asynchronous variant of Algorithm 1 with Hsag schedule we have

By using step size normalized by Δ1/2τ\Delta^{1/2}\tau (similar to Theorem 2) and parameters similar to the ones specified after Theorem 1 we can show speedups similar to the ones obtained in Theorem 2. Please refer to the appendix for more details on the parameters in Theorem 3.

Before ending our discussion on the theoretical analysis, we would like to highlight an important point. Our emphasis throughout the paper was on generality. While the results are presented here in full generality, one can obtain stronger results in specific cases. For example, in the case of Saga, one can obtain per iteration convergence guarantees (see ) rather than those corresponding to per epoch presented in the paper. Also, Saga can be analyzed without any additional synchronization per epoch. However, there is no qualitative difference in these guarantees accumulated over the epoch. Furthermore, in this case, our analysis for both synchronous and asynchronous cases can be easily modified to obtain convergence properties similar to those in .

Experiments

We present our empirical results in this section. For our experiments, we study the problem of binary classification via l2l_{2}-regularized logistic regression. More formally, we are interested in the following optimization problem:

A careful implementation of Svrg is required for sparse gradients since the implementation as stated in Algorithm 1 will lead to dense updates at each iteration. For an efficient implementation, a scheme like the ‘just-in-time’ update scheme, as suggested in , is required. Due to lack of space, we provide the implementation details in the appendix.

We evaluate the following algorithms for our experiments:

Lock-Free Svrg: This is the lock-free asynchronous variant of Algorithm 1 using Svrg schedule; all threads can read and update the parameters with any synchronization. Parameter updates are performed through atomic compare-and-swap instruction . A constant step size that gives the best convergence is chosen for the dataset.

Locked Svrg: This is the locked version of the asynchronous variant of Algorithm 1 using Svrg schedule. In particular, we use a concurrent read exclusive write locking model, where all threads can read the parameters but only one threads can update the parameters at a given time. The step size is chosen similar to Lock-Free Svrg.

Lock-Free Sgd: This is the lock-free asynchronous variant of the Sgd algorithm (see ). We compare two different versions of this algorithm: (i) Sgd with constant step size (referred to as CSgd). (ii) Sgd with decaying step size η0σ0/(t+σ0)\eta_{0}\sqrt{\sigma_{0}/(t+\sigma_{0})} (referred to as DSgd), where constants η0\eta_{0} and σ0\sigma_{0} specify the scale and speed of decay. For each of these versions, step size is tuned for each dataset to give the best convergence progress.

All the algorithms were implemented in C++ All experiments were conducted on a Google Compute Engine n1-highcpu-32 machine with 32 processors and 28.8 GB RAM.. We run our experiments on datasets from LIBSVM websitehttp://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary.html. Similar to , we normalize each example in the dataset so that ∥zi∥2=1\|z_{i}\|_{2}=1 for all i∈[n]i\in[n]. Such a normalization leads to an upper bound of 0.250.25 on the Lipschitz constant of the gradient of fif_{i}. The epoch size mm is chosen as 2n2n (as recommended in ) in all our experiments. In the first experiment, we compare the speedup achieved by our asynchronous algorithm. To this end, for each dataset we first measure the time required for the algorithm to each an accuracy of 10−1010^{-10} (i.e., f(x)−f(x∗)<10−10f(x)-f(x^{*})<10^{-10}). The speedup with PP threads is defined as the ratio of the runtime with a single thread to the runtime with PP threads. Results in Figure 3 show the speedup on various datasets. As seen in the figure, we achieve significant speedups for all the datasets. Not surprisingly, the speedup achieved by Lock-free Svrg is much higher than ones obtained by locking. Furthermore, the lowest speedup is achieved for rcv1 dataset. Similar speedup behavior was reported for this dataset in . It should be noted that this dataset is not sparse and hence, is a bad case for the algorithm (similar to ).

For the second set of experiments we compare the performance of Lock-Free Svrg with stochastic gradient descent. In particular, we compare with the variants of stochastic gradient descent, DSgd and CSgd, described earlier in this section. It is well established that the performance of variance reduced stochastic methods is better than that of Sgd. We would like to empirically verify that such benefits carry over to the asynchronous variants of these algorithms. Figure 4 shows the performance of Lock-Free Svrg, DSgd and CSgd. Since the computation complexity of each epoch of these algorithms is different, we directly plot the objective value versus the runtime for each of these algorithms. We use 10 cores for comparing the algorithms in this experiment. As seen in the figure, Lock-Free Svrg outperforms both DSgd and CSgd. The performance gains are qualitatively similar to those reported in for the synchronous versions of these algorithms. It can also be seen that the DSgd, not surprisingly, outperforms CSgd in all the cases. In our experiments, we observed that Lock-Free Svrg, in comparison to Sgd, is relatively much less sensitive to the step size and more robust to increasing threads.

Discussion & Future Work

In this paper, we presented a unifying framework based on , that captures many popular variance reduction techniques for stochastic gradient descent. We use this framework to develop a simple hybrid variance reduction method. The primary purpose of the framework, however, was to provide a common platform to analyze various variance reduction techniques. To this end, we provided convergence analysis for the framework under certain conditions. More importantly, we propose an asynchronous algorithm for the framework with provable convergence guarantees. The key consequence of our approach is that we obtain asynchronous variants of several algorithms like Svrg, Saga and S2gd. Our asynchronous algorithms exploits sparsity in the data to obtain near linear speedup in settings that are typically encountered in machine learning.

For future work, it would be interesting to perform an empirical comparison of various schedules. In particular, it would be worth exploring the space-time-accuracy tradeoffs of these schedules. We would also like to analyze the effect of these tradeoffs on the asynchronous variants.

References

Appendix A Appendix

Notation: We use DfD_{f} to denote the Bregman divergence (defined below) for function ff.

We would like to clarify the definition of xkmx^{km} here. As noted in the main text, we assume that xkm+mx^{km+m} is replaced with an element chosen randomly from {xkm,…,xkm+m−1}\{x^{km},\dots,x^{km+m-1}\} with probability {p1,⋯ ,pm}\{p_{1},\cdots,p_{m}\} at the end of the (k+1)th(k+1)^{\text{th}} epoch. However, whenever xkmx^{km} appears in the analysis (proofs), it represents the iterate before this replacement.

Implementation Details

Since we are interested in sparse datasets, simply taking fi(x)=log⁡(1+exp⁡(yizi⊤x))+λ∥x∥2f_{i}(x)=\log(1+\exp(y_{i}z_{i}^{\top}x))+\lambda\|x\|^{2} is not efficient as it requires updating the whole vector xx at each iteration. This is due to the regularization term in each of the fif_{i}’s. Instead, similar to , we rewrite problem in (4.1) as follows:

Proof of Theorem 1

We expand function ff as f(x)=g(x)+h(x)f(x)=g(x)+h(x) where g(x)=1n∑i∈Sfi(x)g(x)=\frac{1}{n}\sum_{i\in S}f_{i}(x) and h(x)=1n∑i∉Sfi(x)h(x)=\frac{1}{n}\sum_{i\notin S}f_{i}(x). Let the present epoch be k+1k+1. We define the following:

The last step follows from convexity of ff and the unbiasedness of vtv^{t}. We have the following relationship between Gt+1G_{t+1} and GtG_{t}.

This follows from the definition of the schedule of Hsag for indices in SS. Substituting the above relationship in Equation (A.1) we get the following.

We describe the bounds for btb_{t} (defined below).

The terms T1T_{1} and T2T2 can be bounded in the following fashion:

Substituting these bounds T1T_{1} and T2T_{2} in btb_{t}, we get

The second inequality follows from Lemma 2. In particular, we use the fact that f(x)−f(x∗)=Df(x,x∗)f(x)-f(x^{*})=D_{f}(x,x^{*}) and Df(x,x∗)=Dg(x,x∗)+Dh(x,x∗)≥Dg(x,x∗)D_{f}(x,x^{*})=D_{g}(x,x^{*})+D_{h}(x,x^{*})\geq D_{g}(x,x^{*}). The third inequality follows from the following for the choice of our parameters:

Applying the recursive relationship on Rt+1R_{t+1} for m iterations, we get

Substituting the bound on btb_{t} from Equation (A.3) in the above equation we get the following inequality:

For obtaining the above inequality, we used the strongly convex nature of function ff. Again, using the Bregman divergence based inequality (see Lemma 2)

Using the above notation, we have the following inequality from Equation (A.4).

where θ<1\theta<1 is a constant that depends on the parameters used in the algorithm. ∎

Proof of Theorem 2

Let the present epoch be k+1k+1. Recall that D(t)D(t) denotes the iterate used in the ttht^{\text{th}} iteration of the algorithm. We define the following:

We first bound the last term of the above inequality. We expand the term in the following manner:

The first equality directly follows from the definition of utu^{t} and its property of unbiasedness. The second step follows from simple algebraic calculations. Terms T3T_{3} and T4T_{4} can be bounded in the following way:

This bound directly follows from convexity of function fitf_{i_{t}}.

The first inequality follows from lipschitz continuous nature of the gradient of function fitf_{i_{t}}. The second inequality follows from the definition of Δ\Delta. The last term T5T_{5} can be bounded in the following manner.

The first inequality follows from Cauchy-Schwartz inequality. The second inequality follows from repeated application of triangle inequality. The third step is a simple application of AM-GM inequality and the fact that gradient of the function fitf_{i_{t}} is lipschitz continuous. Finally, the last step can be obtained by using a simple counting argument, the fact that the staleness in gradient is at most τ\tau and the definition of Δ\Delta.

By combining the bounds on T3,T4T_{3},T_{4} and T5T_{5} in Equations (A.7), (A.8) and (A.9) respectively and substituting the sum in Equation (A.6), we get

By substituting the above inequality in Equation (A.5), we get

The first step follows from Lemma 3 for r=2r=2. The third inequality follows from the lipschitz continuous nature of the gradient and simple application of Lemma 3. Adding the above inequalities from t=kmt=km to t=km+m−1t=km+m-1, we get

Here we again used a simple counting argument and the fact that the delay in the gradients is at most τ\tau. From the above inequality, we get

Adding Equation (A.11) from t=kmt=km to t=km+m−1t=km+m-1 and substituting Equation (A.12) in the resultant, we get

Substituting this in the inequality above, we get the following bound:

Proof of Theorem 3

Let the present epoch be k+1k+1. For simplicity, we assume that the iterates xx and AA used in the each iteration are from the same time step (index) i.e., D(t)=D′(t)D(t)=D^{\prime}(t) for all t∈Tt\in T. Recall that D(t)D(t) and D′(t)D^{\prime}(t) denote the index used in the ttht^{\text{th}} iteration of the algorithm. Our analysis can be extended to the case of D(t)≠D′(t)D(t)\neq D^{\prime}(t) in a straightforward manner. We expand function ff as f(x)=g(x)+h(x)f(x)=g(x)+h(x) where g(x)=1n∑i∈Sfi(x)g(x)=\frac{1}{n}\sum_{i\in S}f_{i}(x) and h(x)=1n∑i∉Sfi(x)h(x)=\frac{1}{n}\sum_{i\notin S}f_{i}(x). We define the following:

We use the same Lyapunov function used in Theorem 1. We recall the following definitions:

We bound term T6T_{6} in the following manner:

The first term can be bounded in the following manner:

The second step follows from Lemma 3 for r = 3. The last step follows from simple application of Jensen’s inequality. The first term can be bounded easily in the following manner:

The second and third terms need more delicate analysis. The key insight for our analysis is that at most τ\tau αi\alpha_{i}’s differ from time step D(t)D(t) to tt. This is due to the fact that the delay is bounded by τ\tau and at most one αi\alpha_{i} changes at each iteration. Furthermore, whenever there is a change in αi\alpha_{i}, it changes to one of the iterates xjx^{j} for some j={max⁡{t−τ,km},…,t}j=\{\max\{t-\tau,km\},\dots,t\}. With this intuition we bound the second term in the following fashion.

The first inequality follows from the fact that if αitD(t)\alpha_{i_{t}}^{D(t)} and αitt\alpha_{i_{t}}^{t} differ, then (a) iti_{t} should have been chosen in one of the iteration j∈{D(t),…,t−1}j\in\{D(t),\dots,t-1\} and (b) αit\alpha_{i_{t}} is changed to xjx^{j} in that iteration. The second inequality follows from Lemma 3 for r = 2. The third inequality follows from the fact that the probability P(ij=i)=1/nP(i_{j}=i)=1/n. The last step directly follows from Lemma 1. Note that sum is over indices in SS since αi\alpha_{i}’s for i∉Si\notin S do not change during the epoch.

The third term in Equation (A.15) can be bounded by exactly the same technique we used for the second term. The bound, in fact, turns out to identical to second term since iti_{t} is chosen uniformly random. Combining all the terms we have

The term T7T_{7} can be bounded in a manner similar to one in Theorem 2 to obtain the following (see proof of Theorem 2 for details):

We need the following bound for our analysis:

The above inequality follows directly from the bound on T6T_{6} by adding over all tt in the epoch. Under the condition

The above inequality follows from the fact that

The above relationship is due to the condition on η\eta and the fact that any d∈{D(t),…,t−1}d\in\{D(t),\dots,t-1\} for at most τ\tau values of tt. We have the following:

We bound ete_{t} in the following manner:

The second equality follows from the definition of Gt+1G_{t+1} (see Equation (A.2)).

Applying the recurrence relationship in Equation (A.18) with the derived bound on ete_{t}, we have

where et′e^{\prime}_{t} is defined as follows

The last inequality follows from that fact that the delay is at most τ\tau. In particular, each index j∈{D(t),…,t−1}j\in\{D(t),\dots,t-1\} occurs at most τ\tau times. We use the following notation for ease of exposition:

Substituting the bound in Equation (A.17), we get the following:

We now use the following previously used bound on vtv^{t} (see bound T2T_{2} in the proof of Theorem 1):

Substituting the above bound on vtv^{t} in Equation (A.19), we get the following:

where θa<1\theta_{a}<1 is a constant that depends on the parameters used in the algorithm. ∎

In this section, we briefly remark about the parameters in Theorems 1 & 3. For Theorem 1, suppose we use the following instantiation of the parameters:

In the interesting case of L/λ=nL/\lambda=n (high condition number regime), since κ=Θ(n)\kappa=\Theta(n), one can obtain a constant θ\theta (say θ=0.5\theta=0.5) with m=O(n)m=O(n). This leads to ϵ\epsilon accuracy in the objective function after O(log⁡(1/ϵ))O(\log(1/\epsilon)) epochs of Hsag. When m=O(n)m=O(n), the computational complexity of each epoch of Hsag is O(n)O(n). Hence, the total computational complexity of Hsag is O(nlog⁡(1/ϵ))O(n\log(1/\epsilon)). On the other hand, because L/λ=nL/\lambda=n, batch gradient descent method requires O(nlog⁡(1/ϵ))O(n\log(1/\epsilon)) iterations to achieve ϵ\epsilon accuracy in the objective value. Since the complexity of each iteration of gradient descent is O(n)O(n) (as it passes through the whole dataset for calculating the gradient), the overall computational complexity of batch gradient descent is O(n2log⁡(1/ϵ))O(n^{2}\log(1/\epsilon)). In general, for high condition number regimes (which is typically the case in machine learning applications), Hsag (like Svrg, Saga) will be significantly faster than the batch gradient methods. Furthermore, the convergence rate is strictly better than the sublinear rate obtained for Sgd.

The parameter instantiations for Theorem 3 are much more involved. Suppose Δ1/2τ<1\Delta^{1/2}\tau<1 (this is the sparse regime that is typically of interest to the machine learning community) and m>n>9τm>n>9\tau. The other case ( Δ1/2τ≥1\Delta^{1/2}\tau\geq 1) can be analyzed in a similar fashion. We set the following parameters:

Again, in the case of L/λ=nL/\lambda=n, we can obtain constant θa\theta_{a} (say θa=0.5\theta_{a}=0.5) with m=Θ(n)m=\Theta(n) and κ=Θ(n)\kappa=\Theta(n). The constants in the parameters are not optimized and can be improved by a more careful analysis. Furthermore, sharper constants can be obtained in specific cases. For example, see and Theorem 2 for synchronous and asynchronous convergence rates of Svrg respectively. Similarly, sharper constants for Saga can also be derived by simple modifications of the analysis.

Other Lemmatta

The proof follows trivially from the fact that x∗x^{*} is the optimal solution and linearity and non-negative properties of Bregman divergence. ∎

For random variables z1,…,zrz_{1},\dots,z_{r}, we have