On the Convergence of Asynchronous Parallel Iteration with Unbounded Delays

Recent years have witnessed the surge of asynchronous parallel (async-parallel) iterative algorithms due to problems involving very large-scale data and a large number of decision variables. Because of asynchrony, the iterates are computed with outdated information, and the age of the outdated information, which we call delay, is the number of times it has been updated since its creation. Almost all recent works prove convergence under the assumption of a finite maximum delay and set their stepsize parameters accordingly. However, the maximum delay is practically unknown. This paper presents convergence analysis of an async-parallel method from a probabilistic viewpoint, and it allows for large unbounded delays. An explicit formula of stepsize that guarantees convergence is given depending on delays’ statistics. With p+1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p+1$$\end{document} identical processors, we empirically measured that delays closely follow the Poisson distribution with parameter p, matching our theoretical model, and thus, the stepsize can be set accordingly. Simulations on both convex and nonconvex optimization problems demonstrate the validness of our analysis and also show that the existing maximum-delay-induced stepsize is too conservative, often slows down the convergence of the algorithm.


Introduction
In the "big data" era, the size of the dataset and the number of decision variables involved in many areas such as health care, the Internet, economics, and engineering are becoming tremendously large [34]. It motivates the development of new computational approaches by efficiently utilizing modern multi-core computers or computing clusters.
In this paper, we consider the block-structured optimization problem where x = (x 1 , . . . , x m ) is partitioned into m disjoint blocks, f has a Lipschitz continuous gradient (possibly nonconvex), and r i 's are (possibly nondifferentiable) proper closed convex functions. Note that r i 's can be extended-valued, and thus (1) can have block constraints x i ∈ X i by incorporating the indicator function of X i in r i for all i. Many applications can be formulated in the form of (1), and they include classic machine learning problems: support vector machine (squared hinge loss and its dual formulation) [6], LASSO [30], and logistic regression (linear or multilinear) [37], and also subspace learning problems: sparse principal component analysis [38], nonnegative matrix or tensor factorization [5], just to name a few.
Toward solutions for these problems with extremely large-scale datasets and many variables, first-order methods and also stochastic methods become particularly popular because of their scalability to the problem size, such as FISTA [1], stochastic approximation [21], randomized coordinate descent [22], and their combinations [7,35]. Recently, lots of efforts have been made to the parallelization of these methods, and in particular, asynchronous parallel (async-parallel) methods attract more attention (e.g., [16,24]) over their synchronous counterparts partly due to the better speed-up performance.
This paper focuses on the async-parallel block coordinate update (async-BCU) method (see Algorithm 1) for solving (1). To the best of our knowledge, all works on async-BCU before 2013 consider a deterministic selection of blocks with an exception to [29], and thus they require strong conditions (like a contraction) for convergence. Recent works, e.g., [16,17,24,12], employ randomized block selection and significantly weaken the convergence requirement. However, all of them require bounded delays and/or are restricted to convex problems. The work [12] allows unbounded delays but requires convexity, and [8,3] do not assume convexity but require bounded delays. We consider unbounded delays and deal with nonconvex problems.

Algorithm
We describe the async-BCU method as follows. Assume there are p + 1 processors, and the data and variable x are accessible to all processors. We let all processors continuously and asynchronously update the variable x in parallel. At each time k, one processor reads the variable x asx k from the global memory, randomly picks a block i k ∈ {1, 2, · · · , m}, and renews x i k by a proxlinear update while keeping all the other blocks unchanged. The pseudocode is summarized in Algorithm 1, where the prox operator is defined in (3).
The algorithm first appeared in [16], where the age ofx k relative to x k , which we call the delay of iteration k, was assumed to be bounded by a certain integer τ . For general convex problems, sublinear convergence was established, and for the strongly convex case, linear convergence was shown. However, its convergence for nonconvex problems and/or with unbounded delays was unknown. In addition, numerically, the stepsize is difficult to tune because it depends on τ , which is unknown before the algorithm completes.
Algorithm 1: Async-parallel block coordinate update Input : Any point x 0 ∈ R n in the global memory, maximum number of iterations K, stepsize η > 0 while k < K, each and all processors asynchronously do select i k from [m] uniformly at random; x k ← read x from the global memory; for all i ∈ [m], increase the global counter k ← k + 1; end

Contributions
We summarize our contributions as follows.
-We analyze the convergence of Algorithm 1 and allow for large unbounded delays following a certain distribution. We require the delays to have certain bounded expected quantities (e.g., expected delay, variance of delay). Our results are more general than those requiring bounded delays such as [16,17]. -Both nonconvex and convex problems are analyzed, and those problems include both smooth and nonsmooth functions. For nonconvex problems, we establish the global convergence in terms of first-order optimality conditions and show that any limit point of the iterates is a critical point almost surely. It appears to be the first result of an async-BCU method for general nonconvex problems and allowing unbounded delays. For weakly convex problems, we establish a sublinear convergence result, and for strongly convex problems, we show the linear convergence.
-We show that if all p + 1 processors run at the same speed, the delay follows the Poisson distribution with parameter p. In this case, all the relevant expected quantities can be explicitly computed and are bounded. By setting appropriate stepsizes, we can reach a near-linear speedup if p = o( √ m) for smooth cases and p = o( 4 √ m) for nonsmooth cases. -When the delay follows the Poisson distribution, we can explicitly set the stepsize based on the delay expectation (which equals p). We simulate the async-BCU method on one convex problem: LASSO, and one nonconvex problem: the nonnegative matrix factorization. The results demonstrate that async-BCU performs consistently better with a stepsize set based on the expected delay than on the maximum delay. The number of processors is known while the maximum delay is not. Hence, the setting based on expected delay is practically more useful. Our algorithm updates one (block) coordinate of x in each step and is sharply different from stochastic gradient methods that sample one function in each step to update all coordinates of x. While there are async-parallel algorithms in either classes and how to handle delays is important to both of their convergence, their basic lines of analysis are different with respect to how to absorb the delay-induced errors. The results of the two classes are in general not comparable. That said, for problems with certain proper structures, it is possible to apply both coordinate-wise update and stochastic sampling (e.g., [25,35,20,8]), and our results apply to the coordinate part.

Notation and assumptions
Throughout the paper, bold lowercase letters x, y, . . . , are used for vectors. We denote x i as the i-th block of x and U i as the i-th sampling matrix, i.e., U i x is a vector with x i as its i-th block and 0 for the remaining ones. E i k denotes the expectation with respect to i k conditionally on all previous history, and [m] = {1, . . . , m}.
We consider the Euclidean norm denoted by · , but all our results can be directly extended to problems with general primal and dual norms in a Hilbert space.
The projection to a convex set X is defined as P X (y) = arg min x∈X x − y 2 , and the proximal mapping of a convex function h is defined as Throughout our analysis, we make the following three assumptions to problem (1) and Algorithm 1. Other assumed conditions will be specified if needed.

Assumption 1
The function F is lower bounded. The problem (1) has at least one solution, and the solution set is denoted as X * .
In addition, for each i ∈ [m], fixing all block coordinates but the i-th one, ∇f (x) and ∇ i f (x) are Lipschitz continuous about x i with L r and L c , respectively, i.e., for any x, y, and i, From (6), we have that for any x, y, and i, We denote κ = Lr Lc as the condition number. Assumption 3 For each k ≥ 1, the readingx k is consistent and delayed by j k , namely,x k = x k−j k . The delay j k follows an identical distribution as a random variable j Prob(j = t) = q t , t = 0, 1, 2 . . . , (8) and is independent of i k . We let Remark 1 Although the delay always satisfies 0 ≤ j k ≤ k, the assumption in (8) is without loss of generality if we make negative iterates and regard x k = x 0 , ∀k < 0. For simplicity, we make the identical distribution assumption, which is the same as that in [29]. Our results can still hold for nonidentical distribution; see the analysis for the smooth nonconvex case in the arXiv version of the paper.

Related works
We briefly review block coordinate update (BCU) and async-parallel computing methods. The BCU method is closely related to the Gauss-Seidel method for solving linear equations, which can date back to 1823. In the literature of optimization, BCU method first appeared in [13] as the block coordinate descent method, or more precisely, block minimization (BM), for quadratic programming. The convergence of BM was established early for both convex and nonconvex problems, for example [19,10,31]. However, in general, its convergence rate result was only shown for strongly convex problems (e.g., [19]) until the recent work [14] that shows sublinear convergence for weakly convex cases. [33] proposed a new version of BCU methods, called coordinate gradient descent method, which mimics proximal gradient descent but only updates a block coordinate every time. The block coordinate gradient or block prox-linear update (BPU) becomes popular since [22] proposed to randomly select a block to update. The convergence rate of the randomized BPU is easier to show than the deterministic BPU. It was firstly established for convex smooth problems (both unconstrained and constrained) in [22] and then generalized to nonsmooth cases in [26,18]. Recently, [7,35] incorporated stochastic approximation into the BPU framework to deal with stochastic programming, and both established sublinear convergence for convex problems and also global convergence for nonconvex problems.
The async-parallel computing method (also called chaotic relaxation) first appeared in [28] to solve linear equations arising in electrical network problems. [4] first systematically analyzed (more general) asynchronous iterative methods for solving linear systems. Assuming bounded delays, it gave a necessary and sufficient condition for convergence. [2] proposed an asynchronous distributed iterative method for solving more general fixed-point problems and showed its convergence under a contraction assumption. [32] weakened the contraction assumption to pseudo-nonexpansiveness but made more other assumptions. [9] made a thorough review of asynchronous methods before 2000. It summarized convergence results under nested sets and synchronous convergence conditions, which are satisfied by P-contraction mappings and isotone mappings.
Since it was proposed in 1969, the async-parallel method has not attracted much attention until recent years when the size of data is increasing exponentially in many areas. Motivated by "big data" problems, [16,17] proposed the async-parallel stochastic coordinate descent method (i.e., Algorithm 1) for solving problems in the form of (1). Their analysis focuses on convex problems and assumes bounded delays. Specifically, they established sublinear convergence for weakly convex problems and linear convergence for strongly convex problems. In addition, near-linear speed up was achieved if τ = o( √ m) for unconstrained smooth convex problems and τ = o( 4 √ m) for constrained smooth or nonsmooth cases. For nonconvex problems, [8] introduced an async-parallel coordinate descent method, whose convergence was established under iterate boundedness assumptions and appropriate stepsizes.

Convergence results for the smooth case
Throughout this section, let r i = 0, ∀i, i.e., we consider the smooth optimization problem minimize The general (possibly nonsmooth) case will be analyzed in the next section.
The results for nonsmooth problems of course also hold for smooth ones. How-ever, the smooth case requires weaker conditions for convergence than those required by the nonsmooth case, and their analysis techniques are different. Hence, we consider the two cases separately.

Convergence for the nonconvex case
In this subsection, we establish a subsequence convergence result for the general (possibly nonconvex) case. We begin with some technical lemmas. The first lemma deals with certain infinite sums that will appear later in our analysis.
Lemma 1 For any k and t ≤ k, let Proof To bound ∞ k=0 γ k , we bound the first term where the last equality holds since T : We obtain (11) by combining these two equations. To prove (12), we will use The above inequality yields where the last inequality follows from c k ≤ 1.
The second lemma below bounds the cross term that appears in our analysis.
Lemma 2 (Cross term bound) For any k > 1 and t ≤ k, it holds that By taking expectation, we have Now taking expectation on both sides of (15) and using the above equation, we get Finally, (14) follows from Using the above lemma, we show a result of running one iteration of the algorithm.
Theorem 1 (Fundamental bound) Set γ k , β k and C t,k as in (10). For any k > 1, we have Proof For the first cross term in (19), we write each summand as and we use Young's inequality to bound the second cross term by Now taking expectation over both sides of (19), plugging in (20) and (21), and using Lemma 2, we have the desired result.
We are now ready to show the main result in the following theorem.
Theorem 2 (Convergence for the nonconvex smooth case) Under Assumptions 1 through 3, let {x k } k≥1 be generated from Algorithm 1. Assume T < ∞. Take the stepsize as 0 < η < and any limit point of {x k } k≥1 is almost surely a critical point of (9).
, then η only weakly depends on the delay. The conditions q 0 > 0 or ∇f (x) being bounded can be dropped if S = E[j 2 ] is bounded; see Theorem 5.
Proof Summing up (18) from k = 0 through K and using (85), we have (23) and using the lower boundedness of f , we have from (22) from the above inequality. From the Markov inequality, it follows that ∇f (x k ) converges to zero with probability one. Letx be a limit point of This completes the proof.

Convergence rate for the convex case
In this subsection, we assume the convexity of f and establish convergence rate results of Algorithm 1 for solving (9). Besides Assumptions 1 through 3, we make an additional assumption to the delay as follows. It means the delay follows a sub-exponential distribution.
Assumption 4 There is a constant σ > 1 such that The condition in (24) is stronger than T < ∞, and both of them hold if the delay j k is uniformly bounded by some number τ or follows the Poisson distribution; see the discussions in Section 5. Using this additional assumption and choosing an appropriate stepsize, we are able to control the gradient of f such that it changes not too fast.

Lemma 3 Under Assumptions 2 through 4, for any
with M ρ defined in (24), then for all k, it holds that The proof of Lemma 3 follows an argument similar to [16]. Since it is rather long, it is included in the appendix. Similar to Lemma 2, we can show the following result.
Lemma 4 For any k, it holds that Proof Following an argument similar to how (16) is obtained, we can show Using the above inequalities, we complete the proof by noting (17), Using the above two lemmas, we establish sufficient objective decrease.
Theorem 3 (Sufficient progress) Under Assumptions 1 through 4, we let {x k } k≥1 be the sequence generated from Algorithm 1. For a certain 1 < ρ < σ, define Take the stepsize such that (25) is satisfied and also Then, Proof First note that for any ρ < σ, tρ t is dominated by σ t as t is sufficiently large. Hence, N ρ < ∞ from (24), and it is easy to see T < ∞. Also note that We write the cross terms in (19) to Taking expectation on both sides of (19) and using (27), we have The above inequality together with (26) implies Note that From these relations and (35), we obtain Using (32) and the convexity of f , we establish the following convergence rate.
Theorem 4 (Convergence rate for the convex smooth case) Under the assumptions of Theorem 3, we have where f * denotes the minimum value of (9) and D is given in (31). 2. If f is strongly convex with constant µ, then where D is given in (31).
Remark 3 For the sublinear rate in (36), we assume the boundedness of the iterates. This assumption can be relaxed if we use potentially smaller stepsize; see Theorem 6. For the linear convergence, the assumption on strongly convexity can be weakened to either essential or restrict strong convexity, see [15] and [16].
Substituting (38) into (32) yields Hence, and thus (36) holds. If f is strongly convex with constant µ, then We immediately have (37) from (32) and the above inequality. This completes the proof.

Convergence results for the nonsmooth case
In this section, we analyze the convergence of Algorithm 1 for possibly nonsmooth cases. Throughout this section, we let a virtual full-update iterate, where R is defined in (4), and denote Due to more generality, we will make stronger assumptions on the delay than those made in the previous section. But all these assumptions are satisfied if the delay is uniformly bounded or follows the Poisson distribution, as shown in Section 5.

Convergence for the nonconvex case
We first establish the almost sure global convergence for possibly nonconvex cases starting with the following square summable result.
Lemma 5 (Square summability) Under Assumptions 1 through 3, we let {x k } k≥1 be the sequence generated in Algorithm 1. Assume S < ∞, and the stepsize is taken as 0 < η < 1/Lc Proof By the definition ofx k+1 , we have −∇f (x k−j k ) − 1 η d k ∈ ∂R(x k+1 ), which together with the convexity of R implies that, for any x, To this inequality, take conditional expectation on i k : To bound the right-hand side, we split the cross term as and apply (40) with x = x k , arriving at Following a similar argument in the proof of Lemma 2 and Young's inequality, we get Note that Hence, taking expectation yields Taking expectation on both sides of (41) and substituting (44) yield From Lemma 12, we have that for any K ≥ 0, Summing up (45) from k = 0 through K and substituting (46) and (47), we have Note that Since F is lower bounded, we have (39) from (48) by letting K → ∞.
, the condition S < ∞ implies T < ∞. Equation (39) indicates that E d k → 0 as k → ∞. Together with S < ∞, we are able to show E x k − x k−j k also approaches zero, as summarized in the following.

Lemma 6 Under the assumptions of Lemma 5, we have
Proof Pick any > 0. From (39), there must exist an integer J > 0 such that For the above J, there must exist an integer K > J such that, for any k ≥ K, From Young's inequality, it follows that Hence, for any k ≥ K, using (43) and (83), we have Using Lemmas 5 and 6, we establish the almost sure global convergence of Algorithm 1.
Theorem 5 Under the assumptions of Lemma 5, any limit point x * of {x k } is a critical point of (1) almost surely.
Before proving this theorem, we make two remarks as follows.
Remark 4 From the theorem, we see that if S = E[j 2 ] = o(m), then the stepsize required for convergence only weakly depends on the delay.

Remark 5 (Comparison of stepsize)
The works [8] consider asynchronous coordinate descent for nonconvex problems. To have convergence to critical points, they assume delays bounded by a number τ . Also, they require the boundedness of iterates and the stepsize less than 1/Lc 1+2κτ / √ m . Note that our stepsize in Theorem 5 is larger if κ 2 S ≤ 16m, where S = E[j 2 ] < τ 2 , and that can lead to faster convergence.
Proof Let {x k } k∈K be a subsequence that converges to x * . Since E d k → 0 as K k → ∞, from the Markov inequality, d k converges to zero in probability as K k → ∞. By [11, Theorem 3.4, pp.212], there is a subsubsequence {x k } k∈K such that d k almost surely converges to zero as K k → ∞. Hence,x k+1 almost surely converges to x * as K k → ∞.
Using triangle inequality and the Lipschitz continuity of ∇f , and taking expectation give From Lemmas 5 and 6, it follows that the right-hand side approaches to zero as k → ∞. Hence, Edist(0, ∂F (x k+1 )) → 0 as k → ∞. If necessary, passing to another subsequence, we use Markov inequality and [11, Theorem 3.4, pp.212] again to have dist(0, ∂F (x k+1 )) almost surely converges to zero as K k → ∞. Now use the outer semicontinuity [27] of dist(0, ∂F (x)) to obtain the desired result.

Convergence rate for the convex case
In this subsection, we establish convergence rates of Algorithm 1 for nonsmooth convex cases. Similar to (26), we first show that choosing an appropriate stepsize, the iterate difference does not change too fast.
Lemma 7 (Fundamental bounds) Assume Assumptions 2 through 4. Then for any 1 < ρ < σ, it holds that In addition, if the stepsize is taken such that then, for all k ≥ 1, Proof It is easy to show (51) by noting that tρ t is dominated by σ t as t is sufficiently large. Next we show (53) by induction.
Using the inequality In addition, for all k, ), the nonexpansiveness of prox ηR , and the triangle inequality, we have When k = 1, we have j 0 = 0 and j 1 ∈ {0, 1} because j k ≤ k, ∀k. Hence, from (56), which together with (54) and (55) implies Hence, Assume (53) holds for all k ≤ K − 1. We show it holds for k = K. First, for any d ≤ K − 1, Secondly, we have Note that By the triangle inequality and the Lipschitz of ∇f , it follows that, for any Since d K−1 is independent of j K , we have from the above two equations that Using (58), the definition of γ ρ,1 in (51) and Also, using Young's inequality and (60) with K replaced by K − 1 and t = j K−1 , we have, for any β > 0, Note that Substituting this inequality into (62), noting E x d − x d+1 2 = 1 m E d d 2 , and applying (53) for all k ≤ K − 1, we have and recall the definition of γ ρ,2 in (51). From the above inequality, we have Substituting (61) and (63) into (59) gives and thus Therefore, by induction, it follows that (53) holds for all k, and we complete the proof.
By this lemma, we are able to establish the convergence rate result of Algorithm 1 for solving (1) if the problem is convex.
Theorem 6 (Convergence rate for the nonsmooth convex case) Under Assumptions 1 through 4, let {x k } k≥1 be the sequence generated from Algorithm 1 with stepsize satisfying (52) and also where γ ρ,1 and γ ρ,2 are defined in (51). We have 1. If the function F is convex, then where 2. If F is strongly convex with constant µ, then Before proving this theorem, we make two remarks and present a few lemmas below.
Remark 6 Similar to (37), for the linear convergence result (66), the strong convexity assumption can be weakened to optimal strong convexity. The latter one is strictly weaker than the former one; see [17] for more discussions.

Remark 7 (Comparison of stepsize) For the special case that the delay is bounded by
, we have both γ ρ,1 and γ ρ,2 are O(τ ). Thus we can take stepsize almost 1 Lc , which is larger than the stepsize 1 2Lc given in [17].
Lemma 8 Let γ ρ,2 be defined in (51). We have Proof It is proved via the Cauchy-Schwarz inequality, the bound (63), and Lemma 9 It holds that Lemma 10 Let γ ρ,2 be defined in (51). It holds that Proof Since i k is uniformly distributed and independent of j k , we have We split the term and apply the convexity of f and Lipschitz continuity of ∇f to get Substituting (71) into (70) and taking expectation yield (43) and (53) and using the definition of γ ρ,2 , we complete the proof of (69).

Lemma 11 Under the assumptions of Theorem 6, we have
Proof Taking expectation on both side of (41) and using (67) yield from the condition on η in (64).
Now we are ready to prove Theorem 6.
Proof (of Theorem 6) From the update of x k+1 , we have , and thus for any x i k , it holds from the convexity of r i k that Since From the definition of P X * , it follows that Then using (72) and (73), we have We split the cross term to have From (7), it follows that Plugging the above two equations into (74) gives Substituting (67) through (69) into (75) and rearranging terms yield The above inequality together with (64) implies and thus, with the monotonicity of E[F (x k )] in Lemma 11, Hence, (65) follows. When F is strongly convex with constant µ, we have and thus from (76), it follows that Therefore, (66) follows, and we complete the proof.

Poisson distribution
We can treat the asynchronous reading and writing as a queueing system. Assume the p + 1 processors have the same computing power (i.e., the same speed of reading and writing). At any time k, suppose the update to x i k is performed by the p k -th processor, which can be treated as the server with speed (or service rate) one of reading and writing. All the other p processors can be treated as customers, each with speed (or arrival rate) one, where any update to x from the p processors can be regarded as one customer's arrival. Under this setting, from the p k -th processor starts reading x until it finishes updating x i k , there would be p customers in the queue in average, namely, the delay j k follows the Poisson distribution with parameter p. Summarizing the above discussion, we have the following result.
Claim Suppose Algorithm 1 runs on a system with p + 1 processors, which have the same speed of reading and writing during the iterations. Then the delay j k follows the Poisson distribution with parameter p, i.e., for all k, which implies no delay if p = 0.
In general, if the processors have different computing power, j k would follow Poisson distribution with a parameter being the speed ratio of the other p processors to the p k -th one. However, in a multi-core workstation with shared memory, the processors are usually of the same style and can have the same computing ability. In the following, we assume the distribution in (8) to be Poisson distribution with parameter p and discuss the convergence results we obtained in the previous sections. First we give the values of the expected quantities we used before.
The proof of this proposition is standard. From the quantities in (79) and the theorems we established in the previous sections, we make the following observations: , we can guarantee the convergence of Algorithm 1 for both smooth and nonsmooth problems by setting η

Numerical experiments
In this section, we evaluate the numerical performance of Algorithm 1 on solving two problems: the LASSO problem and the nonnegative matrix factorization (NMF). The tests were carried out on a machine with 64GB of memory and two Intel Xeon E5-2690 v2 processors (20 cores, 40 threads). All of the experiments were coded in C++ and its threading library was used for parallelization. We use the Eigen library for numerical linear algebra operations. To measure the delay, we use an atomic variable to track the number of iterations as defined in the paper. The atomic variable will be incremented by one for each update. For each thread, the delay is calculated based on the difference of the iteration counters before and after the update. For LASSO, two different settings were used. The first one sets the stepsize by the expected delay according to the analysis of this paper, and the other one used the maximum delay from [16,17] and is dubbed as AsySCD. We compared the async-BCU to the serial BCU, which can be regarded as a special case of Algorithm 1 with the delay j k ≡ 0, ∀k. For NMF, we set the stepszie by the expected delay and test its convergence behavior with different numbers of threads.

Parameter settings
According to Theorem 5, the following two stepsizes were used: 1 This paper : η = 1/Lc Max delay : η = 1/Lc where τ equals the maximum number of the generated sequence of delays.

LASSO
We measure the performance of Algorithm 1 on the LASSO problem [30] minimize x∈R n where A ∈ R N ×n , b ∈ R N , and λ is a parameter balancing the fitting term and the regularization term. We randomly generated A and b following the standard normal distribution. The size was fixed to n = 2N and N = 10, 000, and λ = 1 N was used. The Lipschitz constant L c = max{ (A T i A i ) 2 , ∀i}, where A i represents the ith column block of A. Figure 1 shows the delay distribution of Algorithm 1 with different numbers of threads. The blue bars are the normalized histogram so that the bar heights add to 1. Orange curve is the probability density function of Poisson distribution. By using 5 and 10 threads, we observe that the number of delays is concentrated on 4, and 9 respectively. When the number of threads is relatively large, the actual delay distribution closely matches with the theoretical distribution as we discussed in Section 5. For 20 threads, an interesting observation is that, the actual probability density is higher than the theoretical probability density when the number of delays is around 9. We think this is due to the architecture of the testing environment, i.e., the average delay within a CPU is smaller than the average delay across two different CPUs. We observe a similar behavior when 40 threads are used. Figure 2 plots the convergence behavior of Algorithm 1 running on 40 threads with different block sizes. We partition x into m equal-sized blocks with block sizes varying among {10, 50, 100, 500}. The results of the serial randomized coordinate descent method is also plotted for comparison. Here, one epoch is equivalent to updating all coordinates once. Comparing to the serial method, we observe that the delay does affect the convergence speed, and the affect becomes weaker as m increases. Hence, Algorithm 1 can have nearly linear speed-up when the number of blocks is large. In addition, we note that the stepsize setting of AsySCD is too conservative, and Algorithm 1 with stepsize set by the expected delay converges significantly faster. However, we observed that, in general, we could not take larger stepsize than that in (80a). Some divergence behaviors are observed when using stepsizes larger than that in (80a).

Nonnegative matrix factorization (NMF)
This section presents the numerical results of applying Algorithm 1 for solving the NMF problem [23] minimize X,Y where Z ∈ R M ×N + is a given nonnegative matrix. We generated Z = Z L Z R with the elements of Z L and Z R first drawn from the standard normal distribution and then projected into the nonnegative orthant. The size was fixed to M = N = 10, 000 and m = 100.  (80), and also the serial randomized coordinate descent method. The tested problem has 10, 000 samples and 20, 000 coordinates that are evenly partitioned into m blocks. It was simulated as running with 40 threads. We run 100 epochs for each experiments.
We treated one column of X or Y as one block coordinate, and during the iterations, every column of X was kept with unit norm. Therefore, the partial gradient Lipschitz constant equals one if one column of Y is selected to update and y k i k 2 2 if the i k -th column of X is selected. Since y k i k 2 2 could approach to zero, we set the Lipschitz constant to max(0.001, y k i k 2 2 ). This modification can guarantee the whole sequence convergence of the coordinate descent method [36]. Due to nonconvexity, global optimality cannot be guaranteed. Thus, we set the starting point close to Z L and Z R . Specifically, we let X 0 = Z L + 0.5Ξ L and Y 0 = Z R + 0.5Ξ R with the elements of Ξ L and Ξ R following the standard normal distribution. All methods used the same starting point. Figure 3 shows the delay distribution behavior of Algorithm 1 for solving NMF. The observation is similar to Figure 1. Figure 4 plots the convergence results of Algorithm 1 running with 1, 5, 10, 20 and 40 threads. From the results, we observe that Algorithm 1 scales up to 10 threads for the tested problem. Degenerated convergence is observed with 20 and 40 threads. This is mostly due to the following three reasons: (1) since the number of blocks is relatively small (m = 200), as shown in (80a), using more threads leads to smaller stepsize, hence, slower convergence; (2) the gradient used for the current update is more staled when a relative large number of threads are used, which also leads to slow convergence; (3) high cache miss rates and false sharing also downgrade the speedup performance.

Conclusions
We have analyzed the convergence of the async-BCU method for solving both convex and nonconvex problems in a probabilistic way. We showed that the  algorithm is guaranteed to converge for smooth problems if the expected delay is finite and for nonsmooth problems if the variance of the delay is also finite. In addition, we established sublinear convergence of the method for weakly convex problems and linear convergence for strongly convex ones. The stepsize we obtained depends on certain expected quantities. Assuming the given p + 1 processors perform identically, we showed that the delay follows a Poisson distribution with parameter p and thus fully determined the stepsize. We have simulated the performance of the algorithm with our determined stepsize on solving LASSO and the nonnegative matrix factorization, and the numerical results validated our analysis. A.1 Proof of Lemma 3 Proof Following the proof of Theorem 1 in [16], we have and We first show the first inequality in (26). Note that (25) gives us When t = 0, we have from (86) that ∇f (x 0 ) 2 −E ∇f (x 1 ) 2 ≤ 2ηLr √ m ∇f (x 0 ) 2 ≤ (1 + Mρ) ηLr √ m ∇f (x 0 ) 2 . Hence, ∇f (x 0 ) 2 ≤ ρE ∇f (x 1 ) 2 from (88). Now we assume that