ASAGA: Asynchronous Parallel SAGA

Introduction

We consider the unconstrained optimization problem of minimizing a finite sum of smooth convex functions:

where each $f_{i}$ is assumed to be convex with $L$-Lipschitz continuous gradient, $f$ is $\mu$-strongly convex and $n$ is large (for example, the number of data points in a regularized empirical risk minimization setting). We define a condition number for this problem as $\kappa:=\nicefrac{{L}}{{\mu}}$. A flurry of randomized incremental algorithms (which at each iteration select $i$ at random and process only one gradient $f^{\prime}_{i}$) have recently been proposed to solve (1) with a fastTheir complexity in terms of gradient evaluations to reach an accuracy of $\epsilon$ is $O((n+\kappa)\log(\nicefrac{{1}}{{\epsilon}}))$, in contrast to $O(n\kappa\log(\nicefrac{{1}}{{\epsilon}}))$ for batch gradient descent in the worst case. linear convergence rate, such as Sag (Le Roux et al., 2012), Sdca (Shalev-Shwartz and Zhang, 2013), Svrg (Johnson and Zhang, 2013) and Saga (Defazio et al., 2014). These algorithms can be interpreted as variance reduced versions of the stochastic gradient descent (Sgd) algorithm, and they have demonstrated both theoretical and practical improvements over Sgd (for the finite sum optimization problem (1)).

In order to take advantage of the multi-core architecture of modern computers, the aforementioned optimization algorithms need to be adapted to the asynchronous parallel setting, where multiple threads work concurrently. Much work has been devoted recently in proposing and analyzing asynchronous parallel variants of algorithms such as Sgd (Niu et al., 2011), Sdca (Hsieh et al., 2015) and Svrg (Reddi et al., 2015; Mania et al., 2015; Zhao and Li, 2016). Among the incremental gradient algorithms with fast linear convergence rates that can optimize (1) in its general form, only Svrg has had an asynchronous parallel version proposed.We note that Sdca requires the knowledge of an explicit $\mu$-strongly convex regularizer in (1), whereas Sag / Saga are adaptive to any local strong convexity of $f$ (Schmidt et al., 2016; Defazio et al., 2014). This is also true for a variant of Svrg (Hofmann et al., 2015). No such adaptation has been attempted yet for Saga, even though one could argue that it is a more natural candidate as, contrarily to Svrg, it is not epoch-based and thus has no synchronization barriers at all.

In Section 2, we present a novel sparse variant of Saga that is more adapted to the parallel setting than the original Saga algorithm. In Section 3, we present Asaga, a lock-free asynchronous parallel version of Sparse Saga that does not require consistent reads. We propose a simplification of the “perturbed iterate” framework from Mania et al. (2015) as a basis for our convergence analysis. At the same time, through a novel perspective, we revisit and clarify a technical problem present in a large fraction of the literature on randomized asynchronous parallel algorithms (with the exception of Mania et al. (2015), which also highlights this issue): namely, they all assume unbiased gradient estimates, an assumption that is inconsistent with their proof technique without further synchronization assumptions. In Section 3.3, we present a tailored convergence analysis for Asaga. Our main result states that Asaga obtains the same geometric convergence rate per update as Saga when the overlap bound $\tau$ (which scales with the number of cores) satisfies $\tau\leq\mathcal{O}(n)$ and $\tau\leq\mathcal{O}({\scriptstyle\frac{1}{\sqrt{\Delta}}}\max\{1,\frac{n}{\kappa}\})$, where $\Delta\leq 1$ is a measure of the sparsity of the problem, notably implying that a linear speedup is theoretically possible even without sparsity in the well-conditioned regime where $n\gg\kappa$. In Section 4, we provide a practical implementation of Asaga and illustrate its performance on a 40-core architecture, showing improvements compared to asynchronous variants of Svrg and Sgd.

The seminal textbook of Bertsekas and Tsitsiklis (1989) provides most of the foundational work for parallel and distributed optimization algorithms. An asynchronous variant of Sgd with constant step size called Hogwild was presented by Niu et al. (2011); part of their framework of analysis was re-used and inspired most of the recent literature on asynchronous parallel optimization algorithms with convergence rates, including asynchronous variants of coordinate descent (Liu et al., 2015), Sdca (Hsieh et al., 2015), Sgd for non-convex problems (De Sa et al., 2015; Lian et al., 2015), Sgd for stochastic optimization (Duchi et al., 2015) and Svrg (Reddi et al., 2015; Zhao and Li, 2016). These papers make use of an unbiased gradient assumption that is not consistent with the proof technique, and thus suffers from technical problemsExcept Duchi et al. (2015) that can be easily fixed by incrementing their global counter before sampling. that we highlight in Section 3.2.

The “perturbed iterate” framework presented in Mania et al. (2015) is to the best of our knowledge the only one that does not suffer from this problem, and our convergence analysis builds heavily from their approach, while simplifying it. In particular, the authors assumed that $f$ was both strongly convex and had a bound on the gradient, two inconsistent assumptions in the unconstrained setting that they analyzed. We overcome these difficulties by using tighter inequalities that remove the requirement of a bound on the gradient. We also propose a more convenient way to label the iterates (see Section 3.2). The sparse version of Saga that we propose is also inspired from the sparse version of Svrg proposed by Mania et al. (2015). Reddi et al. (2015) presents a hybrid algorithm called Hsag that includes Saga and Svrg as special cases. Their asynchronous analysis is epoch-based though, and thus does not handle a fully asynchronous version of Saga as we do. Moreover, they require consistent reads and do not propose an efficient sparse implementation for Saga, in contrast to Asaga.

Sparse Saga

Borrowing our notation from Hofmann et al. (2015), we first present the original Saga algorithm and then describe a novel sparse variant that is more appropriate for a parallel implementation.

The standard Saga algorithm (Defazio et al., 2014) maintains two moving quantities to optimize (1): the current iterate $x$ and a table (memory) of historical gradients $(\alpha_{i})_{i=1}^{n}$.For linear predictor models, the memory $\alpha_{i}^{0}$ can be stored as a scalar. Following Hofmann et al. (2015), $\alpha_{i}^{0}$ can be initialized to any convenient value (typically ), unlike the prescribed $f^{\prime}_{i}(x_{0})$ analyzed in Defazio et al. (2014). At every iteration, the Saga algorithm samples uniformly at random an index $i\in\{1,\ldots,n\}$, and then executes the following update on $x$ and $\alpha$ (for the unconstrained optimization version):

where $\gamma$ is the step size and $\bar{\alpha}:=\nicefrac{{1}}{{n}}\sum_{i=1}^{n}\alpha_{i}$ can be updated efficiently in an online fashion. Crucially, $\mathbf{E}\alpha_{i}=\bar{\alpha}$ and thus the update direction is unbiased ($\mathbf{E}x^{+}=x-\gamma f^{\prime}(x)$). Furthermore, it can be proven (see Defazio et al. (2014)) that under a reasonable condition on $\gamma$, the update has vanishing variance, which enables the algorithm to converge linearly with a constant step size.

In its current form, every Saga update is dense even if the individual gradients are sparse due to the historical gradient ($\bar{\alpha}$) term. Schmidt et al. (2016) introduced a special implementation with lagged updates where every iteration has a cost proportional to the size of the support of $f_{i}^{\prime}(x)$. However, this subtle technique is not easily adaptable to the parallel setting (see App. F.2). We therefore introduce Sparse Saga, a novel variant which explicitly takes sparsity into account and is easily parallelizable.

As in the Sparse Svrg algorithm proposed in Mania et al. (2015), we obtain Sparse Saga by a simple modification of the parameter update rule in (2) where $\bar{\alpha}$ is replaced by a sparse version equivalent in expectation:

where $D_{i}$ is a diagonal matrix that makes a weighted projection on the support of $f^{\prime}_{i}$. More precisely, let $S_{i}$ be the support of the gradient $f_{i}^{\prime}$ function (i.e., the set of coordinates where $f_{i}^{\prime}$ can be nonzero). Let $D$ be a $d\times d$ diagonal reweighting matrix, with coefficients $\nicefrac{{1}}{{p_{v}}}$ on the diagonal, where $p_{v}$ is the probability that dimension $v$ belongs to $S_{i}$ when $i$ is sampled uniformly at random in $\{1,...,n\}$. We then define $D_{i}:=P_{S_{i}}D$, where $P_{S_{i}}$ is the projection onto $S_{i}$. The normalization from $D$ ensures that $\mathbf{E}D_{i}\bar{\alpha}=\bar{\alpha}$, and thus that the update is still unbiased despite the projection.

For clarity of exposition, we model our convergence result after the simple form of Hofmann et al. (2015, Corollary 3) (note that the rate for Sparse Saga is the same as Saga). The proof is given in Appendix B.

The lagged updates technique in Saga is based on the observation that the updates for component $[x]_{v}$ can be delayed until this coefficient is next accessed. Interestingly, the expected number of iterations between two steps where a given dimension $v$ is involved in the partial gradient is $p_{v}^{-1}$, where $p_{v}$ is the probability that $v$ is involved. $p_{v}^{-1}$ is precisely the term which we use to multiply the update to $[x]_{v}$ in Sparse Saga. Therefore one may view the Sparse Saga updates as anticipated Saga updates, whereas those in the Schmidt et al. (2016) implementation are lagged.

Although Sparse Saga requires the computation of the $p_{v}$ probabilities, this can be done during a first pass through the data (during which constant step size Sgd may be used) at a negligible cost. In our experiments, both Sparse Saga and Saga with lagged updates had similar convergence in terms of number of iterations, with the Sparse Saga scheme being slightly faster in terms of runtime. We refer the reader to Schmidt et al. (2016) and Appendix F for more details.

Asynchronous Parallel Sparse Saga

As most recent parallel optimization contributions, we use a similar hardware model to Niu et al. (2011). We have multiple cores which all have read and write access to a shared memory. They update a central parameter vector in an asynchronous and lock-free fashion. Unlike Niu et al. (2011), we do not assume that the vector reads are consistent: multiple cores can read and write different coordinates of the shared vector at the same time. This means that a full vector read for a core might not correspond to any consistent state in the shared memory at any specific point in time.

We first review the “perturbed iterate” framework recently introduced by Mania et al. (2015) which will form the basis of our analysis. In the sequential setting, stochastic gradient descent and its variants can be characterized by the following update rule:

where $i_{t}$ is a random variable independent from $x_{t}$ and we have the unbiasedness condition $\mathbf{E}g(x_{t},i_{t})=f^{\prime}(x_{t})$ (recall that $\mathbf{E}$ is the relevant-past conditional expectation with respect to $i_{t}$).

Unfortunately, in the parallel setting, we manipulate stale, inconsistent reads of shared parameters and thus we do not have such a straightforward relationship. Instead, Mania et al. (2015) proposed to separate $\hat{x}_{t}$, the actual value read by a core to compute an update, with $x_{t}$, a “virtual iterate” that we can analyze and is defined by the update equation: $x_{t+1}:=x_{t}-\gamma g(\hat{x}_{t},i_{t}).$ We can thus interpret $\hat{x}_{t}$ as a noisy (perturbed) version of $x_{t}$ due to the effect of asynchrony. In the specific case of (Sparse) Saga, we have to add the additional read memory argument $\hat{\alpha}^{t}$ to our update:

We formalize the precise meaning of $x_{t}$ and $\hat{x}_{t}$ in the next section. We first note that all the papers mentioned in the related work section that analyzed asynchronous parallel randomized algorithms assumed that the following unbiasedness condition holds:

This condition is at the heart of most convergence proofs for randomized optimization methods.A notable exception is Sag (Le Roux et al., 2012) which has biased updates, yielding a significantly more complex convergence proof. Making Sag unbiased leads to Saga (Defazio et al., 2014) and a much simpler proof. Mania et al. (2015) correctly pointed out that most of the literature thus made the often implicit assumption that $i_{t}$ is independent of $\hat{x}_{t}$. But as we explain below, this assumption is incompatible with a non-uniform asynchronous model in the analysis approach used in most of the recent literature.

2 On the Difficulty of Labeling the Iterates

The main idea of the perturbed iterate framework is to use this handle on $\hat{x}_{t}-x_{t}$ to analyze the convergence for $x_{t}$. In this paper, we can instead give directly the convergence of $\hat{x}_{t}$, and so unlike in Mania et al. (2015), we do not require that there exists a $T$ such that $x_{T}$ lives in shared memory.

3 Analysis setup

We describe Asaga, a sparse asynchronous parallel implementation of Sparse Saga, in Algorithm 1 in the theoretical form that we analyze, and in Algorithm 2 as its practical implementation. Before stating its convergence, we highlight some properties of Algorithm 1 and make one central assumption.

Given the “after read” global ordering, $i_{r}$ is independent of $\hat{x}_{t}$ $\forall r\geq t$.

We enforce the independence for $r=t$ in Algorithm 1 by having the core read all the shared data parameters and historical gradients before starting their iterations. Although this is too expensive to be practical if the data is sparse, this is required by the theoretical Algorithm 1 that we can analyze. As Mania et al. (2015) stress, this independence property is assumed in most of the parallel optimization literature. The independence for $r>t$ is a consequence of using the “after read” global ordering instead of the “before read” one.

The update, $g_{t}:=g(\hat{x}_{t},\hat{\alpha}^{t},i_{t})$, is an unbiased estimator of the true gradient at $\hat{x}_{t}$ (i.e. (5) yields (6) in conditional expectation).

The shared parameter coordinate update of $[x]_{v}$ on line 11 is atomic.

Since our updates are additions, this means that there are no overwrites, even when several cores compete for the same resources. In practice, this is enforced by using compare-and-swap semantics, which are heavily optimized at the processor level and have minimal overhead. Our experiments with non-thread safe algorithms (i.e. where this property is not verified, see Figure 6 of Appendix G) show that compare-and-swap is necessary to optimize to high accuracy.

Finally, as is standard in the literature, we make an assumption on the maximum delay that asynchrony can cause – this is the partially asynchronous setting as defined in Bertsekas and Tsitsiklis (1989):

We assume that there exists a uniform bound, called $\tau$, on the maximum number of iterations that can overlap together. We say that iterations $r$ and $t$ overlap if at some point they are processed concurrently. One iteration is being processed from the start of the reading of the shared parameters to the end of the writing of its update. The bound $\tau$ means that iterations $r$ cannot overlap with iteration $t$ for $r\geq t+\tau+1$, and thus that every coordinate update from iteration $t$ is successfully written to memory before the iteration $t+\tau+1$ starts.

Our result will give us conditions on $\tau$ subject to which we have linear speedups. $\tau$ is usually seen as a proxy for $p$, the number of cores (which lowerbounds it). However, though $\tau$ appears to depend linearly on $p$, it actually depends on several other factors (notably the data sparsity distribution) and can be orders of magnitude bigger than $p$ in real-life experiments. We can upper bound $\tau$ by $(p-1)R$, where $R$ is the ratio of the maximum over the minimum iteration time (which encompasses theoretical aspects as well as hardware overhead). More details can be found in Appendix E.

By using the overlap Assumption 1 in the expression (8) for the iterates, we obtain the following explicit effect of asynchrony that is crucially used in our proof:

where $G_{u}^{t}$ are $d\times d$ diagonal matrices with terms in $\{0,+1\}$. We know from our definition of $t$ and $x_{t}$ that every update in $\hat{x}_{t}$ is already in $x_{t}$ – this is the case. Conversely, some updates might be late: this is the $+1$ case. $\hat{x}_{t}$ may be lacking some updates from the “past" in some sense, whereas given our global ordering definition, it cannot contain updates from the “future".

4 Convergence and speedup results

We now state our main theoretical results. We give an outline of the proof in Section 3.5 and its full details in Appendix C. We first define a notion of problem sparsity, as it will appear in our results.

As in Niu et al. (2011), we introduce $\Delta_{r}:=\max_{v=1..d}|\{i:v\in S_{i}\}|$. $\Delta_{r}$ is the maximum right-degree in the bipartite graph of the factors and the dimensions, i.e., the maximum number of data points with a specific feature. For succinctness, we also define $\Delta:=\Delta_{r}/n$. We have $1\leq\Delta_{r}\leq n$, and hence $1/n\leq\Delta\leq 1$.

Suppose $\tau<n/10$.Asaga can actually converge for any $\tau$, but the maximum step size then has a term of $\exp(\tau/n)$ in the denominator with much worse constants. See Appendix C.8. Let

This result is very close to Saga’s original convergence theorem, but with the maximum step size divided by an extra $1+\tau\sqrt{\Delta}$ factor. Referring to Hofmann et al. (2015) and our own Theorem 1, the rate factor for Saga is $\min\{1/n,1/\kappa\}$ up to a constant factor. Comparing this rate with Theorem 2 and inferring the conditions on the maximum step size $a^{*}(\tau)$, we get the following conditions on the overlap $\tau$ for Asaga to have the same rate as Saga (comparing upper bounds).

Suppose $\tau\leq\mathcal{O}(n)$ and $\tau\leq\mathcal{O}({\scriptstyle\frac{1}{\sqrt{\Delta}}}\max\{1,\frac{n}{\kappa}\})$. Then using the step size $\gamma=\nicefrac{{a^{*}(\tau)}}{{L}}\,$ from (10), Asaga converges geometrically with rate factor $\Omega(\min\{\frac{1}{n},\frac{1}{\kappa}\})$ (similar to Saga), and is thus linearly faster than its sequential counterpart up to a constant factor. Moreover, if $\tau\leq\mathcal{O}(\frac{1}{\sqrt{\Delta}})$, then a universal step size of $\Theta(\frac{1}{L})$ can be used for Asaga to be adaptive to local strong convexity with a similar rate to Saga (i.e., knowledge of $\kappa$ is not required).

Interestingly, in the well-conditioned regime ($n>\kappa$, where Saga enjoys a range of stepsizes which all give the same contraction ratio), Asaga can get the same rate as Saga even in the non-sparse regime ($\Delta=1$) for $\tau<\mathcal{O}(n/\kappa)$. This is in contrast to the previous work on asynchronous incremental gradient methods which required some kind of sparsity to get a theoretical linear speedup over their sequential counterpart (Niu et al., 2011; Mania et al., 2015). In the ill-conditioned regime ($\kappa>n$), sparsity is required for a linear speedup, with a bound on $\tau$ of $\mathcal{O}(\sqrt{n})$ in the best-case (though degenerate) scenario where $\Delta=1/n$.

We give the first convergence analysis for an asynchronous parallel version of Saga (note that Reddi et al. (2015) only covers an epoch based version of Saga with random stopping times, a fairly different algorithm).

Theorem 2 can be directly extended to a parallel extension of the Svrg version from Hofmann et al. (2015), which is adaptive to the local strong convexity with similar rates (see Appendix C.2).

In contrast to the parallel Svrg analysis from Reddi et al. (2015, Thm. 2), our proof technique handles inconsistent reads and a non-uniform processing speed across $f_{i}$’s. Our bounds are similar (noting that $\Delta$ is equivalent to theirs), except for the adaptivity to local strong convexity: Asaga does not need to know $\kappa$ for optimal performance, contrary to parallel Svrg (see App. C.2 for more details).

In contrast to the Svrg analysis from Mania et al. (2015, Thm. 14), we obtain a better dependence on the condition number in our rate ($1/\kappa$ vs. $1/\kappa^{2}$ for them) and on the sparsity (they get $\tau\leq\mathcal{O}(\Delta^{\nicefrac{{-1}}{{3}}})$), while we remove their gradient bound assumption. We also give our convergence guarantee on $\hat{x}_{t}$ during the algorithm, whereas they only bound the error for the “last” iterate $x_{T}$.

5 Proof outline

We give here the outline of our proof. Its full details can be found in Appendix C.

Let $g_{t}:=g(\hat{x}_{t},\hat{\alpha}^{t},i_{t})$. By expanding the update equation (5) defining the virtual iterate $x_{t+1}$ and introducing $\hat{x}_{t}$ in the inner product term, we get:

With further manipulations on the expectation of (11), including the use of the standard inequality $\|a+b\|^{2}\leq 2\|a\|^{2}+2\|b\|^{2}$ (see Section C.3), we obtain our basic recursive contraction inequality:

The rest of the proof then proceeds as follows:

We bound the first term by $4Le_{t}$ using Hofmann et al. (2015, Equation (8)). To express the second term in terms of past suboptimalities, we note that it can be seen as an expectation of past first terms with an adequate probability distribution which we derive and bound.

By substituting Lemma 2 into Lemma 1, we get a master contraction inequality (28) in terms of $a_{t+1}$, $a_{t}$ and $e_{u},u\leq t$.

We define a novel Lyapunov function $\mathcal{L}_{t}=\sum_{u=0}^{t}(1-\rho)^{t-u}a_{u}$ and manipulate the master inequality to show that $\mathcal{L}_{t}$ is bounded by a contraction, subject to a maximum step size condition on $\gamma$ (given in Lemma 3, see Appendix C.1).

Finally, we unroll the Lyapunov inequality to get the convergence Theorem 2.

Empirical results

We now present the main results of our empirical comparison of asynchronous Saga, Svrg and Hogwild. Additional results, including convergence and speedup figures with respect to the number of iteration and measures on the $\tau$ constant are available in the appendix.

Datasets. We consider two sparse datasets: RCV1 (Lewis et al., 2004) and URL (Ma et al., 2009); and a dense one, Covtype (Collobert et al., 2002), with statistics listed in the table below. As in Le Roux et al. (2012), Covtype is standardized, thus $100\%$ dense. $\Delta$ is $\mathcal{O}(1)$ in all datasets, hence not very insightful when relating it to our theoretical results. Deriving a less coarse sparsity bound remains an open problem.

Hardware and software. Experiments were run on a 40-core machine with 384GB of memory. All algorithms were implemented in Scala. We chose this high-level language despite its typical 20x slowdown compared to C (when using standard libraries, see Appendix G.2) because our primary concern was that the code may easily be reused and extended for research purposes (to this end, we have made all our code available at https://github.com/RemiLeblond/ASAGA).

2 Implementation details

Comparison with the theoretical algorithm. The algorithm we used in the experiments is fully detailed in Algorithm 2. There are two differences with Algorithm 1. First, in the implementation we pick $i_{t}$ at random before we read data. This enables us to only read the necessary data for a given iteration (i.e. $[\hat{x}_{t}]_{S_{i}},[\hat{\alpha}_{i}^{t}],[\bar{\alpha}^{t}]_{S_{i}}$). Although this violates Property 1, it still performs well in practice.

Second, we maintain $\bar{\alpha}^{t}$ in memory. This saves the cost of recomputing it at every iteration (which we can no longer do since we only read a subset data). Again, in practice the implemented algorithm enjoys good performance. But this design choice raises a subtle point: the update is not guaranteed to be unbiased in this setup (see Appendix G.3 for more details).

3 Results

We first compare three different asynchronous variants of stochastic gradient methods on the aforementioned datasets: Asaga, presented in this work, Kromagnon, the asynchronous sparse Svrg method described in Mania et al. (2015) and Hogwild (Niu et al., 2011). Each method had its step size chosen so as to give the fastest convergence (up to $10^{-3}$ in the special case of Hogwild). The results can be seen in Figure 1(a): for each method we consider its asynchronous version with both one (hence sequential) and ten processors. This figure reveals that the asynchronous version offers a significant speedup over its sequential counterpart.

We then examine the speedup relative to the increase in the number of cores. The speedup is measured as time to achieve a suboptimality of $10^{-5}$ ($10^{-3}$ for Hogwild) with one core divided by time to achieve the same suboptimality with several cores, averaged over 3 runs. Again, we choose step size leading to fastest convergence (see Appendix G.2 for information about the step sizes). Results are displayed in Figure 1(b).

As predicted by our theory, we observe linear “theoretical” speedups (i.e. in terms of number of iterations, see Appendix D.2). However, with respect to running time, the speedups seem to taper off after $20$ cores. This phenomenon can be explained by the fact that our hardware model is by necessity a simplification of reality. As noted in Duchi et al. (2015), in a modern machine there is no such thing as shared memory. Each core has its own levels of cache (L1, L2, L3) in addition to RAM. The more cores are used, the lower in the memory stack information goes and the slower it gets. More experimentation is needed to quantify that effect and potentially increase performance.

Conclusions and future work

We have described Asaga, a novel sparse and fully asynchronous variant of the incremental gradient algorithm Saga. Building on the recently proposed “perturbed iterate” framework, we have introduced a novel analysis of the algorithm and proven that under mild conditions Asaga is linearly faster than Saga. Our empirical benchmarks confirm speedups up to 10x.

Our proof technique accommodates more realistic settings than is usually the case in the literature (e.g. inconsistent reads/writes and an unbounded gradient); we obtain tighter conditions than in previous work. In particular, we show that sparsity is not always necessary to get linear speedups. Further, we have proposed a novel perspective to clarify an important technical issue present in most of the recent convergence rate proofs for asynchronous parallel optimization algorithms.

Schmidt et al. (2016) have shown that Sag enjoys much improved performance when combined with non-uniform sampling and line-search. We have also noticed that our $\Delta_{r}$ constant (being essentially a maximum) sometimes fails to accurately represent the full sparsity distribution of our datasets. Finally, while our algorithm can be directly ported to a distributed master-worker architecture, its communication pattern would have to be optimized to avoid prohibitive costs. Limiting communications can be interpreted as artificially increasing the delay, yielding an interesting trade-off between delay influence and communication costs.

A final interesting direction for future analysis is the further exploration of the $\tau$ term, which we have shown encompasses more complexity than previously thought.

We would like to thank Xinghao Pan for sharing with us their implementation of Kromagnon, as well as Alberto Chiappa for spotting a typo in the proof. This work was partially supported by a Google Research Award and the MSR-Inria Joint Center. FP acknowledges financial support from from the chaire Économie des nouvelles données with the data science joint research initiative with the fonds AXA pour la recherche.

References

In Appendix A, we give a simple example illustrating why the “After Write” approach can break the crucial unbiasedness condition (6) needed for standard convergence proofs.

In Appendix B, we adapt the proof from Hofmann et al. (2015) to prove Theorem 1, our convergence result for serial Sparse Saga.

In Appendix C, we first give a detailed outline and then the complete details for the proof of convergence for Asaga (Theorem 2) as well as its linear speedup regimes (Corollary 3).

In Appendix D, we analyze additional experimental results, including a comparison of serial Saga algorithms and a look at “theoretical speedups” for Asaga.

In Appendix E, we take a closer look at the $\tau$ constant. We argue that it encompasses more complexity than is usually implied in the literature, as additional results that we present indicate.

In Appendix F, we compare the lagged updates implementation of Saga with our sparse algorithm, and explain why adapting the former to the asynchronous setting is difficult.

In Appendix G, we give additional details about the datasets and our implementation.

Appendix A Problematic Example for the “After Write” Approach

To understand this subtle point better, note that in this very simple example, $i_{0}$ and $i_{1}$ are not independent. We can show that $P(i_{1}=2\mid i_{0}=2)=1$. They share dependency through the labeling assignment.

The only way we can think to resolve this issue and ensure unbiasedness in the “after write” framework is to assume that the computation time for the algorithm running on a core is independent of the sample $i$ chosen. This assumption seems overly strong in the context of potentially heterogeneous factors $f_{i}$’s, and is thus a fundamental flaw in the “after write” framework that has mostly been ignored in the recent asynchronous optimization literature.

We note that Bertsekas and Tsitsiklis (1989) briefly discussed this issue in Section 7.8.3 of their book, stressing that their analysis for SGD required that the scheduling of computation was independent from the randomness from SGD, but they did not offer any solution if this assumption was not satisfied. Both the “before read” labeling from Mania et al. (2015) and our proposed “after read” labeling resolve this issue.

Appendix B Proof of Theorem 1

As we will heavily reuse the proof technique from Hofmann et al. (2015), we start by giving its sketch.

First, the authors combine classical strong convexity and Lipschitz inequalities to derive the inequality Hofmann et al. (2015, Lemma 1):

This gives a contraction term, as well as two additional terms; $2\gamma^{2}\mathbf{E}\|\alpha_{i}-f^{\prime}_{i}(x^{*})\|^{2}$ is a positive variance term, but (4\gamma^{2}L-2\gamma)\big{(}f(x)-f(x^{*})\big{)} is a negative suboptimality term (provided $\gamma$ is small enough). The suboptimality term can then be used to cancel the variance one.

Second, the authors use a classical smoothness upper bound to control the variance term and relate it to the suboptimality. However, since the $\alpha_{i}$ are partial gradients computed at previous time steps, the upper bounds of the variance involve suboptimality at previous time steps, which are not directly relatable to the current suboptimality.

Third, to circumvent this issue, a Lyapunov function is defined to encompass both current and past terms. To finish the proof, Hofmann et al. (2015) show that the Lyapunov function is a contraction.

Fortunately, we can reuse most of the proof from Hofmann et al. (2015) to show that Sparse Saga converges at the same rate as regular Saga. In fact, once we establish that Hofmann et al. (2015, Lemma 1) is still verified we are done.

To prove this, we derive close variants of equations $(6)$ and $(9)$ in their paper, which we remind the reader of here:

We first show that the update estimator is unbiased. The estimator is unbiased if:

where $e_{v}$ is the vector whose only nonzero component is the $v$ component which is equal to $1$.

By definition, $\sum_{i|v\in S_{i}}1=np_{v},$ which gives us Equation (14).

We define $\bar{\alpha}_{i}:=\alpha_{i}-D_{i}\bar{\alpha}$ (contrary to Hofmann et al. (2015) where the authors define $\bar{\alpha}_{i}:=\alpha_{i}-\bar{\alpha}$ since they do not concern themselves with sparsity). Using the inequality $\|a+b\|^{2}\leq 2\|a\|^{2}+2\|b\|^{2}$, we get:

which is our equivalent to Hofmann et al. (2015, Eq.(6)), where only our definition of $\bar{\alpha}_{i}$ differs.

We want to prove Hofmann et al. (2015, Eq.(9)):

Let $D_{\neg i}:=P_{S_{i}^{c}}D$; we then have the orthogonal decomposition $D\alpha=D_{i}\alpha+D_{\neg i}\alpha$ with $D_{i}\alpha\perp D_{\neg i}\alpha$, as they have disjoint support. We now use the orthogonality of $D_{\neg i}\alpha$ with any vector with support in $S_{i}$ to simplify the expression (17) as follows:

This is our version of Hofmann et al. (2015, Equation (9)), which finishes the proof of Hofmann et al. (2015, Lemma 1). The rest of the proof from Hofmann et al. (2015) can then be reused without modification to obtain Theorem 1. ∎

Appendix C Proof of Theorem 2 and Corollary 3

We first give a detailed outline of the proof. The complete proof is given in the rest of Appendix C.

Let $g_{t}:=g(\hat{x}_{t},\hat{\alpha}^{t},i_{t})$. From the update equation (5) defining the virtual iterate $x_{t+1}$, the perturbed iterate framework (Mania et al., 2015) gives:

Note that we have introduced $\hat{x}_{t}$ in the inner product because $g_{t}$ is a function of $\hat{x}_{t}$, not $x_{t}$.

With further manipulations on the expectation of (C.1), including the use of the standard inequality $\|a+b\|^{2}\leq 2\|a\|^{2}+2\|b\|^{2}$ (see Section C.3), we obtain our basic recursive contraction inequality:

Inequality (23) is a midway point between the one derived in the proof of Lemma 1 in Hofmann et al. (2015) and Equation (2.5) in Mania et al. (2015), because we use the tighter strong convexity bound (22) than in the latter (giving us the important extra term $-2\gamma e_{t}$).

The rest of the proof then proceeds as follows:

By substituting Lemma 2 into Lemma 1, we get a master contraction inequality (28) in terms of $a_{t+1}$, $a_{t}$ and $e_{u},u\leq t$.

We define a novel Lyapunov function $\mathcal{L}_{t}=\sum_{u=0}^{t}(1-\rho)^{t-u}a_{u}$ and manipulate the master inequality to show that $\mathcal{L}_{t}$ is bounded by a contraction, subject to a maximum step size condition on $\gamma$ (given in Lemma 3 below).

Finally, we unroll the Lyapunov inequality to get the convergence Theorem 2.

We list the key lemmas below with their proof sketch, and give the detailed proof in the later sections of Appendix C.

where $C_{1}:=1+\sqrt{\Delta}\tau$ and $C_{2}:=\sqrt{\Delta}+\gamma\mu C_{1}$.

From our Sparse Saga proof we know that (see Appendix B):

From our algorithm, we know that each dimension of the memory vector $[\hat{\alpha}_{i}]_{v}$ contains a partial gradient computed at some point in the past $[f^{\prime}_{i}(\hat{x}_{u_{i,v}^{t}})]_{v}$More precisely: $\forall t,i,v\hskip 5.0pt\exists u_{i,v}^{t}<t$ s.t. $[\hat{\alpha}_{i}^{t}]_{v}=[f^{\prime}_{i}(\hat{x}_{u_{i,v}^{t}})]_{v}$. (unless $u=0$, in which case we replace the partial gradient with $\alpha_{i}^{0}$). We then derive bounds on $P(u_{i,v}^{t}=u)$ and sum on all possible $u$. Together with clever conditioning, we obtain Lemma 2 (see Section C.5).

Let $H_{t}$ be defined as ${H_{t}:=\sum_{u=1}^{t-1}(1-\frac{1}{n})^{(t-2\tau-u-1)_{+}}e_{u}}$. Then, by setting (26) into Lemma 1, we get (see Section C.6):

We now have the beginning of a contraction with additional positive terms which all converge to as we near the optimum, as well as our classical negative suboptimality term. This is not unusual in the variance reduction literature. One successful approach in the sequential case is then to define a Lyapunov function which encompasses all terms and is a true contraction (see Defazio et al. (2014); Hofmann et al. (2015)). We emulate this solution here. However, while all terms in the sequential case only depend on the current iterate, $t$, in the parallel case we have terms “from the past” in our inequality. To resolve this issue, we define a more involved Lyapunov function which also encompasses past iterates:

where $\rho$ is a target contraction rate that we define later.

Using the master inequality (28), we get (see Appendix C.7):

The aim is to prove that $\mathcal{L}_{t}$ is bounded by a contraction. We have two promising terms at the beginning of the inequality, and then we need to handle the last term. Basically, we can rearrange the sums in (28) to expose a simple sum of $e_{u}$ multiplied by factors $r_{u}^{t}$.

Under specific conditions on $\rho$ and $\gamma$, we can prove that $r_{u}^{t}$ is negative for all $u\geq 1$, which coupled with the fact that each $e_{u}$ is positive means that we can safely drop the sum term from the inequality. The $r_{0}^{t}$ term is a bit trickier and is handled separately.

Suppose $\tau<n/10$ and $\rho\leq 1/4n$. If

then for all $u\geq 1$, the $r_{u}^{t}$ from (C.1) verify:

We obtain this result after carefully deriving the $r_{u}^{t}$ terms. We find a second-order polynomial inequality in $\gamma$, which we simplify down to (31) (see Appendix C.8).

We can then finish the argument to bound the suboptimality error $e_{t}$. We have:

We have two linearly contracting terms. The sum contracts linearly with the worst rate between the two (the smallest geometric rate factor). If we define $\rho^{*}:=\nu\min(\rho,\gamma\mu/2)$, with $0<\nu<1$,$\nu$ is introduced to circumvent the problematic case where $\rho$ and $\gamma\mu/2$ are too close together. then we get:

where $\eta:=\frac{1-M}{1-\rho^{*}}$ with $M:=\max(\rho,\gamma\mu/2)$. Our geometric rate factor is thus $\rho^{*}$ (see Appendix C.9).

C.2 Extension to Svrg

Our proof can easily be adapted to accommodate the Svrg variant introduced in Hofmann et al. (2015), which is closer to Saga than the initial Svrg algorithm and which is adaptive to local strong convexity (it does not require the inner loop epoch size $m=\Omega(\kappa)$ as a hyperparameter). In this variant, instead of computing a full gradient every $m$ iterations, a random binary variable $U$ with probability $P(U=1)=1/n$ is sampled at the beginning of every iteration to determine whether a full gradient is computed or a normal Svrg step is made. If $U=1$, then a full gradient is computed. Otherwise the algorithm takes a normal inner Svrg step.Note that the parallel implementation is not very straightforward, as it requires a way to communicate to cores when they should start computing a batch gradient instead of inner steps.

To prove convergence, all one has to do is to modify Lemma 2 very slightly (the only difference is that the $(t-2\tau-u-1)_{+}$ exponent is replaced by $(t-u)$ and the rest of the proof can be used as is). The justification for this small tweak is that the batch steps in Svrg are fully synchronized. More details can be found in Section C.5 (see footnote 16).

By using our “after read” labeling, we were also able to derive a convergence and speedup proof for the original Svrg algorithm, but the proof technique diverges after Lemma 1. This is beyond the scope of this paper, so we omit it here. Using the “after read” labeling and a different proof technique from Mania et al. (2015)), we obtain an epoch size in $\mathcal{O}(\kappa)$ instead of $\mathcal{O}(\kappa^{2})$ and a dependency in our overlap bound in $\mathcal{O}(\Delta^{-1/2})$ instead of $\mathcal{O}(\Delta^{-1/3})$.

There are several scenarios in which Asaga can be practically advantageous over its closely related cousin, asynchronous Svrg (note though that “asynchronous” Svrg still requires a synchronization step to compute the full gradients).

First, while Saga trades memory for less computation, in the case of generalized linear models the memory cost can be reduced to $\mathcal{O}(n)$, which is the same as for Svrg. This is of course also true for their asynchronous counterparts.

Second, as Asaga does not require any synchronization steps, it is better suited to heterogeneous computing environments (where cores have different clock speeds or are shared with other applications).

Finally, Asaga does not require knowing the condition number $\kappa$ for optimal convergence in the sparse regime. It is thus adaptive to local strong convexity, whereas Svrg is not. Indeed, Svrg and its asynchronous variant require setting an additional hyper-parameter – the epoch size $m$ – which needs to be at least $\Omega(\kappa)$ for convergence but yields a slower effective convergence rate than Asaga if it is set much bigger than $\kappa$. Svrg thus requires tuning this additional hyper-parameter or running the risk of either slower convergence (if the epoch size chosen is much bigger than the condition number) or even not converging at all (if $m$ is chosen to be much smaller than $\kappa$).Note that as Saga (and contrary to the original Svrg), the Svrg variant from Hofmann et al. (2015) does not require knowledge of $\kappa$ and is thus adaptive to local strong convexity, which carries over to its asynchronous adaptation.

C.3 Initial recursive inequality derivation

We start by proving Equation (23). Let $g_{t}:=g(\hat{x}_{t},\hat{\alpha}^{t},i_{t})$. From (5), we get:

In order to prove Equation (23), we need to bound the $-2\gamma\langle\hat{x}_{t}-x^{*},g_{t}\rangle$ term. Thanks to Property 1, we can write:

We can now use a classical strong convexity bound as well as a squared triangle inequality to get:

Putting it all together, we get the initial recursive inequality (23), rewritten here explicitly:

C.4 Proof of Lemma 1

We start by proving a relevant property of $\Delta$, which enables us to derive an essential inequality for both these terms, given in Proposition 1 below. We reuse the sparsity constant introduced in Reddi et al. (2015) and relate it to the one we have defined earlier, $\Delta_{r}$:

Let $D$ be the smallest constant such that:

where $\delta_{v}:=\mathbf{card}(i\mid v\in S_{i})$.

Since $D$ is the minimum constant satisfying this inequality, we have:

We need to find $x$ such that it maximizes the right-hand side term. Note that the vector $([x]_{v}^{2}/\|x\|^{2})_{v=1..d}$ is in the unit probability simplex, which means that an equivalent problem is the maximization over all convex combinations of $(\delta_{v})_{v=1..d}$. This maximum is found by putting all the weight on the maximum $\delta_{v}$, which is $\Delta_{r}$ by definition.

This means that $\Delta=\Delta_{r}/n$ is indeed the smallest constant satisfying (39). ∎

Let $u\neq t$. Without loss of generality, $u<t$.One only has to switch $u$ and $t$ if $u>t$. Then:

Thanks to the expansion for $\hat{x}_{t}-x_{t}$ (9), we get:

Using (44) from Proposition 1, we have that for $u\neq v$:

By taking the expectation and using (47), we get:

C.5 Proof of Lemma 2

We now derive our bound on $g_{t}$ with respect to suboptimality. From Appendix B, we know that:

N. B.: In the following, $i_{t}$ is a random variable picked uniformly at random in $\{1,...,n\}$, whereas $i$ is a fixed constant.

Now, with $i$ fixed, let $u_{i,l}^{t}$ be the time of the iterate last used to write the $[\hat{\alpha}_{i}^{t}]_{l}$ quantity, i.e. $[\hat{\alpha}_{i}^{t}]_{l}=[f^{\prime}_{i}(\hat{x}_{u_{i,l}^{t}})]_{l}$. We knowIn the case where $u=0$, one would have to replace the partial gradient with $\alpha_{i}^{0}$. We omit this special case here for clarity of exposition. that $0\leq u_{i,l}^{t}\leq t-1$. To use this information, we first need to split $\hat{\alpha}_{i}$ along its dimensions to handle the possible inconsistencies among them:

We will now rewrite the indicator so as to obtain independent events from the rest of the equality. This will enable us to distribute the expectation. Suppose $u>0$ ($u=0$ is a special case which we will handle afterwards). $\{u_{i,l}^{t}=u\}$ requires two things:

at time $u$, $i$ was picked uniformly at random,

(roughly) $i$ was not picked again between $u$ and $t$.

Note that the third line used the crucial independence assumption $i_{v}\perp\!\!\!\perp\hat{x}_{u},\forall v\geq u$ arising from our “After Read” ordering. Summing over all dimensions $l$, we then get:

Plugging (56) and (C.5) into (50), we get Lemma 2:

C.6 Master inequality derivation

If we define $H_{t}:=\sum_{u=1}^{t-1}(1-\frac{1}{n})^{(t-2\tau-u-1)_{+}}e_{u}$, then we get:

C.7 Lyapunov function and associated recursive inequality

We define $\mathcal{L}_{t}:=\sum_{u=0}^{t}(1-\rho)^{t-u}a_{u}$ for some target contraction rate $\rho<1$ to be defined later. We have:

We now use our new bound on $a_{t+1}$, (60):

We can now rearrange the sums to expose a simple sum of $e_{u}$ multiplied by factors $r_{u}^{t}$:

C.8 Proof of Lemma 3

We want to make explicit what conditions on $\rho$ and $\gamma$ are necessary to ensure that $r_{u}^{t}$ is negative for all $u\geq 1$. Since each $e_{u}$ is positive, we will then be able to safely drop the sum term from the inequality. The $r_{0}^{t}$ term is a bit trickier and is handled separately. Indeed, trying to enforce that $r_{0}^{t}$ is negative results in a significantly worse condition on $\gamma$ and eventually a convergence rate smaller by a factor of $n$ than our final result. Instead, we handle this term directly in the Lyapunov function.

Let’s now make the multiplying factor explicit. We assume $u\geq 1$.

We split $r_{u}^{t}$ into five parts coming from (C.7):

$r_{1}$, the part coming from the $-2\gamma e_{u}$ terms;

$r_{2}$, coming from $4L\gamma^{2}C_{1}e_{u}$;

$r_{3}$, coming from $\frac{4L\gamma^{2}C_{1}}{n}H_{u}$;

$r_{4}$, coming from $4L\gamma^{2}C_{2}\sum_{v=(u-\tau)_{+}}^{u-1}e_{v}$;

$r_{5}$, coming from $\frac{4L\gamma^{2}C_{2}}{n}\sum_{v=(u-\tau)_{+}}^{u-1}H_{v}$.

$r_{1}$ is easy to derive. Each of these terms appears only in one inequality. So for $u$ at time $t$, the term is:

For much the same reasons, $r_{2}$ is also easy to derive and is:

$r_{3}$ is a bit trickier, because for a given $v>0$ there are several $H_{u}$ which contain $e_{v}$. The key insight is that we can rewrite our double sum in the following manner:

Note that we have bounded the $\min(t,v+2\tau)$ term by $v+2\tau$ in the first sub-sum, effectively adding more positive terms.

Finally we compute $r_{5}$ which is the most complicated term. Indeed, to find the factor of $e_{w}$ for a given $w>0$, one has to compute a triple sum, $\sum_{u=0}^{t}(1-\rho)^{t-u}\sum_{v=(u-\tau)_{+}}^{u-1}H_{v}$. We start by computing the factor of $e_{w}$ in the inner double sum, $\sum_{v=(u-\tau)_{+}}^{u-1}H_{v}$.

Now there are at most $\tau$ terms for each $e_{w}$. If $w\leq u-3\tau-1$, then the exponent is positive in every term and it is always bigger than $u-3\tau-1-w$, which means we can bound the sum by $\tau(1-\frac{1}{n})^{u-3\tau-1-w}$. Otherwise we can simply bound the sum by $\tau$. We get:

By combining the five terms together ((64), (65), (68), (70) and (C.8)), we get that $\forall u$ s.t. $1\leq u\leq t$:

$r_{1}$, the part coming from the $-2\gamma e_{u}$ terms;

$r_{2}$, coming from $4L\gamma^{2}C_{1}e_{u}$;

$r_{4}$, coming from $4L\gamma^{2}C_{2}\sum_{v=(u-\tau)_{+}}^{u-1}e_{v}$;

We have $r_{1}=-2\gamma(1-\rho)^{t}$ and $r_{2}=4L\gamma^{2}C_{1}(1-\rho)^{t}$.

We have already computed $r_{4}$ for $u>0$ and the computation is exactly the same for $u=0$. $r_{4}\leq(1-\rho)^{t}4L\gamma^{2}C_{2}\tau(1-\rho)^{-\tau}$ .

Putting it all together, we get that: $\forall t\geq 0$

We need all $r_{u}^{t},u\geq 1$ to be negative so we can safely drop them from (63). Note that for every $u$, this is the same condition. We will reduce that condition to a second-order polynomial sign condition. We also remark that since $\gamma\geq 0$, we can upper bound our terms in $\gamma$ and $\gamma^{2}$ in this upcoming polynomial, which will give us sufficient conditions for convergence.

Now, as $\gamma$ is part of $C_{2}$, we need to expand it once more to find our conditions. We have:

Dividing the bracket in (74) by $\gamma$ and rearranging as a second degree polynomial, we get the condition:

The discriminant of this polynomial is always positive, so $\gamma$ needs to be between its two roots. The smallest is negative, so the condition is not relevant to our case (where $\gamma>0$). By solving analytically for the positive root $\phi$, we get an upper bound condition on $\gamma$ that can be used for any overlap $\tau$ and guarantee convergence. Unfortunately, for large $\tau$, the upper bound becomes exponentially small because of the presence of $\tau$ in the exponent in (80). More specifically, by using the bound $1/(1-\rho)\leq\exp(2\rho)$This bound can be derived from the inequality $(1-x/2)\geq\exp(-x)$ which is valid for $0\leq x\leq 1.59$. and thus $(1-\rho)^{-\tau}\leq\exp(2\tau\rho)$ in (80), we would obtain factors of the form $\exp(\tau/n)$ in the denominator for the root $\phi$ (recall that $\rho<1/n$).

Our Lemma 3 is derived instead under the assumption that $\tau\leq\mathcal{O}(n)$, with the constants chosen in order to make the condition (80) more interpretable and to relate our convergence result with the standard SAGA convergence (see Theorem 1). As explained in Appendix E, the assumption that $\tau\leq\mathcal{O}(n)$ appears reasonable in practice. First, by using Bernoulli’s inequality, we have:

To get manageable constants, we make the following slightly more restrictive assumptions on the target rate $\rho$Note that we already expected $\rho<1/n$. and overlap $\tau$:This bound on $\tau$ is reasonable in practice, see Appendix E.

We can now upper bound loosely the three terms in brackets appearing in (80) as follows:

By plugging (88)–(90) into (80), we get the simpler sufficient condition on $\gamma$:

We simplify it further by using the inequality:This inequality can be derived by using the concavity property $f(y)\leq f(x)+(y-x)f^{\prime}(x)$ on the differentiable concave function $f(x)=\sqrt{x}$ with $y=1$.

Using (93) in (92), and recalling that $\kappa:=L/\mu$, we get:

Since $\frac{\tau}{C_{1}}=\frac{\tau}{1+\sqrt{\Delta}\tau}\leq\min(\tau,\frac{1}{\sqrt{\Delta}})$, we get that a sufficient condition on our stepsize is:

Subject to our conditions on $\gamma$, $\rho$ and $\tau$, we then have that: $r_{u}^{t}\leq 0\ \text{for all}\ u\ \text{s.t.}\ 1\leq u\leq t$. This means we can rewrite (63) as:

We make $e_{t}$ appear on the left side of (96), by adding $\gamma$ to $r_{t}^{t}$ in (63):We could use any multiplier from to $2\gamma$, but choose $\gamma$ for simplicity. For this reason and because our analysis of the $r_{t}^{t}$ term was loose, we could derive a tighter bound, but it does not change the leading terms.

We now require the stronger property that $\gamma+r_{t}^{t}\leq 0$, which translates to replacing $-2\gamma$ with $-\gamma$ in (74):

We can easily derive a new stronger condition on $\gamma$ under which we can drop all the $e_{u},u>0$ terms in (97):

C.9 Proof of Theorem 2

We continue with the assumptions of Lemma 3 which gave us (100). Thanks to (79), we can also rewrite $r_{0}^{t}\leq(1-\rho)^{t+1}A$, where $A$ is a constant which depends on $n$, $\Delta$, $\gamma$ and $L$ but is finite and crucially does not depend on $t$. In fact, by reusing similar arguments as in C.8, we can show the bound $A\leq\gamma n$ under the assumptions of Lemma 3 (including $\gamma\leq\gamma^{*}$).In particular, note that $e_{0}$ does not appear in the definition of $A$ because it turns out that the parenthesis group multiplying $e_{0}$ in (79) is negative. Indeed, it contains less positive terms than (74) which we showed to be negative under the assumptions from Lemma 3. We then have:

We have two linearly contracting terms. The sum contracts linearly with the minimum geometric rate factor between $\gamma\mu/2$ and $\rho$. If we define $m:=\min(\rho,\gamma\mu/2)$, $M:=\max(\rho,\gamma\mu/2)$ and $\rho^{*}:=\nu m$ with $0<\nu<1$,$\nu$ is introduced to circumvent the problematic case where $\rho$ and $\gamma\mu/2$ are too close together, which does not prevent the geometric convergence, but makes the constant $\frac{1}{1-\eta}$ potentially very big (in the case both terms are equal, the sum even becomes an annoying linear term in t). we then get:Note that if $m\neq\rho$, we can perform the index change $t+1-k\rightarrow k$ to get the sum.

where $\eta:=\frac{1-M}{1-\rho^{*}}$. We have $\frac{1}{1-\eta}=\frac{1-\rho^{*}}{M-\rho^{*}}$.

By taking $\nu=\frac{4}{5}$ and setting $\rho=\frac{1}{4n}$ – its maximal value allowed by the assumptions of Lemma 3 – we get $M\geq\frac{1}{4n}$ and $\rho^{*}\leq\frac{1}{5n}$, which means $\frac{1}{1-\eta}\leq 20n$.

C.10 Proof of Corollary 3 (speedup regimes)

Referring to Hofmann et al. (2015) and our own Theorem 1, the geometric rate factor of Saga is $\frac{1}{5}\min\{\frac{1}{n},\frac{a}{\kappa}\}$ for a stepsize of $\gamma=\frac{a}{5L}$. We start by proving the first part of the corollary which considers the step size $\gamma=\frac{a}{L}$ with $a=a^{*}(\tau)$. We distinguish between two regimes to study the parallel speedup our algorithm obtains and to derive a condition on $\tau$ for which we have a linear speedup.

In this regime, $n>\kappa$ and the geometric rate factor of sequential Saga is $\frac{1}{5n}$. To get a linear speedup (up to a constant factor), we need to enforce $\rho^{*}=\Omega(\frac{1}{n})$. We recall that $\rho^{*}=\min\{\frac{1}{5n},a\frac{1}{5\kappa}\}$.

We already have $\frac{1}{5n}=\Omega(\frac{1}{n})$. This means that we need $\tau$ to verify $\frac{a^{*}(\tau)}{5\kappa}=\Omega(\frac{1}{n})$, where $a^{*}(\tau)=\frac{1}{32\left(1+\tau\sqrt{\Delta}\right)\xi(\kappa,\Delta,\tau)}$ according to Theorem 2. Recall that $\xi(\kappa,\Delta,\tau):=\sqrt{1+\frac{1}{8\kappa}\min\{\frac{1}{\sqrt{\Delta}},\tau\}}$. Up to a constant factor, this means we can give the following sufficient condition:

We now consider two alternatives, depending on whether $\kappa$ is bigger than $\frac{1}{\sqrt{\Delta}}$ or not. If $\kappa\geq\frac{1}{\sqrt{\Delta}}$, then $\xi(\kappa,\Delta,\tau)<2$ and we can rewrite the sufficient condition (106) as:

In the alternative case, $\kappa\leq\frac{1}{\sqrt{\Delta}}$. Since $a^{*}(\tau)$ is decreasing in $\tau$, we can suppose $\tau\geq\frac{1}{\sqrt{\Delta}}$ without loss of generality and thus $\xi(\kappa,\Delta,\tau)=\sqrt{1+\frac{1}{8\kappa\sqrt{\Delta}}}$. We can then rewrite the sufficient condition (106) as:

We observe that since we have supposed that $\kappa\leq\frac{1}{\sqrt{\Delta}}$, we have $\sqrt{\kappa\sqrt{\Delta}}\leq\kappa\sqrt{\Delta}\leq 1$, which means that our initial assumption that $\tau<\frac{n}{10}$ is stronger than condition (C.10).

We can now combine both cases to get the following sufficient condition for the geometric rate factor of Asaga to be the same order as sequential Saga when $n>\kappa$:

In this regime, $\kappa>n$ and the geometric rate factor of sequential Saga is $a\frac{1}{\kappa}$. Here, to obtain a linear speedup, we need $\rho^{*}=\mathcal{O}(\frac{1}{\kappa})$. Since $\frac{1}{n}>\frac{1}{\kappa}$, all we require is that $\frac{a^{*}(\tau)}{\kappa}=\Omega(\frac{1}{\kappa})$ where $a^{*}(\tau)=\frac{1}{32\left(1+\tau\sqrt{\Delta}\right)\xi(\kappa,\Delta,\tau)}$, which reduces to $a^{*}(\tau)=\Omega(1)$.

We can give the following sufficient condition:

Using that $\frac{1}{n}\leq\Delta\leq 1$ and that $\kappa>n$, we get that $\xi(\kappa,\Delta,\tau)\leq 2$, which means our sufficient condition becomes:

This finishes the proof for the first part of Corollary 3.

If $\tau=\mathcal{O}(\frac{1}{\sqrt{\Delta}})$, then $\xi(\kappa,\Delta,\tau)=\mathcal{O}(1)$ and $(1+\tau\sqrt{\Delta})=\mathcal{O}(1)$, and thus $a^{*}(\tau)=\Omega(1)$ (for any $n$ and $\kappa$). This means that the universal stepsize $\gamma=\Theta(1/L)$ satisfies $\gamma\leq a^{*}(\tau)$ for any $\kappa$, giving the same rate factor $\Omega(\min\{\frac{1}{n},\frac{1}{\kappa}\})$ that sequential Saga has, completing the proof for the second part of Corollary 3. ∎

Appendix D Additional experimental results

Sparsity plays an important role in our theoretical results, where we find that while it is necessary in the “ill-conditioned” regime to get linear speedups, it is not in the “well-conditioned” regime. We confront this to real-life experiments by comparing the convergence and speedup performance of our three asynchronous algorithms on the Covtype dataset, which is fully dense after standardization. The results appear in Figure 2.

While we still see a significant improvement in speed when increasing the number of cores, this improvement is smaller than the one we observe for sparser datasets. The speedups we observe are consequently smaller, and taper off earlier than on our other datasets. However, since the observed “theoretical” speedup is linear (see Section D.2), we can attribute this worse performance to higher hardware overhead. This is expected because each update is fully dense and thus the shared parameters are much more heavily contended for than in our sparse datasets.

One thing we notice when computing the $\Delta$ variable for our datasets is that it often fails to capture the full sparsity distribution, being essentially a maximum. This means that $\Delta$ can be quite big even for very sparse datasets. Deriving a less coarse bound remains an open problem.

D.2 Theoretical speedups

In the main text of this paper, we show experimental speedup results where suboptimality is a function of the running time. This measure encompasses both theoretical algorithmic properties and hardware overheads (such as contention of shared memory) which are not taken into account in our analysis.

In order to isolate these two effects, we plot our convergence experiments where suboptimality is a function of the number of iterations; thus, we abstract away any potential hardware overhead.To do so, we implement a global counter which is sparsely updated (every $100$ iterations for example) in order not to modify the asynchrony of the system. This counter is used only for plotting purposes and is not needed otherwise. The experimental results can be seen in Figure 3.

For all three algorithms and all three datasets, the curves for $1$ and $10$ cores almost coincide, which means that we are indeed in the “theoretical linear speedup” regime. Indeed, when we plotted the amount of iterations required to converge to a given accuracy as a function of the number of cores, we obtained straight horizontal lines for our three algorithms.

The fact that the speedups we observe in running time are less than linear can thus be attributed to various hardware overheads, including shared variable contention – the compare-and-swap operations are more and more expensive as the number of competing requests increases – and cache effects as mentioned in Section 4.3.

Appendix E A closer look at the τ𝜏\tau constant

In the parallel optimization literature, $\tau$ is often referred to as a proxy for the number of cores. However, intuitively as well as in practice, it appears that there are a number of other factors that can influence this quantity. We will now attempt to give a few qualitative arguments as to what these other factors might be and how they relate to $\tau$.

The first of these factors is indeed the number of cores.

If we have $p$ cores, $\tau\geq p-1$. Indeed, in the best-case scenario where all cores have exactly the same execution speed for a single iteration, $\tau=p-1$.

To get more insight into what $\tau$ really encompasses, let us now try to define the worst-case scenario in the preceding example. Consider $2$ cores. In the worst case scenario, one core runs while the other is stuck. Then the overlap is $t$ for all $t$ and eventually grows to $+\infty$. If we assume that one core runs twice as fast as the other, then $\tau=2$. If both run at the same speed, $\tau=1$.

It appears then that a relevant quantity is $R$, the ratio between the fastest execution time to the slowest execution time for a single iteration. $\tau\leq(p-1)R$, which can be arbitrarily bigger than $p$.

There are several factors at play in $R$ itself.

The first is the speed of execution of the cores themselves (i.e. clock time). The dependency here is quite clear.

The second is the data matrix itself. If one $f_{i}$ has support of size $n$ while all the others have support of size $1$, $r$ may eventually become very big.

The third is the length of the computation itself. The longer our algorithm runs, the more likely it is to explore the potential corner cases of the data matrix.

The overlap is upper bounded by the number of cores times the maximum iteration time over the minimum iteration time (which is linked to the sparsity distribution of the data matrix). This is an upper bound, which means that in some cases it will not really be useful. For example, in the case where one factor has support size $1$ and all others have support size $d$, the probability of the event which corresponds to the upper bound is exponentially small in $d$. We conjecture that a more useful indicator could be the maximum iteration time over the expected iteration time.

To sum up this preliminary theoretical exploration, the $\tau$ term encompasses a lot more complexity than is usually implied in the literature. This is reflected in the experiments we ran, where the constant was orders of magnitude bigger than the number of cores.

E.2 Experimental results

In order to verify our intuition about the $\tau$ variable, we ran several experiments on all three datasets, whose characteristics are reminded in Table 1. $\delta_{l}^{i}$ is the support size of $f_{i}$.

To estimate $\tau$, we compute the average overlap over $100$ iterations, which is a lower bound on the actual overlap (which is a maximum, not an average). We then take the maximum observed quantity. We use an average because computing the overlap requires using a global counter, which we do not want to update every iteration since it would make it a heavily contentious quantity susceptible of artificially changing the asynchrony of our algorithm.

The results we observe are order of magnitude bigger than $p$, indicating that $\tau$ can indeed not be dismissed as a mere proxy for the number of cores, but has to be more carefully analyzed.

First, we plot the maximum observed $\tau$ as a function of the number of cores (see Figure 4). We observe that the relationship does indeed seem to be roughly linear with respect to the number of cores until 30 cores. After 30 cores, we observe what may be a phase transition where the slope increases significantly.

Second, we measured the maximum observed $\tau$ as a function of the number of epochs. We omit the figure since we did not observe any dependency; that is, $\tau$ does not seem to depend on the number of epochs. We know that it must depend on the number of iterations (since it cannot be bigger, and is an increasing function with respect to that number for example), but it appears that a stable value is reached quite quickly (before one full epoch is done).

If we allowed the computations to run forever, we would eventually observe an event such that $\tau$ would reach the upper bound mentioned in the last section, so it may be that $\tau$ is actually a very slowly increasing function of the number of iterations.

Appendix F Lagged updates and Sparsity

The lagged updates technique in Saga is based on the observation that the updates for component $[x]_{v}$ need not be applied until this coefficient needs to be accessed, that is, until the next iteration $t$ such that $v\in S_{i_{t}}$. We refer the reader to Schmidt et al. (2016) for more details.

Interestingly, the expected number of iterations between two steps where a given dimension $v$ is involved in the partial gradient is $p_{v}^{-1}$, where $p_{v}$ is the probability that $v$ is involved in a given step. $p_{v}^{-1}$ is precisely the term which we use to multiply the update to $[x]_{v}$ in Sparse Saga. Therefore one may see the updates in Sparse Saga as anticipated updates, whereas those in the Schmidt et al. (2016) implementation are lagged. The two algorithms appear to be very close, even though Sparse Saga uses an expectation to multiply a given update whereas the lazy implementation uses a random variable (with the same expectation). Sparse Saga therefore uses a slightly more aggressive strategy, which gave faster run-time in our experiments below.

Although Sparse Saga requires the computation of the $p_{v}$ probabilities, this can be done during a first pass throughout the data (during which constant step size Sgd may be used) at a negligible cost.

In our experiments, we compare the Sparse Saga variant proposed in Section 2 to two other approaches: the naive (i.e. dense) update scheme and the lagged updates implementation described in Defazio et al. (2014). Note that we use different datasets from the parallel experiments, including a subset of the RCV1 dataset and the realsim dataset. Figure 5 reveals that sparse and lagged updates have a lower cost per iteration, resulting in faster convergence for sparse datasets. Furthermore, while the two approaches had similar convergence in terms of number of iterations, the Sparse Saga scheme is slightly faster in terms of runtime (and as previously pointed out, sparse updates are better adapted for the asynchronous setting). For the dense dataset (Covtype), the three approaches exhibit a similar performance.

F.2 On the difficulty of parallel lagged updates

In the implementation presented in Schmidt et al. (2016), the dense part ($\bar{\alpha}$) of the updates is deferred. Instead of writing dense updates, counters $c_{d}$ are kept for each coordinate of the parameter vector – which represent the last time these variables were updated – as well as the average gradient $\bar{\alpha}$ for each coordinate. Then, whenever a component $[\hat{x}]_{d}$ is needed (in order to compute a new gradient), we subtract $\gamma(t-c_{d})[\bar{\alpha}]_{d}$ from it and $c_{d}$ is set to $t$. The reason we can do this without modifying the algorithm is that $[\bar{\alpha}]_{d}$ only changes when $[\hat{x}]_{d}$ also does.

In the sequential setting, this is strictly the same as doing the updates in a dense way, since the coordinates are only stale when they’re not used. Note that at the end of an execution all counters have to be subtracted at once to get the true final parameter vector (and to bring every $c_{d}$ counter to the final $t$).

In the parallel setting, several issues arise:

two cores might be attempting to correct the lag at the same time. In which case since updates are done as additions and not replacements (which is necessary to ensure that there are no overwrites), the lag might be corrected multiple times, i.e. overly corrected.

we would have to read and write atomically to each $[\hat{x}_{d}],c_{d},[\bar{\alpha}]_{d}$ triplet, which is highly impractical.

we would need to have an explicit global counter, which we do not in Asaga (our global counter $t$ being used solely for the proof).

in the dense setting, updates happen coordinate by coordinate. So at time $t$ the number of $\bar{\alpha}$ updates a coordinate has received from a fixed past time $c_{d}$ is a random variable, which may differs from coordinate to coordinate. Whereas in the lagged implementation, the multiplier is always $(t-c_{d})$ which is a constant (conditional to $c_{d}$), which means a potentially different $\hat{x}_{t}$.

All these points mean both that the implementation of such a scheme in the parallel setting would be impractical, and that it would actually yields a different algorithm than the dense version, which would be even harder to analyze.

Appendix G Additional empirical details

We run our experiments on four datasets. In every case, we run logistic regression for the purpose of binary classification.

The first is the Reuters Corpus Volume I (RCV1) dataset (Lewis et al., 2004), an archive of over 800,000 manually categorized newswire stories made available by Reuters, Ltd. for research purposes. The associated task is a binary text categorization.

Our second dataset was first introduced in Ma et al. (2009). Its associated task is a binary malicious url detection. This dataset contains more than 2 million URLs obtained at random from Yahoo’s directory listing (for the “benign” URLs) and from a large Web mail provider (for the “malicious” URLs). The benign to malicious ratio is 2. Features include lexical information as well as metadata. This dataset was obtained from the libsvmtools project.http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary.html

On our third dataset, the associated task is a binary classification problem (down from 7 classes originally, following the pre-treatment of Collobert et al. (2002)). The features are cartographic variables. Contrarily to the first two, this is a dense dataset.

We only use our fourth dataset for non-parallel experiments and a specific compare-and-swap test. It constitutes of UseNet articles taken from four discussion groups (simulated auto racing, simulated aviation, real autos, real aviation).

G.2 Implementation details

All experiments were run on a Dell PowerEdge 920 machine with 4 Intel Xeon E7-4830v2 processors with 10 2.2GHz cores each and 384GB 1600 Mhz RAM.

All algorithms were implemented in the Scala language and the software stack consisted of a Linux operating system running Scala 2.11.7 and Java 1.6.

We chose this expressive, high level language for our experimentation despite its typical 20x slower performance compared to C because our primary concern was that the code may easily be reused and extended for research purposes (which is harder to achieve with low level, heavily optimized C code; especially for error prone parallel computing).

As a result our timed experiments exhibit sub-optimal running times, e.g. compared to Konecny and Richtarik (2013). This is as we expected. The observed slowdown is both consistent across datasets (roughly 20x) and with other papers that use Scala code (e.g. Mania et al. (2015), Ma et al. (2015, Fig. 2)).

Despite this slowdown, our experiments show state-of-the-art results in convergence per number of iterations. Furthermore, the speed-up patterns that we observe for our implementation of Hogwild and Kromagnon are similar to the ones given in [MN15], Niu et al. and Reddi et al. (in various languages).

The code we used to run all the experiments is available at https://github.com/RemiLeblond/ASAGA.

Interestingly, we have found necessary to use compare-and-swap instructions in the implementation of Asaga. In Figure 6, we display suboptimality plots using non-thread safe operations and compare-and-swap (CAS) operations. The non-thread safe version starts faster but then fails to converge beyond a specific level of suboptimality, while the compare-and-swap version does converges linearly up to machine precision.

For compare-and-swap instructions we used the AtomicDoubleArray class from the Google library Guava. This class uses an AtomicLongArray under the hood (from package java.util.concurrent.atomic in the standard Java library), which does indeed benefit from lower-level CPU-optimized instructions.

Storing $n$ gradient may seem like an expensive proposition, but for linear predictor models, one can actually store a single scalar per gradient (as proposed in Schmidt et al. (2016)), which is what we do in our implementation of Asaga.

For each algorithm, we picked the best step size among 10 equally spaced values in a grid, and made sure that the best step size was never at the boundary of this interval. For Covtype and RCV1, we used the interval $[\frac{1}{10L},\frac{10}{L}]$, whereas for URL we used the interval $[\frac{1}{L},\frac{100}{L}]$ as it admitted larger step sizes. It turns out that the best step size was fairly constant for different number of cores for both Asaga and Kromagnon, and both algorithms had similar best step sizes.

G.3 Biased update in the implementation

In the implementation detailed in Algorithm 2, $\bar{\alpha}$ is maintained in memory instead of being recomputed for every iteration. This saves both the cost of reading every data point for each iteration and of computing $\bar{\alpha}$ for each iteration.

However, this removes the unbiasedness guarantee. The problem here is the definition of the expectation of $\hat{\alpha}_{i}$. Since we are sampling uniformly at random, the average of the $\hat{\alpha}_{i}$ is taken at the precise moment when we read the $\alpha_{i}^{t}$ components. Without synchronization, between two reads to a single coordinate in $\alpha_{i}$ and in $\bar{\alpha}$, new updates might arrive in $\bar{\alpha}$ that are not yet taken into account in $\alpha_{i}$. Conversely, writes to a component of $\alpha_{i}$ might precede the corresponding write in $\bar{\alpha}$ and induce another source of bias.

In order to alleviate this issue, we can use coordinate-level locks on $\alpha_{i}$ and $\bar{\alpha}$ to make sure they are always synchronized. Such low-level locks are quite inexpensive when $d$ is large, especially when compared to vector-wide locks.

However, as previously noted, experimental results indicate that this fix is not necessary.