Coordinate Descent with Arbitrary Sampling I: Algorithms and Complexity

We study the problem of minimizing the sum of a smooth convex function and a convex block-separable regularizer and propose a new randomized coordinate descent method, which we call ALPHA. Our method at every iteration updates a random subset of coordinates, following an arbitrary distribution. No coordinate descent methods capable to handle an arbitrary sampling have been studied in the literature before for this problem. ALPHA is a remarkably flexible algorithm: in special cases, it reduces to deterministic and randomized methods such as gradient descent, coordinate descent, parallel coordinate descent and distributed coordinate descent -- both in nonaccelerated and accelerated variants. The variants with arbitrary (or importance) sampling are new. We provide a complexity analysis of ALPHA, from which we deduce as a direct corollary complexity bounds for its many variants, all matching or improving best known bounds.


Introduction
With the dawn of the big data age, there has been a growing interest in solving optimization problems of unprecedented sizes.It was soon realized that traditional approaches, which work extremely well for problems of moderate sizes and when solutions of high accuracy are required, are not efficient for modern problems of large enough size and for applications where only rough or moderate accuracy solutions are sufficient.The focus of the optimization, numerical analysis and machine learning communities, and of practitioners in the sciences and industry, shifted to first-order (gradient) algorithms [24].
However, once the size of problems becomes truly big, it is necessary to turn to methods which are able to output a reasonably good solutions after an amount of work roughly equivalent to reading the data describing the problem a few times.For this to be possible, methods need to be able to progress while reading only a small part of the data describing the problem, which often means that a single iteration needs to be based on less information than that contained in the gradient of the objective (loss) function.The most popular methods of this type are stochastic gradient methods [45,23,35,39,47], randomized coordinate descent methods [4,8,33,32,39,41,40,7,29,5,47,37,38,14,27,11] and semi-stochastic gradient descent methods [34,44,9,12,18,19,3,43,10,11].

Randomized coordinate descent
In this paper we focus on randomized coordinate descent methods.After the seminal work of Nesterov [25], which provided an early theoretical justification of these methods for unconstrained convex minimization, the study has been successively extended to L1-regularized [36,31], proximal [33,17], parallel [32,6], distributed [29,5,27] and primal-dual [39,27] variants of coordinate descent.Accelerated coordinate decent-characterized by its O(1/k 2 ) complexity for non-strongly convex problems-was studied in [25,17].However, these methods are of theoretical nature only due to the fact that they rely on the need to perform full-dimensional vector operation at every iteration, which destroys the main advantage of coordinate descent -its ability to reduce the problem into subproblems of smaller sizes.A theoretically and practically efficient accelerated coordinate descent methods were proposed recently by Lee and Sidford [13] and Fercoq and Richtárik [6], the latter work (APPROX algorithm) combining acceleration with parallelism and proximal setup.An accelerated distributed coordinate descent algorithm [5] is obtained by specializing APPROX to a distributed sampling.All above mentioned papers only consider unconstrained or separably constrained problems.Some progress on linearly-coupled constraints has been made by Necoara et al in [20,21,22].Asynchronous variants of parallel coordinate descent methods were developed by Liu, Wright et al [16,15].
Virtually all existing work in stochastic optimization deals with a uniform sampling.In the context of coordinate descent, this means that the random subset (sampling) of coordinates chosen and updated at every iteration has the property that each coordinate is chosen equally likely.The possibility to assign different selection probabilities to different coordinates--also known as importance sampling-was considered in [25,33] and recently in [47,46].However, these works consider the serial case only: a single coordinate is updated in each iteration.Randomized coordinate descent methods updating a subset of coordinates following an arbitrary distribution (i.e., using an arbitrary sampling) was first investigated by Richtárik and Takáč [30] (NSync method) for strongly convex and smooth objective functions, and subsequently by Qu, Richtárik and Zhang [27] (QUARTZ method), for strongly convex and possibly nonsmooth functions, and in a primal-dual framework.

Problem Formulation
In this paper we consider the composite optimization problem minimize F (x) where f : R N → R is convex and differentiable, ψ : R N → R ∪ {+∞} is block-separable: and each ψ i : R N i → R ∪ {+∞} is convex and closed.

Contributions
We now summarize the main contributions of this work.
New algorithm.We propose ALPHA (Algorithm 1) -a new randomized block coordinate descent method for solving the convex composite optimization problem (1).In each iteration, ALPHA picks and updates a random subset of the blocks {1, 2, . . ., n}, using an arbitrary sampling.That is, we allow for the distribution of the random set-valued mapping to be arbitrary.To the best of our knowledge, there are only two randomized coordinate descent methods in the literature with such a property, the NSync method of Richtárik and Takáč [30] (focusing in smooth strongly convex functions) and the QUARTZ method of Qu, Richtárik and Zhang [27] (a primal-dual method; considering strongly convex but possibly nonsmooth functions).Hence, our work is complementary to this development.
Complexity analysis.We study the iteration complexity of ALPHA.That is, we provide bounds on the number of iterations needed to approximately solve the problem, in expectation.Our general bounds are formulated in Section 6: Theorem 6.1 covers the non-accelerated variant with O(1/k) rate and Theorem 6.2 covers the accelerated variant with O(1/k 2 ) rate, where k is the iteration counter.To the best of our knowledge, these are the first complexity results for a randomized coordinate descent methods utilizing an arbitrary sampling for problem (1).
Expected separable overapproximation.Besides the dependence of the complexity bound on the iteration counter k, it is important to study its dependence on the sampling Ŝ and the objective function.Our results make this dependence explicit: they hold under the assumption that f admits an expected separable overapproximation (ESO) with respect to the sampling Ŝ (Assumption 2.1).This is an inequality involving f and Ŝ which determines certain important parameters v 1 , . . ., v n which are needed to run the method (they determine the stepsizes) and which also appear in the complexity bounds.In some cases it is possible to design a sampling which optimizes the complexity bound.
In the case of a serial sampling, which is by far the most common type of sampling studied in conjunction with randomized coordinate descent methods ( Ŝ is serial if | Ŝ| = 1 with probability 1), the parameter v i can simply be set to the Lipschitz constant of the block-derivative of f corresponding to block i.In particular, if n = 1, then v 1 is the Lipschitz constant of the gradient of f [25,33].The situation is more complicated in the case of a parallel sampling ( Ŝ is parallel if it is not serial; that is, if we allow for multiple blocks to be updated at every iteration) -and this why there is a need for the ESO inequality.Intuitively speaking, the parameters v 1 , . . ., v n capture certain smoothness properties of the gradient of f in a random subspace spanned by the blocks selected by the sampling Ŝ.
The ESO concept is of key importance in the design and analysis of randomized coordinate descent methods [32,7,29,30,6,5,27,26].We provide a systematic study of ESO inequalities in a companion paper [26].
Simple complexity analysis in the smooth case.In order to make the exposition more accessible, we first focus on ALPHA applied to problem (1) with ψ ≡ 0 (we call this the "smooth case").In this simpler setting, it is possible to provide a simplified complexity analysis -we do this in Section 3; see Theorem 3.1 (non-accelerated variant with O(1/k) rate) and Theorem 3.2 (accelerated variant with O(1/k2 ) rate).For convenience, ALPHA specialized to the smooth case is formulated as Algorithm 2. Our analysis in this case is different from the one we give in Section 6, where we analyze the method in the general proximal setup.
Flexibility.ALPHA is a remarkably flexible1 algorithm, encoding a number of classical, recent and new algorithms in special cases, depending on the choice of the parameters of the method: sampling Ŝ and "stepsize sequence" {θ k }.We devote Section 4 to highlighting several of the many algorithms ALPHA reduces to in special cases, focusing on the smooth case for simplicity (special cases in the general proximal setting are discussed in the appendix).In particular, if Ŝ = {1, 2, . . ., n} with probability 1, ALPHA reduces to a deterministic method: gradient descent (GD; Algorithm 3) or accelerated gradient descent (AGD; Algorithm 4), depending on the choice of the sequence {θ k }.For a non-deterministic sampling, we obtain parallel coordinate descent (PCD; Algorithm 5) and accelerated parallel coordinate descent (APCD; Algorithm 6) with arbitrary sampling -which is new.If a uniform sampling is used, PCD reduces to the PCDM algorithm [32].If a distributed 2 sampling is used instead, PCD reduces to Hydra [29].Similarly, if a uniform sampling is used, APCD reduces to APPROX [6] (in fact, our version of APPROX is a bit more flexible with respect to choice of θ 0 , which leads to a better complexity result).APCD specialized to a distributed sampling reduces to Hydra 2 [5].

Improved complexity results
. In all cases where ALPHA reduces to an existing method, our complexity bound either matches the best known bound for that method or improves upon the best known bound.For instance, while the complexity of PCDM [32] (which coincides with PCD specialized to a uniform sampling; Algorithm 5) depends on the size of a certain level-set of f (and in particular, requires the level set to be bounded), our bound does not involve this quantity (see Section A.3). Another example is APPROX [6] (which is closely related to APCD specialized to a uniform sampling): we obtain a more compact and improved result.
Two in one.We provide a unified complexity analysis covering the nonaccelerated and accelerated variants of ALPHA.This is achieved by establishing a certain key technical recursion (Lemma 6.4) for an arbitrary choice of the parameters {θ k }.Since the two variants of ALPHA differ in the choice of this sequence only, the analysis of both is identical up to this point.The recursion is then analyzed in two different ways, depending on the sequence {θ k }, which leads to the final complexity result.
Efficient implementation.As formulated in Algorithm 1, ALPHA seems to require that two vectors in R N be added at each iteration (unless the three sequences coincide, which happens in some important special cases).Motivated by [13,6], in Section 5 we give an equivalent form of Algorithm 1, which under some structural assumptions on f (see (51)) does not require such fulldimensional operations.This is important as the efficiency of coordinate descent methods largely stems from their ability to decompose the problem into subproblems, in an iterative fashion, of much smaller size than is the size of the original problem.

Outline of the paper
The paper is organized as follows.In Section 2 we establish notation, describe the ALPHA algorithm and comment on the key assumption: Expected Separable Overapproximation (ESO).We defer the in-depth study ESO inequalities to a companion paper [26].In Section 3 we give a simple complexity proof of ALPHA in the smooth case (ψ = 0).Subsequently, in Section 4 we present four algorithms that ALPHA reduces to in special cases, and state the corresponding complexity results, which follow from our general result, in a simplified form.In Section 5 we provide an equivalent form of writing ALPHA-one leading to an efficient implementation avoiding full dimensional operations.In Section 6 we state and prove the convergence result of ALPHA when applied to the general proximal minimization problem (1).Finally, in Section 7 we conclude and in the appendix we comment on several special cases ALPHA reduces to in the general proximal setup.

The Algorithm
In this section we first formalize the block structure of R N and establish necessary notation (Section 2.3), then proceed to describing the ALPHA algorithm (Section 2.2) and finally comment in the key assumption needed for our complexity results (Section 2.3).

Preliminaries
Blocks.We first describe the block setup which has become standard in the analysis of block coordinate descent methods [25,33,32,6].The space R N is decomposed into n subspaces: Let U be the N × N identity matrix and U = [U 1 , . . ., U n ] be its decomposition into column submatrices U i ∈ R N ×N i .For x ∈ R N , let x i ∈ R N i be the block of coordinates corresponding to the columns of U i , i.e., x i = U ⊤ i x.For any h ∈ R N and S ⊆ [n] def = {1, . . ., n} we define: For function f , we denote by ∇f (x) the gradient of f at point x ∈ R N and by Norms.The standard Euclidean inner product (with respect to the standard basis) in spaces R N and R N i , i ∈ [n], will be denoted by •, • .That is, for vectors x, y of equal size, we have x, y = x ⊤ y.Each space R N i is equipped with a Euclidean norm: where B i is an N i -by-N i positive definite matrix.For w ∈ R n ++ and x, y ∈ R N we further define and For x ∈ R N , by Bx we mean the vector Bx = n i=1 U i B i x i .That is, Bx is the vector in R N whose ith block is equal to B i x i .For vectors x, y ∈ R N we have Vectors.For any two vectors x and y of the same size, we denote by x • y their Hadamard (i.e., elementwise) product.By abuse of notation, we denote by u 2 the elementwise square of the vector u, by u −1 the elementwise inverse of vector u and by u −2 the elementwise square of u −1 .For vector v ∈ R n and x ∈ R N we will write That is, vx is the vector obtained from x by multiplying its block i by v i for each i ∈ [n].If all blocks are of size one (N i = 1 for all i), then vx = Diag(v)x where Diag(v) is the diagonal matrix with diagonal vector v.

ALPHA
In this section we describe ALPHA (Algorithm 1) -a randomized block coordinate descent method for solving (1).
We denote by dom ψ the domain of the proximal term ψ.

Algorithm 1 ALPHA
1: Parameters: proper sampling Ŝ with probability vector p Generate a random set of blocks S k ∼ Ŝ 6: for i ∈ S k do 8: end for 10: To facilitate the presentation, we first recall some basic facts and terminology related to samplings (for a broader coverage see [32,6,27,26].A sampling Ŝ is a random set valued mapping with values in 2 [n] .We say that sampling Ŝ is and serial if P(| Ŝ| = 1) = 1.The probability that the block i is chosen is denoted by: It is easy to see that if Ŝ is uniform, then Algorithm 1 starts from an initial vector x 0 ∈ R N and generates three sequences {x k , y k , z k } k 0 .At iteration k, a random subset of blocks, S k ⊆ [n], is generated according to the distribution of sampling Ŝ (a parameter to enter at the beginning of the algorithm).In order to guarantee that each block has a nonzero probability to be selected, it is necessary to assume that Ŝ is proper.
To move from z k to z k+1 , we only need to evaluate |S k | partial derivatives of f at point y k and update only the blocks of z k belonging to S k to the solutions of |S k | proximal problems (Steps 6 to 8).The vector x k+1 is obtained from y k by changing only the blocks of y k belonging to S k (Step 10).The vector y k is a convex combination of x k and z k (Step 4) with coefficient θ k , a parameter to be chosen between (0, 1].Note that unless θ k p i = 1 for all i ∈ [n] and k ∈ N, in which case the three sequences {x k , y k , z k } k 0 reduce to one same sequence, the update in Step 4 requires a full-dimensional vector operation, as previously remarked in [25,6].In Section 5 we will provide an equivalent form of Algorithm 1, which avoids full-dimensional vector operations for special forms of f .
Let us extract the relations between the three sequences.Define Then and hence Note also that from the definition of y k in Algorithm 1, we have:

Expected separable overapproximation
To guarantee the convergence of Algorithm 1, we shall require that f admits an expected separable overapproximation (ESO) with respect to the sampling Ŝ with parameter v ∈ R n ++ : Assumption 2.1 (ESO assumption).Let Ŝ be a sampling and v ∈ R n ++ a vector of positive weights.We say that the function f admits an ESO with respect to Ŝ with parameter v, denoted as (f, Ŝ) ∼ ESO(v), if the following inequality holds for all x, h ∈ R N , where p = (p 1 , . . ., p n ) ∈ R n is the vector of probabilities associated with Ŝ defined in (8).
Looking behind the compact notation in which the assumption is formulated, observe that the upper bound is a quadratic function of h = (h 1 , . . ., h n ), separable in the blocks h i : As a tool for the design and analysis of randomized coordinate descent methods, ESO was first formulated in [32] for the complexity study of parallel coordinate descent method (PCDM).It is a powerful technical tool which provides a generic approach to establishing the convergence of randomized coordinate descent methods of many flavours [32,39,7,40,30,29,5,27].As shown in the listed papers as well as our results which follow, the convergence of Algorithm 1 can be established for arbitrary sampling Ŝ as long as the parameter vector v is chosen such that (f, Ŝ) ∼ ESO(v).Moreover, the vector v appears in the convergence result and directly influences the complexity of the method.
Since [32], the problem of computing efficiently a vector v such that the ESO assumption 2.1 holds has been addressed in many papers for special uniform samplings relevant to practical implementation including serial sampling, τ -nice sampling [32,6] and distributed sampling [29,5], and also for a particular example of nonuniform parallel sampling [30].In this paper we focus on the complexity analysis of Algorithm 1 and refer the reader to the companion paper [26] for a systematic study of the computation of admissible vector parameter v for arbitrary sampling Ŝ.

Simple complexity analysis in the smooth case
In this section we give a brief complexity analysis of ALPHA in the case when ψ ≡ 0; that is, when applied to the following unconstrained smooth convex minimization problem: While our general theory, which we develop in Section 6, covers also this special case, the analysis we present here is different and simpler.When applied to problem 14, Step 8 in ALPHA has an explicit solution and the method reduces to Algorithm 2.
Algorithm 2 ALPHA specialized to the smooth minimization problem (14) 1: Parameters: proper sampling Ŝ with probability vector p Generate a random set of blocks S k ∼ Ŝ 6: for i ∈ S k do 8: end for 10: We now state the complexity result for ALPHA (Algorithm 2) in its nonaccelerated variant.Theorem 3.1 (ALPHA -smooth & nonaccelerated).Let Ŝ be an arbitrary proper sampling and v ∈ R n ++ be such that (f, Ŝ) ∼ ESO(v).Choose θ k = θ 0 ∈ (0, 1] for all k 0. Then for any y ∈ R N , the iterates {x k } k 1 of Algorithm 2 satisfy: where In particular, if we choose The next result gives a complexity bound for ALPHA (Algorithm 2) in its accelerated variant.
Theorem 3.2 (ALPHA -smooth and accelerated).Let Ŝ be an arbitrary proper sampling and v ∈ R n ++ be such that (f, Ŝ) ∼ ESO(v).Choose θ 0 ∈ (0, 1] and define the sequence {θ k } k 0 by Then for any y ∈ R N such that C 0, the iterates {x k } k 1 of Algorithm 2 satisfy: where In particular, if we choose θ 0 = 1, then for all k 1, In the rest of this section, we provide a short proof of Theorems 3.1 and 3.2.In Section 6 we shall present complexity bounds (Theorem 6.1 and 6.2) for ALPHA (Algorithm 1) as applied to the general regularized problem (1).The proof in the general case is more involved, which is why we prefer to present the smooth case first and also provide a separate briefer proof.

Two Lemmas
We first establish two lemmas and then proceed directly to the proofs of the theorems.Lemma 3.1.For any sampling Ŝ and any x, a ∈ R N and w ∈ R n ++ , the following identity holds: Lemma 3.2.Let Ŝ be an arbitrary proper sampling and v ∈ R n ++ be such that (f, Ŝ) ∼ ESO(v).Let {θ k } k 0 be an arbitrary sequence of positive numbers in (0, 1] and fix y ∈ R N .Then for the sequence of iterates produced by Algorithm 2 and all k 0, the following recursion holds: Proof.Based on line 8 of Algorithm 2, we can write or equivalently, −∇f (y k ) = θ k (v•p −1 )Ba.Using this notation, the update on line 10 of Algorithm 2 can be written as Letting b = zk+1 − y and t = θ 2 k (v • p −1 ), we apply the ESO assumption and rearrange the result: Substituting these expressions to (23), we obtain the recursion: It now only remains to apply Lemma 3.1 to the last two terms in (24), with x ← z k − y, w ← v • p −2 and Ŝ ← S k , and rearrange the resulting inequality.

Proof of Theorem 3.1
Using the fact that θ k = θ 0 , for all k and taking expectation in both sides of (20), we obtain the recursion where . Combining these inequalities, we get Let Finally, subtracting f (y) from both sides and taking expectations, we obtain

Proof of Theorem 3.2
If θ 0 ∈ (0, 1], then the sequence {θ k } k 0 has the following properties (see [42]): After dividing both sides of (20) by θ 2 k , using (27) and taking expectations, we obtain: where φ k and r k are as in the proof of Theorem 3.1.Finally, Note that in the last inequality we used the assumption that C 0.

The many variants of ALPHA (in the smooth case)
The purpose of this section is to demonstrate that ALPHA is a very flexible method, encoding several classical as well as modern optimization methods for special choices of the parameters of the method.In order to achieve this goal, it is enough to focus on the smooth case, i.e., on Algorithm 2. Similar reasoning can be applied to the proximal case.
Note that in ALPHA we have the liberty to choose the sampling Ŝ and the sequence {θ k } k 0 .As we have already seen, by modifying the sequence we can obtain simple (i.e., nonaccelerated) and accelerated variants of the method.By the choice of the sampling, we can force the method to be deterministic or randomized.In the latter case, there are many ways of choosing the distribution of the sampling.Here we will constrain ourselves to a basic classification between uniform samplings (samplings for which p i = p i ′ for all i ∈ [n]) and non-uniform or importance samplings.This is summarized in Table 1.
The deterministic variants of ALPHA (Algorithm 3 and 4) are obtained by choosing the sampling which always selects all blocks: Ŝ = [n] with probability 1.The ESO assumption in this special case has the form which simply requires the gradient of f to be 1-Lipschitz with respect to the norm Likewise, if there is just a single block in our block setup (i.e., if n = 1) then it is natural to only consider a sampling which picks this block with probability 1, which again results in a deterministic method.However, in this case the norm • v can be an arbitrary Euclidean norm (that is, it does not need to be block diagonal).
In the randomized variants of ALPHA (Algorithm 5 and 6) we allow for the sampling to have an arbitrary distribution.

Special case 1: gradient descent
By specializing Algorithm 2 to the choice Ŝ = [n] and θ k = 1 for all k, we obtain classical gradient descent (with fixed stepsize).Indeed, note that in this special case we have Recall that the ESO assumption reduces to (29) when Ŝ = [n].
Algorithm 3 Gradient Descent (GD) for solving ( 14) end for 7: end for The complexity of the method is a corollary of Theorem 3.1.
Corollary 4.1.For any optimal solution x * of (14), the output of Algorithm 3 for all k 1 satisfies: In particular, for ǫ > 0, if Proof.By letting y = x k in (20) we know that: Note that for this special case where the second inequality follows from applying Theorem 3.1 to θ 0 = 1 and Ŝ = [n].
Corollary 4.1 is a basic result and can be found in many textbooks on convex optimization; see for example [24].

Special case 2: accelerated gradient descent
Let us still keep Ŝ = [n], but assume now the sequence {θ k } k 0 is chosen accoridng to (17).In this case, Algorithm 2 reduces to accelerated gradient descent.

Special case 3: Parallel coordinate descent
We now allow the method (Algorithm 2) to use an arbitrary sampling Ŝ, but keep θ k = θ 0 for all k 1.This leads to Algorithm 5.
Algorithm 5 Parallel Coordinate Descent (PCD) for solving (14) 1: Parameters: proper sampling Ŝ with probability vector p = (p 1 , . . ., p n ), vector v ∈ R n ++ for which (f, Ŝ) ∼ ESO(v) 2: Initialization: choose x 0 ∈ R N , set z 0 = x 0 and θ 0 = min i p i 3: for k 0 do 4: Generate a random set of blocks S k ∼ Ŝ 6: for i ∈ S k do 8: end for 10: Note that in classical non-accelerated coordinate descent methods, only a single sequence of iterates is needed.This is indeed the case for our method as well, in the special case when the sampling Ŝ is uniform and θ 0 = E[| Ŝ|]/n, so that the three sequences are equal to each other.We now state a direct corollary of Theorem 3.1.
Corollary 4.3.For any optimal solution x * of (14), the output of Algorithm 5 for all k 1 satisfies: In particular, for ǫ > 0, if then In the special case when Ŝ is the serial uniform sampling, the three sequences {x k , y k , z k } k 0 coincide, and one can show that the following bound holds: Randomized coordinate descent with serial and importance sampling (in a form different from Algorithm 5) was considered by Nesterov [25].In the special case of serial uniform sampling when Algorithm 5 is the same as in [25], the following convergence rate was proved in [25]: where R(x 0 ) is a weighted level-set distance to the set of optimal points X * : Our result does not require the level sets of f to be bounded.

Special case 4: Accelerated parallel coordinate descent
To obtain the accelerated coordinate descent method, as a special case of Algorithm 2, we only need to let the sequence {θ k } k 0 satisfy (17).
Corollary 4.4.For any optimal solution x * of (14), the output of Algorithm 6 for all k 1 satisfies: In particular, for ǫ > 0, if When specialized to the serial uniform sampling (sampling Ŝ for which P( Ŝ = {i}) = 1/n for i ∈ [n]), the bound 36 simplifies to: An accelerated coordinate descent method for unconstrained minimization in the special case of serial uniform sampling was first proposed and analyzed by Nesterov [25], where the following bound was proved: Comparing ( 38) and (39), it is clear that we obtain a better bound.An accelerated coordinate descent method (APPROX) utilizing an arbitrary uniform sampling was studied by Fercoq and Richtárik [6].Algorithm 6, when restricted to a uniform sampling, is similar to this method.The main difference is in the value of θ 0 .Indeed, in [6], θ 0 is chosen to be E[| Ŝ|]/n, while our analysis allows θ 0 to be chosen as large as 1.This lead to larger stepsizes and a simplified and improved convergence bound.Each serial sampling Ŝ is uniquely characterized by the vector of probabilities p = (p 1 , . . ., p n ) where p i is defined by (8).Suppose that the function f has block-Lipschitz gradient with constants L 1 , . . ., L n : If Ŝ is a serial sampling, then which means that (f, Ŝ) ∼ ESO(L), where L = (L 1 , . . ., L n ) ∈ R n ++ .We can now find a sampling Ŝ for which the complexity bound (37) is minimized.This leads to the choice: The optimal serial sampling given by ( 41) is not very useful without (at least some) knowledge of x * , which is not known.However, note that the formula (41) confirms the intuition that blocks with larger L i and larger distance to the optimal block x i * − x i 0 i should be picked (and hence updated) more often.

Efficient implementation
As mentioned in Section 2.2, Algorithm 1 requires full-dimensional operations at each iteration unless θ k p i = 1 for all i ∈ [n] and k ∈ N. In this section we provide an equivalent form of Algorithm 1 which is suitable for efficient implementation under some additional assumptions on the computation of the gradient ∇f .

Equivalent form
Focusing on the iterates x k , y k , z k in Algorithm 1 only, the general algorithm can schematically be written as follows: Consider the change of variables from {x k , y k , z k } to {z k , g k } where and {α k } k 0 is a sequence defined by: Note that, in all the special cases presented in Section 4 and 6, either θ k < 1 for all k 1 or θ k = 1 for all k 1.The latter case does not require the full-dimensional operation y k = (1 − θ k )x k + θ k z k because the three sequences equal to each other.We thus only address the case when θ k < 1 for all k 1 which implies that α k = 0 for all k 1.Then from {z k , g k } and {α k } we can recover {x k , y k } as follows: +( 45) Moreover, g k+1 can be computed recursively as follows: where e ∈ R n is the vector of all ones.Therefore the updating scheme ( 42)-( 44) can thus be written in the form: Hence Algorithm 1 can be written in the following equivalent form.

Cost of a single iteration
In order to perform Step 7, it is important that we have access to without actually computing y k .In [6], the authors show that this is possible for problems (1) where f can be written as: for some matrix A ∈ R m×N .Let us write For i ∈ [n] and j ∈ [m] denote by A ji the ith block of the jth row vector of the matrix A, i.e., A ji = U ⊤ i A ⊤ e j .For each i ∈ [n], denote by I i the number of row containing a non-zero jth block, i.e., Then for f taking the form of (51) we have where by abuse of notation u j k and w j k denote respectively the jth element of the vectors u k and w k .With the knowledge of the vectors u k and w k , computing ∇ i f (y k ) requires O(|I i N i |) operations.Now in order to keep record of the vectors u k and w k , we use the following equality: and Since AU i is a matrix with |I i N i | nonzero elements, the updating schemes (52) and (53) require then operations.Note that this is also the total complexity of gradient computation ∇ i f (y k ) at kth iteration.Denote by nnz(A) the total number of nonzero blocks of the matrix A, i.e., nnz(A) Let us consider the special case when | Ŝ| = τ and N i = N/n for all i ∈ [n].In this case, the expected one iteration computational complexity is: To make it more direct to understand, let us consider the case when each block contains only one coordinate, i.e., N = n.Then the latter expected one iteration complexity becomes Hence, in this case the one iteration complexity in expectation of Algorithm 7 is of order O(τ ω) where ω is the average number of nonzero elements of the columns of A.Not considering the time spent on synchronization and handling read/write conflicts, the average processing time would be O(ω) if we use a parallel implementation with τ processors.

Proximal minimization
In this section we present and prove complexity results for ALPHA (Algorithm 1) as applied to the general problem (1) involving the proximal term.We leave the discussion concerning special cases to the appendix.

Complexity results
In the presence of the proximal term ψ, the same complexity bounds as those given in Theorems 3.1 and 3.2 hold for the output of Algorithm 1, with the exception that θ 0 is only allowed to be chosen between (0, min i p i ].We now state the formal complexity theorems, first in the nonaccelerated and then in the accelerated case.Theorem 6.1 (ALPHA -proximal & nonaccelerated).Let Ŝ be arbitrary proper sampling and v ∈ R n ++ be such that (f, Ŝ) ∼ ESO(v).Choose θ k = θ 0 ∈ (0, min i p i ] for all k 0. Then for any y ∈ R N , the iterates {x k } k 1 of Algorithm 1 satisfy: where Theorem 6.2 (ALPHA -proximal & accelerated).Let Ŝ be arbitrary proper sampling and v ∈ R n ++ be such that (f, Ŝ) ∼ ESO(v).Choose θ 0 ∈ (0, min i p i ] and define the sequence {θ k } k 0 by Then for any y ∈ R N such that C 0, the iterates {x k } k 1 of Algorithm 1 satisfy: where In the remainder of the section we will provide the complexity analysis.Our approach is similar to that presented in [6], but with many modifications required because we allow for an arbitrary sampling.We begin with some technical lemmas.

Technical lemmas
Lemma 6.1 shows that each individual block x i k of the variable x k is a convex combination of all the history blocks z i 0 , . . ., z i k .Note that due to the importance sampling, the combination coefficients γ i k,0 , . . ., γ i k,k is now block-dependent, in contrast with the block-independent coefficients proved in [6].Lemma 6.1.Let {x k , z k } k 0 be the iterates of Algorithm 1. Then for all k ∈ N and i ∈ [n] we have where for each i, the coefficients {γ i k,l } l=0,...,k are defined recursively by setting and for k 1, so that the following identity holds, Moreover, if θ 0 ∈ (0, min i p i ] and {θ k } k 0 is a decreasing positive sequence, then for all k ∈ N and i ∈ [n], the coefficients {γ i k,l } l=0,...,k are all positive and sum to 1. Proof.Fix any i ∈ [n].We proceed by induction on k.It is clear from x 0 = z 0 that γ i 0,0 = 1.Since x 1 = y 0 + θ 0 p −1 (z 1 − z 0 ) and y 0 = x 0 , we get that Assuming that (56) holds for some k 1, then Therefore the recursive equation ( 57) holds.The identity (58) can then be verified by direct substitution.Next we assume that θ 0 ∈ (0, min i p i ] and {θ k } k 0 is a decreasing positive sequence and show that the linear combination in (56) is a convex combination.Let k 1.Since θ k 1, we deduce from (57) that {γ i k+1,l } l=0,...,k−1 are positive if {γ i k,l } l=0,...,k−1 are positive.Moreover, Then using θ k 0, we conclude that γ i k+1,k 1. Besides, we have: We deduce from the above facts that the coefficients {γ i k+1,l } l=0,...,k+1 are all positive and sum to 1 if the same holds for {γ i k,l } l=0,...,k .Since θ 0 min i p i , we know that {γ i 1,0 , γ i 1,1 } are positive and sum to 1.It follows that the same property holds for all k ∈ N. Lemma 6.2.For k ∈ N and i ∈ Moreover, Proof.
The next result was previously stated and used in [6].

Recursion
For all k 0 define: We next prove an inequality similar to the one we established in the smooth case (20).
Lemma 6.4.Let Ŝ be an arbitrary proper sampling and v ∈ R n ++ be such that (f, Ŝ) ∼ ESO(v).Let {θ k } k 0 be arbitrary sequence of positive numbers in (0, 1] and fix y ∈ R N .Then for the sequence of iterates produced by Algorithm 1 and all k 0, the following recursion holds: We first write and then bound the expectation of Fk+1 as follows: Therefore, for all y ∈ R N , 6.4 Proof of Theorems 6.1 and 6.2 Using the same reasoning as that in the proof of Theorems 3.1 and 3.2, we analyze recursion of Lemma 6.4 and obtain: , ∀k 1. for all k 0, θ 0 ∈ (0, 1] and F0 F (y), then Finally, by Lemma 6.1, for the latter two choices of {θ k }, if in addition θ 0 ∈ (0, min i p i ], then for k 1, each block of the vector x k is a convex combination of the corresponding blocks of the vectors z 0 , . . ., z k .By the convexity of each function ψ i , we get: Hence, Theorem 6.1 and 6.2 hold by the fact that Fk F (x k ) for all k ∈ N and F0 = F (x 0 ).Note that the condition θ 0 ∈ (0, min i p i ] is only needed to prove (67).Thus it can be relaxed to θ 0 ∈ (0, 1] if ψ ≡ 0, in which case Fk = F (x k ) and Theorem 3.1 and 3.2 follows.

Conclusion
In this paper we propose a general randomized coordinate descent method which can be specialized to serial or parallel and accelerated or non-accelerated variants, with or without importance sampling.Based on the technical assumption which captures in a compact way certain smoothness properties of the function in a random subspace spanned by the sampled coordinates, we provide a unified complexity analysis which allows to derive as direct corollary the convergence results for the multiple variants of the general algorithm.We focused on the minimization of non-strongly convex function.Further study on a unified algorithm and complexity analysis for both strongly and non-strongly convex objective functions can be investigated.

A Special cases of ALPHA in the proximal setup
Extending the discussion presented in Section 4 which focused on the smooth case, we now present four special cases of Algorithm 1 for solving problem (1).
A.1 Special case 1: proximal gradient descent Specializing Algorithm 1 to Ŝ = [n] and θ k = 1 for all k 1, we obtain the classical proximal gradient descent algorithm.Note that in this case the three sequences {x k , y k , z k } k 0 reduce to one sequence.
Algorithm 8 Proximal Gradient Descent for solving 1 1: Parameters: vector v ∈ R n ++ such that (29) holds 2: Initialization: choose x 0 ∈ dom ψ 3: for k 0 do 4: end for 7: end for Corollary A.1.For any optimal solution x * of (1), the output of Algorithm 8 for all k 1 satisfies: In particular, for ǫ > 0, if The proof follows as a corollary of Theorem 6.1, with additional remark (30) which holds in this special case.The reader can refer to the proof of Corollary 4.1.
Corollary A.1 can be found in classical textbooks on convex optimization, see for example [25].
As discussed for the unconstrained case in Section 4.2, Algorithm 9 and Corollary A.2 can be attributed to Tseng [42].
Algorithm 10 Parallel Coordinate Descent Method (PCDM) [32] 1: Parameters: uniform proper random sampling Ŝ, positive vector v ∈ R n ++ such that (f, Ŝ) ∼ ESO(v) Generate S k ∼ Ŝ 5: x k+1 ← x k 6: for i ∈ S k do 7: x i k+1 = arg min end for 9: end for By applying Theorem 6.1 and using the same reasoning as in the proof of Corollary 4.1 we obtain the following new convergence result for PCDM.
Corollary A.3.For any optimal solution x * of (1), the output of Algorithm 10 for all k 1 satisfies: In particular, for 0 < ǫ < (1 A high-probability result involving the level-set distance was provided in [32] for PCDM (Algorithm 10).Although not explicitly stated in the paper, it is apparent from the proof that their approach yields the following rate: , it is clear that for sufficiently large k, our rate (71) is better than (73) .

A.4 Special case 4: APPROX with importance sampling
Finally, let {θ k } k 0 be chosen in accordance with (17).In this case, ALPHA (Algorithm 1) reduces to an accelerated coordinate descent method with arbitrary sampling Ŝ for solving the proximal minimization problem 1.
If Ŝ is a uniform sampling with τ = E[| Ŝ|] and θ 0 = τ /n, then we recover the convergence result established for APPROX [6,Theorem 3], as well as all the special cases that APPROX can recover, including the fast distributed coordinate descent method Hydra 2 [5].

k 2 9
: end for Note that only two sequences {x k , z k } k 0 are explicitly used in Algorithm 4. This is achieved by replacing y k in Algorithm 2 by (1 − θ k )x k + θ k z k .The following result follows directly from Theorem 3.2 by letting θ 0 = 1 and p i = 1 for all i ∈ [n].

4 :
Generate a random set of blocks S k ∼ Ŝ 5:

k 2 12 :
end forAPPROXis is a generalization of APPROX[6] from a uniform sampling to an arbitrary sampling: we recover APPROX if Ŝ is uniform with τ = E[| Ŝ|] and θ 0 = τ /n.