On optimal probabilities in stochastic coordinate descent methods

We propose and analyze a new parallel coordinate descent method—NSync—in which at each iteration a random subset of coordinates is updated, in parallel, allowing for the subsets to be chosen using an arbitrary probability law. This is the first method of this type. We derive convergence rates under a strong convexity assumption, and comment on how to assign probabilities to the sets to optimize the bound. The complexity and practical performance of the method can outperform its uniform variant by an order of magnitude. Surprisingly, the strategy of updating a single randomly selected coordinate per iteration—with optimal probabilities—may require less iterations, both in theory and practice, than the strategy of updating all coordinates at every iteration.


Introduction
In this work we consider the unconstrained minimization problem where φ is strongly convex and differentiable. We propose a new randomized algorithm for solving this problem-NSync (Nonuniform SYNchronous Coordinate descent)and analyze its iteration complexity. The main novelty of this paper is the algorithm itself. In particular, NSync is the first method which in each iteration updates a random subset of coordinates, allowing for an arbitrary probability law (sampling) to be used for this.

Algorithm 1 (NSync)
Input: Initial point x 0 ∈ R n , subset probabilities { p S } and stepsize parameters w 1 , . . . , w n > 0 for k = 0, 1, 2, . . . do Select a random set of coordinatesŜ ⊆ {1, . . . , n} such that Prob(Ŝ = S) = p S Update selected coordinates: At every iteration, a random setŜ is generated, independently from previous iterations, following the law Prob(Ŝ = S) = p S , S ⊆ [n], and then coordinates i ∈Ŝ are updated in parallel by moving in the direction of the negative partial derivative with stepsize 1/w i . By ∇ i φ(x) we mean ∇φ(x), e i , where e i ∈ R n is the ith unit coordinate vector. The updates are synchronized: no processor/thread is allowed to proceed before all updates are applied, generating the new iterate x k+1 . We study the complexity of NSync for arbitrary samplingŜ. In particular,Ŝ can be non-uniform in the sense that the probability that coordinate i is chosen, is allowed to vary with i.
NSync is the first randomized method in the literature which is capable of updating a subset of the coordinates without any restrictions, i.e., according to an arbitrary probability law, except for the necessary requirement that p i > 0 for all i. In particular, NSync is the first nonuniform parallel coordinate descent method.
In the time between the first online appearance of this work on arXiv (October 2013; arXiv:1310.3438), and the time this paper went to press, this work led to a number of extensions [3,7,[16][17][18]. All of these papers share the defining feature of NSync, namely, its ability to work with an arbitrary probability law defining the selection of the active coordinates in each iteration. These works also utilize the nonuniform ESO assumption introduced here (Assumption 1), as it appears to be key in the study of such methods.

Analysis
In this section we provide a complexity analysis of NSync.

Assumptions
Our analysis of NSync is based on two assumptions. The first assumption generalizes the ESO concept introduced in [21] and later used in [5,6,22,27,28] to nonuniform samplings. The second assumption requires that φ be strongly convex.
Notation For x, y, u ∈ R n we write x 2 u := i u i x 2 i , x, y u := n i=1 u i y i x i , x • y := (x 1 y 1 , . . . , x n y n ) and u −1 := (1/u 1 , . . . , 1/u n ). For S ⊆ [n] and h ∈ R n , let h [S] := i∈S h i e i .

Assumption 1 (Nonuniform ESO: Expected Separable Overapproximation)
Assume that p = ( p 1 , . . . , p n ) T > 0 and that for some positive vector w ∈ R n and all x, h ∈ R n , the following inequality holds: As soon as φ has a Lipschitz continuous gradient, then for every random samplinĝ S there exist positive weights w 1 , . . . , w n such that Assumption 1 holds. In this sense, the assumption is not restrictive. Inequalities of the type (2), in the uniform case ( p i = p j for all i, j), were studied in [6,21,22,27]. Motivated by the introduction of the nonuniform ESO assumption in this paper, and the development in Sect. 3 of our work, an entire paper was recently written, dedicated to the study of nonuniform ESO inequalities [16]. 1 We now turn to the second and final assumption.

Complexity
We can now establish a bound on the number of iterations sufficient for NSync to approximately solve (1) with high probability. We believe it is remarkable that the proof is very concise.

Theorem 3 Let Assumptions 1 and 2 be satisfied. Choose x
If {x k } are the random iterates generated by NSync, then Moreover, we have the lower bound Proof We first claim that φ is μ-strongly convex with respect to the norm 1 A clarifying comment answering a question raised by the reviewer: The authors of [16] give explicit formulas for w for which (2) holds, under an assumption that is slightly weaker than Lipschitz continuity of the gradient of φ. In particular, they study functions φ admitting the global quadratic upper bound for all x, h ∈ R n , where A ∈ R m×n . One of the consequence of their work is that the parameters w 1 , . . . , w n must necessarily satisfy the inequalities: w i ≥ A :i 2 , where A :i is the ith column of A. Moreover, as long as Prob(|Ŝ| ≤ τ ) = 1 for some τ , then (2) holds for w i = τ A :i 2 . However, this choice of parameters is rather conservative. The goal of [16] is to give explicit and tight formulas for w, where hopefully w i will be much smaller than τ A :i 2 , utilizing specific properties of the samplingŜ and data matrix A.
, and utilizing Assumption 1, we get Taking expectations in the last inequality and rearranging the terms, we obtain Using this, Markov inequality, and the definition of K , we finally get Let us now establish the last claim. Letting where the last equality follows since optimal p i is proportional to v i /w i .
Theorem 3 is generic in the sense that we do not say when Assumption 1 is satisfied and how should one go about choosing the stepsizes {w i } and probabilities { p S }. In the next section we address these issues. On the other hand, this abstract setting allowed us to write a brief complexity proof.
The quantity , defined in (4), can be interpreted as a condition number associated with the problem and our method. Hence, as we vary the distribution ofŜ, will vary.
It is clear intuitively that can be arbitrarily bad. Indeed, by choosing a samplingŜ which "nearly" ignores one or more of the coordinates (by setting p i ≈ 0 for some i), we should expect the number of iterations to grow as the method will necessarily be very slow in updating these coordinates.
In the light of this, inequality (6) is useful as it gives a useful expression for bounding from below.

Change of variables
Consider

Nonuniform samplings and ESO
In this section we consider a problem with standard assumptions and show that the (admittedly nonstandard) ESO assumption, Assumption 1, is satisfied.
Consider now problem (1) with φ of the form where v > 0. Note that Assumption 2 is satisfied. We further make the following two assumptions.

Assumption 4 (Smoothness)
Function f has Lipschitz gradient with respect to the coordinates, with positive constants L 1 , . . . , L n . That is, for all x ∈ R n and t ∈ R.

Assumption 5 (Partial separability) Function f has the form
where J is a finite collection of nonempty subsets of [n] and f J are differentiable convex functions such that f J depends on coordinates i ∈ J only. Let ω := max J |J |.
We say that f is separable of degree ω.
Uniform parallel coordinate descent methods for regularized problems with f of the above structure were analyzed in [21].
where A ∈ R m×n . Then L i = A :i 2 and where A :i is the ith column of A, A j: is the jth row of A and · is the standard L2 norm. Then ω is the maximum # of nonzeros in a row of A.
Nonuniform sampling Instead of considering the general case of arbitrary p S assigned to all subsets of [n], here we consider a special kind of sampling having two advantages: (i) sets can be generated easily, (ii) it leads to larger stepsizes 1/w i and hence improved convergence rate.
Fix τ ∈ [n] and c ≥ 1 and let S 1 , . . . , S c be a collection of (possibly overlapping) subsets of [n] such that Moreover, let q = (q 1 , . . . , q c ) > 0 be a probability vector. LetŜ j be τ -nice sampling from S j ; that is,Ŝ j picks subsets of S j having cardinality τ , uniformly at random. We assume these samplings are independent. Now,Ŝ is defined as follows: We first pick j ∈ {1, . . . , c} with probability q j , and then drawŜ j .
Note that we do not need to compute the quantities p S , S ⊆ [n], to execute NSync. In fact, it is much easier to implement the sampling via the two-tier procedure explained above. SamplingŜ is a nonuniform variant of the τ -nice sampling studied in [21], which here arises as a special case for c = 1.
Note that where δ i j = 1 if i ∈ S j , and 0 otherwise. In our next result we show that Assumption 1 is satisfied for f and the sampling described above.
Theorem 6 Let Assumptions 4 and 5 be satisfied, and letŜ be the sampling described above. Then Assumption 1 is satisfied with p given by (11) and any w = (w 1 , . . . , w n ) T for which for all i ∈ [n], where where the last inequality follows from the ESO for τ -nice samplings established in [21,Theorem 15]. The claim now follows by comparing the above expression and (2).

Optimal probabilities
Observe that the formula (12) can be used to design a sampling (characterized by the sets S j and probabilities q j ) that maximizes μ, which in view of Theorem 3 optimizes the convergence rate of the method.

Serial setting
Consider the serial version of NSync (Prob(|Ŝ| = 1) = 1). We can model this via c = n, with S i = {i} and p i = q i for all i ∈ [n]. In this case, using (11) and (12), we get w i = w * i = L i + v i . Minimizing in (4) over the probability vector p gives the optimal probabilities (we refer to this as the optimal serial method) and optimal complexity Note that the uniform sampling, defined by p i = 1/n for all i ∈ [n], leads to Note that this can be much larger than O S . We refer to NSync utilizing this sampling as the uniform serial method.
Moreover, the condition numbers L i /v i can not be improved via such a change of variables. Indeed, under the change of variables y = Diag(d)x, the gradient of f d (y) := f (Diag(d −1 )y) has coordinate Lipschitz constants L d i = L i /d 2 i , while the weights in (10)

Optimal serial method can be faster than the fully parallel method
To model the "fully parallel" setting (i.e., the variant of NSync updating all coordinates at every iteration), we can set c = 1 and τ = n, which yields Since ω ≤ n, it is clear that U S ≥ F P . However, for large enough ω it will be the case that F P ≥ O S , implying, surprisingly, that the optimal serial method can be faster than the fully parallel method.

Parallel setting
Fix τ and sets S j , j = 1, 2, . . . , c, and define Consider running NSync with stepsizes w i = θ(L i + v i ) (note that w i ≥ w * i , so we are fine). From (4), (11) and (12) we see that the complexity of NSync is determined by The probability vector q minimizing this quantity can be computed by solving a linear program with c + 1 variables (q 1 , . . . , q c , α), 2n linear inequality constraints and a single linear equality constraint: where b i ∈ R c , i ∈ [n], are given by  Fig. 1 Left optimal sampling (OS) is better than uniform sampling (US). Right nonuniform serial method (NS), updating a single coordinate in each iteration, can be faster than the fully parallel (FP) method, which updates all coordinates in each iteration

Experiments
We now conduct two preliminary small scale experiments to illustrate the theory; the results are depicted in Fig. 1. All experiments are with problems of the form (10) with f chosen as in Example 1.
In the left plot we chose A ∈ R 2×30 , γ = 1, v 1 = 0.05, v i = 1 for i = 1 and L i = 1 for all i. We compare the US method ( p i = 1/n, blue) with the OS method [ p i given by (13), red]. The dashed lines show 95 % confidence intervals (we run the methods 100 times, the line in the middle is the average behavior). While OS can be faster, it is sensitive to over/under-estimation of the constants L i , v i . In the right plot we show that a nonuniform serial (NS) method can be faster than the fully parallel (FP) variant (we have chosen m = 8, n = 10 and three values of ω). On the horizontal axis we display the number of epochs, where one epoch corresponds to updating n coordinates (for FP this is a single iteration, whereas for NS it corresponds to n iterations).