Parallel coordinate descent methods for big data optimization

In this work we show that randomized (block) coordinate descent methods can be accelerated by parallelization when applied to the problem of minimizing the sum of a partially separable smooth convex function and a simple separable convex function. The theoretical speedup, as compared to the serial method, and referring to the number of iterations needed to approximately solve the problem with high probability, is a simple expression depending on the number of parallel processors and a natural and easily computable measure of separability of the smooth component of the objective function. In the worst case, when no degree of separability is present, there may be no speedup; in the best case, when the problem is separable, the speedup is equal to the number of processors. Our analysis also works in the mode when the number of blocks being updated at each iteration is random, which allows for modeling situations with busy or unreliable processors. We show that our algorithm is able to solve a LASSO problem involving a matrix with 20 billion nonzeros in 2 h on a large memory node with 24 cores.


Introduction
Big data optimization.Recently there has been a surge in interest in the design of algorithms suitable for solving convex optimization problems with a huge number of variables [16,11].Indeed, the size of problems arising in fields such as machine learning [1], network analysis [29], PDEs [27], truss topology design [15] and compressed sensing [5] usually grows with our capacity to solve them, and is projected to grow dramatically in the next decade.In fact, much of computational science Our starting point was the following simple observation.Assume that we wish to minimize a separable function F of n variables (i.e., a function that can be written as a sum of n functions each of which depends on a single variable only).For simplicity, in this thought experiment, assume that there are no constraints.Clearly, the problem of minimizing F can be trivially decomposed into n independent univariate problems.Now, if we have n processors/threads/cores, each assigned with the task of solving one of these problems, the number of parallel iterations should not depend on the dimension of the problem 1 .In other words, we get an n-times speedup compared to the situation with a single processor only.Note that any parallel algorithm of this type can be viewed as a parallel coordinate descent method.Hence, a PCDM with n processors should be n-times faster than a serial one.If τ processors are used instead, where 1 ≤ τ ≤ n, one would expect a τ -times speedup.
By extension, one would perhaps expect that optimization problems with objective functions which are "close to being separable" would also be amenable to acceleration by parallelization, where the acceleration factor τ would be reduced with the reduction of the "degree of separability".One of the main messages of this paper is an affirmative answer to this.Moreover, we give explicit and simple formulae for the speedup factors.
As it turns out, and as we discuss later in this section, many real-world big data optimization problems are, quite naturally, "close to being separable".We believe that this means that PCDMs is a very promising class of algorithms when it comes to solving structured big data optimization problems.
Minimizing a partially separable composite objective.In this paper we study the problem minimize F (x) where f is a (block) partially separable smooth convex function and Ω is a simple (block) separable convex function.We allow Ω to have values in R ∪ {∞}, and for regularization purposes we assume Ω is proper and closed.While (1) is seemingly an unconstrained problem, Ω can be chosen to model simple convex constraints on individual blocks of variables.Alternatively, this function can be used to enforce a certain structure (e.g., sparsity) in the solution.For a more detailed account we refer the reader to [16].Further, we assume that this problem has a minimum (F * > −∞).What we mean by "smoothness" and "simplicity" will be made precise in the next section.
Let us now describe the key concept of partial separability.Let x ∈ R N be decomposed into n non-overlapping blocks of variables x (1) , . . ., x (n) (this will be made precise in Section 2).We assume throughout the paper that f : R N → R is partially separable of degree ω, i.e., that it can be written in the form where J is a finite collection of nonempty subsets of [n] def = {1, 2, . . ., n} (possibly containing identical sets multiple times), f J are differentiable convex functions such that f J depends on blocks x (i) for i ∈ J only, and |J| ≤ ω for all J ∈ J .
Clearly, 1 ≤ ω ≤ n.The PCDM algorithms we develop and analyze in this paper only need to know ω, they do not need to know the decomposition of f giving rise to this ω.
Examples of partially separable functions.Many objective functions naturally encountered in the big data setting are partially separable.Here we give examples of three loss/objective functions frequently used in the machine learning literature and also elsewhere.For simplicity, we assume all blocks are of size 1 (i.e., N = n).Let where m is the number of examples, x ∈ R n is the vector of features, (A j , y j ) ∈ R n × R are labeled examples and L is one of the three loss functions listed in Table 1.Let A ∈ R m×n with row j equal Square Loss 1 2 (A T j x − y j ) 2 Logistic Loss log(1 + e −y j A T j x ) Hinge Square Loss 1 2 max{0, 1 − y j A T j x} 2 Table 1: Three examples of loss of functions to A T j .Often, each example depends on a few features only; the maximum over all features is the degree of partial separability ω.More formally, note that the j-th function in the sum (4) in all cases depends on A j 0 coordinates of x (the number of nonzeros in the j-th row of A) and hence f is partially separable of degree ω = max j A j 0 .
All three functions of Table 1 are smooth (based on the definition of smoothness in the next section).We refer the reader to [13] for more examples of interesting (but nonsmooth) partially separable functions arising in graph cuts and matrix completion.
Brief literature review.Several papers were written recently studying the iteration complexity of serial CDMs of various flavours and in various settings.We will only provide a brief summary here, for a more detailed account we refer the reader to [16].Classical CDMs update the coordinates in a cyclic order; the first attempt at analyzing the complexity of such a method is due to [21].Stochastic/randomized CDMs, that is, methods where the coordinate to be updated is chosen randomly, were first analyzed for quadratic objectives [24,4], later independently generalized to L 1 -regularized problems [23] and smooth block-structured problems [10], and finally unified and refined in [19,16].The problems considered in the above papers are either unconstrained or have (block) separable constraints.Recently, randomized CDMs were developed for problems with linearly coupled constraints [7,8].
A greedy CDM for L 1 -regularized problems was first analyzed in [15]; more work on this topic include [5,2].A CDM with inexact updates was first proposed and analyzed in [26].Partially separable problems were independently studied in [13], where an asynchronous parallel stochastic gradient algorithm was developed to solve them.
When writing this paper, the authors were aware only of the parallel CDM proposed and analyzed in [1].Several papers on the topic appeared around the time this paper was finalized or after [6,28,22,22,14].Further papers on various aspects of the topic of parallel CDMs, building on the work in this paper, include [25,17,3,18].
Contents.We start in Section 2 by describing the block structure of the problem, establishing notation and detailing assumptions.Subsequently we propose and comment in detail on two parallel coordinate descent methods.In Section 3 we summarize the main contributions of this paper.In Section 4 we deal with issues related to the selection of the blocks to be updated in each iteration.It will involve the development of some elementary random set theory.Sections 5-6 deal with issues related to the computation of the update to the selected blocks and develop a theory of Expected Separable Overapproximation (ESO), which is a novel tool we propose for the analysis of our algorithms.In Section 7 we analyze the iteration complexity of our methods and finally, Section 8 reports on promising computational results.For instance, we conduct an experiment with a big data (cca 350GB) LASSO problem with a billion variables.We are able to solve the problem using one of our methods on a large memory machine with 24 cores in 2 hours, pushing the difference between the objective value at the starting iterate and the optimal point from 10 22 down to 10 −14 .We also conduct experiments on real data problems coming from machine learning.

Parallel Block Coordinate Descent Methods
In Section 2.1 we formalize the block structure of the problem, establish notation2 that will be used in the rest of the paper and list assumptions.In Section 2.2 we propose two parallel block coordinate descent methods and comment in some detail on the steps.

Block structure, notation and assumptions
Some elements of the setup described in this section was initially used in the analysis of block coordinate descent methods by Nesterov [10] (e.g., block structure, weighted norms and block Lipschitz constants).
The block structure of ( 1) is given by a decomposition of R N into n subspaces as follows.Let U ∈ R N ×N be a column permutation3 of the N × N identity matrix and further let U = [U 1 , U 2 , . . ., U n ] be a decomposition of U into n submatrices, with U i being of size N × N i , where i N i = N .Proposition 1 (Block decomposition4 ).Any vector x ∈ R N can be written uniquely as where for every j we get In view of the above proposition, from now on we write x (i) def = U T i x ∈ R N i , and refer to x (i) as the i-th block of x.The definition of partial separability in the introduction is with respect to these blocks.For simplicity, we will sometimes write x = (x (1) , . . ., x (n) ).
Projection onto a set of blocks.For S ⊂ [n] and x ∈ R N we write That is, given x ∈ R N , x [S] is the vector in R N whose blocks i ∈ S are identical to those of x, but whose other blocks are zeroed out.In view of Proposition 1, we can equivalently define x [S] block-by-block as follows Inner products.The standard Euclidean inner product in spaces R N and R N i , i ∈ [n], will be denoted by •, • .Letting x, y ∈ R N , the relationship between these inner products is given by x, y = n j=1 x (i) , y (i) .
For any w ∈ R n and x, y ∈ R N we further define For are equipped with a pair of conjugate norms: where B i is an N i ×N i positive definite matrix and t * ++ , define a pair of conjugate norms in R N by Note that these norms are induced by the inner product (9) and the matrices B 1 , . . ., B n .Often we where the constants L i are defined below.
Smoothness of f .We assume throughout the paper that the gradient of f is block Lipschitz, uniformly in x, with positive constants L 1 , . . ., L n , i.e., that for all where (11) is the following standard inequality [9]: Separability of Ω.We assume that5 Ω : R N → R ∪ {+∞} is (block) separable, i.e., that it can be decomposed as follows: where the functions Ω i : R N i → R ∪ {+∞} are convex and closed.
Strong convexity.In one of our two complexity results (Theorem 20) we will assume that either f or Ω (or both) is strongly convex.A function φ : R N → R ∪ {+∞} is strongly convex with respect to the norm • w with convexity parameter µ φ (w) ≥ 0 if for all x, y ∈ dom φ, where φ (x) is any subgradient of φ at x.The case with µ φ (w) = 0 reduces to convexity.Strong convexity of F may come from f or Ω (or both); we write µ f (w) (resp.µ Ω (w)) for the (strong) convexity parameter of f (resp.Ω).It follows from ( 14) that The following characterization of strong convexity will be useful.For all x, y ∈ dom φ and It can be shown using ( 12) and ( 14) that µ f (w) ≤ L i w i .

Algorithms
In this paper we develop and study two generic parallel coordinate descent methods.The main method is PCDM1; PCDM2 is its "regularized" version which explicitly enforces monotonicity.As we will see, both of these methods come in many variations, depending on how Step 3 is performed.
Algorithm 1 Parallel Coordinate Descent Method 1 (PCDM1) Randomly generate a set of blocks S k ⊆ {1, 2, . . ., n} 4: Let us comment on the individual steps of the two methods.
Step 3. At the beginning of iteration k we pick a random set (S k ) of blocks to be updated (in parallel) during that iteration.The set S k is a realization of a random set-valued mapping Ŝ with values in 2 [n] or, more precisely, it the sets S k are iid random sets with the distribution of Ŝ.For brevity, in this paper we refer to such a mapping by the name sampling.We limit our attention to Algorithm 2 Parallel Coordinate Descent Method 2 (PCDM2) Randomly generate a set of blocks S k ⊆ {1, 2, . . ., n} 4: If F (x k+1 ) > F (x k ), then x k+1 ← x k 6: end for uniform samplings, i.e., random sets having the following property: That is, the probability that a block gets selected is the same for all blocks.Although we give an iteration complexity result covering all such samplings (provided that each block has a chance to be updated, i.e., P(i ∈ Ŝ) > 0), there are interesting subclasses of uniform samplings (such as doubly uniform and nonoverlapping uniform samplings; see Section 4) for which we give better results.
Step 4. For x ∈ R N we define 6 . where and β > 0, w = (w 1 , . . ., w n ) T ∈ R n ++ are parameters of the method that we will comment on later.Note that in view of ( 5), ( 10) and ( 13), H β,w (x, •) is block separable; Consequently, we have h(x) = (h (1) where We mentioned in the introduction that besides (block) separability, we require Ω to be "simple".By this we mean that the above optimization problem leading to h (i) (x) is "simple" (e.g., it has a closed-form solution).Recall from (8) that (h(x k )) [S k ] is the vector in R N identical to h(x k ) except for blocks i / ∈ S k , which are zeroed out.Hence, Step 4 of both methods can be written as follows: In parallel for i ∈ S k do: Parameters β and w depend on f and Ŝ and stay constant throughout the algorithm.We are not ready yet to explain why the update is computed via (17) and ( 18) because we need technical tools, which will be developed in Section 4, to do so.Here it suffices to say that the parameters β and w come from a separable quadratic overapproximation of E[f (x + h [ Ŝ] )], viewed as a function of h ∈ R N .Since expectation is involved, we refer to this by the name Expected Separable Overapproximation (ESO).This novel concept, developed in this paper, is one of the main tools of our complexity analysis.Section 5 motivates and formalizes the concept, answers the why question, and develops some basic ESO theory.
Section 6 is devoted to the computation of β and w for partially separable f and various special classes of uniform samplings Ŝ. Typically we will have w i = L i , while β will depend on easily computable properties of f and Ŝ.For example, if Ŝ is chosen as a subset of [n] of cardinality τ , with each subset chosen with the same probability (we say that Ŝ is τ -nice) then, assuming n > 1, we may choose w = L and β = 1 + (ω−1)(τ −1) , where ω is the degree of partial separability of f .More generally, if Ŝ is any uniform sampling with the property | Ŝ| = τ with probability 1, then we may choose w = L and β = min{ω, τ }.Note that in both cases w = L and that the latter β is always larger than (or equal to) the former one.This means, as we will see in Section 7, that we can give better complexity results for the former, more specialized, sampling.We analyze several more options for Ŝ than the two just described, and compute parameters β and w that should be used with them (for a summary, see Table 4).
Step 5.The reason why, besides PCDM1, we also consider PCDM2, is the following: in some situations we are not able to analyze the iteration complexity of PCDM1 (non-strongly-convex F where monotonicity of the method is not guaranteed by other means than by directly enforcing it by inclusion of Step 5).Let us remark that this issue arises for general Ω only.It does not exist for Ω = 0, Ω(•) = λ • 1 and for Ω encoding simple constraints on individual blocks; in these cases one does not need to consider PCDM2.Even in the case of general Ω we sometimes get monotonicity for free, in which case there is no need to enforce it.Let us stress, however, that we do not recommend implementing PCDM2 as this would introduce too much overhead; in our experience PCDM1 works well even in cases when we can only analyze PCDM2.

Smmary of Contributions
In this section we summarize the main contributions of this paper (not in order of significance).
1. Problem generality.We give the first complexity analysis for parallel coordinate descent methods for problem (1) in its full generality.

2.
Complexity.We show theoretically (Section 7) and numerically (Section 8) that PCDM accelerates on its serial counterpart for partially separable problems.In particular, we establish two complexity theorems giving lower bounds on the number of iterations k sufficient for one or both of the PCDM variants (for details, see the precise statements in Section 7) to produce a random iterate x k for which the problem is approximately solved with high probability, i.e., P(F The results, summarized in Table 2, hold under the standard assumptions listed in Section 2.1 and the additional assumption that f, Ŝ, β and w satisfy the following inequality for all x, h ∈ R N : This inequality, which we call Expected Separable Overapproximation (ESO), is the main new theoretical tool that we develop in this paper for the analysis of our methods (Sections 4-6 are devoted to the development of this theory).

Setting Complexity Theorem
Convex f O βn The main observation here is that as the average number of block updates per iteration increases (say, τ = E[| Ŝ|]), enabled by the utilization of more processors, the leading term in the complexity estimate, n/τ , decreases in proportion.However, β will generally grow with τ , which has an adverse effect on the speedup.Much of the theory in this paper goes towards producing formulas for β (and w), for partially separable f and various classes of uniform samplings Ŝ. Naturally, the ideal situation is when β does not grow with τ at all, or if it only grows very slowly.We show that this is the case for partially separable functions f with small ω.For instance, in the extreme case when f is separable (ω = 1), we have β = 1 and we obtain linear speedup in τ .As ω increases, so does β, depending on the law governing Ŝ. Formulas for β and ω for various samplings Ŝ are summarized in Table 4. 3. Algorithm unification.Depending on the choice of the block structure (as implied by the choice of n and the matrices U 1 , . . ., U n ) and the way blocks are selected at every iteration (as given by the choice of Ŝ), our framework encodes a family of known and new algorithms7 (see Table 3).In particular, PCDM is the first method which "continuously" interpolates between serial coordinate descent and gradient (by manipulating n and/or E[| Ŝ|]).

Method
4. Partial separability.We give the first analysis of a coordinate descent type method dealing with a partially separable loss / objective.In order to run the method, we need to know the Lipschitz constants L i and the degree of partial separability ω.It is crucial that these quantities are often easily computable/predictable in the huge-scale setting.For example, if and we choose all blocks to be of size 1, then L i is equal to the squared Euclidean norm of the i-th column of A and ω is equal to the maximum number of nonzeros in a row of A. Many problems in the big data setting have small ω compared to n.

Choice of blocks.
To the best of our knowledge, existing randomized strategies for paralleling gradient-type methods (e.g., [1]) assume that Ŝ (or an equivalent thereof, based on the method) is chosen as a subset of [n] of a fixed cardinality, uniformly at random.We refer to such Ŝ by the name nice sampling in this paper.We relax this assumption and our treatment is hence much more general.In fact, we allow for Ŝ to be any uniform sampling.It is possible to further consider nonuniform samplings9 , but this is beyond the scope of this paper.
In particular, as a special case, our method allows for a variable number of blocks to be updated throughout the iterations (this is achieved by the introduction of doubly uniform samplings).This may be useful in some settings such as when the problem is being solved in parallel by τ unreliable processors each of which computes its update h (i) (x k ) with probability p b and is busy/down with probability 1 − p b (binomial sampling).
Uniform, doubly uniform, nice, binomial and other samplings are defined, and their properties studied, in Section 4.
6. ESO and formulas for β and w.In Table 4 we list parameters β and w for which ESO inequality (19) holds.Each row corresponds to a specific sampling Ŝ (see Section 4 for the definitions).The last 5 samplings are special cases of one or more of the first three samplings.
Details such as what is ν, γ and "monotonic" ESO are explained in appropriate sections later in the text.When a specific sampling Ŝ is used in the algorithm to select blocks in each iteration, the corresponding parameters β and w are to be used in the method for the computation of the update (see (17) and ( 18)). sampling doubly uniform  En route to proving the iteration complexity results for our algorithms, we develop a theory of deterministic and expected separable overapproximation (Sections 5 and 6) which we believe is of independent interest, too.For instance, methods based on ESO can be compared favorably to the Diagonal Quadratic Approximation (DQA) approach used in the decomposition of stochastic optimization programs [20].
7. Parallelization speedup.Our complexity results can be used to derive theoretical parallelization speedup factors.For several variants of our method, in case of a non-strongly convex objective, these are given in Section 7.1 (Table 5).For instance, in the case when all block are updated at each iteration (we later refer to Ŝ having this property by the name fully parallel sampling), the speedup factor is equal to n ω .If the problem is separable (ω = 1), the speedup is equal to n; if the problem is not separable (ω = n), there may be no speedup.For strongly convex F the situation is even better; the details are given in Section 7.2.

8.
Relationship to existing results.To the best of our knowledge, there are just two papers analyzing a parallel coordinate descent algorithm for convex optimization problems [1,6].In the first paper all blocks are of size 1, Ŝ corresponds to what we call in this paper a τ -nice sampling (i.e., all sets of τ coordinates are updated at each iteration with equal probability) and hence their algorithm is somewhat comparable to one of the many variants of our general method.While the analysis in [1] works for a restricted range of values of τ , our results hold for all τ ∈ [n].Moreover, the authors consider a more restricted class of functions f and the special case Ω = λ x 1 , which is simpler to analyze.Lastly, the theoretical speedups obtained in [1], when compared to the serial CDM method, depend on a quantity σ that is hard to compute in big data settings (it involves the computation of an eigenvalue of a huge-scale matrix).Our speedups are expressed in terms of natural and easily computable quantity: the degree ω of partial separability of f .In the setting considered by [1], in which more structure is available, it turns out that ω is an upper bound10 on σ.Hence, we show that one can develop the theory in a more general setting, and that it is not necessary to compute σ (which may be complicated in the big data setting).The parallel CDM method of the second paper [6] only allows all blocks to be updated at each iteration.Unfortunately, the analysis (and the method) is too coarse as it does not offer any theoretical speedup when compared to its serial counterpart.In the special case when only a single block is updated in each iteration, uniformly at random, our theoretical results specialize to those established in [16].

9.
Computations.We demonstrate that our method is able to solve a LASSO problem involving a matrix with a billion columns and 2 billion rows on a large memory node with 24 cores in 2 hours (Section 8), achieving a 20× speedup compared to the serial variant and pushing the residual by more than 30 degrees of magnitude.While this is done on an artificial problem under ideal conditions (controlling for small ω), large speedups are possible in real data with ω small relative to n.We also perform additional experiments on real data sets from machine learning (e.g., training linear SVMs) to illustrate that the predictions of our theory match reality.

10.
Code.The open source code with an efficient implementation of the algorithm(s) developed in this paper is published here: http://code.google.com/p/ac-dc/.

Block Samplings
In Step 3 of both PCDM1 and PCDM2 we choose a random set of blocks S k to be updated at the current iteration.Formally, S k is a realization of a random set-valued mapping Ŝ with values in 2 [n] , the collection of subsets of [n].For brevity, in this paper we refer to Ŝ by the name sampling.
A sampling Ŝ is uniquely characterized by the probability mass function that is, by assigning probabilities to all subsets of [n].Further, we let p = (p 1 , . . ., p n ) T , where In Section 4.1 we describe those samplings for which we analyze our methods and in Section 4.2 we prove several technical results, which will be useful in the rest of the paper.

Uniform, Doubly Uniform and Nonoverlapping Uniform Samplings
A sampling is proper if p i > 0 for all blocks i.That is, from the perspective of PCDM, under a proper sampling each block gets updated with a positive probability at each iteration.Clearly, PCDM can not converge for a sampling that is not proper.A sampling Ŝ is uniform if all blocks get updated with the same probability, i.e., if p i = p j for all i, j.We show in (33) that, necessarily, Note that a uniform sampling is proper if and only if it is not nil.
All our iteration complexity results in this paper are for PCDM used with a proper uniform sampling (see Theorems 19 and 20) for which we can compute β and w giving rise to an inequality (we we call "expected separable overapproximation") of the form (43).We derive such inequalities for all proper uniform samplings (Theorem 12) as well as refined results for two special subclasses thereof: doubly uniform samplings (Theorem 15) and nonoverlapping uniform samplings (Theorem 13).We will now give the definitions: 1. Doubly Uniform (DU) samplings.A DU sampling is one which generates all sets of equal cardinality with equal probability.That is, P(S ) = P(S ) whenever |S | = |S |.The name comes from the fact that this definition postulates a different uniformity property, "standard" uniformity is a consequence.Indeed, let us show that a DU sampling is necessarily uniform.Let q j = P(| Ŝ| = j) for j = 0, 1, . . ., n and note that from the definition we know that whenever S is of cardinality j, we have P(S) = q j / n j .Finally, using this we obtain It is clear that each DU sampling is uniquely characterized by the vector of probabilities q; its density function is given by 2. Nonoverlapping Uniform (NU) samplings.A NU sampling is one which is uniform and which assigns positive probabilities only to sets forming a partition of [n].Let S 1 , S 2 , . . ., S l be a partition of [n], with |S j | > 0 for all j.The density function of a NU sampling corresponding to this partition is given by Note that Let us now describe several interesting special cases of DU and NU samplings: 3. Nice sampling.Fix 1 ≤ τ ≤ n.A τ -nice sampling is a DU sampling with q τ = 1.
Interpretation: There are τ processors/threads/cores available.At the beginning of each iteration we choose a set of blocks using a τ -nice sampling (i.e., each subset of τ blocks is chosen with the same probability), and assign each block to a dedicated processor/thread/core.Processor assigned with block i would compute and apply the update h (i) (x k ).This is the sampling we use in our computational experiments.
4. Independent sampling.Fix 1 ≤ τ ≤ n.A τ -independent sampling is a DU sampling with where Interpretation: There are τ processors/threads/cores available.Each processor chooses one of the n blocks, uniformly at random and independently of the other processors.It turns out that the set Ŝ of blocks selected this way is DU with q as given above.Since in one parallel iteration of our methods each block in Ŝ is updated exactly once, this means that if two or more processors pick the same block, all but one will be idle.On the other hand, this sampling can be generated extremely easily and in parallel!For τ n this sampling is a good (and fast) approximation of the τ -nice sampling.For instance, for n = 10 3 and τ = 8 we have q 8 = 0.9723, q 7 = 0.0274, q 6 = 0.0003 and q k ≈ 0 for k = 1, . . ., 5.
5. Binomial sampling.Fix 1 ≤ τ ≤ n and 0 < p b ≤ 1.A (τ, p b )-binomial sampling is defined as a DU sampling with Notice that Interpretation: Consider the following situation with independent equally unreliable processors.
We have τ processors, each of which is at any given moment available with probability p b and busy with probability 1 − p b , independently of the availability of the other processors.Hence, the number of available processors (and hence blocks that can be updated in parallel) at each iteration is a binomial random variable with parameters τ and p b .That is, the number of available processors is equal to k with probability q k .
-Case 1 (explicit selection of blocks): We learn that k processors are available at the beginning of each iteration.Subsequently, we choose k blocks using a k-nice sampling and "assign one block" to each of the k available processors.
-Case 2 (implicit selection of blocks): We choose τ blocks using a τ -nice sampling and assign one to each of the τ processors (we do not know which will be available at the beginning of the iteration).With probability q k , k of these will send their updates.It is easy to check that the resulting effective sampling of blocks is (τ, p b )-binomial.
6. Serial sampling.This is a DU sampling with q 1 = 1.Also, this is a NU sampling with l = n and S j = {j} for j = 1, 2, . . ., l.That is, at each iteration we update a single block, uniformly at random.This was studied in [16].
7. Fully parallel sampling.This is a DU sampling with q n = 1.Also, this is a NU sampling with l = 1 and That is, at each iteration we update all blocks.
The following simple result says that the intersection between the class of DU and NU samplings is very thin.A sampling is called vacuous if P(∅) > 0.
Proposition 2. There are precisely two nonvacuous samplings which are both DU and NU: i) the serial sampling and ii) the fully parallel sampling.
Proof.Assume Ŝ is nonvacuous, NU and DU.Since Ŝ is nonvacuous, P( Ŝ = ∅) = 0. Let S ⊂ [n] be any set for which P( Ŝ = S) > 0. If 1 < |S| < n, then there exists S = S of the same cardinality as S having a nonempty intersection with S. Since Ŝ is doubly uniform, we must have P( Ŝ = S ) = P( Ŝ = S ) > 0. However, this contradicts the fact that Ŝ is non-overlapping.Hence, Ŝ can only generate sets of cardinalities 1 or n with positive probability, but not both.One option leads to the fully parallel sampling, the other one leads to the serial sampling.

Technical results
For a given sampling Ŝ and i, j ∈ [n] we let The following simple result has several consequences which will be used throughout the paper.
Lemma 3 (Sum over a random index set).Let ∅ = J ⊂ [n] and Ŝ be any sampling. and 11 Sum over an empty index set will, for convenience, be defined to be zero.
Proof.We prove the first statement, proof of the remaining statements is essentially identical: The consequences are summarized in the next theorem and the discussion that follows.
Theorem 4. Let ∅ = J ⊂ [n] and Ŝ be an arbitrary sampling.Further, let a, h ∈ R N , w ∈ R n + and let g be a block separable function, i.e., g(x) = i g i (x (i) ).Then Moreover, the matrix all five identities follow directly by applying Lemma 3. Finally, for any θ = (θ 1 , . . ., θ n ) T ∈ R n , The above results hold for arbitrary samplings.Let us specialize them, in order of decreasing generality, to uniform, doubly uniform and nice samplings.
• Uniform samplings.If Ŝ is uniform, then from (28) using J = [n] we get Plugging (33) into ( 28), (30), ( 31) and (32) yields • Doubly uniform samplings.Consider the case n > 1; the case n = 1 is trivial.For doubly uniform Ŝ, p ij is constant for i = j: Indeed, this follows from Substituting (38) and ( 33) into (29) then gives Then for all i ∈ J, Substituting this into (26) yields 5 Expected Separable Overapproximation Recall that given x k , in PCDM1 the next iterate is the random vector for a particular choice of h ∈ R N .Further recall that in PCDM2, x k , otherwise, again for a particular choice of h.While in Section 2 we mentioned how h is computed, i.e., that h is the minimizer of H β,w (x, •) (see (17) and ( 18)), we did not explain why is h computed this way.
The reason for this is that the tools needed for this were not yet developed at that point (as we will see, some results from Section 4 are needed).In this section we give an answer to this why question.Given x k ∈ R N , after one step of PCDM1 performed with update h we get On the the other hand, after one step of PCDM2 we have So, for both PCDM1 and PCDM2 the following estimate holds, A good choice for h to be used in the algorithms would be one minimizing the right hand side of inequality (42).At the same time, we would like the minimization process to be decomposable so that the updates h (i) , i ∈ Ŝ, could be computed in parallel.However, the problem of finding such h is intractable in general even if we do not require parallelizability.Instead, we propose to construct/compute a "simple" separable overapproximation of the right-hand side of (42).Since the overapproximation will be separable, parallelizability is guaranteed; "simplicity" means that the updates h (i) can be computed easily (e.g., in closed form).
From now on we replace, for simplicity and w.l.o.g., the random vector x k by a fixed deterministic vector x ∈ R N .We can thus remove conditioning in (42) and instead study the quantity we indeed find a simple separable overapproximation of where we recall from ( 18) that H β,w (x, h) = f (x) + ∇f (x), h + β 2 h 2 w + Ω(x + h).That is, (44) says that the expected objective value after one parallel step of our methods, if block i ∈ Ŝ is updated by h (i) , is bounded above by a convex combination of F (x) and H β,w (x, h).The natural choice of h is to set Note that this is precisely the choice we make in our methods.Since H β,w (x, 0) = F (x), both PCDM1 and PCDM2 are monotonic in expectation.
The above discussion leads to the following definition.

2.
Reshuffling.Since for any c > 0 we have h 2 cw = c h 2 w , one can "shuffle" constants between β and w as follows: Indeed, it suffices to take expectation in ( 14) with y replaced by x + h [ Ŝ] and compare the resulting inequality with (43) (this gives β h 2 w ≥ µ f (w) h 2 w , which must hold for all h).Recall that Step 5 of PCDM2 was introduced so as to explicitly enforce monotonicity into the method as in some situations, as we will see in Section 7, we can only analyze a monotonic algorithm.However, sometimes even PCDM1 behaves monotonically (without enforcing this behavior externally as in PCDM2).The following definition captures this.Definition 6 (Monotonic ESO).Assume (f, Ŝ) ∼ ESO(β, w) and let h(x) be as in (45).We say that the ESO is monotonic if F (x + (h(x)) [ Ŝ] ) ≤ F (x), with probability 1, for all x ∈ dom F .

Deterministic Separable Overapproximation (DSO) of Partially Separable Functions
The following theorem will be useful in deriving ESO for uniform samplings (Section 6.1) and nonoverlapping uniform samplings (Section 6.2).It will also be useful in establishing monotonicity of some ESOs (Theorems 12 and 13).
Theorem 7 (DSO).Assume f is partially separable (i.e., it can be written in the form (2)). Letting Proof.Let us fix x and define φ(h) One can define functions φ J in an analogous fashion from the constituent functions f J , which satisfy φ J (0) = 0, Note that ( 12) can be written as Now, since φ J depends on the intersection of J and the support of its argument only, we have The argument in the last expression can be written as a convex combination of 1 + θ J vectors: the zero vector (with weight θ−θ J θ ) and the θ J vectors {θU i h (i) : i ∈ J ∩ Supp(h)} (with weights 1 θ ): Finally, we plug (53) into (52) and use convexity and some simple algebra: Besides the usefulness of the above result in deriving ESO inequalities, it is interesting on its own for the following reasons.

2.
Global Lipschitz continuity of ∇f .The DSO inequality also says that the gradient of f is Lipschitz with Lipschitz constant ω with respect to the norm • L : Indeed, this follows from (48) via max J∈J |J ∩ Supp(h)| ≤ max J∈J |J| = ω.For ω = n this has been shown in [10]; our result for partially separable functions appears to be new.
3. Tightness of the global Lipschitz constant.The Lipschitz constant ω is "tight" in the following sense: there are functions for which ω cannot be replaced in (54) by any smaller number.We will show this on a simple example.Let f (x) = 1 2 Ax 2 with A ∈ R m×n (blocks are of size 1).Note that we can write f (x + h) = f (x) + ∇f (x), h + 1 2 h T A T Ah, and that L = (L 1 , . . ., L n ) = diag(A T A).Let D = Diag(L).We need to argue that there exists A for which = ω.Since we know that σ ≤ ω (otherwise (54) would not hold), all we need to show is that there is A and h for which , where A j is the j-th row of A, we assume that each row of A has at most ω nonzeros (i.e., f is partially separable of degree ω).Let us pick A with the following further properties: a) A is a 0-1 matrix, b) all rows of A have exactly ω ones, c) all columns of A have exactly the same number (k) of ones.Immediate consequences: L i = k for all i, D = kI n and ωm = kn.If we let e m be the m × 1 vector of all ones and e n be the n × 1 vector of all ones, and set h = k −1/2 e n , then . Using similar techniques one can easily prove the following more general result: Tightness also occurs for matrices A which in each row contain ω identical nonzero elements (but which can vary from row to row).

ESO for a convex combination of samplings
Let Ŝ1 , Ŝ2 , . . ., Ŝm be a collection of samplings and let q ∈ R m be a probability vector.By j q j Ŝj we denote the sampling Ŝ given by This procedure allows us to build new samplings from existing ones.A natural interpretation of Ŝ is that it arises from a two stage process as follows.Generating a set via Ŝ is equivalent to first choosing j with probability q j , and then generating a set via Ŝj .Lemma 8. Let Ŝ1 , Ŝ2 , . . ., Ŝm be arbitrary samplings, q ∈ R m a probability vector and κ : 2 [n] → R any function mapping subsets of [n] to reals.If we let Ŝ = j q j Ŝj , then (iii) P(i ∈ Ŝ) = m j=1 q j P(i ∈ Ŝj ), for any i = 1, 2, . . ., n, (iv) If Ŝ1 , . . ., Ŝm are uniform (resp.doubly uniform), so is Ŝ. .

Proof. Statement (i) follows by writing E[κ( Ŝ)] as
As the last expression depends on S via |S| only, Ŝ is doubly uniform. Remarks: 1.If we fix S ⊂ [n] and define k(S ) = 1 if S = S and k(S ) = 0 otherwise, then statement (i) of Lemma8 reduces to (56).
2. All samplings arise as a combination of elementary samplings, i.e., samplings whose all weight is on one set only.Indeed, let Ŝ be an arbitrary sampling.For all subsets S j of [n] define Ŝj by P( Ŝj = S j ) = 1 and let q j = P( Ŝ = S j ).Then clearly, Ŝ = j q j Ŝj .

All doubly uniform samplings arise as convex combinations of nice samplings.
Often it is easier to establish ESO for a simple class of samplings (e.g., nice samplings) and then use it to obtain an ESO for a more complicated class (e.g., doubly uniform samplings as they arise as convex combinations of nice samplings).The following result is helpful in this regard.
Theorem 9 (Convex Combination of Uniform Samplings).Let Ŝ1 , . . ., Ŝm be uniform samplings satisfying (f, Ŝj ) ∼ ESO(β j , w j ) and let q ∈ R m be a probability vector.If j q j Ŝj is not nil, then Proof.First note that from part (iv) of Lemma 8 we know that Ŝ def = j q j Ŝj is uniform and hence it makes sense to speak about ESO involving this sampling.Next, we can write It now remains to use (43) and part (ii) of Lemma 8: where w = j q j E[| Ŝj |]β j w j .In the third step we have also used the fact that E[| Ŝ|] > 0 which follows from the assumption that Ŝ is not nil.

ESO for a conic combination of functions
We now establish an ESO for a conic combination of functions each of which is already equipped with an ESO.It offers a complementary result to Theorem 9.

Uniform samplings
Consider an arbitrary proper sampling Ŝ and let ν = (ν 1 , . . ., ν n ) T be defined by Proof.Let us use Theorem 7 with h replaced by h Taking expectations of both sides of (48) we therefore get It remains to bound the last term in the expression above.Letting The above lemma will now be used to establish ESO for arbitrary (proper) uniform samplings.
Proof.First, (60) follows from (57) since for a uniform sampling one has p i = E[| Ŝ|]/n for all i.
If P(| Ŝ| = τ ) = 1, we get ν i = min{ω, τ } for all i; (61) therefore follows from (60).Let us now establish monotonicity.Using the deterministic separable overapproximation (48 Now let h(x) = arg min h H β,w (x, h) and recall that So, by definition, (h(x)) (i) minimizes κ i (t) and hence, (h(x)) [ Ŝ] (recall ( 7)) minimizes the upper bound (63).In particular, (h(x)) [ Ŝ] is better than a nil update, which immediately gives Besides establishing an ESO result, we have just shown that, in the case of τ -uniform samplings with a conservative estimate for β, PCDM1 is monotonic, i.e., F (x k+1 ) ≤ F (x k ).In particular, PCDM1 and PCDM2 coincide.We call the estimate β = min{ω, τ } "conservative" because it can be improved (made smaller) in special cases; e.g., for the τ -nice sampling.Indeed, Theorem 14 establishes an ESO for the τ -nice sampling with the same w (w = L), but with β = 1 + (ω−1)(τ −1) n−1 , which is better (and can be much better than) min{ω, τ }.Other things equal, smaller β directly translates into better complexity.The price for the small β in the case of the τ -nice sampling is the loss of monotonicity.This is not a problem for strongly convex objective, but for merely convex objective this is an issue as the analysis techniques we developed are only applicable to the monotonic method PCDM2 (see Theorem 19).
Theorem 13.If Ŝ a nonoverlapping uniform sampling, then Moreover, this ESO is monotonic.
Proof.By Theorem 7, used with h replaced by h [S j ] for j = 1, 2, . . ., l, we get Since Ŝ = S j with probability 1 l , which establishes (65).It now only remains to establish monotonicity.Adding Ω(x + h [ Ŝ] ) to (66) with S j replaced by Ŝ, we get

Nice samplings
In this section we establish an ESO for nice samplings.
Theorem 14.If Ŝ is the τ -nice sampling and τ = 0, then Proof.Let us fix x and define φ and φ J as in the proof of Theorem 7. Since max(1,n−1) Let us now adopt the convention that expectation conditional on an event which happens with probability 0 is equal to 0. Letting η J def = |J ∩ Ŝ|, and using this convention, we can write Note that the last identity follows if we assume, without loss of generality, that all sets J have the same cardinality ω (this can be achieved by introducing "dummy" dependencies).Indeed, in such a case P(η J = k) does not depend on J. Now, for any k ≥ 1 for which P(η J = k) > 0 (for some J and hence for all), using convexity of φ J , we can now estimate If we now sum the inequalities (70) for all J ∈ J , we get Finally, (68) follows after plugging (71) into (69): max(1,n−1)

Doubly uniform samplings
We are now ready, using a bootstrapping argument, to formulate and prove a result covering all doubly uniform samplings.
Theorem 15.If Ŝ is a (proper) doubly uniform sampling, then Proof.Letting q k = P(| Ŝ| = k) and d = max{1, n − 1}, we have This theorem could have alternatively been proved by writing Ŝ as a convex combination of nice samplings and applying Theorem 9.
Note that Theorem 15 reduces to that of Theorem 14 in the special case of a nice sampling, and gives the same result as Theorem 13 in the case of the serial and fully parallel samplings.

Iteration Complexity
In this section we prove two iteration complexity theorems 12 .The first result (Theorem 19) is for non-strongly-convex F and covers PCDM2 with no restrictions and PCDM1 only in the case when a monotonic ESO is used.The second result (Theorem 20) is for strongly convex F and covers PCDM1 without any monotonicity restrictions.Let us first establish two auxiliary results. Proof.
= min ≤ min Lemma 17. (i) Let x * be an optimal solution of (1), x ∈ dom F and let R = x − x * w .Then (ii) If µ f (w) + µ Ω (w) > 0 and β ≥ µ f (w), then for all x ∈ dom F , Proof.Part (i): Since we do not assume strong convexity, we have µ f (w) = 0, and hence Minimizing the last expression in λ gives λ * = min 1, (F (x) − F * )/(βR 2 ) ; the result follows.Part (ii): Letting ≤ min The last inequality follows from the identity We could have formulated part (ii) of the above result using the weaker assumption µ F (w) > 0, leading to a slightly stronger result.However, we prefer the above treatment as it gives more insight.

Iteration complexity: convex case
The following lemma will be used to finish off the proof of the complexity result of this section.
Lemma 18 (Theorem 1 in [16]).Fix x 0 ∈ R N and let {x k } k≥0 be a sequence of random vectors in R N with x k+1 depending on x k only.Let φ : R N → R be a nonnegative function and define ξ k = φ(x k ).Lastly, choose accuracy level 0 < < ξ 0 , confidence level 0 < ρ < 1, and assume that the sequence of random variables {ξ k } k≥0 is nonincreasing and has one of the following properties: If property (i) holds and we choose K ≥ 2 + c 1 (1 − ξ 0 + log( 1 ρ )), or if property (ii) holds, and we choose K ≥ c 2 log( ξ 0 ρ ), then P(ξ K ≤ ) ≥ 1 − ρ.This lemma was recently extended in [26] so as to aid the analysis of a serial coordinate descent method with inexact updates, i.e., with h(x) chosen as an approximate rather than exact minimizer of H 1,L (x, •) (see (17)).While in this paper we deal with exact updates only, the results can be extended to the inexact case.
Theorem 19.Assume that (f, Ŝ) ∼ ESO(β, w), where Ŝ is a proper uniform sampling, and let where x * is an optimal point of (1).Further, choose target confidence level 0 < ρ < 1, target accuracy level > 0 and iteration counter K in any of the following two ways: If {x k }, k ≥ 0, are the random iterates of PCDM (use PCDM1 if the ESO is monotonic, otherwise use PCDM2), then P(F Proof.Since either PCDM2 is used (which is monotonic) or otherwise the ESO is monotonic, we must have F (x k ) ≤ F (x 0 ) for all k.In particular, in view of (75) this implies that Consider case (i) and let c 1 = 2 β α max{R 2 w (x 0 , x * ), ξ 0 β }.Continuing with (78), we then get it suffices to apply Lemma 18(i).Consider now case (ii) and let By assumption, c 2 > 1, and hence it remains to apply Lemma 18(ii).
The important message of the above theorem is that the iteration complexity of our methods in the convex case is O( β α 1 ).Note that for the serial method (PCDM1 used with Ŝ being the serial sampling) we have α = 1 n and β = 1 (see Table 4), and hence β α = n.It will be interesting to study the parallelization speedup factor defined by parallelization speedup factor = Table 5, computed from the data in Table 4, gives expressions for the parallelization speedup factors for PCDM based on a DU sampling (expressions for 4 special cases are given as well).

Ŝ
Parallelization speedup factor doubly uniform fully parallel The speedup of the serial sampling (i.e., of the algorithm based on it) is 1 as we are comparing it to itself.On the other end of the spectrum is the fully parallel sampling with a speedup of n ω .If the degree of partial separability is small, then this factor will be high -especially so if n is huge, which is the domain we are interested in.This provides an affirmative answer to the research question stated in italics in the introduction.
Let us now look at the speedup factor in the case of a τ -nice sampling.Letting r = ω−1 max(1,n−1) ∈ [0, 1] (degree of partial separability normalized), the speedup factor can be written as .
Note that as long as r ≤ k−1 τ −1 ≈ k τ , the speedup factor will be at least τ k .Also note that max{1, τ ω } ≤ s(r) ≤ min{τ, n ω }.Finally, if a speedup of at least s is desired, where s ∈ [0, n ω ], one needs to use at least 1−r 1/s−r processors.For illustration, in Figure 1 we plotted s(r) for a few values of τ .Note that for small values of τ , the speedup is significant and can be as large as the number of processors (in the separable case).We wish to stress that in many applications ω will be a constant independent of n, which means that r will indeed be very small in the huge-scale optimization setting.

Iteration complexity: strongly convex case
In this section we assume that F is strongly convex with respect to the norm • w and show that F (x k ) converges to F * linearly, with high probability.
Theorem 20.Assume F is strongly convex with µ f (w) + µ Ω (w) > 0. Further, assume (f, Ŝ) ∼ ESO(β, w), where Ŝ is a proper uniform sampling and let α = E[| Ŝ|] n .Choose initial point x 0 ∈ dom F , target confidence level 0 < ρ < 1, target accuracy level 0 < < F (x 0 ) − F * and If {x k } are the random points generated by PCDM1 or PCDM2, then P(F Note that 0 < γ ≤ 1 since 0 < α ≤ 1 and β ≥ µ f (w) by (47).By taking expectation in x k , we obtain E[ξ k ] ≤ (1 − γ) k ξ 0 .Finally, it remains to use Markov inequality: Instead of doing a direct calculation, we could have finished the proof of Theorem 20 by applying Lemma 18(ii) to the inequality E[ξ k+1 | x k ] ≤ (1 − γ)ξ k .However, in order to be able to use Lemma 18, we would have to first establish monotonicity of the sequence {ξ k }, k ≥ 0. This is not necessary using the direct approach of Theorem 20.Hence, in the strongly convex case we can analyze PCDM1 and are not forced to resort to PCDM2.Consider now the following situations: 1. µ f (w) = 0. Then the leading term in (80) is 1+β/µ Ω (w) α .
In a similar way as in the non-strongly convex case, define the parallelization speedup factor as the ratio of the leading term in (80) for the serial method (which has α = 1 n and β = 1) and the leading term for a parallel method: . (81) First, note that the speedup factor is independent of µ f .Further, note that as µ Ω (w) → 0, the speedup factor approaches the factor we obtained in the non-strongly convex case (see (79) and also Table 5).That is, for large values of µ Ω (w), the speedup factor is approximately equal αn = E[| Ŝ|], which is the average number of blocks updated in a single parallel iteration.Note that thuis quantity does not depend on the degree of partial separability of f .

Numerical Experiments
In Section 8.1 we present preliminary but very encouraging results showing that PCDM1 run on a system with 24 cores can solve huge-scale partially-separable LASSO problems with a billion variables in 2 hours, compared with 41 hours on a single core.In Section 8.2 we demonstrate that our analysis is in some sense tight.In particular, we show that the speedup predicted by the theory can be matched almost exactly by actual wall time speedup for a particular problem.

A LASSO problem with 1 billion variables
In this experiment we solve a single randomly generated huge-scale LASSO instance, i.e., (1) with where A = [a 1 , . . ., a n ] has 2 × 10 9 rows and N = n = 10 9 columns.We generated the problem using a modified primal-dual generator [16] enabling us to choose the optimal solution x * (and hence, indirectly, F * ) and thus to control its cardinality x * 0 , as well as the sparsity level of A. In particular, we made the following choices: x * 0 = 10 5 , each column of A has exactly 20 nonzeros and the maximum cardinality of a row of A is ω = 35 (the degree of partial separability of f ).The histogram of cardinalities is displayed in Figure 2.
We solved the problem using PCDM1 with τ -nice sampling Ŝ, β = 1 + (ω−1)(τ −1) n−1 2 ), for τ = 1, 2, 4, 8, 16, 24, on a single large-memory computer utilizing τ of its 24 cores.The problem description took around 350GB of memory space.In fact, in our implementation we departed from the just described setup in two ways.First, we implemented an asynchronous version of the method; i.e., one in which cores do not wait for others to update the current iterate within an iteration before reading x k+1 and proceeding to another update step.Instead, each core reads the current iterate whenever it is ready with the previous update step and applies the new update as soon as it is computed.Second, as mentioned in Section 4, the τ -independent sampling is for τ n a very good approximation of the τ -nice sampling.We therefore allowed each processor to pick a block uniformly at random, independently from the other processors.
Choice of the first column of Table 6.In Table 6 we show the development of the gap F (x k ) − F * as well as the elapsed time.The choice and meaning of the first column of the table, τ k n , needs some commentary.Note that exactly τ k coordinate updates are performed after k iterations.Hence, the first column denotes the total number of coordinate updates normalized by the number of coordinates n.As an example, let τ 1 = 1 and τ 2 = 24.Then if the serial method is run for k 1 = 24 iterations and the parallel one for k 2 = 1 iteration, both methods would have updated the same number (τ 1 k 1 = τ 2 k 2 = 24) of coordinates; that is, they would "be" in the same row of Table 6.In summary, each row of the table represents, in the sense described above, the "same amount of work done" for each choice of τ .
Progress to solving the problem.One can conjecture that the above meaning of the phrase "same amount of work done" would perhaps be roughly equivalent to a different one: "same progress to solving the problem".Indeed, it turns out, as can be seen from the table and also from Figure 3(a), that in each row for all algorithms the value of F (x k ) − F * is roughly of the same order of magnitude.This is not a trivial finding since, with increasing τ , older information is used to update the coordinates, and hence one would expect that convergence would be slower.It does seem to be slower-the gap F (x k ) − F * is generally higher if more processors are used-but the slowdown is limited.Looking at Table 6 and/or Figure 3(a), we see that for all choices of τ , PCDM1 managed to push the gap below 10 −13 after 34n to 37n coordinate updates.
The progress to solving the problem during the final 1 billion coordinate updates (i.e., when moving from the last-but-one to the last nonempty line in each of the columns of Table 6 showing F (x k ) − F * ) is remarkable.The method managed to push the optimality gap by 9-12 degrees of magnitude.We do not have an explanation for this phenomenon; we do not give local convergence  estimates in this paper.It is certainly the case though that once the method managed to find the nonzero places of x * , fast local convergence comes in.
Parallelization speedup.Since a parallel method utilizing τ cores manages to do the same number of coordinate updates as the serial one τ times faster, a direct consequence of the above observation is that doubling the number of cores corresponds to roughly halving the number of iterations (see Figure 3(b).This is due to the fact that ω n and τ n.It turns out that the number of iterations is an excellent predictor of wall time; this can be seen by comparing Figures 3(b) and 3(c).Finally, it follows from the above, and can be seen in Figure 3(d), that the speedup of PCDM1 utilizing τ cores is roughly equal to τ .Note that this is caused by the fact that the problem is, relative to its dimension, partially separable to a very high degree.

Theory versus reality
In our second experiment we demonstrate numerically that our parallelization speedup estimates are in some sense tight.For this purpose it is not necessary to reach for complicated problems and high dimensions; we hence minimize the function 1  2 Ax − b 2 2 with A ∈ R 3000×1000 .Matrix A was generated so that its every row contains exactly ω non-zero values all of which are equal (recall the construction in point 3 at the end of Section 5.1).We generated 4 matrices with ω = 5, 10, 50 and 100 and measured the number of iterations needed for PCDM1 used with τ -nice sampling to get within = 10 −6 of the optimal value.The experiment was done for a range of values of τ (between 1 core and 1000 cores).
The solid lines in Figure 4 present the theoretical speedup factor for the τ -nice sampling, as presented in Table 5.The markers in each case correspond to empirical speedup factor defined as # of iterations till -solution is found by PCDM1 used with serial sampling # of iterations till -solution is found by PCDM1 used with τ -nice sampling .
As can be seen in Figure 4, the match between theoretical prediction and reality is remarkable!A partial explanation of this phenomenon lies in the fact that we have carefully designed the problem so as to ensure that the degree of partial separability is equal to the Lipschitz constant σ of ∇f (i.e., that it is not a gross overestimation of it; see Section 5.1).This fact is useful since it is possible to prove complexity results with ω replaced by σ.However, this answer is far from satisfying, and a deeper understanding of the phenomenon remains an open problem.x (i) , where Z ∈ R n×n with Z ii = y i y j A i , A j .It is a standard practice to apply serial coordinate descent to the dual.Here we apply parallel coordinate descent (PCDM; with τ -nice sampling of coordinates) to the dual; i.e., minimize the convex function f subject to box constraints.In this setting all blocks are of size N i = 1.The dual can be written in the form (1), i.e., min x∈R n {F (x) = f (x) + Ω(x)}, where Ω(x) = 0 whenever x (i) ∈ [0, 1] for all i = 1, 2, . . ., n, and Ω(x) = +∞ otherwise.
We consider the rcv1.binarydataset 13 .The training data has n = 677, 399 examples, d = 47, 236 features, 49, 556, 258 nonzero elements and requires cca 1GB of RAM for storage.Hence, this is a small-scale problem.The degree of partial separability of f is ω = 291, 516 (i.e., the maximum number of examples sharing a given feature).This is a very large number relative to n, and hence our theory would predict rather bad behavior for PCDM.We use PCDM1 with τnice sampling ( approximating it by τ -independent sampling for added efficiency) with β following Theorem 14: β = 1 + (τ −1)(ω−1) n−1 .The results of our experiments are summarized in Figure 5.Each column corresponds to a different level of regularization: λ ∈ {1, 10 −3 , 10 −5 }.The rows show the 1) duality gap, 2) dual suboptimality, 3) train error and 4) test error; each for 1,4 and 16 processors (τ = 1, 4, 16).Observe that the plots in the first two rows are nearly identical; which means that the method is able to solve the primal problem at about the same speed as it can solve the dual problem 14 .
Observe also that in all cases, duality gap of around 0.01 is sufficient for training as training error (classification performance of the SVM on the train data) does not decrease further after this point.Also observe the effect of λ on training accuracy: accuracy increases from about 92% for λ = 1, through 95.3% for λ = 10 −3 to above 97.8% with λ = 10 −5 .In our case, choosing smaller λ does not lead to overfitting; the test error on test dataset (# features =677,399, # examples = 20,242) increases as λ decreases, quickly reaching about 95% (after 2 seconds of training) for λ = 0.001 and for the smallest λ going beyond 97%.
Note that PCDM with τ = 16 is about 2.5× faster than PCDM with τ = 1.This is much less than linear speedup, but is fully in line with our theoretical predictions.Indeed, for τ = 16 we get β = 7.46.Consulting Table 5, we see that the theory says that with τ = 16 processors we should expect the parallelization speedup to be P SF = τ /β = 2.15.This training dataset is good for PCDM as each example depends on at most 75 features.That is, ω = 75, which is much smaller than n.As before, we will use PCDM1 with τ -nice sampling (approximated by τ -independent sampling) for τ = 1, 2, 4, 8 and set λ = 1.
Figure 6 depicts the evolution of the regularized loss F (x k ) throughout the run of the 4 versions of PCDM (starting with x 0 for which F (x 0 ) = 13, 352, 855).Each marker corresponds to approximately n/3 coordinate updates (n coordinate updates will be referred to as an "epoch").Observe that as more processors are used, it takes less time to achieve any given level of loss; nearly in exact proportion to the increase in the number of processors.Table 7 offers an alternative view of the same experiment.In the first 4 columns (F (x 0 )/F (x k )) we can see that no matter how many processors are used, the methods produce similar loss values after working through the same number of coordinates.However, since the method utilizing τ = 8 processors updates 8 coordinates in parallel, it does the job approximately 8 times faster.Indeed, we can see this speedup in the table.

n ω serial 1 Table 5 :
Convex F : Parallelization speedup factors for DU samplings.The factors below the line are special cases of the general expression.Maximum speedup is naturally obtained by the fully parallel sampling: n ω .

Figure 1 :
Figure 1: Parallelization speedup factor of PCDM1/PCDM2 used with τ -nice sampling as a function of the normalized/relative degree of partial separability r.

Figure 2 :
Figure 2: Histogram of the cardinalities of the rows of A.
For each τ , PCDM1 needs roughly the same number of coordinate updates to solve the problem.Doubling the number of cores corresponds to roughly halving the number of iterations.Doubling the number of cores corresponds to roughly halving the wall time.Parallelization speedup is essentially equal to the number of cores.

Figure 3 :
Figure 3: Four computational insights into the workings of PCDM1.

Figure 4 :
Figure 4: Theoretical speedup factor predicts the actual speedup almost exactly for a carefully constructed problem.

8. 3 [ 1 −λ 2 w 2 2 ,
Training linear SVMs with bad data for PCDM In this experiment we test PCDM on the problem of training a linear Support Vector Machine (SVM) based on n labeled training examples:(y i , A i ) ∈ {+1, −1} × R d , i = 1, 2, . . ., n.In particular, we consider the primal problem of minimizing L2-regularized average hinge-loss, y i w, a i ] + + and the dual problem of maximizing a concave quadratic subject to zero-one box constraints,max x∈R n , 0≤x (i) ≤1 −f (x)

Figure 5 :
Figure5: The performance of PCDM on the rcv1 dataset (this dataset is not good for the method.).

8. 4
L2-regularized logistic regression with good data for PCDMIn our last experiment we solve a problem of the form (1) with f being a sum of logistic losses and Ω being an L2 regularizer, j , A j ) ∈ {+1, −1} × R n , j = 1, 2, . . ., d, are labeled examples.We have used the the KDDB dataset from the same source as the rcv1.binarydataset considered in the previous experiment.The data contains n = 29, 890, 095 features and is divided into two parts: a training set with d = 19, 264, 097 examples (and 566, 345, 888 nonzeros; cca 8.5 GB) and a testing with d = 748, 401 examples (and 21, 965, 075 nonzeros; cca 0.32 GB).

Figure 6 :
Figure 6: PCDM accelerates well with more processors on a dataset with small ω.

Table 2 :
Summary of the main complexity results for PCDM established in this paper.

Table 3 :
New and known gradient methods obtained as special cases of our general framework.

Table 4 :
Values of parameters β and w for various samplings Ŝ.