Large-scale randomized-coordinate descent methods with non-separable linear
constraints

Sashank J. Reddi ⇤ Ahmed Hefny⇤ Carlton Downey
Machine Learning Department
Carnegie Mellon University

Avinava Dubey

Suvrit Sra †Massachusetts Institute of Technology

Abstract
We develop randomized block coordinate descent (CD) methods for linearly constrained convex optimization. Unlike other large-scale CDmethods, we do not assume the constraints to be
separable, but allow them be coupled linearly.
To our knowledge, ours is the first CD methodthat allows linear coupling constraints, without
making the global iteration complexity have an
exponential dependence on the number of constraints. We present algorithms and theoretical analysis for four key (convex) scenarios: (i)smooth; (ii) smooth + separable nonsmooth; (iii)
asynchronous parallel; and (iv) stochastic. We
discuss some architectural details of our methods
and present preliminary results to illustrate thebehavior of our algorithms.

1

INTRODUCTION

Coordinate descent (CD) methods are conceptually among
the simplest schemes for unconstrained optimization—they
have been studied for a long time (see e.g., [1, 4, 28]), andare now enjoying greatly renewed interest. Their resurgence is rooted in successful applications in machine learning [15, 16], statistics [8, 17], and many other areas—see [31, 32, 35] and references therein for more examples.
A catalyst to the theoretical as well as practical success of CD methods has been randomization. (The ideaof randomized algorithms for optimization methods is ofcourse much older, see e.g., [29].) Indeed, generic nonrandomized CD has resisted complexity analysis, though
there is promising recent work [14, 34, 40]; remarkably
for randomized CD for smooth convex optimization, Nesterov [25, 26] presented an analysis of global iteration complexity. This work triggered several improvements, such as⇤

Indicates equal contribution.
A part of this work was performed when the author was at
Carnegie Mellon University†[32, 33], who simplified and extended the analysis to include separable nonsmooth terms. Randomization has also
been crucial to a host of other CD algorithms and analyses [5, 16, 20, 23, 30, 31, 33, 35–37].
Almost all of the aforementioned CD methods assume essentially unconstrained problems, which at best allow separable constraints. In contrast, we develop, analyze, andimplement randomized CD methods for the following composite objective convex problem with non-separable linear
constraints
minx F (x) := f (x) + h(x)s.t.

Ax = 0.

(1)Here, f : R ! R is assumed to be continuously differentiable and convex, while h : Rn ! R [ {1} islower semi-continuous, convex, coordinate-wise separable,
but not necessarily smooth; the linear constraints (LC) are
specified by a matrix A 2 Rm⇥n , for which m ⌧ n, anda certain structure (see §4) is assumed. The reader may
wonder whether one cannot simply rewrite Problem (1) in
the form of f + h (without additional constraints) using
suitable indicator functions. However, the resulting regularized problem then no longer fits the known efficient CD
frameworks [32], since the nonsmooth part is not blockseparable.
nProblem (1) subsumes the usual regularized optimization
problems pervasive in machine learning for the simplest(m = 0) case. In the presence of linear constraints (m >0), Problem (1) assumes a form used in the classic Alternating Direction Method of Multipliers (ADMM) [9, 10].
The principal difference between our approach and ADMM
is that the latter treats the entire variable x 2 Rn as asingle block, whereas we use the structure of A to split
x into b smaller blocks. Familiar special cases of Problem (1) include SVM (with bias) dual, fused Lasso andgroup Lasso [38], and linearly constrained least-squares regression [11, 19].
Recently, Necoara et al. [23] studied a special case of Problem (1) that sets h ⌘ 0 and assumes a single sum constraint. They presented a randomized CD method that starts
with a feasible solution and at each iteration updates a pair

Paper
[23]
[39]
[22]
[7]
[5]
[33]
Ours

of coordinates to ensure descent on the objective while
maintaining feasibility. This scheme is reminiscent of thewell-known SMO procedure for SVM optimization [27].
For smooth convex problems with n variables, Necoara
et al. [23] prove an O(1/✏) rate of convergence. More recently, in [22] considered a generalization to the generalcase Ax = 0 (assuming h is coordinatewise separable).
Unfortunately, the analysis in [22] yields an extremely pessimistic complexity result:
Theorem 1 ([22]). Consider Problem (1) with h being coordinatewise separable, and A 2 Rm⇥n with b blocks.
Then, the CD algorithm in [22] takes no more thanO(bm /✏) iterations to obtain a solution of ✏-accuracy.
This result is exponential in the number of constraints andtoo severe even for small-scale problems!
We present randomized CD methods, and prove that forimportant special cases (mainly h ⌘ 0 or A is a sum constraint) we can obtain global iteration complexity that does
not have an intractable dependence on either the number ofcoordinate blocks (b), or on the number of linear constraints(m). Previously, Tseng and Yun [39] also studied a linearly
coupled block-CD method based on the Gauss-Southwell
choice; however, their complexity analysis applies only tothe special m = 0 and m = 1 cases.
To our knowledge, ours is the first work on CD for problems with more than one (m > 1) linear constraints that
presents such results.
Contributions. In light of the above background, the primary contributions of this paper are as follows:
Convergence rate analysis of a randomized block-CD
method for the smooth case (h ⌘ 0) with m
1
general linear constraints.
A tighter convergence analysis for the composite function optimization (h 6= 0) than [22] in the case of sum
constraint.
An asynchronous CD algorithm for Problem (1).
A stochastic CD method with convergence analysis
for solving problems with a separable loss f (x) =PN(1/N ) i=1 fi (x).

Table 1 summarizes our contributions and compares it withexisting state-of-the-art coordinate descent methods. The
detailed proofs of all our theoretical claims are available in
the appendix.
Additional related work. As noted, CD methods have along history in optimization and they have gained tremendous recent interest. We cannot hope to do full justice to all
the related work, but refer the reader to [32, 33] and [20] formore thorough coverage. Classically, local linear convergence was analyzed in [21]. Global rates for randomized

LC
YES
YES
YES
⇥
⇥
⇥
YES

Prox
⇥
YES
YES
YES
`1
YES
YES

Parallel
⇥
⇥
⇥
YES
YES
YES
YES

Stochastic
⇥
⇥
⇥
⇥
YES
⇥
YES

Table 1: Summary comparison of our method with other CD

methods; LC denotes ‘linear constraints’; Prox signifies an extension using proximal operators (to handle h 6= 0).

block coordinate descent (BCD) were pioneered by Nesterov [25], and have since then been extended by various authors [2, 32, 33, 40]. The related family of GaussSeidel like analyses for ADMM have also recently gained
prominence [13]. A combination of randomized blockcoordinate ideas with Frank-Wolfe methods was recently
presented in [18], though algorithmically the Frank-Wolfe
approach is very different as it relies on non projection
based oracles.

2

PRELIMINARIESIn this section, we further explain our model and assumptions. We assume that the entire space Rn is decomposed>
into b blocks, i.e., x = [x>, x>where x 2 Rn ,1 ,···Pb ]ni
xi 2 R for all i 2 [b], and n = i ni . For any x 2 Rn ,we use xi to denote the ith block of x. We model communication constraints in our algorithms by viewing variables
as nodes in a connected graph G := (V, E). Specifically,
node i 2 V ⌘ [b] corresponds to variable xi , while an edge(i, j) 2 E ⇢ V ⇥V is present if nodes i and j can exchange
information. We use “pair” and “edge” interchangeably.
For a differentiable function f , we use fi1 ···ip andri1 ···ip f (x) (or rxi1 ···xip f (x)) to denote the restriction of
the function and its partial gradient to coordinate blocks(xi1 , · · · , xip ). For any matrix B with n columns, we use
Bi to denote the columns of B corresponding to xi and
Bij to denote the columns of B corresponding to xi andxj . We use U to denote the n ⇥ n identity matrix and hence
Ui is a matrix that places an ni dimensional vector into thecorresponding block of an n dimensional vector.
We make the following standard assumption on the partial
gradients of f .
Assumption 1. The function has block-coordinate Lipschitz continuous gradient, i.e.,
kri f (x)ri f (x + Ui h)k  Li khi k for all x 2 Rn , .

Assumption 1 is similar to the typical Lipschitz continuous
gradients assumed in first-order methods and it is necessary to ensure convergence of block-coordinate methods.

When functions fi and fj have Lipschitz continuous gradients with constants Li and Lj respectively, one can showthat the function fij has a Lipschitz continuous gradientwith Lij = Li + Lj [22; Lemma 1]. The following result
is standard.
Lemma 2. For any function g : Rn ! R with L-Lipschitz
continuous gradient rg, we haveyi + L2 kx

g(x)  g(y) + hrg(y), x

yk2 x, y 2 Rn .

Following [32, 39], we also make the following assumption
on the structure of h.
Assumption 2. The
P nonsmooth function h is block separable, i.e., h(x) = i hi (xi ).

This assumption is critical to composite optimization using CD methods. We also assume access to an oracle that
returns function values and partial gradients at any points
and iterates of the optimization algorithm.

3

ALGORITHM

We are now ready to present our randomized CD methods
for Problem (1) in various settings. We first study composite minimization (§3.1) and later look at asynchronous(§3.2) and stochastic (§3.3) variants. The main idea underlying our algorithms is to pick a random pair (i, j) 2 E ofvariables (blocks) at each iteration, and to update them ina manner which maintains feasibility and ensures progressin optimization.
3.1

Composite Minimization

We begin with the nonsmooth setting, where h 6⌘ 0. We
start with a feasible point x0 . Then, at each iteration wepick a random pair (i, j) 2 E of variables and minimize
the first-order Taylor expansion of the loss f around thecurrent iterate while maintaining feasibility. Formally, thisinvolves performing the update
Z(f, x, (i, j), ↵) := arg min f (x) + hrij f (x), dij i (2)Aij dij =0

+(2↵)

1

kdij k2 + h(x + Uij dij ),where ↵ > 0 is a stepsize parameter and dij is the update.
The right hand side of Equation (2) upper bounds f at x +Uij dij , as seen by using Assumption 1 and Lemma 2. Ifh(x) ⌘ 0, minimizing Equation (2) yields> +↵(Ai A>i + Aj Aj ) (Ai ri f (x) + Aj rj f (x))

di↵ri f (x) + A>idj

↵rj f (x) + A>j

(3)

Algorithm 1 presents the resulting method.
Note that since we start with a feasible point x0 and the
update dk satisfies Adk = 0, the iterate xk is always feasible. However, it can be shown that a necessary condition for Equation (2) to result in a non-zero update is thatAi and Aj span the same column space. If the constraints
are not block separable (i.e. for any partitioning of blocks
x1 , . . . , xb into two groups, there is a constraint that involves blocks from both groups), a typical way to satisfy
the aforementioned condition is to require Ai to be fullrow-rank for all i 2 [b]. This constraints the minimum
block size to be chosen in order to apply randomized CD.
Theorem 3 describes convergence of Algorithm 1 for thesmooth case (h ⌘ 0), while Theorem 6 considers the nonsmooth case under a suitable assumption on the structure of
the interdependency graph G—both results are presented
in Section 4.
1: x0 2 Rn such that Ax0 = 0
2: for k
0 do
3:Select a random edge (ik , jk ) 2 E with probability

p ik j k4:
dk
Uik jk Z(f, xk , (ik , jk ), ↵k /Lik jk )k+1
5:x
xk + dk
6:kk+1
7: end for
Algorithm 1: Composite Minimization with Linear Constraints

3.2

Asynchronous Parallel Algorithm for Smooth
Minimization

Although the algorithm described in the previous section
solves a simple subproblem at each iteration, it is inherently sequential. This can be a disadvantage when addressing large-scale problems. To overcome this concern, wedevelop an asynchronous parallel method that solves Problem (1) for the smooth case.
Our parallel algorithm is similar to Algorithm 1, except fora crucial difference: now we may have multiple processors, and each of these executes the loop 2–6 independently
without the need for coordination. This way, we can solve
subproblems (i.e., multiple pairs) simultaneously in parallel, and due to the asynchronous nature of our algorithm,
we can execute updates as they complete, without requiring any locking.
The critical issue, however, with implementing an asynchronous algorithm in the presence of non-separable constraints is ensuring feasibility throughout the course of thealgorithm. This requires the operation xi
xi + to be
executed in an atomic (i.e., sequentially consistent) fashion. Modern processors facilitate that without an additional
locking structure through the “compare-and-swap” instruction [30]. Since the updates use atomic increments andkeach
PTupdate ksatisfies Ad = 0, the net effect of T updates
is k=1 Ad = 0, which is feasible despite asynchronicityof the algorithm.
The next key issue is that of convergence. In an asynchronous setting, the updates are based on stale gradientsthat are computed using values of x read many iterations
earlier. But provided that gradient staleness is bounded,
we can establish a sublinear convergence rate of the asynchronous parallel algorithm (Theorem 4). More formally,
we assume that in iteration k, stale gradients are computed
based on xD(k) such that k
D(k)  ⌧ . The boundon staleness, denoted by ⌧ , captures the degree of parallelism in the method: such parameters are typical in asynchronous systems and provides a bound on the delay of the
updates [20].
Before concluding the discussion on our asynchronous algorithm, it is important to note the difficulty of extending our algorithm to nonsmooth problems. For example,
consider the case where h = IC (indicator function ofsome convex set). Although a pairwise update as suggested
above maintains feasibility with respect to the linear constraint Ax = 0, it may violate the feasibility of being in theconvex set C. This complication can be circumvented byusing a convex combination of the current iterate with theupdate, as this would retain overall feasibility. However,
it would complicate the convergence analysis. We plan toinvestigate this direction in future work.
3.3

Stochastic Minimization

An important subclass
PN of Problem (1) assumes separable
losses f (x) = N1 i=1 fi (x). This class arises naturally
in many machine learning applications where the loss separates over training examples. To take advantage of thisadded separability of f , we can derive a stochastic blockCD procedure.
Our key innovation here is the following: in addition to
randomly picking an edge (i, j), we also pick a function
randomly from {f1 , · · · , fN } and perform our update using this function. This choice substantially reduces the costof each iteration when N is large, since now the gradientcalculations involve only the randomly selected functionfi (i.e., we now use a stochastic-gradient). Pseudocode isgiven in Algorithm 2.
1: Choose x0 2 Rn such that Ax0 = 0.
2: for k
0 do
3:Select a random edge (ik , jk ) 2 E with probability

p ik j k

4:
Select random integer l 2 [N ]5:xk+1
xk + Uik jk Z(fl , xk , (ik , jk ), ↵k /Lik jk )6:
kk+1
7: end for
Algorithm 2: Stochastic Minimization with Linear Constraints

Notice that the per iteration cost of Algorithm 2 is lower
than Algorithm 1 by a factor of N . However, as we will
see later, this speedup comes at a price of slower convergence rate (Theorem 5). Moreover, to ensure convergence,

decaying step sizes {↵k }k 0 are generally chosen.

4

CONVERGENCE ANALYSIS

In this section, we outline convergence results for the algorithms described above. The proofs are somewhat technical, and hence left in the appendix due to lack of space;here we present only the key ideas.
For simplicity, we present our analysis for the following
reformulation of the main problem:minf (y, z) +subject to

Xb

y,z

i=1

Xb

i=1

h(yi , zi )(4)

yi = 0,where yi 2 Rny and zi 2 Rnz . Let y = [y1> · · · yb> ]> andz = [z1> · · · zb> ]> . We use x to denote the concatenatedvector [y > z > ]> and hence we assume (unless otherwise
mentioned) that the constraint matrix A is defined as follows✓◆ ✓ Pb◆y
i=1 yi
A
=
.
(5)
z
0
It is worth emphasizing that this analysis does not resultin any loss of generality. This is due to the fact that Problem (1) with a general constraint matrix Ã having full rowrank submatrices Ãi ’s can be rewritten in the form of Problem (4) by using the transformation specified in Section Eof the appendix. It is important to note that this reduction
is presented only for the ease of exposition. For our experiments, we directly solve the problem in Equation 2.
Let ⌘k = {(i0 , j0 ), . . . , (ik 1 , jk 1 )} denote the pairs selected up to iteration k 1. To simplify notation, assume(without loss of generality) that the Lipschitz constant forthe partial gradient ri f (x) and rij f (x) is L for all i 2 [n]and (i, j) 2 E.
Similar to [23], we introduce a Laplacian matrix L 2 Rb⇥bthat represents the communication graph G. Since we also
have unconstrained variables zi , we introduce a diagonal
matrix D 2 Rb⇥b .
⇢ P
⇢ pi
pir
i=j
i=j
r6
=
i
2L
L
Dij =
Lij =
pij
0
i 6= ji=6
j2L
We use K to denote the concatenation of the Laplacian L
and the diagonal matrix D. More formally,L ⌦ I ny0K=.
0D ⌦ I nzpThis matrix induces a norm kxkK = x> Kx on the feasible subspace, with a corresponding dual norms ✓◆L + ⌦ I ny
0⇤
kxkK = x>x0
D 1 ⌦ I nz

Let X ⇤ denote the set of optimal solutions and let x0 denotethe initial point. We define the following distance, whichquantifies how far the initial point is from the optimal, taking into account the graph layout and edge selection probabilities
R(x0 ) :=max kxmax

x:f (x)f (x0 ) x⇤ 2X ⇤

⇤

x⇤ kK

(6)

Note. Before delving into the details of the convergence results, we would like to draw the reader’s attention to the impact of the communication network G on convergence. In
general, the convergence results depend on R(x0 ), which
in turn depends on the Laplacian L of the graph G. As arule of thumb, the larger the connectivity of the graph, thesmaller the value of R(x0 ), and hence, faster the convergence.
4.1

convergence rates improve upon those of [22], since they
do not involve an exponential dependence of the form bm
on the number of constraints.
We now turn our attention towards the convergence analysis of our asynchronous algorithm under a consistent reading model [20]. In this context we would like to emphasize that while our theoretical analysis assumes consistentreads, we do not enforce this assumption in our experiments.
Theorem 4. Let ⇢ > 1 and ↵k =p ↵ be such that ↵ <2/(1 + ⌧ + ⌧ ⇢⌧ ) and ↵ < (⇢ 1)/( 2(⌧ + 2)(⇢⌧ +1 + ⇢)).
Let {xk }k 0 be the sequence generated by asynchronousalgorithm using step size ↵k and let f ⇤ denote the optimalvalue. Then, we have the following rate of convergence forthe expected values of the objective function
E[f (xk )]Convergence results for the smooth caseWe first consider the case when h = 0. Here the subproblem at k th iteration has a very simple update dik jk =kk
Uik dk Ujk dk where dk = ↵rik f (xk )).
2L (rjk f (x )We now prove that Algorithm 1 attains an O(1/k) convergence rate.
Theorem 3. Let ↵k = 1 for k0, and let {xk }k 0 be
the sequence generated by Algorithm 1; let f ⇤ denote theoptimal value. Then, we have the following rate of convergence:f⇤ 

E[f (xk )]2R2 (x0 )kwhere R(x0 ) is as defined in Equation 6.
Proof Sketch. We first prove that each iteration leads to descent in expectation. More formally, we get
Eik jk [f (xk+1 )|⌘k ]  f (xk )k >k1
2 rf (x ) Krf (x ).

The above step can be proved using Lemma 2. Let
E[f (xk )] f ⇤ . It can be proved that
1
k



1

1

k+1

2R2 (x0 )k =This follows from the fact thatf (xk+1 )f ⇤  kxkx⇤ k⇤K krf (xk )kK R(x0 )krf (xk )kK

8k

0

Telescoping the sum, we get the desired result.
Note that Theorem 3 is a strict generalization of the analysis in [23] and [22] due to: (i) the presence of unconstrained
variables z; and (ii) the presence of a non-decomposable
objective function. it is also worth emphasizing that ourf⇤ R2 (x0 )
µk

0
where
⇣ R(x ) is as⌘ defined in Equation 6 and µ =
↵2k
1+⌧
+⌧ ⇢⌧
1
.
2
↵k
2

Proof Sketch. For ease of exposition, we describe the casewhere the unconstrained variables z are absent. The analysis of case with z variables can be carried out in a similar manner. Let D(k) denote the iterate of the variables used in the k th iteration (the existence of D(k) follows
Let dk =⇣ from the consistent reading assumption).
⌘
↵k
D(k)
) ryik f (xD(k) ) and dkik jk = xk+1
2L ryjk f (xxk = Uik dk Ujk dk . Using Lemma 2 and the assumption that staleness in the variables is bounded by ⌧ , i.e.,k D(k)  ⌧ and definition of dkij , we can derive thefollowing bound:
✓
◆
1
1+⌧
k+1
k
E[f (x)]  E[f (x )] L
E[kdkik jk k2 ]
↵k
2
" ⌧
#
X
L
k t
2
+ E
kdik t jk t k .
2
t=1
In order to obtain an upper bound on the norms of dkik jk ,we prove thathi⇥⇤
E kdikk 11 jk 1 k2  ⇢E kdkik jk k2

This can proven using mathematical induction. Using theabove bound on kdkik jk k2 , we get
E[f (xk+1 )] 
E[f (xk )]

L

✓

1
↵k

1 + ⌧ + ⌧ ⇢⌧2◆

E[kdkik jk k2 ]This proves that the method is a descent method in expectation. Following similar analysis as Theorem 3, we get therequired result.

Note the dependence of convergence rate on the stalenessbound ⌧ . For larger values of ⌧ , the stepsize ↵k needs tobe decreased to ensure convergence, which in turn slows
down the convergence rate of the algorithm. Nevertheless,the convergence rate remains O(1/k).
The last smooth case we analyze is our stochastic algorithm.
ppTheorem 5. Let ↵i =0 in
0 L/(M i + 1) for i
Algorithm 2. Let {xk }k 0 be the sequence generated by
Algorithm 2 and let f ⇤ denote the optimal value. We denotex̄k = arg min0ik f (xk ). Then, we have the following
rate of convergence for the expected values of the objective:✓◆
1
E[f (x̄k )] f ⇤  O p
4
k
where

0
0 = f (x )f ⇤.

The convergence rate is O(1/k 1/4 ) as opposed to O(1/k)
of Theorem 3. On the other hand, the iteration complexity is lower by a factor of N ; this kind of tradeoff is typical in stochastic algorithms, where the slower rate is theprice we pay for a lower iteration complexity. Wepbelievethat the convergence rate can be improved to O(1/ k), therate generally observed in stochastic algorithms, by a more
careful analysis.
4.2

Nonsmooth caseWe finally state the convergence rate for the nonsmoothcase (h 6⌘ 0) in the case of a sum constraint. Similar to[22], we assume h is coordinatewise separable (i.e. we can
Pb Pthwrite h(x) =coori=1
j xij ), where xij is the jthdinate in the i block. For this analysis, we assume thatthe graph G is a clique 1 with uniform probability, i.e.,= pij = 2/b(b 1).
P
Theorem 6. Assume Ax = i Ai xi . Let {xk }k 0 be thesequence generated by Algorithm 1 and let F ⇤ denote theoptimal value. Assume that the graph G is a clique withuniform probability. Then we have the following:
E[F (xk )F ⇤] 

b2 LR2 (x0 )2202k + b LR 0(x ),where R(x ) is as defined in Equation 6.
0This convergence rate is a generalization of the convergence rate obtained in Necoara and Patrascu [22] for a single linear constraint (see Theorem 1 in [22]). It is also an
1
We believe our results also easily extend to the general casealong the lines of [31–33], using the concept of Expected Separable Overapproximation (ESO). Moreover, the assumption is nottotally impractical, e.g., in a multicore setting with a zero-sum
constraint (i.e. Ai = I), the clique-assumption introduces little
cost.

improvement of the rate obtained in Necoara and Patrascu[22] for general linear constraints (see Theorem 4 in [22])when applied to the special case of a sum constraint. Our
improvement comes in the form of a tractable constant,as opposed to the exponential dependence O(bm ) shownin [22].

5

APPLICATIONS

To gain a better understanding of our approach, we statesome applications of interest, while discussing details of
Algorithm 1 and Algorithm 2 for them. While there aremany applications of problem (1), due to lack of space weonly mention a few prominent ones here.
Support Vector Machines: The SVM dual (with bias
term) assumes the form (1); specifically,X
Xn
min 12↵i ↵j yi yj zi> zj
↵i
↵
i,j
i=1
X
s.t.
↵i yi = 0, 0  ↵i  C 8 i 2 [n].
(7)iHere, zi denotes the feature vector of the ith training example and yi 2 {1, 1}
laPdenotes the correspondingP
1>bel. By letting
f
(↵)
=
↵
↵
y
y
z
z
↵
and
i
j
i
j
j
i
i
i,j
i
2
P
h(↵) =i I(0  ↵i  C) and A = [y1 , . . . , yn ] thisproblem can be written in form of Problem (1). Using Algorithm 1 for SVM involves solving a sub-problem similar
to one used in SMO in the scalar case (i.e., ↵i 2 R) andcan be solved in linear time in the block case (see [3]).
Generalized Lasso: The objective is to solve the following
optimization problem.
1
2 kY

min

X k22 + kD k1

where Y 2 RN denotes the output, X 2 RN ⇥n is the input
and D 2 Rq⇥n represents a specified penalty matrix. This
problem can also be seen as a specific case of Problem (1)by introducing an auxiliary variable t andPslack variablesu, v. Then, f ( , t) = 12 kY
X k22 + i ti , h(u, v) =I(u 0)+I(v 0) and, t D u = 0 and t+D v = 0
are the linear constraints. To solve this problem, we can use
either Algorithm 1 or Algorithm 2. In general, optimizationof convex functions on a structured convex polytope can besolved in a similar manner.
Unconstrained Separable Optimization: Another interesting application is for unconstrained
separable optimizaPtion. For any problem minx i fi (x)—a form generally
encountered across machine learning—can be rewritten using variable-splitting as min{xi =x,8i2[N ]} fi (xi ). Solving
the problem in distributed environment requires considerable synchronization (for the consensus constraint), which
can slow down the algorithm significantly. However, thedual of the problem isXXNminfi⇤ ( i ) s.t
i = 0.
ii=1

where fi⇤ is the Fenchel conjugate of fi . This reformulation perfectly fits our framework and can be solved in
an asynchronous manner using the procedure described in
Section 3.2.

Objective Value

1000

Other interesting application include constrained leastsquare problem, multi-agent planning problems, resource
allocation—see [22, 23] and references therein for more
examples.

6

In this section, we present our empirical results. In particular, we examine the behavior of random coordinate descentalgorithms analyzed in this paper under different communication constraints and concurrency conditions. 2
6.1 Effect of Communication Constraints
Our first set of experiments test the affect of the connectivity of the graph on the convergence rate. In particular, recall that the convergence analysis established in Theorem 3
depends on the Laplacian of the communication graph. In
this experiment we demonstrate how communication constraints affect convergence in practice. We experiment withthe following graph topologies of graph G: Ring, Clique,Star + Ring (i.e., the union of edges of a star and a ring)and Tree + Ring. On each layout we run the sequential
Algorithm 1 on the following quadratic problem
C
s.t.

XN

i=1
XN

i=1

kxi

(i mod 10)1k2

Ai xi = 0,(8)Note the decomposable structure of the problem. For thisexperiment, we use N = 1000 and xi 2 R50 . We have
10 constraints whose coefficients are randomly generated
from U [0, 1] and we choose C such that the objective evaluates to 1000 when x = 0.
The results for Algorithm 1 on each topology for 10000 iterations are shown in Figure 1. The results clearly showthat better connectivity implies better convergence rate.
Note that while the clique topology has significantly better
convergence than other topologies, acceptable long-term
performance can be achieved by much sparser topologies
such as Star + Ring and Tree + Ring.
Having a sparse communication graph is important to lower
the cost of a distributed system. Furthermore, it is worth
mentioning that the sparsity of the communication graphis also important in a multicore setting; since Algorithm 1
2

900
800
700
600
0

2000

4000

6000

8000

IterationhFigure 1: Objective value vs. number of iterations for different

EXPERIMENTS

min

Tree + Ring
Star + Ring
Ring
CliqueAll experiments were conducted on a Google Compute Engine virtual machine of type “n1-highcpu-16”, which comprises
16 virtual CPUs and 14.4 GB of memory. For more details, pleaserefer to https://cloud.google.com/compute/docs/
machine-types#highcpu.

graph topologies. Note that larger the connectivity of the graph,faster is the convergence.
> +requires computing (Ai A>
i + Aj Aj ) for each communicating pair of nodes (i, j). Our analysis shows that thiscomputation takes a significant portion of the running timeand hence it is essential to minimize the number of variable
pairs that are allowed to be updated.

6.2

Concurrency and Synchronization

As seen earlier, compared to Tree + Ring, Star + Ring isa low diameter layout (diameter = 2). Hence, in a sequential setting, it indeed results in a faster convergence. However, Star + Ring requires a node to be connected to allother nodes. This high-degree node could be a contention
point in a parallel setting. We test the performance of ourasynchronous algorithm in this setting. To assess how theperformance would be affected with such contention andhow asynchronous updates would increase performance,
we conduct another experiment on the synthetic problem(8) but on a larger scale (N = 10000, xi 2 R100 , 100constraints).
Our concurrent update follows a master/slave scheme.
Each thread performs a loop where in each iteration it elects
a master i and slave j and then applies the following sequence of actions:
1. Obtain the information required for the update fromthe master (i.e., information for calculating the gradients used for solving the subproblem).
2. Send the master information to the slave, update theslave variable and get back the information needed to
update the master.
3. Update the master based (only) on the information obtained from steps 1 and 2.
We emphasize that the master is not allowed to read its own
state at step 3 except to apply an increment, which is computed based on steps 1 and 2. This ensures that the master’s increment is consistent with that of the slave, even if
one or both of them was being concurrently overwritten byanother thread. More details on the implementation can befound in [12].

Following [30], we use spinlocks instead of mutex locks
to implement locking. Spinlocks are preferred over mutex locks when the resource is locked for a short period oftime, which is the case in our algorithm. For each locking
mechanism, we vary the number of threads from 1 to 15.
We stop when f0 ft > 0.99(f0 f ⇤ ), where f ⇤ is computed beforehand up to three significant digits. Similar to[30], we add artificial delay to steps 1 and 2 in the update
scheme to model complicated gradient calculations and/or
network latency in a distributed setting.
Figure 2 shows the speedup for Tree + Ring and Star + Ring
layouts. The figure clearly shows that a fully synchronous
method suffers from contention in the Star + Ring topology
whereas asynchronous method does not suffer from thisproblem and hence, achieves higher speedups. Although
the Tree + Ring layouts achieves higher speedup than Star +Ring, the latter topology results in much less running time(⇠ 67 seconds vs 91 seconds using 15 threads).
6.3

Practical Case Study: Parallel Training of Linear
SVM

In this section, we explore the effect of parallelism on randomized CD for training a linear SVM based on the dual
formulation stated in (7). Necoara et. al. [22] have shown
that, in terms of CPU time, a sequential randomized CD
outperforms coordinate descent using Gauss-Southwell selection rule. It was also observed that randomized CD outperforms LIBSVM [6] for large datasets while maintaining
reasonable performance for small datasets.
In this experiment we use a clique layout. For SVM training in a multicore setting, using a clique layout does not introduce additional cost compared to a more sparse layout.
To maintain the box constraint, we use the double-locking
scheme described in Section 6.2 for updating a pair of dual
variables.
One advantage of coordinate descent algorithms is that they
do not require the storage of the Gram matrix; instead they
can compute its elements on the fly. That comes, however,at the expense of CPU time. Similar to [22], to speed up
gradient computations without increasing memory require-

Speedup10

Optimal
Double lockSingle lock
Lock free

5

0
0

15

Speedup

Given this update scheme, we experiment with three levelsof synchronization: (a) Double Locking: Locks the masterand the slave through the entire update. Because the objective function is decomposable, a more conservative locking(e.g. locking all nodes) is not needed. (b) Single Locking:Locks the master during steps 1 and 3 (the master is unlocked during step 2 and locks the slave during step 2). (c)Lock-free: No locks are used. Master and slave variables
are updated through atomic increments similar to Hogwild!
method.

15

10

5

P

10

15

10

15

Optimal
Double lockSingle lock
Lock free

5

0
0

5

P

Figure 2: Speedup for Tree + Ring (top) and Star + Ring(bottom) topologies and different levels of synchronization.
Note for Star + Ring topology, speedup of asynchronousalgorithm is significantly higher than that of synchronous
version.

ments, we maintain the primal weight vector of the linear
SVM and use it to compute gradients. Basically, if weincrement ↵i by i and ↵j by j , then we increment theweight vector by i yi xi + j yj xj . This increment is accomplished using atomic additions. However, this implies
that all threads will be concurrently updating the primal
weight vector. Similar to [30], we require these updates tobe sparse with small overlap between non-zero coordinatesin order to ensure convergence. In other words, we require
training examples to have sparse features with small overlap between non-zero features.
We report speedups on two datasets used in [22].3 Table 2 provides a description of both the datasets. For
each dataset, we train the SVM model until f0ft >0.9999(f0 f ⇤ ), where f ⇤ is the objective reported in [22].
In Figure 3, we report speedup for both the datasets. The
figure shows that parallelism indeed increases the performance of randomized CD training of linear SVM.

3
Datasets can be downloaded from http://www.csie.
ntu.edu.tw/˜cjlin/libsvmtools/datasets.

Dataset# of
instances# offeatures

a7a
w8a

16100
49749

122
300

Avg # ofnon-zerofeatures14
12

Table 2: Datasets used for linear SVM Speedup experiment

3.5

Speedup

3
2.5
2
1.5
1
0.5
0

5

P

10

15

10

15

5

Speedup

4
3
2
1
0

5

P

Figure 3: Speedup for linear SVM training on a7a (top)and w8a datasets.

7

DISCUSSION AND FUTURE WORK

We presented randomized coordinate descent methods forsolving convex optimization problems with linear constraints that couple the variables. Moreover, we also presented composite objective, stochastic, and asynchronous
versions of our basic method and provided their convergence analysis. We demonstrated the empirical performance of the algorithms. The experimental results of asynchronous algorithm look very promising.
There are interesting open problems for our problem
in consideration: First, we would like to obtain highprobability results not just in expectation; another interesting direction is to extend the asynchronous algorithm tothe non-smooth setting. Finally, while we obtain O(1/k)
for general convex functions, obtaining an accelerated
O(1/k 2 ) rate is a natural question.

Acknowledgments
SS is partly supported by NSF grant: IIS-1409802. We
thanks the anonymous reviewers for the helpful comments.

References[1] A. Auslender. Optimisation Méthodes Numériques.
Masson, 1976.
[2] A. Beck and L. Tetruashvili. On the convergence ofblock coordinate descent type methods. SIAM Journal on Optimization, 23(4):2037–2060, 2013. doi:
10.1137/120887679.
[3] P. Berman, N. Kovoor, and P. M. Pardalos. A lineartime algorithm for the least-distance problem. Technical report, Pennsylvania State University, Departmentof Computer Science, 1992.
[4] D. P. Bertsekas. Nonlinear Programming. Athena
Scientific, Belmont, MA, second edition, 1999.
[5] J. Bradley, A. Kyrola, D. Bickson, and C. Guestrin.
Parallel coordinate descent for L1-regularized loss
minimization. In L. Getoor and T. Scheffer, editors,
Proceedings of the 28th International Conference on
Machine Learning, pages 321–328. Omnipress, 2011.
[6] C.-C. Chang and C.-J. Lin. Libsvm: A library for support vector machines. ACM Trans. Intell. Syst. Technol., 2(3):27:1–27:27, May 2011. ISSN 2157-6904.
doi: 10.1145/1961189.1961199. URL http://
doi.acm.org/10.1145/1961189.1961199.
[7] O. Fercoq and P. Richtárik. Accelerated, parallel andproximal coordinate descent. CoRR, abs/1312.5799,
2013.
[8] J. Friedman, T. Hastie, H. Höfling, R. Tibshirani, et al.
Pathwise coordinate optimization. The Annals of Applied Statistics, 1(2):302–332, 2007.
[9] D. Gabay and B. Mercier. A dual algorithm for the solution of nonlinear variational problems via finite element approximation. Computers & Mathematics with
Applications, 2(1):17–40, 1976. doi: http://dx.doi.
org/10.1016/0898-1221(76)90003-1. URL http://www.sciencedirect.com/science/
article/pii/0898122176900031.
[10] R. Glowinski and A. Marrocco. Sur l’approximation,
par éléments finis d’ordre un, et la résolution, parpénalisation-dualité d’une classe de problèmes
de dirichlet non linéares.
Revue Françaised’Automatique,
Informatique,et RechercheOpérationelle, 1975.
[11] G. H. Golub and C. F. Van Loan. Matrix Computations. John Hopkins University Press, Baltimore,
MD, 3rd edition, 1996.
[12] A. Hefny, S. Reddi, and S. Sra. Coordinate descent algorithms with coupling constraints: Lessons learned.
In NIPS Workshop on Software Engineering For Machine Learning, 2014.

[13] M. Hong and Z.-Q. Luo. On the linear convergence of
the alternating direction method of multipliers. arXiv
preprint arXiv:1208.3922, 2012.
[14] M. Hong, X. Wang, M. Razaviyayn, and Z.-Q. Luo.
Iteration Complexity Analysis of Block Coordinate
Descent Methods. arXiv:1310.6957, 2013.
[15] C. J. Hsieh and I. S. Dhillon. Fast coordinate descentmethods with variable selection for non-negative
matrix factorization. In Proceedings of the 17th
ACM SIGKDD International Conference on Knowledge Discovery and Data Mining(KDD), pages 1064–
1072, August 2011.
[16] C. J. Hsieh, K. W. Chang, C. J. Lin, S. S. Keerthi, and
S. Sundararajan. A dual coordinate descent method
for large-scale linear SVM. In W. Cohen, A. McCallum, and S. Roweis, editors, ICML, pages 408–415.
ACM, 2008.
[17] C.-J. Hsieh, M. A. Sustik, I. S. Dhillon, and P. D.
Ravikumar. Sparse inverse covariance matrix estimation using quadratic approximation. In NIPS, pages
2330–2338, 2011.
[18] S. Lacoste-Julien, M. Jaggi, M. Schmidt, and
P. Pletscher.
Block-coordinate frank-wolfe optimization for structural svms.
arXiv preprint
arXiv:1207.4747, 2012.
[19] C. L. Lawson and R. J. Hanson. Solving Least
Squares Problems. Prentice–Hall, Englewood Cliffs,
NJ, 1974. Reissued with a survey on recent developments by SIAM, Philadelphia, 1995.
[20] J. Liu, S. J. Wright, C. Ré, V. Bittorf, and S. Sridhar. An asynchronous parallel stochastic coordinate
descent algorithm. arXiv:1311.1873, 2013.
[21] Z.-Q. Luo and P. Tseng. On the convergence of thecoordinate descent method for convex differentiable
minimization. Journal of Optimization Theory andApplications, 72(1):7–35, 1992.
[22] I. Necoara and A. Patrascu. A random coordinate
descent algorithm for optimization problems withcomposite objective function and linear coupled constraints. Comp. Opt. and Appl., 57(2):307–337, 2014.
[23] I. Necoara, Y. Nesterov, and F. Glineur. A random coordinate descent method on large optimization problems with linear constraints. Technical report, Technical Report, University Politehnica Bucharest, 2011,2011.
[24] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro.
Robust stochastic approximation approach to stochastic programming. SIAM J. on Optimization, 19(4):1574–1609, Jan. 2009. ISSN 1052-6234. doi: 10.
1137/070704277.
[25] Y. Nesterov. Efficiency of coordinate descent methods
on huge-scale optimization problems. Core discussion papers, Université catholique de Louvain, Center
for Operations Research and Econometrics (CORE),

2010.
[26] Y. Nesterov. Efficiency of coordinate descent methods
on huge-scale optimization problems. SIAM Journal
on Optimization, 22(2):341–362, 2012.
[27] J. C. Platt. Sequential minimal optimization: A fastalgorithm for training support vector machines. Technical report, ADVANCES IN KERNEL METHODS- SUPPORT VECTOR LEARNING, 1998.
[28] B. T. Polyak. Introduction to Optimization. Optimization Software Inc., 1987. Nov 2010 revision.
[29] L. A. Rastrigin. Statisticheskie Metody Poiska Ekstremuma (Statistical Extremum Seeking Methods).
Nauka, Moscow, 1968.
[30] B. Recht, C. Re, S. J. Wright, and F. Niu. Hogwild: A
Lock-Free Approach to Parallelizing Stochastic Gradient Descent. In J. Shawe-Taylor, R. S. Zemel, P. L.
Bartlett, F. C. N. Pereira, and K. Q. Weinberger, editors, NIPS, pages 693–701, 2011.
[31] P. Richtárik and M. Takáč. Distributed coordinate
descent method for learning with big data. ArXiv eprints, Oct. 2013.
[32] P. Richtárik and M. Takáč.
Iteration complexity of randomized block-coordinate descentmethods for minimizing a composite function.
arXiv:1107.2848v1, July 2011.
[33] P. Richtárik and M. Takáč.
Parallel coordinate descent methods for big data optimization.
arXiv:1212.0873v1, Dec 2012.
[34] A. Saha and A. Tewari. On the nonasymptotic convergence of cyclic coordinate descent methods. SIAM
Journal on Optimization, 23(1):576–601, 2013.
[35] S. Shalev-Shwartz and T. Zhang. Stochastic dual
coordinate ascent methods for regularized loss minimization. JMLR, 14, 2013.
[36] S. Shalev-Shwartz and T. Zhang.
Accelerated mini-batch stochastic dual coordinate ascent.
arXiv:1305.2581v1, 2013.
[37] R. Tappenden, P. Richtárik, and J. Gondzio. Inexact
coordinate descent: complexity and preconditioning.
arXiv:1304.5530, 2013.
[38] R. Tibshirani, M. Saunders, S. Rosset, J. Zhu, andK. Knight. Sparsity and smoothness via the fusedlasso. Journal of the Royal Statistical Society: Series
B (Statistical Methodology), 67(1):91–108, 2005.
[39] P. Tseng and S. Yun. A block-coordinate gradient descent method for linearly constrained nonsmooth separable optimization. Journal of Optimization Theory
and Applications, 2009.
[40] P.-W. Wang and C.-J. Lin. Iteration complexity offeasible descent methods for convex optimization.
JMLR, 15:1523–1548, 2014.