Title: A2BCD: Asynchronous Acceleration with Optimal Complexity

Abstract: In this paper, we propose the Asynchronous Accelerated Nonuniform Randomized Block Coordinate Descent algorithm (A2BCD). We prove A2BCD converges linearly to a solution of the convex minimization problem at the same rate as NU_ACDM, so long as the maximum delay is not too large. This is the first asynchronous Nesterov-accelerated algorithm that attains any provable speedup. Moreover, we then prove that these algorithms both have optimal complexity. Asynchronous algorithms complete much faster iterations, and A2BCD has optimal complexity. Hence we observe in experiments that A2BCD is the top-performing coordinate descent algorithm, converging up to 4 -5× faster than NU_ACDM on some data sets in terms of wall-clock time. To motivate our theory and proof techniques, we also derive and analyze a continuous-time analogue of our algorithm and prove it converges at the same rate.

1. Introduction:

In this paper, we propose and prove the convergence of the Asynchronous Accelerated Nonuniform Randomized Block Coordinate Descent algorithm (A2BCD), the first asynchronous Nesterovaccelerated algorithm that achieves optimal complexity. No previous attempts have been able to prove a speedup for asynchronous Nesterov acceleration. We aim to find the minimizer x * of the unconstrained minimization problem:

min x∈R d f (x) = f x (1) , . . . , x (n) (1.1)

where f is σ-strongly convex for σ > 0 with L-Lipschitz gradient ∇f = (∇ 1 f, . . . , ∇ n f ). x ∈ R d is composed of coordinate blocks x (1) , . . . , x (n) . The coordinate blocks of the gradient ∇ i f are assumed L i -Lipschitz with respect to the ith block. That is, ∀x, h ∈ R d :

∇ i f (x + P i h) -∇ i f (x) ≤ L i h (1.2)

where P i is the projection onto the ith block of R d . Let L 1 n n i=1 L i be the average block Lipschitz constant. These conditions on f are assumed throughout this whole paper. Our algorithm can also be applied to non-strongly convex objectives (σ = 0) or non-smooth objectives using the black box reduction techniques proposed in Allen-Zhu & Hazan (2016) . Hence we consider only the coordinate smooth, strongly-convex case. Our algorithm can also be applied to the convex regularized ERM problem via the standard dual transformation (see for instance Lin et al. (2014) ):

f (x) = 1 n n i=1 f i ( a i , x ) + λ 2 x 2 (1.3)

Hence A2BCD can be used as an asynchronous Nesterov-accelerated finite-sum algorithm.

Coordinate descent methods, in which a chosen coordinate block i k is updated at every iteration, are a popular way to solve equation 1.1. Randomized block coordinate descent (RBCD, Nesterov (2012) ) updates a uniformly randomly chosen coordinate block i k with a gradient-descent-like step:

x k+1 = x k -(1/L i k )∇ i k f (x k ).

The complexity K( ) of an algorithm is defined as the number of iterations required to decrease the error E(f (x k ) -f (x * )) to less than (f (x 0 ) -f (x * )). Randomized coordinate descent has a complexity of K( ) = O(n( L/σ) ln(1/ )).

Using a series of averaging and extrapolation steps, accelerated RBCD Nesterov (2012) improves RBCD's iteration complexity K( ) to O(n L/σ ln(1/ )), which leads to much faster convergence when L σ is large. This rate is optimal when all L i are equal Lan & Zhou (2015) . Finally, using a special probability distribution for the random block index i k , the non-uniform accelerated coordinate descent method Allen-Zhu et al. (2015) (NU_ACDM) can further decrease the complexity to O( n i=1 L i /σ ln(1/ )), which can be up to √ n times faster than accelerated RBCD, since some L i can be significantly smaller than L. NU_ACDM is the current state-of-the-art coordinate descent algorithm for solving equation 1.1.

Our A2BCD algorithm generalizes NU_ACDM to the asynchronous-parallel case. We solve equation 1.1 with a collection of p computing nodes that continually read a shared-access solution vector y into local memory then compute a block gradient ∇ i f , which is used to update shared solution vectors (x, y, v) . Proving convergence in the asynchronous case requires extensive new technical machinery.

A traditional synchronous-parallel implementation is organized into rounds of computation: Every computing node must complete an update in order for the next iteration to begin. However, this synchronization process can be extremely costly, since the lateness of a single node can halt the entire system. This becomes increasingly problematic with scale, as differences in node computing speeds, load balancing, random network delays, and bandwidth constraints mean that a synchronous-parallel solver may spend more time waiting than computing a solution.

Computing nodes in an asynchronous solver do not wait for others to complete and share their updates before starting the next iteration. They simply continue to update the solution vectors with the most recent information available, without any central coordination. This eliminates costly idle time, meaning that asynchronous algorithms can be much faster than traditional ones, since they have much faster iterations. For instance, random network delays cause asynchronous algorithms to complete iterations Ω(ln(p)) time faster than synchronous algorithms at scale. This and other factors that influence the speed of iterations are discussed in Hannah & Yin (2017a) . However, since many iterations may occur between the time that a node reads the solution vector, and the time that its computed update is applied, effectively the solution vector is being updated with outdated information. At iteration k, the block gradient ∇ i k f is computed at a delayed iterate ŷk defined as1 : ŷk = y k-j(k,1) (1) , . . . , y k-j(k,n) (n)

(1.4)

1.1 Summary of Contributions for delay parameters j(k, 1), . . . , j(k, n) ∈ N. Here j(k, i) denotes how many iterations out of date coordinate block i is at iteration k. Different blocks may be out of date by different amounts, which is known as an inconsistent read. We assume2 that j(k, i) ≤ τ for some constant τ < ∞.

Asynchronous algorithms were proposed in Chazan & Miranker (1969) to solve linear systems. General convergence results and theory were developed later in Bertsekas (1983) ; Bertsekas & Tsitsiklis (1997) ; Tseng et al. (1990) ; Luo & Tseng (1992; 1993) ; Tseng (1991) for partially and totally asynchronous systems, with essentially-cyclic block sequence i k . More recently, there has been renewed interest in asynchronous algorithms with random block coordinate updates. Linear and sublinear convergence results were proven for asynchronous RBCD Liu & Wright (2015) ; Liu et al. (2014) ; Avron et al. (2014) , and similar was proven for asynchronous SGD Recht et al. (2011) , and variance reduction algorithms Reddi et al. (2015) ; Leblond et al. (2017) ; Mania et al. (2015) ; Huo & Huang (2016) , and primal-dual algorithms Combettes & Eckstein (2018) .

There is also a rich body of work on asynchronous SGD. In the distributed setting, Zhou et al. (2018) showed global convergence for stochastic variationally coherent problems even when the delays grow at a polynomial rate. In Lian et al. (2018) , an asynchronous decentralized SGD was proposed with the same optimal sublinear convergence rate as SGD and linear speedup with respect to the number of workers. In Liu et al. (2018) , authors obtained an asymptotic rate of convergence for asynchronous momentum SGD on streaming PCA, which provides insight into the tradeoff between asynchrony and momentum. In Dutta et al. (2018) , authors prove convergence results for asynchronous SGD that highlight the tradeoff between faster iterations and iteration complexity. Further related work is discussed in Section 4.

1.1. Summary of Contributions:

In this paper, we prove that A2BCD attains NU_ACDM's state-of-the-art iteration complexity to highest order for solving equation 1.1, so long as delays are not too large (see Section 2). The proof is very different from that of Allen-Zhu et al. (2015) , and involves significant technical innovations and complexity related to the analysis of asynchronicity.

We also prove that A2BCD (and hence NU_ACDM) has optimal complexity to within a constant factor over a fairly general class of randomized block coordinate descent algorithms (see Section 2.1). This extends results in Lan & Zhou (2015) to asynchronous algorithms with L i not all equal. Since asynchronous algorithms complete faster iterations, and A2BCD has optimal complexity, we expect A2BCD to be faster than all existing coordinate descent algorithms. We confirm with numerical experiments that A2BCD is the current fastest coordinate descent algorithm (see Section 5).

We are only aware of one previous and one contemporaneous attempt at proving convergence results for asynchronous Nesterov-accelerated algorithms. However, the first is not accelerated and relies on extreme assumptions, and the second obtains no speedup. Therefore, we claim that our results are the first-ever analysis of asynchronous Nesterov-accelerated algorithms that attains a speedup. Moreover, our speedup is optimal for delays not too large3 .

The work of Meng et al. claims to obtain square-root speedup for an asynchronous accelerated SVRG.

In the case where all component functions have the same Lipschitz constant L, the complexity they obtain reduces to (n + κ) ln(1/ ) for κ = O τ n 2 (Corollary 4.4). Hence authors do not even obtain accelerated rates. Their convergence condition is τ < In a contemporaneous preprint, authors in Fang et al. (2018) skillfully devised accelerated schemes for asynchronous coordinate descent and SVRG using momentum compensation techniques. Although their complexity results have the improved √ κ dependence on the condition number, they do not prove any speedup. Their complexity is τ times larger than the serial complexity. Since τ is necessarily greater than p, their results imply that adding more computing nodes will increase running time. The authors claim that they can extend their results to linear speedup for asynchronous, accelerated SVRG under sparsity assumptions. And while we think this is quite likely, they have not yet provided proof.

We also derive a second-order ordinary differential equation (ODE), which is the continuous-time limit of A2BCD (see Section 3). This extends the ODE found in Su et al. (2014) to an asynchronous accelerated algorithm minimizing a strongly convex function. We prove this ODE linearly converges to a solution with the same rate as A2BCD's, without needing to resort to the restarting techniques. The ODE analysis motivates and clarifies the our proof strategy of the main result.

2. Main results:

We should consider functions f where it is efficient to calculate blocks of the gradient, so that coordinate-wise parallelization is efficient. That is, the function should be "coordinate friendly" Peng et al. (2016b) . This is a very wide class that includes regularized linear regression, logistic regression, etc. The L 2 -regularized empirical risk minimization problem is not coordinate friendly in general, however the equivalent dual problem is, and hence can be solved efficiently by A2BCD (see Lin et al. (2014) , and Section 5).

To calculate the k + 1'th iteration of the algorithm from iteration k, we use only one block of the gradient ∇ i k f . We assume that the delays j(k, i) are independent of the block sequence i k , but otherwise arbitrary (This is a standard assumption found in the vast majority of papers, but can be relaxed Sun et al. (2017) ; Leblond et al. (2017) ; Cannelli et al. (2017) ). Definition 1. Asynchronous Accelerated Randomized Block Coordinate Descent (A2BCD). Let f be σ-strongly convex, and let its gradient ∇f be L-Lipschitz with block coordinate Lipschitz parameters L i as in equation 1.2. We define the condition number κ = L/σ, and let L = min i L i . Using these parameters, we sample i k in an independent and identically distributed (IID) fashion according to

P[i k = j] = L 1/2 j /S, j ∈ {1, . . . , n}, for S = n i=1 L 1/2 i .

(2.1)

Let τ be the maximum asynchronous delay. We define the dimensionless asynchronicity parameter ψ, which is proportional to τ , and quantifies how strongly asynchronicity will affect convergence:

ψ = 9 S -1/2 L -1/2 L 3/4 κ 1/4 × τ (2.2)

We use the above system parameters and ψ to define the coefficients α, β, and γ via eqs. (2.3) to (2.5). Hence A2BCD algorithm is defined via the iterations: eqs. (2.6) to (2.8).

α 1 + (1 + ψ)σ -1/2 S -1

(2.3)

β 1 -(1 -ψ)σ 1/2 S -1 (2.4) h 1 - 1 2 σ 1/2 L -1/2 ψ.

(2.5)

y k = αv k + (1 -α)x k , (2.6) x k+1 = y k -hL -1 i k ∇i k f (ŷ k ), (2.7) v k+1 = βv k + (1 -β)y k -σ -1/2 L -1/2 i k ∇i k f (ŷ k ).

(2.8)

See Section A for a discussion of why it is practical and natural to have the gradient ∇ i k f (ŷ k ) to be outdated, while the actual variables x k , y k , v k can be efficiently kept up to date. Essentially it is because most of the computation lies in computing ∇ i k f (ŷ k ). After this is computed, x k , y k , v k can be updated more-or-less atomically with minimal overhead, meaning that they will always be up to date. However our main results still hold for more general asynchronicity.

A natural quantity to consider in asynchronous convergence analysis is the asynchronicity error, a powerful tool for analyzing asynchronous algorithms used in several recent works Peng et al. (2016a) ; Hannah & Yin (2017b) ; Sun et al. (2017) ; Hannah & Yin (2017a) . We adapt it and use a weighted sum of the history of the algorithm with decreasing weight as you go further back in time. Definition 2. Asynchronicity error. Using the above parameters, we define:

A k = τ j=1 c j y k+1-j -y k-j 2 (2.9) for c i = 6 S L 1/2 κ 3/2 τ τ j=i 1 -σ 1/2 S -1 i-j-1 ψ -1 . (2.10)

Here we define y k = y 0 for all k < 0. The determination of the coefficients c i is in general a very involved process of trial and error, intuition, and balancing competing requirements. The algorithm doesn't depend on the coefficients, however; they are only an analytical tool.

We define E k [X] as the expectation of X conditioned on (x 0 , . . . , x k ), (y 0 , . . . , y k ), (v 0 , . . . , v k ), and (i 0 , . . . , i k-1 ). To simplify notation4 , we assume that the minimizer x * = 0, and that f (x * ) = 0 with no loss in generality. We define the Lyapunov function:

ρ k = v k 2 + A k + cf (x k ) (2.11) for c = 2σ -1/2 S -1 βα -1 (1 -α) + 1 .

(2.12)

We now present this paper's first main contribution. Theorem 1. Let f be σ-strongly convex with a gradient ∇f that is L-Lipschitz with block Lipschitz constants

{L i } n i=1 . Let ψ defined in equation 2.2 satisfy ψ ≤ 3 7 (i.e. τ ≤ 1 21 S 1/2 L 1/2 L -3/4 κ -1/4

). Then for A2BCD we have:

E k [ρ k+1 ] ≤ 1 -(1 -ψ)σ 1/2 S -1 ρ k .

To obtain E[ρ k ] ≤ ρ 0 , it takes K A2BCD ( ) iterations for:

KA2BCD( ) = σ -1/2 S + O(1) ln(1/ ) 1 -ψ , (2.13)

where O(•) is asymptotic with respect to σ -1/2 S → ∞, and uniformly bounded.

This result is proven in Section B. A stronger result for L i ≡ L can be proven, but this adds to the complexity of the proof; see Section E for a discussion. In practice, asynchronous algorithms are far more resilient to delays than the theory predicts. τ can be much larger without negatively affecting the convergence rate and complexity. This is perhaps because we are limited to a worst-case analysis, which is not representative of the average-case performance.

Allen-Zhu et al. (2015) (Theorem 5.1) shows a linear convergence rate of 1 -2/ 1 + 2σ -1/2 S for NU_ACDM, which leads to the corresponding iteration complexity of

K NU_ACDM ( ) = σ -1/2 S + O(1) ln(1/ )

. Hence, we have:

K A2BCD ( ) = 1 1 -ψ (1 + o(1))K NU_ACDM ( )

Published as a conference paper at ICLR 2019

2.1. Optimality:

When 0 ≤ ψ 1, or equivalently, when τ S 1/2 L 1/2 L -3/4 κ -1/4 , the complexity of A2BCD asymptotically matches that of NU_ACDM. Hence A2BCD combines state-of-the-art complexity with the faster iterations and superior scaling that asynchronous iterations allow. We now present some special cases of the conditions on the maximum delay τ required for good complexity.

Corollary 3. Let the conditions of Theorem 1 hold. If all coordinate-wise Lipschitz constants L i are equal (i.e.

L i = L 1 , ∀i), then we have K A2BCD ( ) ∼ K NU_ACDM ( ) when τ n 1/2 κ -1/4 (L 1 /L) 3/4 . If we further assume all coordinate-wise Lipschitz constants L i equal L. Then K A2BCD ( ) ∼ K NU_ACDM ( ) = K ACDM ( ), when τ n 1/2 κ -1/4 .

Remark 1. Reduction to synchronous case. Notice that when τ = 0, we have ψ = 0, c i ≡ 0 and hence A k ≡ 0. Thus A2BCD becomes equivalent to NU_ACDM, the Lyapunov function5 ρ k becomes equivalent to one found in Allen-Zhu et al. (2015) (pg. 9), and Theorem 1 yields the same complexity.

The maximum delay τ will be a function τ (p) of p, number of computing nodes. Clearly τ ≥ p, and experimentally it has been observed that τ = O(p) Leblond et al. (2017) . Let gradient complexity K( , τ ) be the number of gradients required for an asynchronous algorithm with maximum delay τ to attain suboptimality . τ (1) = 0, since with only 1 computing node there can be no delay. This corresponds to the serial complexity. We say that an asynchronous algorithm attains a complexity speedup if pK( ,τ (0)) K( ,τ (p) is increasing in p. We say it attains linear complexity speedup if pK( ,τ (0)) K( ,τ (p) = Ω(p). In Theorem 1, we obtain a linear complexity speedup (for p not too large), whereas no other prior attempt can attain even a complexity speedup with Nesterov acceleration.

In the ideal scenario where the rate at which gradients are calculated increases linearly with p, algorithms that have linear complexity speedup will have a linear decrease in wall-clock time. However in practice, when the number of computing nodes is sufficiently large, the rate at which gradients are calculated will no longer be linear. This is due to many parallel overhead factors including too many nodes sharing the same memory read/write bandwidth, and network bandwidth. However we note that even with these issues, we obtain much faster convergence than the synchronous counterpart experimentally.

NU_ACDM and hence A2BCD are in fact optimal in some sense. That is, among a fairly wide class of coordinate descent algorithms A, they have the best-possible worst-case complexity to highest order. We extend the work in Lan & Zhou (2015) to encompass algorithms are asynchronous and have unequal L i . For a subset S ∈ R d , we let IC(S) (inconsistent read) denote the set of vectors v whose components are a combination of components of vectors in the set

S. That is, v = (v 1,1 , v 2,2 , . . . , v d,d ) for some vectors v 1 , v 2 , . . . , v d ∈ S.

Here v i,j denotes the jth component of vector v i . Definition 4. Asynchronous Randomized Incremental Algorithms. Consider the unconstrained minimization problem equation 1.1 for function f satisfying the conditions stated in Section 1. We define the class A as algorithms G on this problem such that:

1. For each parameter set (σ, L 1 , . . . , L n , n), G has an associated IID random variable i k with some fixed distribution P[i k ] = p i for n i=1 p i = 1. 2. The iterates of A satisfy: x k+1 ∈ span{IC(X k ), ∇i 0 f (IC(X0)), ∇i 1 f (IC(X1)), . . . , ∇i k f (IC(X k ))}

This is a rather general class: x k+1 can be constructed from any inconsistent reading of past iterates IC(X k ), and any past gradient of an inconsistent read ∇ ij f (IC(X j )).

Theorem 2. For any algorithm G ∈ A that solves eq. (1.1), and parameter set (σ, L 1 , . . . , L n , n), there is a dimension d, a corresponding function f on R d , and a starting point x 0 , such that

E x k -x * 2 / x 0 -x * 2 ≥ 1 2 1 -4/ n j=1 L i /σ + 2n k Hence A has a complexity lower bound: K( ) ≥ 1 4 (1 + o(1)) n j=1 L i /σ + 2n ln(1/2 )

Our proof in Section D follows very similar lines to Lan & Zhou (2015) ; Nesterov (2013) .

3. ODE Analysis:

In this section we present and analyze an ODE which is the continuous-time limit of A2BCD. This ODE is a strongly convex, and asynchronous version of the ODE found in Su et al. (2014) . For simplicity, assume L i = L, ∀i. We rescale (I.e. we replace f (x) with 1 σ f .) f so that σ = 1, and hence κ = L/σ = L. Taking the discrete limit of synchronous A2BCD (i.e. accelerated RBCD), we can derive the following ODE6 (see Section equation C.1):

Ÿ + 2n -1 κ -1/2 Ẏ + 2n -2 κ -1 ∇f (Y ) = 0 (3.1)

We define the parameter η nκ 1/2 , and the energy:

E(t) = e n -1 κ -1/2 t (f (Y ) + 1 4 Y + η Ẏ 2 )

. This is very similar to the Lyapunov function discussed in equation 2.11, with 1 4 Y (t) + η Ẏ (t) 2 fulfilling the role of v k 2 , and A k = 0 (since there is no delay yet). Much like the traditional analysis in the proof of Theorem 1, we can derive a linear convergence result with a similar rate. See Section C.2. Lemma 5. If Y satisfies equation 3.1, the energy satisfies E (t) ≤ 0, E(t) ≤ E(0), and hence:

f (Y (t)) + 1 4 Y (t) + nκ 1/2 Ẏ (t) 2 ≤ f (Y (0)) + 1 4 Y (0) + η Ẏ (0) 2 e -n -1 κ -1/2 t

We may also analyze an asynchronous version of equation 3.1 to motivate the proof of our main theorem. Here Ŷ (t) is a delayed version of Y (t) with the delay bounded by τ .

Ÿ + 2n -1 κ -1/2 Ẏ + 2n -2 κ -1 ∇f Ŷ = 0, (3.2)

Unfortunately, this energy satisfies (see Section equation C.4, equation C.7):

e -η -1 t E (t) ≤ - 1 8 η Ẏ 2 + 3κ 2 η -1 τ D(t), for D(t) t t-τ Ẏ (s) 2 ds.

Hence this energy E(t) may not be decreasing in general. But, we may add a continuous-time asynchronicity error (see Sun et al. (2017) ), much like in Definition 2, to create a decreasing energy. Let c 0 ≥ 0 and r > 0 be arbitrary constants that will be set later. Define:

A(t) = t t-τ c(t -s) Ẏ (s) 2 ds, for c(t) c 0 e -rt + e -rτ 1 -e -rτ e -rt -1 . Lemma 6. When rτ ≤ 1

2 , the asynchronicity error A(t) satisfies:

e -rt d dt e rt A(t) ≤ c 0 Ẏ (t) 2 - 1 2 τ -1 c 0 D(t).

See Section C.3 for the proof. Adding this error to the Lyapunov function serves a similar purpose in the continuous-time case as in the proof of Theorem 1 (see Lemma 11). It allows us to negate 1 2 τ -1 c 0 units of D(t) for the cost of creating c 0 units of Ẏ (t) 2 . This restores monotonicity.

Theorem 3. Let c 0 = 6κ 2 η -1 τ 2 , and r = η -1 . If τ ≤ 1 √ 48 nκ -1/2 then we have:

e -η -1 t d dt E(t) + e η -1 t A(t) ≤ 0. (3.3) Hence f (Y (t)) convergence linearly to f (x * ) with rate O exp -t/(nκ 1/2 )

Notice how this convergence condition is similar to Corollary 3, but a little looser. The convergence condition in Theorem 1 can actually be improved to approximately match this (see Section E).

Proof.

e -η -1 t d dt E(t) + e η -1 t A(t) ≤ c 0 - 1 8 η Ẏ 2 + 3κ 2 η -1 τ - 1 2 τ -1 c 0 D(t) = 6η -1 κ 2 τ 2 - 1 48 n 2 κ -1 Ẏ 2 ≤ 0

The preceding should hopefully elucidate the logic and general strategy of the proof of Theorem 1.

4. Related work:

We now discuss related work that was not addressed in Section 1. Nesterov acceleration is a method for improving an algorithm's iteration complexity's dependence the condition number κ.

Nesterov-accelerated methods have been proposed and discovered in many settings Nesterov (1983) f i (x) that use 1 function gradient ∇f i (x) per iteration). Such algorithms can often be augmented to solve composite minimization problems (minimization for objective of the form f (x) + g(x), especially for nonsomooth g), or include constraints.

In Peng et al. (2016a) , authors proposed and analyzed an asynchronous fixed-point algorithm called ARock, that takes proximal algorithms, forward-backward, ADMM, etc. as special cases. Work has also been done on asynchronous algorithms for finite sums in the operator setting Davis (2016) ; Johnstone & Eckstein (2018) . In Hannah & Yin (2017b) ; Sun et al. (2017) ; Peng et al. (2016c) ; Cannelli et al. (2017) showed that many of the assumptions used in prior work (such as bounded delay τ < ∞) were unrealistic and unnecessary in general. In Hannah & Yin (2017a) the authors showed that asynchronous iterations will complete far more iterations per second, and that a wide class of asynchronous algorithms, including asynchronous RBCD, have the same iteration complexity as their synchronous counterparts. Hence certain asynchronous algorithms can be expected to significantly outperform traditional ones.

In Xiao et al. (2017) authors propose a novel asynchronous catalyst-accelerated Lin et al. (2015) primal-dual algorithmic framework to solve regularized ERM problems. They structure the parallel updates so that the data that an update depends on is up to date (though the rest of the data may not be). However catalyst acceleration incurs a log(κ) penalty over Nesterov acceleration in general. In Allen-Zhu (2017) , the author argues that the inner iterations of catalyst acceleration are hard to tune, making it less practical than Nesterov acceleration.

5. Numerical experiments:

To investigate the performance of A2BCD, we solve the ridge regression problem. Consider the following primal and corresponding dual objective (see for instance Lin et al. (2014) ): min

w∈R d P (w) = 1 2n A T w -l 2 + λ 2 w 2 , min α∈R n D(α) = 1 2d 2 λ Aα 2 + 1 2d α + l 2 (5.1)

where A ∈ R d×n is a matrix of n samples and d features, and l is a label vector. We let A = [A 1 , . . . , A m ] where A i are the column blocks of A. We compare A2BCD (which is asynchronous accelerated), synchronous NU_ACDM (which is synchronous accelerated), and asynchronous RBCD (which is asynchronous non-accelerated). Nodes randomly select a coordinate block according to equation 2.1, calculate the corresponding block gradient, and use it to apply an update to the shared solution vectors. synchronous NU_ACDM is implemented in a batch fashion, with batch size p (1 block per computing node). Nodes in synchronous NU_ACDM implementation must wait until all nodes apply their computed gradients before they can start the next iteration, but the asynchronous algorithms simply compute with the most up-to-date information available.

We use the datasets w1a (47272 samples, 300 features), wxa which combines the data from from w1a to w8a (293201 samples, 300 features), and aloi (108000 samples, 128 features) from LIBSVM Chang & Lin (2011) . The algorithm is implemented in a multi-threaded fashion using C++11 and GNU Scientific Library with a shared memory architecture. We use 40 threads on two 2.5GHz 10-core Intel Xeon E5-2670v2 processors. See Section A.1 for a discussion of parameter tuning and estimation. The parameters for each algorithm are tuned to give the fastest performance, so that a fair comparison is possible.

A critical ingredient in the efficient implementation of A2BCD and NU_ACDM for this problem is the efficient update scheme discussed in Lee & Sidford (2013b; a) . In linear regression applications such as this, it is essential to be able to efficiently maintain or recover Ay. This is because calculating block gradients requires the vector A T i Ay, and without an efficient way to recover Ay, block gradient evaluations are essentially 50% as expensive as full-gradient calculations. Unfortunately, every accelerated iteration results in dense updates to y k because of the averaging step in equation 2.6. Hence Ay must be recalculated from scratch.

However Lee & Sidford (2013a) introduces a linear transformation that allows for an equivalent iteration that results in sparse updates to new iteration variables p and q. The original purpose of this transformation was to ensure that the averaging steps (e.g. equation 2.6) do not dominate the computational cost for sparse problems. However we find a more important secondary use which applies to both sparse and dense problems. Since the updates to p and q are sparse coordinate-block updates, the vectors Ap, and Aq can be efficiently maintained, and therefore block gradients can be efficiently calculated. The specifics of this efficient implementation are discussed in Section A.2.

In Table 5 , we plot the sub-optimality vs. time for decreasing values of λ, which corresponds to increasingly large condition numbers κ. When κ is small, acceleration doesn't result in a significantly better convergence rate, and hence A2BCD and async-RBCD both outperform sync-NU_ACDM since they complete faster iterations at similar complexity. Acceleration for low κ has unnecessary overhead, which means async-RBCD can be quite competitive. When κ becomes large, async-RBCD is no longer competitive, since it has a poor convergence rate. We observe that A2BCD and sync-NU_ACDM have essentially the same convergence rate, but A2BCD is up to 4 -5× faster than sync-NU_ACDM because it completes much faster iterations. We observe this advantage despite the fact that we are in an ideal environment for synchronous computation: A small, homogeneous, high-bandwidth, low-latency cluster. In large-scale heterogeneous systems with greater synchronization overhead, bandwidth constraints, and latency, we expect A2BCD's advantage to be much larger.

Table 1 : Sub-optimality f (y k ) -f (x * ) (y-axis) vs time in seconds (x-axis) for A2BCD, synchronous NU_ACDM, and asynchronous RBCD for data sets w1a, wxa and aloi for various values of λ.

6. Acknowledgement:



A Efficient Implementation:

An efficient implementation will have coordinate blocks of size greater than 1. This to ensure the efficiency of linear algebra subroutines. Especially because of this, the bulk of the computation for each iteration is computing ∇ i k f (ŷ k ), and not the averaging steps. Hence the computing nodes only need a local copy of y k in order to do the bulk of an iteration's computation. Given this gradient ∇ i k f (ŷ k ), updating y k and v k is extremely fast (x k can simply be eliminated). Hence it is natural to simply store y k and v k centrally, and update them when the delayed gradients ∇ i k f (ŷ k ). Given the above, a write mutex over (y, v) has minuscule overhead (which we confirm with experiments), and makes the labeling of iterates unambiguous. This also ensures that v k and y k are always up to date when (y, v) are being updated. Whereas the gradient ∇ i k f (ŷ k ) may at the same time be out of date, since it has been calculated with an outdated version of y k . However a write mutex is not necessary in practice, and does not appear to affect convergence rates or computation time. Also it is possible to prove convergence under more general asynchronicity.

A.1 Parameter selection and tuning:

When defining the coefficients, σ may be underestimated, and L, L 1 , . . . , L n may be overestimated if exact values are unavailable. Notice that x k can be eliminated from the above iteration, and the block gradient ∇ i k f (ŷ k ) only needs to be calculated once per iteration. A larger (or overestimated) maximum delay τ will cause a larger asynchronicity parameter ψ, which leads to more conservative step sizes to compensate.

To estimate ψ, one can first performed a dry run with all coefficient set to 0 to estimate τ . All function parameters can be calculated exactly for this problem in terms of the data matrix and λ. We can then use these parameters and this tau to calculate ψ. ψ and τ merely change the parameters, and do not change execution patterns of the processors. Hence their parameter specification doesn't affect the observed delay. Through simple tuning though, we found that ψ = 0.25 resulted in good performance.

In tuning for general problems, there are theoretical reasons why it is difficult to attain acceleration without some prior knowledge of σ, the strong convexity modulus Arjevani (2017) . Ideally σ is pre-specified for instance in a regularization term. If the Lipschitz constants L i cannot be calculated directly (which is rarely the case for the classic dual problem of empirical risk minimization objectives), the line-search method discussed in Roux et al. (2012) Section 4 can be used.

A.2 Sparse update formulation:

As mentioned in Section 5, authors in Lee & Sidford (2013a) proposed a linear transformation of an accelerated RBCD scheme that results in sparse coordinate updates. Our proposed algorithm can be given a similar efficient implementation. We may eliminate x k from A2BCD, and derive the equivalent iteration below:

y k+1 v k+1 = 1 -αβ, αβ 1 -β, β y k v k -   ασ -1/2 L -1/2 i k + h(1 -α)L -1 i k ∇ i k f ŷk σ -1/2 L -1/2 i k ∇ i k f ŷk   C y k v k -Q k

where C and Q k are defined in the obvious way. Hence we define auxiliary variables p k , q k defined via:

y k v k = C k p k q k (A.1)

These clearly follow the iteration:

p k+1 q k+1 = p k q k -C -(k+1) Q k (A.2)

Since the vector Q k is sparse, we can evolve variables p k , and q k in a sparse manner, and recover the original iteration variables at the end of the algorithm via A.1.

The gradient of the dual function is given by:

∇D(y) = 1 λd 1 d A T Ay + λ(y + l)

As mentioned before, it is necessary to maintain or recover Ay k to calculate block gradients. Since Ay k can be recovered via the linear relation in equation A.1, and the gradient is an affine function, we maintain the auxiliary vectors Ap k and Aq k instead.

Hence we propose the following efficient implementation in Algorithm 1. We used this to generate the results in Table 5 . We also note also that it can improve performance to periodically recover v k and y k , reset the values of p k , q k , and C to v k , y k , and I respectively, and restarting the scheme (which can be done cheaply in time O(d)).

We let B ∈ R 2×2 represent C k , and b represent B -1 . ⊗ is the Kronecker product. Each computing node has local outdated versions of p, q, Ap, Aq which we denote p, q, Âp, Âq respectively. We also find it convenient to define:

D k 1 D k 2 = ασ -1/2 L -1/2 i k + h(1 -α)L -1 i k σ -1/2 L -1/2 i k (A.3)

Algorithm 1 Shared-memory implementation of A2BCD Randomly select block i via equation 2.1.

8::

Read shared data into local memory: p ← p, q ← q, Âp ← Ap, Âq ← Aq, B ← B.

9::

Compute block gradient:

∇ i f (ŷ) = 1 nλ 1 n A T i B1,1 Âp + B1,2 Âq + λ B1,1 p + B1,2 q 10: Compute quantity g i = A T i ∇ i f (ŷ) Shared memory updates: 11: Update B ← 1 -αβ αβ 1 -β β × B, calculate inverse b ← B -1 .

12:

p q -= b D k 1 D k 2 ⊗ ∇ i f (ŷ) , Ap Aq -= b D k 1 D k 2 ⊗ g i 13:

Increase iteration count: k ← k + 1 14: end while 15: Recover original iteration variables:

y v ← B p q . Output y.

B Proof of the main result:

We first recall a couple of inequalities for convex functions.

Lemma 7. Let f be σ-strongly convex with L-Lipschitz gradient. Then we have:

f (y) ≤ f (x) + y -x, ∇f (x) + 1 2 L y -x 2 , ∀x, y (B.1) f (y) ≥ f (x) + y -x, ∇f (x) + 1 2 σ y -x 2 , ∀x, y (B.2)

We also find it convenient to define the norm:

s * = n i=1 L -1/2 i s i 2 (B.3)

Published as a conference paper at ICLR 2019

B.1 Starting point

B.1 Starting point:

First notice that using the definition equation 2.8 of v k+1 we have:

v k+1 2 = βv k + (1 -β)y k 2 -2σ -1/2 L -1/2 i k βv k + (1 -β)y k , ∇ i k f (ŷ k ) + σ -1 L -1 i k ∇ i k f (ŷ k ) 2 E k v k+1 2 = βv k + (1 -β)y k 2 A -2σ -1/2 S -1 βv k + (1 -β)y k , ∇f (ŷ k ) B (B.4) + S -1 σ -1 n i=1 L -1/2 i ∇ i f (ŷ k ) 2 C

We have the following general identity:

βx + (1 -β)y 2 = β x 2 + (1 -β) y 2 -β(1 -β) x -y 2 , ∀x, y (B.5)

It can also easily be verified from equation 2.6 that we have: 

v k = y k + α -1 (1 -α)(y k -x k ) (B.

E k v k+1 2 = β v k 2 + (1 -β) y k 2 -β(1 -β) v k -y k 2 + S -1 σ -1/2 ∇f (ŷ k ) 2 * (B.7) -2σ -1/2 S -1 βα -1 (1 -α) y k -x k , ∇f (ŷ k ) -2σ -1/2 S -1 y k , ∇f (ŷ k )

This inequality is our starting point. We analyze the terms on the second line in the next section.

B.2 The Cross Term:

To analyze these terms, we need a small lemma. This lemma is fundamental in allowing us to deal with asynchronicity. Lemma 8. Let χ, A > 0. Let the delay be bounded by τ . Then:

A ŷk -y k ≤ 1 2 χ -1 A 2 + 1 2 χτ τ j=1 y k+1-j -y k-j 2

Proof. See Hannah & Yin (2017a) .

Lemma 9. We have:

-∇f (ŷ k ), y k ≤ -f (y k ) - 1 2 σ(1 -ψ) y k 2 + 1 2 Lκψ -1 τ τ j=1 y k+1-j -y k-j 2 (B.8) ∇f (ŷ k ), x k -y k ≤ f (x k ) -f (y k ) (B.9) + 1 2 Lα(1 -α) -1   κ -1 ψβ v k -y k 2 + κψ -1 β -1 τ τ j=1 y k+1-j -y k-j 2  

The terms in bold in equation B.8 and equation B.9 are a result of the asynchronicity, and are identically 0 in its absence.

Published as a conference paper at ICLR 2019

B.3 Function-value term:

Proof. Our strategy is to separately analyze terms that appear in the traditional analysis of Nesterov (2012) , and the terms that result from asynchronicity. We first prove equation B.8: 

-∇f (ŷ k ), y k = -∇f (y k ), y k -∇f (ŷ k ) -∇f (y k ), y k ≤ -f (y k ) - 1 2 σ y k 2 + L ŷk -y k y k (B.

L ŷk -y k y k ≤ 1 2 Lκ -1 ψ y k 2 + 1 2 Lκψ -1 τ τ j=1 y k+1-j -y k-j 2

using Lemma 8 with χ = κψ -1 , A = y k . Combining this with equation B.10 completes the proof of equation B.8.

We now prove equation B.9:

∇f (ŷ k ), x k -y k = ∇f (y k ), x k -y k + ∇f (ŷ k ) -∇f (y k ), x k -y k ≤ f (x k ) -f (y k ) + L ŷk -y k x k -y k ≤ f (x k ) -f (y k ) + 1 2 L   κ -1 ψβα -1 (1 -α) x k -y k 2 + κψ -1 β -1 α(1 -α) -1 τ τ j=1 y k+1-j -y k-j 2  

Here the last line follows from Lemma 8 with

χ = κψ -1 β -1 α(1 -α) -1 , A = nx k -y k .

We can complete the proof using the following identity that can be easily obtained from equation 2.6:

y k -x k = α(1 -α) -1 (v k -y k )

Much like Nesterov (2012) , we need a f (x k ) term in the Lyapunov function (see the middle of page 357). However we additionally need to consider asynchronicity when analyzing the growth of this term. Again terms due to asynchronicity are emboldened.

Lemma 10. We have:

E k f (x k+1 ) ≤ f (y k ) - 1 2 h 2 -h 1 + 1 2 σ 1/2 L -1/2 ψ S -1 ∇f (ŷ k ) 2 * + S -1 Lσ 1/2 κψ -1 τ τ j=1 y k+1-j -y k-j 2

Proof. From the definition equation 2.7 of x k+1 , we can see that x k+1 -y k is supported on block i k . Since each gradient block ∇ i f is L i Lipschitz with respect to changes to block i, we can use equation B.1 to obtain:

f (x k+1 ) ≤ f (y k ) + ∇f (y k ), x k+1 -y k + 1 2 L i k x k+1 -y k 2 (from equation 2.7) = f (y k ) -hL -1 i k ∇ i k f (y k ), ∇ i k f (ŷ k ) + 1 2 h 2 L -1 i k ∇ i k f (ŷ k ) 2 = f (y k ) -hL -1 i k ∇ i k f (y k ) -∇ i k f (ŷ k ), ∇ i k f (ŷ k ) - 1 2 h(2 -h)L -1 i k ∇ i k f (ŷ k ) 2 E k f (x k+1 ) ≤ f (y k ) -hS -1 n i=1 L -1/2 i ∇ i f (y k ) -∇ i f (ŷ k ), ∇ i f (ŷ k ) - 1 2 h(2 -h)S -1 ∇f (ŷ k ) 2 * (B.11)

Here the last line followed from the definition equation B.3 of the norm • * 1/2 . We now analyze the middle term:

-

n i=1 L -1/2 i ∇ i f (y k ) -∇ i f (ŷ k ), ∇ i f (ŷ k ) = - n i=1 L -1/4 i (∇ i f (y k ) -∇ i f (ŷ k )), n i=1 L -1/4 i ∇ i f (ŷ k ) (Cauchy Schwarz) ≤ n i=1 L -1/4 i (∇ i f (y k ) -∇ i f (ŷ k )) n i=1 L -1/4 i ∇ i f (ŷ k ) = n i=1 L -1/2 i ∇ i f (y k ) -∇ i f (ŷ k ) 2 1/2 n i=1 L -1/2 i ∇ i f (ŷ k ) 2 1/2 (L ≤ L i , ∀i and equation B.3) ≤ L -1/4 ∇f (y k ) -∇f (ŷ k ) ∇f (ŷ k ) * (∇f is L-Lipschitz) ≤ L -1/4 L y k -ŷk ∇f (ŷ k ) *

We then apply Lemma 8 to this with χ = 2h -1 σ 1/2 L 1/4 κψ -1 , A = ∇f (ŷ k ) * to yield:

-

n i=1 L -1/2 i ∇ i f (y k ) -∇ i f (ŷ k ), ∇ i f (ŷ k ) ≤ h -1 Lσ 1/2 κψ -1 τ τ j=1 y k+1-j -y k-j 2 (B.12) + 1 4 hσ 1/2 L -1/2 ψ ∇f (ŷ k ) 2 *

Finally to complete the proof, we combine equation B.11, with equation B.12.

B.4 Asynchronicity error:

The previous inequalities produced difference terms of the form y k+1-j -y k-j 2 . The following lemma shows how these errors can be incorporated into a Lyapunov function. Lemma 11. Let 0 < r < 1 and consider the asynchronicity error and corresponding coefficients:

A k = ∞ j=1 c j y k+1-j -y k-j 2 c i = ∞ j=i r i-j-1 s j This sum satisfies: E k [A k+1 -rA k ] = c 1 E k y k+1 -y k 2 - ∞ j=1 s j y k+1-j -y k-j 2

Remark 2. Interpretation. This result means that an asynchronicity error term A k can negate a series of difference terms -∞ j=1 s j y k+1-j -y k-j 2 at the cost of producing an additional error

c 1 E k y k+1 -y k 2

, while maintaining a convergence rate of r. This essentially converts difference terms, which are hard to deal with, into a y k+1 -y k 2 term which can be negated by other terms in the Lyapunov function. The proof is straightforward.

Proof.

E k [A k+1 -rA k ] = E k ∞ j=0 c j+1 y k+1-j -y k-j 2 -rE k ∞ j=1 c j y k+1-j -y k-j 2 = c 1 E k y k+1 -y k 2 + E k ∞ j=1 (c j+1 -rc j ) y k+1-j -y k-j 2

Noting the following completes the proof:

c i+1 -rc i = ∞ j=i+1 r i+1-j-1 s j -r ∞ j=i r i-j-1 s j = -s i

Given that A k allows us to negate difference terms, we now analyze the cost c 1 E k y k+1 -y k 2 of this negation.

Lemma 12. We have:

E k y k+1 -y k 2 ≤ 2α 2 β 2 v k -y k 2 + 2S -1 L -1 ∇f (ŷ k ) 2 Proof. y k+1 -y k = (αv k+1 + (1 -α)x k+1 ) -y k = α βv k + (1 -β)y k -σ -1/2 L -1/2 i k ∇ i k f (ŷ k ) + (1 -α) y k -hL -1 i k ∇ i k f (ŷ k ) -y k (B.13) = αβv k + α(1 -β)y k -ασ -1/2 L -1/2 i k ∇ i k f (ŷ k ) -αy k -(1 -α)hL -1 i k ∇ i k f (ŷ k ) = αβ(v k -y k ) -ασ -1/2 L -1/2 i k + h(1 -α)L -1 i k ∇ i k f (ŷ k ) y k+1 -y k 2 ≤ 2α 2 β 2 v k -y k 2 + 2 ασ -1/2 L -1/2 i k + h(1 -α)L -1 i k 2 ∇ i k f (ŷ k ) 2 (B.14)

Here equation B.13 following from equation 2.8, the definition of v k+1 . equation B.14 follows from the inequality x + y 2 ≤ 2 x 2 + 2 y 2 . The rest is simple algebraic manipulation.

y k+1 -y k 2 ≤ 2α 2 β 2 v k -y k 2 + 2L -1 i k ασ -1/2 + h(1 -α)L -1/2 i k 2 ∇ i k f (ŷ k ) 2 (L ≤ L i , ∀i) ≤ 2α 2 β 2 v k -y k 2 + 2L -1 i k ασ -1/2 + h(1 -α)L -1/2 2 ∇ i k f (ŷ k ) 2 = 2α 2 β 2 v k -y k 2 + 2L -1 i k L -1 L 1/2 σ -1/2 α + h(1 -α) 2 ∇ i k f (ŷ k ) 2 E y k+1 -y k 2 ≤ 2α 2 β 2 v k -y k 2 + 2S -1 L -1 L 1/2 σ -1/2 α + h(1 -α) 2 ∇f (ŷ k ) 2 *

Finally, to complete the proof, we prove

L 1/2 σ -1/2 α + h(1 -α) ≤ 1. L 1/2 σ -1/2 α + h(1 -α) = h + α L 1/2 σ -1/2 -h

(definitions of h and α: equation 2.3, and equation 2.

5) = 1 - 1 2 σ 1/2 L -1/2 ψ + σ 1/2 S -1 L 1/2 σ -1/2 ≤ 1 -σ 1/2 L -1/2 1 2 ψ -σ -1/2 S -1 L 1 (B.15)

Rearranging the definition of ψ, we have:

S -1 = 1 9 2 ψ 2 L 1 L -3/2 κ -1/2 τ -2 (τ ≥1 and ψ ≤ 1 2 ) ≤ 1 18 2 L 1 L -3/2 κ -1/2

Using this on equation B.15, we have:

L 1/2 ασ -1/2 + h(1 -α) ≤ 1 -σ 1/2 L -1/2 1 2 ψ - 1 18 2 L 1 L -3/2 κ -1/2 σ -1/2 L 1 = 1 -σ 1/2 L -1/2 1 2 ψ - 1 18 2 (L/L) 2 (ψ ≤ 1 2 ) = 1 -σ 1/2 L -1/2 1 24 - 1 18 2 ≤ 1.

This completes the proof.

B.5 Master inequality:

We are finally in a position to bring together all the all the previous results together into a master inequality for the Lyapunov function ρ k (defined in equation 2.11). After this lemma is proven, we will prove that the right hand size is negative, which will imply that ρ k linearly converges to 0 with rate β.

Lemma 13. Master inequality. We have:

E k [ρ k+1 -βρ k ] ≤ + y k 2 × 1 -β -σ -1/2 S -1 σ(1 -ψ) (B.16) + v k -y k 2 ×β 2α 2 βc 1 + S -1 βL 1/2 κ -1/2 ψ -(1 -β) + f (y k ) × c -2σ -1/2 S -1 βα -1 (1 -α) + 1 + f (x k ) ×β 2σ -1/2 S -1 α -1 (1 -α) -c + τ j=1 y k+1-j -y k-j 2 ×S -1 Lκψ -1 τ σ 1/2 2σ -1 + c -s + ∇f (ŷ k ) 2 * ×S -1 σ -1 + 2L -1 c 1 - 1 2 ch 2 -h 1 + 1 2 σ 1/2 L -1/2 ψ

Proof.

E k v k+1 2 -β v k 2 (B.7) = (1 -β) y k 2 -β(1 -β) v k -y k 2 + S -1 σ -1 ∇f (ŷ k ) 2 * -2σ -1/2 S -1 y k , ∇f (ŷ k ) -2σ -1/2 S -1 βα -1 (1 -α) y k -x k , ∇f (ŷ k ) ≤ (1 -β) y k 2 -β(1 -β) v k -y k 2 + S -1 σ -1 ∇f (ŷ k ) 2 * (B.17) (B.8) + 2σ -1/2 S -1   -f (y k ) - 1 2 σ(1 -ψ) y k 2 + 1 2 Lκψ -1 τ τ j=1 y k+1-j -y k-j 2   (B.9) -2σ -1/2 S -1 βα -1 (1 -α)(f (x k ) -f (y k )) + σ -1/2 S -1 βL   κ -1 ψβ v k -y k 2 + κψ -1 β -1 τ τ j=1 y k+1-j -y k-j 2  

We now collect and organize the similar terms of this inequality.

≤ + y k 2 × 1 -β -σ -1/2 S -1 σ(1 -ψ) + v k -y k 2 ×β σ -1/2 S -1 βLκ -1 ψ -(1 -β) -f (y k ) ×2σ -1/2 S -1 βα -1 (1 -α) + 1 + f (x k ) ×2σ -1/2 S -1 βα -1 (1 -α) + τ j=1 y k+1-j -y k-j 2 ×2σ -1/2 S -1 Lκψ -1 τ + ∇f (ŷ k ) 2 * ×σ -1 S -1

Now finally, we add the function-value and asynchronicity terms to our analysis. We use Lemma 11 is with r = 1 -σ 1/2 S -1 , and

s i = s = 6S -1 L 1/2 κ 3/2 ψ -1 τ, 1 ≤ i ≤ τ 0, i > τ (B.18)

Notice that this choice of s i will recover the coefficient formula given in equation 2.9. Hence we have:

E k [cf (x k+1 ) + A k+1 -β(cf (x k ) + A k )] (Lemma 10) ≤ cf (y k ) - 1 2 ch 2 -h 1 + 1 2 σ 1/2 L -1/2 ψ S -1 ∇f (ŷ k ) 2 * -βcf (x k ) (B.19) + S -1 Lσ 1/2 κψ -1 τ τ j=1 y k+1-j -y k-j 2 (Lemmas 11 and 12) + c 1 2α 2 β 2 v k -y k 2 + 2S -1 L -1 ∇f (ŷ k ) 2 (B.20) - ∞ j=1 s j y k+1-j -y k-j 2 + A k (r -β)

Notice A k (r -β) ≤ 0. Finally, combining equation B.17 and equation B.19 completes the proof.

In the next section, we will prove that every coefficient on the right hand side of equation B.16 is 0 or less, which will complete the proof of Theorem 1.

B.6 Proof of main theorem:

Lemma 14.

The coefficients of y k 2 , f (y k ), and τ j=1 y k+1-j -y k-j 2 in Lemma 13 are non-positive.

Proof. The coefficient 1 -(1 -ψ)σ 1/2 S -1 -β of y k 2 is identically 0 via the definition equation 2.4 of β. The coefficient c -2σ -1/2 S -1 βα -1 (1 -α) + 1 of f (y k ) is identically 0 via the definition equation 2.12 of c.

First notice from the definition equation 2.12 of c:

c = 2σ -1/2 S -1 βα -1 (1 -α) + 1 (definitions of α, β) = 2σ -1/2 S -1 1 -σ 1/2 S -1 (1 -ψ) (1 + ψ)σ -1/2 S + 1 = 2σ -1/2 S -1 (1 + ψ)σ -1/2 S + ψ 2 = 2σ -1 (1 + ψ) + ψ 2 σ 1/2 S -1 (B.21) c ≤ 4σ -1 (B.22)

Here the last line followed since ψ ≤ 1 2 and σ 1/2 S -1 ≤ 1. We now analyze the coefficient of

τ j=1 y k+1-j -y k-j 2 . S -1 Lκψ -1 τ σ 1/2 2σ -1 + c -s (B.22) ≤ 6L 1/2 κ 3/2 ψ -1 τ -s (definition equation B.18 of s) ≤ 0 Lemma 15. The coefficient β 2σ -1/2 S -1 α -1 (1 -α) -c of f (x k ) in Lemma 13 is non-positive. Proof. 2σ -1/2 S -1 α -1 (1 -α) -c (B.21) = 2σ -1/2 S -1 (1 + ψ)σ -1/2 S -2σ -1 (1 + ψ) + ψ 2 σ 1/2 S -1 = 2σ -1 (1 + ψ) -(1 + ψ) + ψ 2 σ 1/2 S -1 = -2ψ 2 σ -1/2 S -1 ≤ 0 Lemma 16. The coefficient S -1 σ -1 + 2L -1 c 1 -1 2 ch 2 -h 1 + 1 2 σ 1/2 L -1/2 ψ of ∇f (ŷ k ) 2 *

in Lemma 13 is non-positive.

Proof. We first need to bound c 1 .

(equation B.18 and equation 2.9)

c 1 = s τ j=1 1 -σ 1/2 S -1 -j equation B.18 ≤ 6S -1 L 1/2 κ 3/2 ψ -1 τ τ j=1 1 -σ 1/2 S -1 -j ≤ 6S -1 L 1/2 κ 3/2 ψ -1 τ 2 1 -σ 1/2 S -1 -τ

It can be easily verified that if x ≤ 1 2 and y ≥ 0, then (1 -x) -y ≤ exp(2xy). Using this fact with x = σ 1/2 S -1 and y = τ , we have:

≤ 6S -1 L 1/2 κ 3/2 ψ -1 τ 2 exp τ σ 1/2 S -1 (since ψ ≤ 3/7 and hence τ σ 1/2 S -1 ≤ 1 7 ) ≤ S -1 L 1/2 κ 3/2 ψ -1 τ 2 × 6 exp 1 7 c 1 ≤ 7S -1 L 1/2 κ 3/2 ψ -1 τ 2 (B.23)

We now analyze the coefficient of ∇f (ŷ k )

2 * σ -1 + 2L -1 c 1 - 1 2 ch 2 -h 1 + 1 2 σ 1/2 L -1/2 ψ (B.23 and 2.5) ≤ σ -1 + 14S -1 L -1 L 1/2 κ 3/2 ψ -1 τ 2 - 1 2 ch 1 + 1 4 σ 1 L -1 ψ 2 ≤ σ -1 + 14S -1 L -1 L 1/2 κ 3/2 ψ -1 τ 2 - 1 2 ch (definition 2.2 of ψ) = σ -1 + 14 81 σ -1 ψ - 1 2 ch (B.21, definition 2.5 of h) = σ -1 1 + 14 81 ψ -(1 + ψ) + ψ 2 σ 1/2 S -1 1 - 1 2 σ 1/2 L -1/2 ψ (σ 1/2 L -1/2 ≤ 0 and σ 1/2 S -1 ≤ 1) ≤ σ -1 1 + 14 81 ψ -(1 + ψ) 1 - 1 2 ψ = σ -1 ψ 14 81 + 1 2 ψ - 1 2 (ψ ≤ 1 2 ) ≤ 0 Lemma 17. The coefficient β 2α 2 βc 1 + S -1 βL 1/2 κ -1/2 ψ -(1 -β) of v k -y k 2 in 13 is non- positive.

Proof.

2α 2 βc 1 + σ 1/2 S -1 βψ -(1 -ψ)σ 1/2 S -1 (B.23) ≤ 14α 2 βS -1 L 1/2 κ 3/2 ψ -1 τ 2 + σ 1/2 S -1 βψ -(1 -ψ)σ 1/2 S -1 ≤ 14σS -3 L 1/2 κ 3/2 ψ -1 τ 2 + σ 1/2 S -1 ψ -(1 -ψ)σ 1/2 S -1 = σ 1/2 S -1 14S -2 Lκτ 2 ψ -1 + 2ψ -1

Here the last inequality follows since β ≤ 1 and α ≤ σ 1/2 S -1 . We now rearrange the definition of ψ to yield the identity:

S -2 κ = 1 9 4 L 2 L -3 τ -4 ψ 4

Using this, we have:

14S -2 Lκτ 2 ψ -1 + 2ψ -1 = 14 9 4 L 2 L -2 ψ 3 τ -2 + 2ψ -1 ≤ 14 9 4 3 7 3 1 -2 + 6 7 -1 ≤ 0

Here the last line followed since L ≤ L, ψ ≤ 3 7 , and τ ≥ 1. Hence the proof is complete.

Proof of Theorem 1. Using the master inequality 13 in combination with the previous Lemmas 14, 15, 16, and 17, we have:

E k [ρ k+1 ] ≤ βρ k = 1 -(1 -ψ)σ 1/2 S -1 ρ k

Published as a conference paper at ICLR 2019 When we have:

1

-(1 -ψ)σ 1/2 S -1 k ≤

then the Lyapunov function ρ k has decreased below ρ 0 in expectation. Hence the complexity K( ) satisfies:

K( ) ln 1 -(1 -ψ)σ 1/2 S -1 = ln( ) K( ) = -1 ln 1 -(1 -ψ)σ 1/2 S -1 ln(1/ )

Now it can be shown that for 0 < x ≤ 1 2 , we have:

1 x -1 ≤ -1 ln(1 -x) ≤ 1 x - 1 2 -1 ln(1 -x) = 1 x + O(1)

Since n ≥ 2, we have σ 1/2 S -1 ≤ 1 2 . Hence:

K( ) = 1 1 -ψ σ -1/2 S + O(1) ln(1/ )

An expression for K NU_ACDM ( ), the complexity of NU_ACDM follows by similar reasoning.

K NU_ACDM ( ) = σ -1/2 S + O(1) ln(1/ ) (B.24)

Finally we have:

K( ) = 1 1 -ψ σ -1/2 S + O(1) σ -1/2 S + O(1) K NU_ACDM ( ) = 1 1 -ψ (1 + o(1))K NU_ACDM ( )

which completes the proof.

C Ordinary Differential Equation Analysis:

C.1 Derivation of ODE for synchronous A2BCD

If we take expectations with respect to E k , then synchronous (no delay) A2BCD becomes:

y k = αv k + (1 -α)x k E k x k+1 = y k -n -1 κ -1 ∇f (y k ) E k v k+1 = βv k + (1 -β)y k -n -1 κ -1/2 ∇f (y k )

We find it convenient to define η = nκ 1/2 . Inspired by this, we consider the following iteration: We now look at some of the terms in this equation to find the highest-order dependence on s.

y k = αv k + (1 -α)x k (C.1) x k+1 = y k -s 1/2 κ -1/2 η -1 ∇f (y k ) (C.2) v k+1 = βv k + (1 -β)y k -s 1/2 η -1 ∇f (y k ) (C.3)

β(1 -α) = 1 -s 1/2 η -1 1 - 1 1 + s -1/2 η = 1 -s 1/2 η -1 s -1/2 η 1 + s -1/2 η = s -1/2 η -1 s -1/2 η + 1 = 1 -s 1/2 η -1 1 + s 1/2 η -1 = 1 -2s 1/2 η -1 + O(s)

We also have:

(β -1)(1 -α) -α = β(1 -α) -1 = -2s 1/2 η -1 + O(s)

Hence using these facts on equation C.4, we have:

0 = -s 1/2 Ẏ - 1 2 s Ÿ + 1 -2s 1/2 η -1 + O(s) s 1/2 Ẏ - 1 2 s Ÿ + s 1/2 η -1 -2s 1/2 η -1 + O(s) ∇f (y k ) + O s 3/2 0 = -s 1/2 Ẏ - 1 2 s Ÿ + s 1/2 Ẏ - 1 2 s Ÿ -2s 1 η -1 Ẏ + O s 3/2 -2s 1 η -2 + O s 3/2 ∇f (y k ) + O s 3/2 0 = -s Ÿ -2sη -1 Ẏ -2sη -2 ∇f (y k ) + O s 3/2 0 = -Ÿ -2η -1 Ẏ -2η -2 ∇f (y k ) + O s 1/2

Taking the limit as s → 0, we obtain the ODE:

Ÿ (t) + 2η -1 Ẏ + 2η -2 ∇f (Y ) = 0 C.2 Convergence proof for synchronous ODE e -η -1 t E (t) = ∇f (Y (t)), Ẏ (t) + η -1 f (Y (t)) + 1 2 Y (t) + η Ẏ (t), Ẏ (t) + η Ÿ (t) + η -1 1 4 Y (t) + η Ẏ (t) 2 (strong convexity equation B.2) ≤ ∇f (Y ), Ẏ + η -1 ∇f (Y ), Y - 1 2 η -1 Y 2 + 1 2 Y + η Ẏ , -Ẏ -2η -1 ∇f (Y ) + η -1 1 4 Y (t) + η Ẏ (t) 2 = - 1 4 η -1 Y 2 - 1 4 η Ẏ 2 ≤ 0

Published as a conference paper at ICLR 2019

C.3 Asynchronicity error lemma:

Hence we have E (t) ≤ 0. Therefore E(t) ≤ E(0). That is:

E(t) = e n -1 κ -1/2 t f (Y ) + 1 4 Y + η Ẏ 2 ≤ E(0) = f (Y (0)) + 1 4 Y (0) + η Ẏ (0) 2 (C.5) which implies: f (Y (t)) + 1 4 Y (t) + η Ẏ (t) 2 ≤ e -n -1 κ -1/2 t f (Y (0)) + 1 4 Y (0) + η Ẏ (0) 2 (C.6)

This result is the continuous-time analogue of Lemma 11. First notice that c(0) = c 0 and c(τ ) = 0. We also have:

c (t)/c 0 = -re -rt -re -rt e -rτ 1 -e -rτ

= -r e -rt + e -rt e -rτ 1 -e -rτ = -r e -rt + e -rt -1 e -rτ 1 -e -rτ + e -rτ 1 -e -rτ c (t) = -rc(t) -rc 0 e -rτ 1 -e -rτ Hence using c(τ ) = 0:

A (t) = c 0 Ẏ (t) 2 + t t-τ c (t -s) Ẏ (s) 2 ds = c 0 Ẏ (t) 2 -rA(t) -rc 0 e -rτ 1 -e -rτ D(t) Now when x ≤ 1 2 , we have e -x

1-e -x ≥ 1 2 x -1 . Hence when rτ ≤ 1 2 , we have:

A (t) ≤ c 0 Ẏ (t) 2 -rA(t) - 1 2 τ -1 c 0 D(t)

and the result easily follows.

C.4 Convergence analysis for the asynchronous ODE:

We consider the same energy as in the synchronous case (that is, the ODE in equation 3.1). Similar to before, we have:

e -η -1 t E (t) ≤ ∇f (Y ), Ẏ + η -1 ∇f (Y ), Y - 1 2 η -1 Y 2 + 1 2 Y + η Ẏ , -Ẏ -2η -1 ∇f Ŷ + η -1 1 4 Y (t) + η Ẏ (t) 2 = ∇f (Y ), Ẏ + η -1 ∇f (Y ), Y - 1 2 η -1 Y 2 + 1 2 Y + η Ẏ , -Ẏ -2η -1 ∇f (Y ) + η -1 1 4 Y (t) + η Ẏ (t) 2 -η -1 Y + η Ẏ , ∇f Ŷ -∇f (Y ) = - 1 4 η -1 Y 2 - 1 4 η Ẏ 2 -η -1 Y + η Ẏ , ∇f Ŷ -∇f (Y )

where the final equality follows from the proof in Section C.2. Hence

e -η -1 t E (t) ≤ - 1 4 η -1 Y 2 - 1 4 η Ẏ 2 + Lη -1 Y Ŷ -Y + L Ẏ Ŷ -Y (C.7)

Now we present an inequality that is similar to equation 8.

Lemma 18. Let A, χ > 0. Then:

Y (t) -Ŷ (t) A ≤ 1 2 χτ D(t) + 1 2 χ -1 A 2 Proof. Since Ŷ (t) is a delayed version of Y (t), we have: Ŷ (t) = Y (t -j(t))

for some function j(t) ≥ 0 (though this can be easily generalized to an inconsistent read). Recall that for χ > 0, we have ab

≤ 1 2 χa 2 + χ -1 b 2 . Hence X(t) -X(t) = t s=t-j(t) X (s)ds X(t) -X(t) A = t s=t-j(t) X (s)ds A ≤ 1 2 χ t s=t-j(t) X (s)ds 2 + 1 2 χ -1 A 2 (Holder's inequality) ≤ 1 2 χ t s=t-j(t) X (s) 2 ds t s=t-j(t) 1ds + 1 2 χ -1 A 2 ≤ 1 2 χτ t s=t-j(t) X (s) 2 ds + 1 2 χ -1 A 2

We use this lemma twice on Y Ŷ -Y and Ẏ Ŷ -Y in equation C.7 with χ = 2L, A = Y and χ = 4Lη -1 , A = Ẏ respectively, to yield:

e -η -1 t E (t) ≤ - 1 4 η -1 Y 2 - 1 4 η Ẏ 2 + Lη -1 Lτ D(t) + 1 4 L -1 Y 2 + L 2Lη -1 τ D(t) + 1 8 L -1 η Ẏ 2 = - 1 8 η Ẏ 2 + 3L 2 η -1 τ D(t)

The proof of convergence is completed in Section 3.

D Optimality proof:

For parameter set σ, L 1 , . . . , L n , n, we construct a block-separable function f on the space R bn (separated into n blocks of size b), which will imply this lower bound. Define κ i = L i /σ. We define Published as a conference paper at ICLR 2019 the matrix A i ∈ R b×b via:

A i          2 -1 0 -1 2 . . . . . . 0 . . . . . . -1 0 . . . -1 2 -1 0 -1 θ i          , for θ i = κ 1/2 i + 3 κ 1/2 i + 1 .

Hence we define f i on R b via:

f i = L i -σ 4 1 2 x, A i x -e 1 , x + σ 2 x 2

which is clearly σ-strongly convex and L i -Lipschitz on R b . From Lemma 8 of Lan & Zhou (2015) , we know that this function has unique minimizer

x * ,(i) q i , q 2 i , . . . , q b i , for q = κ 1/2 i -1

κ 1/2 i + 1

.

Finally, we define f via:

f (x) n i=1 f i x (i) .

Now let e(i, j) be the jth unit vector of the ith block of size b in R bn . For I 1 , . . . , I n ∈ N, we define the subspaces V i (I) = span{e(i, 1), . . . , e(i, I)}, V (I 1 , . . . , I n ) = V 1 (I 1 ) ⊕ . . . ⊕ V n (I n ).

V (I 1 , . . . , I n ) is the subspace with the first I 1 components of block 1 nonzero, the first I 2 components of block 2 nonzero, etc. First notice that IC(V (I 1 , . . . , I n )) = V (I 1 , . . . , I n ). Also, clearly, we have: ∇ i f (V (I 1 , . . . , I n )) ⊂ V (0, . . . , 0, min{I i + 1, b}, 0, . . . , 0).

(D.1)

∇ i f is supported on the ith block, hence why all the other indices are 0. The patten of nonzeros in A means that the gradient will have at most 1 more nonzero on the ith block (see Nesterov (2013) ).

Let the initial point x 0 belong to V Ī1 , . . . , Īn . Let I K,i be the number of times we have had i k = i for k = 0, . . . , K -1. By induction on condition 2 of Definition 4 using equation D.1, we have: Hence if x 0,(i) ∈ V i (0) and k ≤ b, then

x k ∈ V min Ī1 + I k,

x k,(i) -x * ,(i) 2 ≥ min x∈Vi(I k,i ) x -x * ,(i) 2 = b j=I k,i +1 q 2j i = q 2I k,i +2 i -q 2b+2 i / 1 -q 2 i

Therefore for all i, we have:

E x k -x * 2 ≥ E x k,(i) -x * ,(i) 2 ≥ E q 2I k,i +2 i -q 2b+2 i / 1 -q 2 i

To evaluate this expectation, we note:

Eq 2I k,i i = k j=0 k j p j i (1 -p i ) k-j q 2j i = (1 -p i ) k k j=0 k j q 2 i p i (1 -p i ) -1 j = (1 -p i ) k 1 + q 2 i p i (1 -p i ) -1 k = 1 -1 -q 2 i p i k Hence E x k -x * 2 ≥ 1 -1 -q 2 i p i k -q 2b i q 2 i / 1 -q 2 i .

For any i, we may select the starting iterate x 0 by defining its block j = 1, . . . , n via:

x 0,(j) = (1 -δ ij )x * ,(j)

For such a choice of x 0 , we have

x 0 -x * 2 = x * ,(i) 2 = q 2 i + . . . + q 2b i = q 2 i 1 -q 2b i 1 -q 2 i

Hence for this choice of x 0 :

E x k -x * 2 / x 0 -x * 2 ≥ 1 -1 -q 2 i p i k -q 2b i / 1 -q 2b i

Now notice:

1 -1 -q 2 i p i k = q -2 i -q -2 i -1 p i k q 2k i ≥ q 2k

i Hence

E x k -x * 2 / x 0 -x * 2 ≥ 1 -1 -q 2 i p i k 1 -q 2b-2k i / 1 -q 2b

i Now if we let b = 2k, then we have:

E x k -x * 2 / x 0 -x * 2 ≥ 1 -1 -q 2 i p i k 1 -q 2k i / 1 -q 4k i = 1 -1 -q 2 i p i k / 1 + q 2k i E x k -x * 2 / x 0 -x * 2 ≥ 1 2 max i 1 -1 -q 2 i p i k

Now let us take the minimum of the right-hand side over the parameters p i , subject to n i=1 p i = 1. The solution to this minimization is clearly: 

p i = 1 -q 2 i -1 /   n j=1 1 -q 2 j -1   Hence E x k -x * 2 / x 0 -x * 2 ≥ 1 2   1 -   n j=1 1 -q 2 j -1   -1    k n j=1 1 -q 2 j -1 = 1 4 n j=1 κ 1/2 i + 2 + κ -1/2 i ≥ 1 4   n j=1 κ 1/2 i + 2n   E x k -x * 2 / x 0 -x * 2 ≥ 1 2 1 - 4 n

E Extensions:

As mentioned, a stronger result than Theorem 1 is possible. In the case when L i = L for all i, we can consider a slight modification of the coefficients:

α 1 + (1 + ψ)σ -1/2 S -1 (E.1) β 1 -(1 + ψ) -1 σ 1/2 S -1 (E.2) h 1 + 1 2 σ 1/2 L -1/2 ψ -1

.

(E.3)

y k = αv k + (1 -α)x k , (E.4)

x k+1 = y k -hL -1 ∇i k f (ŷ k ), (E.5)

v k+1 = βv k + (1 -β)y k -σ -1/2 L -1/2 ∇i k f (ŷ k ). (E.6)

for the asynchronicity parameter:

ψ = 6κ 1/2 n -1 × τ (E.7)

This leads to complexity:

K( ) = (1 + ψ)nκ 1/2 ln(1/ ) (E.8)

Here there is no restriction on ψ as in Theorem 1, and hence there is no restriction on τ . Assuming ψ ≤ 1 gives optimal complexity to within a constant factor. Notice then that the resulting condition of τ τ ≤ 1 6 nκ -1/2 (E.9)

:

Every coordinate can be outdated by a different amount without significantly changing the proofs.

This condition can be relaxed however by techniques inHannah & Yin (2017b);Sun et al. (2017);Peng et al. (2016c);Hannah & Yin (2017a)

Speedup is defined precisely in Section 2

We can assume x * = 0 with no loss in generality since we may translate the coordinate system so that x * is at the origin. We can assume f (x * ) = 0 with no loss in generality, since we can replace f (x) with f (x)-f (x * ). Without this assumption, the Lyapunov function simply becomes: v k -x * 2 + A k + c(f (x k ) -f (x * )).

Their Lyapunov function is in fact a generalization of the one found inNesterov (2012).

For compactness, we have omitted the (t) from time-varying functions Y (t), Ẏ (t), ∇Y (t), etc.

Figures and Tables:

6) Using equation B.5 on term A, equation B.6 on term B, and recalling the definition equation B.3 on term C, we have from equation B.4:

10) equation B.10 follows from strong convexity (equation B.2 with x = y k and y = x * ), and the fact that ∇f is L-Lipschitz. The term due to asynchronicity becomes:

1 , b , . . . , min Īn + I k,m , b

Solution vectors p, q; auxiliary vectors Ap, Aq; sparsifying matrix B 3: Node local data: Solution vectors p, q, auxiliary vectors Âp, Âq, sparsifying matrix B. 4: Calculate parameters ψ, α, β, h via 1. Set k = 0. 5: Initializations: p ← y, q ← v, Ap ← Ay, Aq ← Av, B ← I. 6: while not converged, each computing node asynchronous do

Table content: <table><tr><td>7:</td></tr></table>

Peer review (iclr2019):
 This paper studies the combination of the asynchronous parallelization and the accelerated stochastic coordinate descent method. The proved convergence rate is claimed to be consistent with the non parallel counterpart. The linear speedup is achievable when the maximal staleness is bounded by n^{1/2} roughly, that sounds very interesting result to me. However, I have a few questions about the correctness of the results: - Theorem 1 essentially shows that every single step is guaranteed to improve the last step in the expectation sense. However, this violates my my experiences to study Nesterov*s accelerated methods. To my knowledge, Nesterov*s accelerated methods generally do not guarantee improvement over each single step, because accelerate methods essentially constructs a sequence z_{t+1} = A z_t where A is a nonsymmetric matrix with spectral norm greater than 1. - The actual implemented algorithm is using the sparse update other than the analyzed version, since the analyzed version is not efficient or suitable for parallelization. However, the sparse updating rule is equivalent to the original version only for the non asynchronous version. Therefore, the analysis does not apply the actual implementation. minors: - pp2 line 8, K(epsilon) is not defined - Eq. (1.4), the index is missing. - missing reference: An Asynchronous Parallel Stochastic Coordinate Descent Algorithm, ICML 2014.