arXiv:1601.00863v3 [math.OC] 14 Aug 2016

Coordinate friendly structures, algorithms and
applications∗
Zhimin Peng, Tianyu Wu, Yangyang Xu, Ming Yan, and
Wotao Yin

This paper focuses on coordinate update methods, which are useful
for solving problems involving large or high-dimensional datasets.
They decompose a problem into simple subproblems, where each
updates one, or a small block of, variables while fixing others. These
methods can deal with linear and nonlinear mappings, smooth and
nonsmooth functions, as well as convex and nonconvex problems.
In addition, they are easy to parallelize.
The great performance of coordinate update methods depends
on solving simple subproblems. To derive simple subproblems for
several new classes of applications, this paper systematically studies coordinate friendly operators that perform low-cost coordinate
updates.
Based on the discovered coordinate friendly operators, as well
as operator splitting techniques, we obtain new coordinate update
algorithms for a variety of problems in machine learning, image
processing, as well as sub-areas of optimization. Several problems
are treated with coordinate update for the first time in history.
The obtained algorithms are scalable to large instances through
parallel and even asynchronous computing. We present numerical
examples to illustrate how effective these algorithms are.
Keywords and phrases: coordinate update, fixed point, operator splitting, primal-dual splitting, parallel, asynchronous.

∗

This work is supported by NSF Grants DMS-1317602 and ECCS-1462398.

1

2

1. Introduction
This paper studies coordinate update methods, which reduce a large problem to smaller subproblems and are useful for solving large-sized problems.
These methods handle both linear and nonlinear maps, smooth and nonsmooth functions, and convex and nonconvex problems. The common special
examples of these methods are the Jacobian and Gauss-Seidel algorithms for
solving a linear system of equations, and they are also commonly used for
solving differential equations (e.g., domain decomposition) and optimization
problems (e.g., coordinate descent).
After coordinate update methods were initially introduced in each topic
area, their evolution had been slow until recently, when data-driven applications (e.g., in signal processing, image processing, and statistical and
machine learning) impose strong demand for scalable numerical solutions;
consequently, numerical methods of small footprints, including coordinate
update methods, become increasingly popular. These methods are generally
applicable to many problems involving large or high-dimensional datasets.
Coordinate update methods generate simple subproblems that update
one variable, or a small block of variables, while fixing others. The variables
can be updated in the cyclic, random, or greedy orders, which can be selected to adapt to the problem. The subproblems that perform coordinate
updates also have different forms. Coordinate updates can be applied either
sequentially on a single thread or concurrently on multiple threads, or even
in an asynchronous parallel fashion. They have been demonstrated to give
rise to very powerful and scalable algorithms.
Clearly, the strong performance of coordinate update methods relies on
solving simple subproblems. The cost of each subproblem must be proportional to how many coordinates it updates. When there are totally m coordinates, the cost of updating one coordinate should not exceed the average per-coordinate cost of the full update (made to all the coordinates
at once). Otherwise, coordinate update is not computationally worthy. For
example, let f : Rm → Rbe a C 2 function, and consider the Newton update
−1
xk+1 ← xk − ∇2 f (xk ) ∇f (xk ). Since updating each xi (keeping others fixed) still requires forming the Hessian matrix ∇2 f (x) (at least O(m2 )
operations) and factorizing it (O(m3 ) operations), there is little to save in
computation compared to updating all the components of x at once; hence,
the Netwon’s method is generally not amenable to coordinate update.
The recent coordinate-update literature has introduced new algorithms.
However, they are primarily applied to a few, albeit important, classes of
problems that arise in machine learning. For many complicated problems,

Coordinate friendly structures, algorithms, and applications

3

it remains open whether simple subproblems can be obtained. We provide
positive answers to several new classes of applications and introduce their
coordinate update algorithms. Therefore, the focus of this paper is to build
a set of tools for deriving simple subproblems and extending coordinate
updates to new territories of applications.
We will frame each application into an equivalent fixed-point problem
(1)

x=Tx

by specifying the operator T : H → H, where x = (x1 , . . . , xm ) ∈ H, and
H = H1 × · · · × Hm is a Hilbert space. In many cases, the operator T itself
represents an iteration:
(2)

xk+1 = T xk

such that the limit of the sequence {xk } exists and is a fixed point of T ,
which is also a solution to the application or from which a solution to the
application can be obtained. We call the scheme (2) a full update, as opposed
to updating one xi at a time. The scheme (2) has a number of interesting
special cases including methods of gradient descent, gradient projection,
proximal gradient, operator splitting, and many others.
We study the structures of T that make the following coordinate update
algorithm computationally worthy
(3)

xk+1
= xki − ηk (xk − T xk )i ,
i

where ηk is a step size and i ∈ [m] := {1, . . . , m} is arbitrary. Specifically, the
1
, or lower, of that of performing (2). We
cost of performing (3) is roughly m
call such T a Coordinate Friendly (CF) operator, which we will formally
define.
This paper will explore a variety of CF operators. Single CF operators
include linear maps, projections to certain simple sets, proximal maps and
gradients of (nearly) separable functions, as well as gradients of sparsely
supported functions. There are many more composite CF operators, which
are built from single CF and non-CF operators under a set of rules. The fact
that some of these operators are CF is not obvious.
These CF operators let us derive powerful coordinate update algorithms
for a variety of applications including, but not limited to, linear and secondorder cone programming, variational image processing, support vector machine, empirical risk minimization, portfolio optimization, distributed computing, and nonnegative matrix factorization. For each application, we present

4
an algorithm in the form of (2) so that its coordinate update (3) is efficient.
In this way we obtain new coordinate update algorithms for these applications, some of which are treated with coordinate update for the first time.
The developed coordinate update algorithms are easy to parallelize. In
addition, the work in this paper gives rise to parallel and asynchronous extensions to existing algorithms including the Alternating Direction Method
of Multipliers (ADMM), primal-dual splitting algorithms, and others.
The paper is organized as follows. §1.1 reviews the existing frameworks
of coordinate update algorithms. §2 defines the CF operator and discusses
different classes of CF operators. §3 introduces a set of rules to obtain composite CF operators and applies the results to operator splitting methods.
§4 is dedicated to primal-dual splitting methods with CF operators, where
existing ones are reviewed and a new one is introduced. Applying the results
of previous sections, §5 obtains novel coordinate update algorithms for a variety of applications, some of which have been tested with their numerical
results presented in §6.
Throughout this paper, all functions f, g, h are proper closed convex and
can take the extended value ∞, and all sets X, Y, Z are nonempty closed
convex. The indicator function ιX (x) returns 0 if x ∈ X, and ∞ elsewhere.
For a positive integer m, we let [m] := {1, . . . , m}.
1.1. Coordinate Update Algorithmic Frameworks
This subsection reviews the sequential and parallel algorithmic frameworks
for coordinate updates, as well as the relevant literature.
The general framework of coordinate update is
1. set k ← 0 and initialize x0 ∈ H = H1 × · · · × Hm
2. while not converged do
3.
select an index ik ∈ [m];
4.
update xk+1
for i = ik while keeping xk+1
= xki , ∀ i 6= ik ;
i
i
5. k ← k + 1;
Next we review the index rules and the methods to update xi .
1.1.1. Sequential Update. In this framework, there is a sequence of
coordinate indices i1 , i2 , . . . chosen according to one of the following rules:
cyclic, cyclic permutation, random, and greedy rules. At iteration k, only
the ik th coordinate is updated:
(
xk+1
= xki − ηk (xk − T xk )i , i = ik ,
i
k+1
xi = xki ,
for all i 6= ik .

Coordinate friendly structures, algorithms, and applications

5

Sequential updates have been applied to many problems such as the GaussSeidel iteration for solving a linear system of equations, alternating projection [75, 4] for finding a point in the intersection of two sets, ADMM [31, 30]
for solving monotropic programs, and Douglas-Rachford Splitting (DRS) [26]
for finding a zero to the sum of two operators.
In optimization, coordinate descent algorithms, at each iteration, minimize the function f (x1 , . . . , xm ) by fixing all but one variable xi . Let
xi− := (x1 , . . . , xi−1 ),

xi+ = (xi+1 , . . . , xm )

collect all but the ith coordinate of x. Coordinate descent solves one of the
following subproblems:
(4a)

(T xk )i = arg min f (xki− , xi , xki+ ),
xi

(4b)
(4c)
(4d)

1
kxi − xki k2 ,
2ηk
xi
1
(T xk )i = arg min h∇i f (xk ), xi i +
kxi − xki k2 ,
2ηk
xi
1
(T xk )i = arg min h∇i f diff (xk ), xi i + fiprox (xi ) +
kxi − xki k2 ,
2ηk
xi

(T xk )i = arg min f (xki− , xi , xki+ ) +

which are called direct update, proximal update, gradient update, and proxgradient update, respectively. The last update applies to the function
f (x) = f diff (x) +

m
X

fiprox (xi ),

i=1
prox
where f diff is differentiable and each
is proximable (its proximal map
 fi
takes O dim(xi ) polylog(dim(xi )) operations to compute).
Sequential-update literature. Coordinate descent algorithms date
back to the 1950s [35], when the cyclic index rule was used. Its convergence
has been established under a variety of cases, for both convex and nonconvex
objective functions; see [77, 85, 57, 33, 45, 70, 32, 72, 58, 8, 36, 78]. Proximal
updates are studied in [32, 1] and developed into prox-gradient updates
in [74, 73, 13] and mixed updates in [81].
The random index rule first appeared in [48] and then [61, 44]. Recently, [82, 80] compared the convergence speeds of cyclic and stochastic
update-orders. The gradient update has been relaxed to stochastic gradient
update for large-scale problems in [21, 83].

6
The greedy index rule leads to fewer iterations but is often impractical
since it requires a lot of effort to calculate scores for all the coordinates.
However, there are cases where calculating the scores is inexpensive [11, 41,
79] and the save in the total number of iterations significantly outweighs the
extra calculation [74, 25, 55, 49].
A simple example. We present the coordinate update algorithms under different index rules for solving a simple least squares problem:
1
minimize f (x) := kAx − bk2 ,
x
2
where A ∈ Rp×m and b ∈ Rp are Gaussian random. Our goal is to numerically demonstrate the advantages of coordinate updates over the full update
of gradient descent:
xk+1 = xk − ηk A> (Axk − b).
The four tested index rules are: cyclic, cyclic permutation, random, and
greedy under the Gauss-Southwell1 rule. Note that because this example is
very special, the comparisons of different index rules are far from conclusive.
In the full update, the step size ηk is set to the theoretical upper bound
2
2 , where kAk2 denotes the matrix operator norm and equals the largest
kAk2
singular value of A. For each coordinate update to xi , the step size ηk is set
to (A>1A)ii . All of the full and coordinate updates have the same per-epoch
complexity, so we plot the objective errors in Figure 1.
1.1.2. Parallel Update. As one of their main advantages, coordinate
update algorithms are easy to parallelize. In this subsection, we discuss both
synchronous (sync) and asynchronous (async) parallel updates.
Sync-parallel (Jacobi) update specifies a sequence of index subsets
I1 , I2 , . . . ⊆ [m], and at each iteration k, the coordinates in Ik are updated
in parallel by multiple agents:
(
xk+1
= xki − ηk (xk − T xk )i , i ∈ Ik ,
i
xk+1
= xki ,
i 6∈ Ik .
i
Synchronization across all agents ensures that all xi in Ik are updated and
also written to the memory before the next iteration starts. Note that, if
1

it selects ik = arg maxi k∇i f (xk )k.

Coordinate friendly structures, algorithms, and applications

7

5

10

full update
greedy
random
cyclic
cyclic permutation

0

fk − f∗

10

−5

10

−10

10

−15

10

0

20

40
60
Epochs

80

100

Figure 1: Gradient descent: the coordinate updates are faster than the full
update since the former can take larger steps at each step.

Ik = [m] for all k, then all the coordinates are updated and, thus, each
iteration reduces to the full update: xk+1 = xk − ηk (xk − T xk ).
Async-parallel update. In this setting, a set of agents still perform
parallel updates, but synchronization is eliminated or weakened. Hence, each
agent continuously applies (5), which reads x from and writes xi back to
the shared memory (or through communicating with other agents without
shared memory):
(

xk+1
= xki − ηk (I − T )xk−dk i , i = ik ,
i
(5)
xk+1
= xki ,
for all i 6= ik .
i
Unlike before, k increases whenever any agent completes an update.
The lack of synchronization often results in computation with out-ofdate information. During the computation of the kth update, other agents
make dk updates to x in the shared memory; when the kth update is written,
its input is already dk iterations out of date. This number is referred to as the
asynchronous delay. In (5), the agent reads xk−dk and commits the update
to xkik . Here we have assumed consistent reading, i.e., xk−dk lying in the set
{xj }kj=1 . This requires implementing a memory lock. Removing the lock can
lead to inconsistent reading, which still has convergence guarantees; see [54,
Section 1.2] for more details.
Synchronization across all agents means that all agents will wait for the
last (slowest) agent to complete. Async-parallel updates eliminate such idle
time, spread out memory access and communication, and thus often run

8
Agent 1
Agent 2

idle

Agent 2

idle

Agent 3
k: 0

Agent 1

idle

Agent 3

idle

1

(a) sync-parallel computing

2

k: 0

1

2

3 4 5 6 7 8 9 10

(b) async-parallel computing

Figure 2: Sync-parallel computing (left) versus async-parallel computing
(right). On the left, all the agents must wait at idle (white boxes) until
the slowest agent has finished.

much faster. However, async-parallel is more difficult to analyze because of
the asynchronous delay.
Parallel-update literature. Async-parallel methods can be traced
back to [17] for systems of linear equations. For function minimization, [12]
introduced an async-parallel gradient projection method. Convergence rates
are obtained in [69]. Recently, [14, 62] developed parallel randomized methods.
For fixed-point problems, async-parallel methods date back to [3] in
1978. In the pre-2010 methods [2, 10, 6, 27] and the review [29], each agent
updates its own subset of coordinates. Convergence is established under the
P -contraction condition and its variants [10]. Papers [6, 7] show convergence
for async-parallel iterations with simultaneous reading and writing to the
same set of components. Unbounded but stochastic delays are considered
in [67].
Recently, random coordinate selection appeared in [19] for fixed-point
problems. The works [47, 59, 43, 42, 37] introduced async-parallel stochastic
methods for function minimization. For fixed-point problems, [54] introduced
async-parallel stochastic methods, as well as several applications.
1.2. Contributions of This Paper
The paper systematically discusses the CF properties found in both single and composite operators underlying many interesting applications. We
introduce approaches to recognize CF operators and develop coordinateupdate algorithms based on them. We provide a variety of applications to
illustrate our approaches. In particular, we obtain new coordinate-update
algorithms for image deblurring, portfolio optimization, second-order cone
programming, as well as matrix decomposition. Our analysis also provides

Coordinate friendly structures, algorithms, and applications

9

guidance to the implementation of coordinate-update algorithms by specifying how to compute certain operators and maintain certain quantities in
memory. We also provide numerical results to illustrate the efficiency of the
proposed coordinate update algorithms.
This paper does not focus on the convergence perspective of coordinate
update algorithms, though a convergence proof is provided in the appendix
for a new primal-dual coordinate update algorithm. In general, in fixed-point
algorithms, the iterate convergence is ensured by the monotonic decrease of
the distance between the iterates and the solution set, while in minimization problems, the objective value convergence is ensured by the monotonic
decrease of a certain energy function. The reader is referred to the existing
literature for details.
The structural properties of operators discussed in this paper are irrelevant to the convergence-related properties such as nonexpansiveness (for an
operator) or convexity (for a set or function). Hence, the algorithms developed can be still applied to nonconvex problems.

2. Coordinate Friendly Operators
2.1. Notation
For convenience, we do not distinguish a coordinate from a block of coordinates throughout this paper. We assume our variable x consists of m
coordinates:
x = (x1 , . . . , xm ) ∈ H := H1 × · · · × Hm

and xi ∈ Hi , i = 1, . . . , m.

For simplicity, we assume that H1 , . . . , Hm are finite-dimensional real Hilbert
spaces, though most results hold for general Hilbert spaces. A function maps
from H to R, the set of real numbers, and an operator maps from H to G,
where the definition of G depends on the context.
Our discussion often involves two points x, x+ ∈ H that differ over one
coordinate: there exists an index i ∈ [m] and a point δ ∈ H supported on
Hi , such that
(6)

x+ = x + δ.

+
Note that x+
j = xj for all j 6= i. Hence, x = (x1 , . . . , xi + δi , . . . , xm ).

Definition 1 (number of operations). We let M [a 7→ b] denote the number
of basic operations that it takes to compute the quantity b from the input a.

10
For example, M [x 7→ (T x)i ] denotes the number of operations to compute the ith component of T x given x. We explore the possibility to compute
(T x)i with much fewer operations than what is needed to first compute T x
and then take its ith component.
2.2. Single Coordinate Friendly Operators
This subsection studies a few classes of CF operators and then formally
defines the CF operator. We motivate the first class through an example.
In the example below, we let Ai,: and A:,j be the ith row and jth column
>
of a matrix A, respectively. Let A> be the transpose of A and A>
i,: be (A )i,: ,
i.e., the ith row of the transpose of A.
Example 1 (least squares I). Consider the least squares problem
(7)

1
minimize f (x) := kAx − bk2 ,
x
2

where A ∈ Rp×m and b ∈ Rp . In this example, assume that m = Θ(p),
namely, m and p are of the same order. We compare the full update of
gradient descent to its coordinate update.2 The full update is referred to as
the iteration xk+1 = T xk where T is given by
(8)

T x := x − η∇f (x) = x − ηA> Ax + ηA> b.

Assuming that A> A and A> b are already computed, we have M [x 7→ T x] =
O(m2 ). The coordinate update at the kth iteration performs
= (T xk )ik = xkik − η∇ik f (xk ),
xk+1
ik
and xk+1
= xkj , ∀j 6= ik , where ik is some selected coordinate.
j

Since for all i, ∇i f (xk ) = A> (Ax − b) i = (A> A)i,: · x − (A> b)i , we
1
have M [x 7→ (T x)i ] = O(m) and thus M [x 7→ (T x)i ] = O( m
M [x 7→ T x]).
Therefore, the coordinate gradient descent is computationally worthy.
The operator T in the above example is a special Type-I CF operator.
Definition 2 (Type-I CF). For an operator T : H → H, let M [x 7→ (T x)i ]
be the number of operations for computing the ith coordinate of T x given x
2
Although gradient descent is seldom used to solve least squares, it often appears
as a part in first-order algorithms for problems involving a least squares term.

Coordinate friendly structures, algorithms, and applications

11

and M [x 7→ T x] the number of operations for computing T x given x. We
say T is Type-I CF (denoted as F1 ) if for any x ∈ H and i ∈ [m], it holds


1
M [x 7→ (T x)i ] = O
M [x 7→ T x] .
m
Example 2 (least squares II). We can implement the coordinate update in
Example 1 in a different manner by maintaining the result T xk in the memory. This approach works when m = Θ(p) or p  m. The full update (8) is
unchanged. At each coordinate update, from the maintained quantity T xk , we
immediately obtain xk+1
= (T xk )ik . But we need to update T xk to T xk+1 .
ik
k+1
k
Since x
and x differ only over the coordinate ik , this update can be
computed as
k
>
T xk+1 = T xk + xk+1 − xk − η(xk+1
ik − xik )(A A):,ik ,

which is a scalar-vector multiplication followed by vector addition, taking
only O(m) operations. Computing T xk+1 from scratch involves a matrixvector multiplication, taking O(M [x 7→ T (x)]) = O(m2 ) operations. Therefore,

h
i
i
h
1
k
k k+1
k+1
k+1
k+1
=O
M x
7→ T x
.
M {x , T x , x } 7→ T x
m
The operator T in the above example is a special Type-II CF operator.
Definition 3 (Type-II CF). An operator T is called Type-II
CF (denoted

as F2 ) if, for any i, x and x+ := x1 , . . . , (T x)i , . . . , xm , the following holds
(9)



M {x, T x, x+ } 7→ T x+ = O






1
M x+ 7→ T x+ .
m

The next example illustrates an efficient coordinate update by maintaining certain quantity other than T x.
Example 3 (least squares III). For the case p  m, we should avoid
pre-computing the relative large matrix A> A, and it is cheaper to compute
A> (Ax) than (A> A)x. Therefore, we change the implementations of both
the full and coordinate updates in Example 1. In particular, the full update
xk+1 = T xk = xk − η∇f (xk ) = xk − ηA> (Axk − b),


pre-multiplies xk by A and then A> . Hence, M xk 7→ T (xk ) = O(mp).

12
We change the coordinate update to maintain the intermediate quantity
In the first step, the coordinate update computes

Axk .

(T xk )ik = xkik − η(A> (Axk ) − A> b)ik ,
k
k+1
by pre-multiplying Axk by A>
ik ,: . Then, the second step updates Ax to Ax
k+1
k
k
by adding (xik − xik )A:,ik to Ax . Both steps take O(p) operations, so

h
i
h
i
1
k
k
k
k
k+1
k+1
.
M x 7→ T x
M {x , Ax } 7→ {x , Ax } = O(p) = O
m

Combining Type-I and Type-II CF operators with the last example, we
arrive at the following CF definition.
Definition 4 (CF operator). We say that anoperator T : H → H is CF if,
for any i, x and x+ := x1 , . . . , (T x)i , . . . , xm , the following holds




1
+
+
(10)
M {x, M(x)} 7→ {x , M(x )} = O
M [x 7→ T x] ,
m
where M(x) is some quantity maintained in the memory to facilitate each
coordinate update and refreshed to M(x+ ). M(x) can be empty, i.e., except
x, no other varying quantity is maintained.
The left-hand side of (10) measures the cost of performing one coordinate update (including the cost of updating M(x) to M(x+ )) while the
right-hand side measures the average per-coordinate cost of updating all
the coordinates together. When (10) holds, T is amenable to coordinate
updates.
By definition, a Type-I CF operator T is CF without maintaining any
quantity, i.e., M(x) = ∅.
A Type-II CF operator T satisfies (10) with M(x) = T x, so it is also
CF. Indeed, given any x and i, we can compute x+ by immediately letting
+
x+
i = (T x)i (at O(1) cost) and keeping xj = xj , ∀j 6= i; then, by (9), we
+
update T x to T x at a low cost. Formally, letting M(x) = T x,


M {x, M(x)} 7→ {x+ , M(x+ )}




≤M {x, T x} 7→ x+ + M {x, T x, x+ } 7→ T x+


 +

1
(9)
+
= O(1) + O
M x 7→ T x
m


1
=O
M [x 7→ T x] .
m

Coordinate friendly structures, algorithms, and applications

13

In general, the set of CF operators is much larger than the union of
Type-I and Type-II CF operators.
Another important subclass of CF operators are operators T : H → H
where (T x)i only depends on one, or a few, entries among x1 , . . . , xm . Based
on how many input coordinates they depend on, we partition them into three
subclasses.
Definition 5 (separable operator). Consider T := {T | T : H → H}. We
have the partition T = C1 ∪ C2 ∪ C3 , where
• separable operator: T ∈ C1 if, for any index i, there exists Ti : Hi → Hi
such that (T x)i = Ti xi , that is, (T x)i only depends on xi .
• nearly-separable operator: T ∈ C2 if, for any index i, there exists Ti
and index set Ii such that (T x)i = Ti ({xj }j∈Ii ) with |Ii |  m, that is,
each (T x)i depends on a few coordinates of x.
• non-separable operator: C3 := T \ (C1 ∪ C2 ). If T ∈ C3 , there exists
some i such that (T x)i depends on many coordinates of x.
Throughout the paper, we assume the coordinate update of a (nearly-)
separable operator costs roughly the same for all coordinates. Under this
assumption, separable operators are both Type-I CF and Type-II CF, and
nearly-separable operators are Type-I CF.3
2.3. Examples of CF Operators
In this subsection, we give examples of CF operators arising in different
areas including linear algebra, optimization, and machine learning.
Example 4 ((block) diagonal matrix). Consider the diagonal matrix


a1,1

A=
0

0
..



m×m
.
∈R

.
am,m

Clearly T : x 7→ Ax is separable.
3

Not all nearly-separable operators are Type-II CF. Indeed, consider a sparse
matrix A ∈ Rm×m whose non-zero entries are only located in the last column. Let
T x = Ax and x+ = x+δm . As x+ and x differ over the last entry, T x+ = T x+(x+
m−
xm )A:,m takes m operations. Therefore, we have M [{x, T x, x+ } 7→ T x+ ] = O(m).
+
+
Since T x+ = x+
m A:,m takes m operations, we also have M [x 7→ T x ] = O(m).
Therefore, (9) is violated, and there is no benefit from maintaining T x.

14
Example 5 (gradient and proximal maps of a separable function). Consider
a separable function
m
X
f (x) =
fi (xi ).
i=1

Then, both ∇f and proxγf are separable, in particular,
(∇f (x))i = ∇fi (xi )

and

(proxγf (x))i = proxγfi (xi ).

Here, proxγf (x) (γ > 0) is the proximal operator that we define in Definition 10 in Appendix A.
Example 6 (projection to box constraints). Consider the “box” set B :=
{x : ai ≤ xi ≤ bi , i ∈ [m]} ⊂ Rm . Then, the projection operator projB is
separable. Indeed,


projB (x) i = max(bi , min(ai , xi )).
Example 7 (sparse matrices). If every row of the matrix A ∈ Rm×m is
sparse, T : x 7→ Ax is nearly-separable.
Examples of sparse matrices arise from various finite difference schemes
for differential equations, problems defined on sparse graphs. When most
pairs of a set of random variables are conditionally independent, their inverse
covariance matrix is sparse.
Example 8 (sum of sparsely supported functions). Let E be a class of
index sets and every e ∈ E be a small subset of [m], |e|  m. In addition
#{e : i ∈ e}  #{e} for all i ∈ [m]. Let xe := (xi )i∈e , and
f (x) =

X

fe (xe ).

e∈E

The gradient map ∇f is nearly-separable.
An application of this example arises in wireless communication over a
graph of m nodes. Let each xi be the spectrum assignment to node i, each e be
a neighborhood of nodes, and each fe be a utility function. The input of fe is
xe since the utility depends on the spectra assignments in the neighborhood.
In machine learning, if each observation only involves a few features,
then each function of the optimization objective will depend on a small
number of components of x. This is the case when graphical models are
used [64, 9].

Coordinate friendly structures, algorithms, and applications

15

Example 9 (squared hinge loss function). Consider for a, x ∈ Rm ,
f (x) :=

2
1
max(0, 1 − βa> x) ,
2

which is known as the squared hinge loss function. Consider the operator
(11)

T x := ∇f (x) = −β max(0, 1 − βa> x)a.

Let us maintain M(x) = a> x. For arbitrary x and i, let
>
x+
i := (T x)i = −β max(0, 1 − βa x)ai
+
>
and x+
j := xj , ∀j 6= i. Then, computing xi from x and a x takes O(1) (as
>
a> x is maintained), and computing a> x+ from x+
i − xi and a x costs O(1).
Formally, we have
h
i
M {x, a> x} 7→ {x+ , a> x+ }
i
i
h
h
> +
−
x
}
→
7
a
x
≤M {x, a> x} 7→ x+ + M {a> x, x+
i
i

=O(1) + O(1) = O(1).
On the other hand, M [x 7→ T x] = O(m). Therefore, (10) holds, and T
defined in (11) is CF.

3. Composite Coordinate Friendly Operators
Compositions of two or more operators arise in algorithms for problems that
have composite functions, as well as algorithms that are derived from operator splitting methods. To update the variable xk to xk+1 , two or more
operators are sequentially applied, and therefore the structures of all operators determine whether the update is CF. This is where CF structures
become less trivial but more interesting. This section studies composite CF
operators. The exposition leads to the recovery of existing algorithms, as
well as powerful new algorithms.
3.1. Combinations of Operators
We start by an example with numerous applications. It is a generalization
of Example 9.

16
Example 10 (scalar map pre-composing affine function). Let aj ∈ Rm , bj ∈
R, and φj : R → R be differentiable functions, j ∈ [p]. Let
f (x) =

p
X

φj (a>
j x + bj ).

j=1

Assume that evaluating φ0j costs O(1) for each j. Then, ∇f is CF. Indeed,
let
T1 y := A> y, T2 y := [φ01 (y1 ); . . . ; φ0p (yp )], T3 x := Ax + b,
>
>
p×m and b = [b ; b ; . . . ; b ] ∈ Rp×1 . Then we
where A = [a>
1 2
p
1 ; a2 ; . . . ; ap ] ∈ R
have ∇f (x) = T1 ◦ T2 ◦ T3 x. For any x and i ∈ [m], let x+
i = ∇i f (x) and
x+
=
x
,
∀j
=
6
i,
and
let
M(x)
:=
T
x.
We
can
first
compute
T2 ◦ T3 x from
j
3
j
+
T3 x for O(p) operations, then compute ∇i f (x) and thus x from {x, T2 ◦T3 x}
for O(p) operations, and finally update the maintained T3 x to T3 x+ from
{x, x+ , T3 x} for another O(p) operations. Formally,


M {x, T3 x} 7→ {x+ , T3 x+ }




≤M [T3 x 7→ T2 ◦ T3 x] + M {x, T2 ◦ T3 x} 7→ x+ + M {x, T3 x, x+ } 7→ {T3 x+ }
=O(p) + O(p) + O(p) = O(p).

Since M [x 7→ ∇f (x)] = O(pm), therefore ∇f = T1 ◦ T2 ◦ T3 is CF.
If p = m, T1 , T2 , T3 all map from Rm to Rm . Then, it is easy to check
that T1 is Type-I CF, T2 is separable, and T3 is Type-II CF. The last one is
crucial since not maintaining T3 x would disqualify T from CF. Indeed, to
obtain (T x)i , we must multiply A>
i to all the entries of T2 ◦ T3 x, which in
turn needs all the entries of T3 x, computing which from scratch would cost
O(pm).
There are general rules to preserve Type-I and Type-II CF. For example, T1 ◦ T2 is still Type-I CF, and T2 ◦ T3 is still CF, but there are counter
examples where T2 ◦ T3 can be neither Type-I nor Type-II CF. Such properties are important for developing efficient coordinate update algorithms for
complicated problems; we will formalize them in the following.
The operators T2 and T3 in the above example are prototypes of cheap
and easy-to-maintain operators from H to G that arise in operator compositions.
Definition 6 (cheap operator). For a composite operator T = T1 ◦ · · · ◦ Tp ,
an operator Ti : H → G is cheap if M [x 7→ Ti x] is less than or equal to the
number of remaining coordinate-update operations, in order of magnitude.

Coordinate friendly structures, algorithms, and applications
Case
1
2
3
4

T1
C1
C2
C2
C3

∈
(separable)
(nearly-sep.)
(non-sep.)

T2 ∈
C1 , C2 , C3
C1 , C3
C2
C1 ∪ C 2 ∪ C 3

17

(T1 ◦ T2 ) ∈
C1 , C2 , C3 , respectively
C2 , C3 , resp.
C2 or C3 , case by case
C3

Table 1: T1 ◦ T2 inherits the weaker separability property from those of T1
and T2 .
Case
5
6
7
8
9

T1 ∈
C1 ∪ C 2
F, F2
F1
cheap
F1

T2 ∈
F, F1
C1
F2
F2
cheap

(T1 ◦ T2 ) ∈
F, F1 , resp.
F, F2 , resp.
F
F
F1

Example
Examples 11 and 13
Example 10
Example 12
Example 13
Examples 10 and 13

Table 2: Summary of how T1 ◦ T2 inherits CF properties from those of T1
and T2 .
Definition 7 (easy-to-maintain operator). For a composite operator T =
T1 ◦ · · · ◦ Tp , the operator Tp : H → G is easy-to-maintain, if for any x, i, x+
satisfying (6), M [{x, Tp x, x+ } 7→ Tp x+ ] is less than or equal to the number
of remaining coordinate-update operations, in order of magnitude, or belongs
1
+
+
to O( dim
G M [x 7→ T x ]).
The splitting schemes in §3.2 below will be based on T1 + T2 or T1 ◦ T2 ,
as well as a sequence of such combinations. If T1 and T2 are both CF, T1 + T2
remains CF, but T1 ◦ T2 is not necessarily so. This subsection discusses how
T1 ◦ T2 inherits the properties from T1 and T2 . Our results are summarized
in Tables 1 and 2 and explained in detail below.
The combination T1 ◦ T2 generally inherits the weaker property from T1
and T2 .
The separability (C1 ) property is preserved by composition. If T1 , . . . , Tn
are separable, then T1 ◦ · · · ◦ Tn is separable. However, combining nearlyseparable (C2 ) operators may not yield a nearly-separable operator since
composition introduces more dependence among the input entries. Therefore, composition of nearly-separable operators can be either nearly-separable
or non-separable.
Next, we discuss how T1 ◦ T2 inherits the CF properties from T1 and
T2 . For simplicity, we only use matrix-vector multiplication as examples to
illustrate the ideas; more interesting examples will be given later.

18
• If T1 is separable or nearly-separable (C1 ∪ C2 ), then as long as T2 is
CF (F), T1 ◦ T2 remains CF. In addition, if T2 is Type-I CF (F1 ), so
is T1 ◦ T2 .
Example 11. Let A ∈ Rm×m be sparse and B ∈ Rm×m dense. Then
T1 x = Ax is nearly-separable and T2 x = Bx is Type-I CF4 . For any
i, let Ii index the set of nonzeros on the ith row of A. We first compute (Bx)Ii , which costs O(|Ii |m), and then ai,Ii (Bx)Ii , which costs
O(|Ii |), where ai,Ii is formed by the nonzero entries on the ith row of
A. Assume O(|Ii |) = O(1), ∀i. We have, from the above discussion,
that M [x 7→ (T1 ◦ T2 x)i ] = O(m), while M [x 7→ T1 ◦ T2 x] = O(m2 ).
Hence, T1 ◦ T2 is Type-I CF.
• Assume that T2 is separable (C1 ). It is easy to see that if T1 is CF (F),
then T1 ◦ T2 remains CF. In addition if T1 is Type-II CF (F2 ), so is
T1 ◦ T2 ; see Example 10.
Note that, if T2 is nearly-separable, we do not always have CF properties for T1 ◦ T2 . This is because T2 x and T2 x+ can be totally different
(so updating T2 x is expensive) even if x and x+ only differ over one
coordinate; see the footnote 3 on Page 13.
• Assume that T1 is Type-I CF (F1 ). If T2 is Type-II CF (F2 ), then
T1 ◦ T2 is CF (F).
Example 12. Let A, B ∈ Rm×m be dense. Then T1 x = Ax is Type-I
CF and T2 x = Bx Type-II CF (by maintaining Bx; see Example 2).
For any x and i, let x+ satisfy (6). Maintaining T2 x, we can compute
(T1 ◦ T2 x)j for O(m) operations for any j and update T2 x+ for O(m)
operations. On the other hand, computing T1 ◦ T2 x+ without maintaining T2 x takes O(m2 ) operations.
• Assume that one of T1 and T2 is cheap. If T2 is cheap, then as long as
T1 is Type-I CF (F1 ), T1 ◦ T2 is Type-I CF. If T1 is cheap, then as long
as T2 is Type-II CF (F2 ), T1 ◦ T2 is CF (F); see Example 13.
We will see more examples of the above cases in the rest of the paper.
3.2. Operator Splitting Schemes
We will apply our discussions above to operator splitting and obtain new
algorithms. But first, we review several major operator splitting schemes and
4

For this example, one can of course pre-compute AB and claim that (T1 ◦ T2 )
is Type-I CF. Our arguments keep A and B separate and only use the nearlyseparability of T1 and Type-I CF property of T2 , so our result holds for any such
composition even when T1 and T2 are nonlinear.

Coordinate friendly structures, algorithms, and applications

19

discuss their CF properties. We will encounter important concepts such as
(maximum) monotonicity and cocoercivity, which are given in Appendix A.
For a monotone operator A, the resolvent operator JA and the reflectiveresolvent operator RA are also defined there, in (67) and (68), respectively.
Consider the following problem: given three operators A, B, C, possibly
set-valued,
(12)

find x ∈ H

such that

0 ∈ Ax + Bx + Cx,

where “+” is the Minkowski sum. This is a high-level abstraction of many
problems or their optimality conditions. The study began in the 1960s, followed by a large number of algorithms and applications over the last fifty
years. Next, we review a few basic methods for solving (12).
When A, B are maximally monotone (think it as the subdifferential ∂f
of a proper convex function f ) and C is β-cocoercive (think it as the gradient
∇f of a 1/β-Lipschitz differentiable function f ), a solution can be found by
the iteration (2) with T = T3S , introduced recently in [24], where
(13)

T3S := I − JγB + JγA ◦ (2JγB − I − γC ◦ JγB ).

2β
Indeed, by setting γ ∈ (0, 2β), T3S is ( 4β−γ
)-averaged (think it as a property
weaker than the Picard contraction; in particular, T may not have a fixed
point). Following the standard convergence result (cf. textbook [5]), provided
that T has a fixed point, the sequence from (2) converges to a fixed-point
x∗ of T . Note that, instead of x∗ , JγB (x∗ ) is a solution to (12).
Following §3.1, T3S is CF if JγA is separable (C1 ), JγB is Type-II CF
(F2 ), and C is Type-I CF (F1 ).
We give a few special cases of T3S below, which have much longer history.
They all converge to a fixed point x∗ whenever a solution exists and γ is
properly chosen. If B 6= 0, then JγB (x∗ ), instead of x∗ , is a solution to (12).
Forward-Backward Splitting (FBS): Letting B = 0 yields JγB = I.
Then, T3S reduces to FBS [52]:

(14)

TFBS := JγA ◦ (I − γC)

for solving the problem 0 ∈ Ax + Cx.
Backward-Forward Splitting (BFS): Letting A = 0 yields JγA = I.
Then, T3S reduces to BFS:
(15)

TBFS := (I − γC) ◦ JγB

20
for solving the problem 0 ∈ Bx + Cx. When A = B, TFBS and TBFS apply
the same pair of operators in the opposite orders, and they solve the same
problem. Iterations based on TBFS are rarely used in the literature because
they need an extra application of JγB to return the solution, so TBFS is
seemingly an unnecessary variant of TFBS . However, they become different
for coordinate update; in particular, TBFS is CF (but TFBS is generally not)
when JγB is Type-II CF (F2 ) and C is Type-I CF (F1 ). Therefore, TBFS is
worth discussing alone.
Douglas-Rachford Splitting (DRS): Letting C = 0, T3S reduces to
(16)

1
TDRS := I − JγB + JγA ◦ (2JγB − I) = (I + RγA ◦ RγB )
2

introduced in [26] for solving the problem 0 ∈ Ax+Bx. A more general splitting is the Relaxed Peaceman-Rachford Splitting (RPRS) with λ ∈ [0, 1]:
(17)

TRPRS = (1 − λ) I + λ RγA ◦ RγB ,

which recovers TDRS by setting λ = 12 and Peaceman-Rachford Splitting
(PRS) [53] by letting λ = 1.
Forward-Douglas-Rachford Splitting (FDRS): Let V be a linear
subspace, and NV and PV be its normal cone and projection operator, respectively. The FDRS [15]
TFDRS = I − PV + JγA ◦ (2PV − I − γPV ◦ C˜ ◦ PV ),
aims at finding a point x such that 0 ∈ A x+C˜ x+NV x. If an optimal x exists,
we have x ∈ V and NV x is the orthogonal complement of V . Therefore, the
problem is equivalent to finding x such that 0 ∈ A x + PV ◦ C˜ ◦ PV x + NV x.
Thus, T3S recovers TFDRS by letting B = NV and C = PV ◦ C˜ ◦ PV .
Forward-Backward-Forward Splitting (FBFS): Composing TFBS
with one more forward step gives TFBFS introduced in [71]:
(18)

TFBFS = −γC + (I − γC)JγA (I − γC).

TFBFS is not a special case of T3S . At the expense of one more application
of (I − γC), TFBFS relaxes the convergence condition of TFBS from the cocoercivity of C to its monotonicity. (For example, a nonzero skew symmetric
matrix is monotonic but not cocoercive.) From Table 2, we know that TFBFS
is CF if both C and JγA are separable.

Coordinate friendly structures, algorithms, and applications
3.2.1. Examples in Optimization.
(19)

21

Consider the optimization problem

minimize f (x) + g(x),
x∈X

where X is the feasible set and f and g are objective functions. We present
examples of operator splitting methods discussed above.
Example 13 (proximal gradient method). Let X = Rm , f be differentiable,
and g be proximable in (19). Setting A = ∂g and C = ∇f in (14) gives
JγA = proxγg and reduces xk+1 = TFBS (xk ) to prox-gradient iteration:
xk+1 = proxγg (xk − γ∇f (xk )).

(20)

A special case of (20) with g = ιX is the projected gradient iteration:
xk+1 = PX (xk − γ∇f (xk )).

(21)

If ∇f is CF and proxγg is (nearly-)separable (e.g., g(x) = kxk1 or the
indicator function of a box constraint) or if ∇f is Type-II CF and proxγg
is cheap (e.g., ∇f (x) = Ax − b and g = kxk2 ), then the FBS iteration (20)
is CF. In the latter case, we can also apply the BFS iteration (15) (i.e,
compute proxγg and then perform the gradient update), which is also CF.
Example 14 (ADMM). Setting X = Rm simplifies (19) to
(22)

minimize f (x) + g(y),
x,y

subject to x − y = 0.

The ADMM method iterates:
(23a)

xk+1 = proxγf (y k − γsk ),

(23b)

y k+1 = proxγg (xk+1 + γsk ),
1
sk+1 = sk + (xk+1 − y k+1 ).
γ

(23c)

(The iteration can be generalized to handle the constraint Ax − By = b.)
The dual problem of (22) is mins f ∗ (−s) + g ∗ (s), where f ∗ is the convex
conjugate of f . Letting A = −∂f ∗ (−·) and B = ∂g ∗ in (16) recovers the
iteration (23) through (see the derivation in Appendix B)
tk+1 = TDRS (tk ) = tk − JγB (tk ) + JγA ◦ (2JγB − I)(tk ).
From the results in §3.1, a sufficient condition for the above iteration to be
CF is that JγA is (nearly-)separable and JγB being CF.

22
The above abstract operators and their CF properties will be applied
in §5 to give interesting algorithms for several applications.

4. Primal-dual Coordinate Friendly Operators
We study how to solve the problem
(24)

minimize f (x) + g(x) + h(Ax),
x∈H

with primal-dual splitting algorithms, as well as their coordinate update
versions. Here, f is differentiable and A is a “p-by-m” linear operator from
H = H1 × · · · × Hm to G = G1 × · · · × Gp . Problem (24) abstracts many
applications in image processing and machine learning.
Example 15 (image deblurring/denoising). Let u0 be an image, where
u0i ∈ [0, 255], and B be the blurring linear operator. Let k∇uk1 be the
anisotropic5 total variation of u (see (49) for definition). Suppose that b
is a noisy observation of Bu0 . Then, we can try to recover u0 by solving
(25)

minimize
u

1
kBu − bk2 + ι[0,255] (u) + λk∇uk1 ,
2

which can be written in the form of (24) with f = 21 kB · −bk2 , g = ι[0,255] ,
A = ∇, and h = λk · k1 .
More examples with the formulation (24) will be given in §4.2. In general,
primal-dual methods are capable of solving complicated problems involving
constraints and the compositions of proximable and linear maps like k∇uk1 .
In many applications, although h is proximable, h ◦ A is generally nonproximable and non-differentiable. To avoid using slow subgradient methods,
we can consider the primal-dual splitting approaches to separate h and A so
that proxh can be applied. We derive that the equivalent form (for convex
cases) of (24) is to find x such that
(26)

0 ∈ (∇f + ∂g + A> ◦ ∂h ◦ A)(x).

Introducing the dual variable s ∈ G and applying the biconjugation prop5

Generalization to the isotropic case is straightforward by grouping variables
properly.

Coordinate friendly structures, algorithms, and applications

23

erty: s ∈ ∂h(Ax) ⇔ Ax ∈ ∂h∗ (s), yields the equivalent condition
 
  
 
 
x
∂g 0
∇f 0
0 A>
+
+
(27)
0∈
,
∗
0 ∂h
0 0
s
−A 0
| {z } |
{z
} |{z}
z
operator A
operator B
which we shorten as 0 ∈ Az + Bz, with z ∈ H × G =: F.
Problem (27) can be solved by the Condat-Vũ algorithm [20, 76]:
 k+1
s
= proxγh∗ (sk + γAxk ),
(28)
k+1
x
= proxηg (xk − η(∇f (xk ) + A> (2sk+1 − sk ))),
which explicitly applies A and A> and updates s, x in a Gauss-Seidel style 6 .
We introduce an operator TCV : F → F and write
iteration (28)

⇐⇒

z k+1 = TCV (z k ).

Switching the orders of x and s yields the following algorithm:
 k+1
x
= proxηg (xk − η(∇f (xk ) + A> sk )),
0
(29)
as z k+1 = TCV
zk .
k+1
s
= proxγh∗ (sk + γA(2xk+1 − xk )),
It is known from [18, 23] that both (28) and (29) reduce to iterations of
nonexpansive operators (under a special metric), i.e., TCV is nonexpansive;
see Appendix C for the reasoning.
Remark 1. Similar primal-dual algorithms can be used to solve other problems such as saddle point problems [40, 46, 16] and variational inequalities [68]. Our coordinate update algorithms below apply to these problems as
well.
4.1. Primal-dual Coordinate Update Algorithms
In this subsection, we make the following assumption.
Assumption 1. Functions g and h∗ in the problem (24) are separable and
proximable. Specifically,
g(x) =

m
X
i=1

gi (xi )

and

∗

h (y) =

p
X

h∗i (yi ).

j=1

By the Moreau identity: proxγh∗ = I − γprox γ1 h ( γ· ), one can compute prox γ1 h
instead of proxγh∗ , which inherits the same separability properties from prox γ1 h .
6

24
Furthermore, ∇f is CF.
Proposition 1. Under Assumption 1, the followings hold:
(a) when p = O(m), the Condat-Vu operator TCV in (28) is CF, more specifically,

h

k

+

+

i



M {z , Ax} 7→ {z , Ax } = O

h
i
1
k
k
;
M z 7→ TCV z
m+p

0
(b) when m  p and M [x 7→ ∇f (x)] = O(m), the Condat-Vu operator TCV
in (29) is CF, more specifically,


h
i
M {z k , A> s} 7→ {z + , A> s+ } = O

h
i
1
0
M z k 7→ TCV
zk .
m+p

Proof. Computing z k+1 = TCV z k involves evaluating ∇f , proxg , and proxh∗ ,


applying A and A> , and adding vectors. It is easy to see M z k 7→ TCV z k =
 k

0 z k is the same.
O(mp + m + p) + M[x → ∇f (x)], and M z 7→ TCV
(a) We assume ∇f ∈ F1 for simplicity, and other cases are similar.
1. If (TCV z k )j = sk+1
, computing it involves: adding ski and γ(Axk )i , and
i


evaluating proxγh∗i . In this case M {z k , Ax} 7→ {z + , Ax+ } = O(1).
2. If (TCV z k )j = xk+1
, computing it involves evaluating: the entire sk+1
i
for O(p) operations, (A> (2sk+1 − sk ))i for O(p) operations, proxηgi
1
for O(1) operations, ∇i f (xk ) for O( m
M [x 7→ ∇f (x)]) operations, as
+
well as updating Ax for O(p) operations. In this case


1
M {z k , Ax} 7→ {z + , Ax+ } = O(p + m
M [x 7→ ∇f (x)]).



Therefore, M {z k , Ax} 7→ {z + , Ax+ } = O

1
m+p M

 k

z 7→ TCV z k .

(b) When m  p and M [x 7→ ∇f (x)] = O(m), following arguments similar
to the above, we have


0 z k ) = xk+1 ;
M {z k , A> s} 7→ {z + , A> s+ } = O(1)+ M [x 7→ ∇i f (x)] if (TCV
j
i
 k >

+
>
+
0 zk ) =
and M {z , A s} 7→ {z , A s } = O(m) + M [x 7→ ∇f (x)] if (TCV
j
sk+1
.
i




1
0 z k ).
M z k 7→ TCV
In both cases M {z k , A> s} 7→ {z + , A> s+ } = O( m+p

Coordinate friendly structures, algorithms, and applications

25

4.2. Extended Monotropic Programming
We develop a primal-dual coordinate update algorithm for the extended
monotropic program:
(30)

minimize

g1 (x1 ) + g2 (x2 ) + · · · + gm (xm ) + f (x),

subject to

A1 x1 + A2 x2 + · · · + Am xm = b,

x∈H

where x = (x1 , . . . , xm ) ∈ H = H1 ×. . .× Hm with Hi being Euclidean spaces.
It generalizes linear, quadratic, second-order cone, semi-definite programs by
allowing extended-valued objective functions gi and f . It is a special case
m
X
of (24) by letting g(x) =
gi (xi ), A = [A1 , · · · , Am ] and h = ι{b} .
i=1

Example 16 (quadratic programming). Consider the quadratic program
(31)

minimize
m
x∈R

1 >
x U x + c> x, subject to Ax = b, x ∈ X,
2

where U is a symmetric positive semidefinite matrix and X = {x : xi ≥
0 ∀i}. Then, (31) is a special case of (30) with gi (xi ) = ι·≥0 (xi ), f (x) =
1 >
>
2 x U x + c x and h = ι{b} .
Example 17 (Second Order Cone Programming (SOCP)). The SOCP
minimize
c> x,
m
x∈R

subject to Ax = b,
x ∈ X = Q1 × · · · × Qn ,

(where the number of cones n may not be equal to the number of blocks m,)
can be written in the form of (30): minimizex∈Rm ιX (x) + c> x + ι{b} (Ax).
Applying iteration (28) to problem (30) and eliminating sk+1 from the
second row yield the Jacobi-style update (denoted as Temp ):


(32)

sk+1 = sk + γ(Axk − b),
xk+1 = proxηg (xk − η(∇f (xk ) + A> sk + 2γA> Axk − 2γA> b)).

To the best of our knowledge, this update is never found in the literature.
Note that xk+1 no longer depends on sk+1 , making it more convenient to
perform coordinate updates.

26
Remark 2. In general, when the s update is affine, we can decouple sk+1
and xk+1 by plugging the s update into the x update. It is the case when h
is affine or quadratic in problem (24).
A sufficient condition for Temp to be CF is proxg ∈ C1 i.e., separable.
Indeed, we have Temp = T1 ◦ T2 , where


  

I
0
s + γ(Ax − b)
s
T1 =
, T2
=
.
0 proxηg
x
x − η(∇f (x) + A> s + 2γA> Ax − 2γA> b)
Following Case 5 of Table 2, Temp is CF. When m = Θ(p), the separability
condition on proxg can be relaxed to proxg ∈ F1 since in this case T2 ∈ F2 ,
and we can apply Case 7 of Table 2 (by maintaining ∇f (x), A> s, Ax and
A> Ax.)
4.3. Overlapping-Block Coordinate Updates
In the coordinate update scheme based on (28), if we select xi to update then
we must first compute sk+1 , because the variables xi ’s and sj ’s are coupled
through the matrix A. However, once xk+1
is obtained, sk+1 is discarded. It
i
is not used to update s or cached for further use. This subsection introduces
ways to utilize the otherwise wasted computation.
We define, for each i, J(i) ⊂ [p] as the set of indices j such that A>
i,j 6= 0,
and, for each j, I(j) ⊂ [m] as the set of indices of i such that A>
=
6
0. We
i,j
also let mj := |I(j)|, and assume mj 6= 0, ∀j ∈ [p] without loss of generality.
We arrange the coordinates of z = [x; s] into m overlapping blocks. The
ith block consists of the coordinate xi and all sj ’s for j ∈ J(i). This way, each
sj may appear in more than one block. We propose a block coordinate update
scheme based on (28). Because the blocks overlap, each sj may be updated
in multiple blocks, so the sj P
update is relaxed with parameters ρi,j ≥ 0
(see (33) below) that satisfy
i∈I(j) ρi,j = 1, ∀j ∈ [p]. The aggregated
effect is to update sj without scaling. (Following the KM iteration [39], we
can also assign a relaxation parameter ηk for the xi update; then, the sj
update should be relaxed with ρi,j ηk .)
We propose the following update scheme:

select i ∈ [m], and then compute



k+1
k
k

all j ∈ J(i),

 s̃j = proxγh∗j (sj + γ(Ax )j ), forP
k+1
k+1
k
k
x̃i = proxηgi (xi − η(∇i f (x ) + j∈J(i) A>
− skj ))),
(33)
i,j (2s̃j



update xk+1
= xki + (x̃k+1
− xki ),

i
i


k+1
k+1
k
update sj = sj + ρi,j (s̃j − skj ), for all j ∈ J(i).

Coordinate friendly structures, algorithms, and applications

27

Remark 3. The use of relaxation parameters ρi,j makes our scheme different from that in [56].
Following the assumptions and arguments in §4.1, if we maintain Ax,
the cost for each block coordinate update is O(p) + M [x 7→ ∇i f (x)], which
1
is O( m
M [z 7→ TCV z]). Therefore the coordinate update scheme (33) is computationally worthy.
Typical choices of ρi,j include: (1) one of the ρi,j ’s is 1 for each j, others
all equal to 0. This can be viewed as assigning the update of sj solely to a
block containing xi . (2) ρi,j = m1j for all i ∈ I(j). This approach spreads the
update of sj over all the related blocks.
Remark 4. The recent paper [28] proposes a different primal-dual coordinate update algorithm. The authors produce a new matrix Ā based on A,
with only one nonzero entry in each row, i.e. mj = 1 for each j. They also
modify h to h̄ so that the problem
(34)

minimize f (x) + g(x) + h̄(Āx)
x∈H

has the same solution as (24). Then they solve (34) by the scheme (33).
Because they have mj = 1, every dual variable coordinate is only associated
with one primal variable coordinate. They create non-overlapping blocks of
z by duplicating each dual variable coordinate sj multiple times. The computation cost for each block coordinate update of their algorithm is the same
as (33), but more memory is needed for the duplicated copies of each sj .

4.4. Async-Parallel Primal-Dual Coordinate Update Algorithms
and Their Convergence

In this subsection, we propose two async-parallel primal-dual coordinate
update algorithms using the algorithmic framework of [54] and state their
convergence results. When there is only one agent, all algorithms proposed in
this section reduce to stochastic coordinate update algorithms [19], and their
convergence is a direct consequence of Theorem 1. Moreover, our convergence
analysis also applies to sync-parallel algorithms.

28
The two algorithms are based on §4.1 and §4.3, respectively.
Algorithm 1: Async-parallel primal-dual coordinate update algorithm using TCV
Input :P
z 0 ∈ F, K > 0, a discrete distribution (q1 , . . . , qm+p ) with
m+p
i=1 qi = 1 and qi > 0, ∀i,
set global iteration counter k = 0;
while k < K, every agent asynchronously and continuously do
select ik ∈ [m + p] with Prob(ik = i) = qi ;
perform an update to zik according to (35);
update the global counter k ← k + 1;

Whenever an agent updates a coordinate, the global iteration number k
increases by one. The kth update is applied to zik , with ik being independent
random variables: zi = xi when i ≤ m and zi = si−m when i > m. Each
coordinate update has the form:

(
(35)

= zikk −
zik+1
k
zik+1 = zik ,

ηk
(m+p)qik

(ẑikk − (TCV ẑ k )ik ),

∀i 6= ik ,

where ηk is the step size, z k denotes the state of z in global memory just
before the update (35) is applied, and ẑ k is the result that z in global memory
is read by an agent to its local cache (see [54, §1.2] for both consistent and
inconsistent cases). While (ẑikk −(TCV ẑ k )ik ) is being computed, asynchronous
parallel computing allows other agents to make updates to z, introducing
so-called asynchronous delays. Therefore, ẑ k can be different from z k . We
refer the reader to [54, §1.2] for more details.
The async-parallel algorithm using the overlapping-block coordinate update (33) is in Algorithm 2 (recall that the overlapping-block coordinate

Coordinate friendly structures, algorithms, and applications

29

update is introduced to save computation).
Algorithm 2: Async-parallel primal-dual overlapping-block coordinate update algorithm using TCV
Input :P
z 0 ∈ F, K > 0, a discrete distribution (q1 , . . . , qm ) with
m
i=1 qi = 1 and qi > 0, ∀i,
set global iteration counter k = 0;
while k < K, every agent asynchronously and continuously do
select ik ∈ [m] with Prob(ik = i) = qi ;
Compute s̃k+1
, ∀j ∈ J(ik ) and x̃k+1
according to (33);
j
ik
ηk
k+1
k+1
k
k
update xik = xik + mqi (x̃ik − xik );
k

let xk+1
= xki for i 6= ik ;
i
i,j ηk
update sk+1
= skj + ρmq
(s̃k+1
− skj ), for all j ∈ J(ik );
j
j
i
k

/ J(ik );
let sk+1
= skj , for all j ∈
j
update the global counter k ← k + 1;
Here we still allow asynchronous delays, so x̃ik and s̃k+1
are computed
j
k
using some ẑ .
Remark 5. If shared memory is used, it is recommended to set all but one
ρi,j ’s to 0 for each i.
Theorem 1. Let Z ∗ be the set of solutions to problem (24) and (z k )k≥0 ⊂ F
be the sequence generated by Algorithm 1 or Algorithm 2 under the following
conditions:
(i) f, g, h∗ are closed proper convex functions, f is differentiable, and ∇f
is Lipschitz continuous with constant β;
(ii) the delay for every coordinate is bounded by a positive number τ , i.e.
k−dk
for every 1 ≤ i ≤ m + p, ẑik = zi i for some 0 ≤ dki ≤ τ ;
(iii) ηk ∈ [ηmin , ηmax ] for certain ηmin , ηmax > 0.
Then (z k )k≥0 converges to a Z ∗ -valued random variable with probability 1.
The formulas for ηmin and ηmax , as well as the proof of Theorem 1, are
given in Appendix D along with additional remarks. The algorithms can be
applied to solve problem (24). A variety of examples are provided in §5.1
and §5.2.

30

5. Applications
In this section, we provide examples to illustrate how to develop coordinate
update algorithms based on CF operators. The applications are categorized
into five different areas. The first subsection discusses three well-known machine learning problems: empirical risk minimization, Support Vector Machine (SVM), and group Lasso. The second subsection discusses image processing problems including image deblurring, image denoising, and Computed Tomography (CT) image recovery. The remaining subsections provide
applications in finance, distributed computing as well as certain stylized optimization models. Several applications are treated with coordinate update
algorithms for the first time.
For each problem, we describe the operator T and how to efficiently
calculate (T x)i . The final algorithm is obtained after plugging the update
in a coordinate update framework in §1.1 along with parameter initialization,
an index selection rule, as well as some termination criteria.
5.1. Machine Learning
5.1.1. Empirical Risk Minimization (ERM). We consider the following regularized empirical risk minimization problem
p

(36)

minimize
m
x∈R

1X
φj (a>
j x) + f (x) + g(x),
p
j=1

where aj ’s are sample vectors, φj ’s are loss functions, and f + g is a regularization function. We assume that f is differentiable and g is proximable.
Examples of (36) include linear SVM, regularized logistic regression, ridge
regression, and Lasso. Further information on ERM can be found in [34]. The
need for coordinate update algorithms arises in many applications of (36)
where the number of samples or the dimension of x is large.
1 Pp
>
>
We define A = [a>
1 ; a2 ; . . . ; ap ] and h(y) := p
j=1 φj (yj ). Hence,
1 Pp
>
h(Ax) = p j=1 φj (aj x), and problem (36) reduces to form (24). We can
apply the primal-dual update scheme to solve this problem, for which we introduce the dual variable s = (s1 , ..., sp )> . We use p + 1 coordinates, where
the 0th coordinate is x ∈ Rm and the jth coordinate is sj , j ∈ [p]. The

Coordinate friendly structures, algorithms, and applications

31

operator T is given in (29). At each iteration, a coordinate is updated:

(37)


if x is chosen (the index 0), then compute




xk+1 = proxηg (xk − η(∇f (xk ) + A> sk )),


 if s is chosen (an index j ∈ [p]), then compute
j
x̃k+1 = proxηg (xk − η(∇f (xk ) + A> sk )),




and



k+1 − xk )).
sk+1
= p1 proxpγφ∗j (pskj + pγa>
j (2x̃
j

We maintain A> s ∈ Rm in the memory. Depending on the structure of
∇f , we can compute it each time or maintain it. When proxg ∈ F1 , we
can consider breaking x into coordinates xi ’s and also select an index i to
update xi at each time.
5.1.2. Support Vector Machine. Given the training data {(ai , βi )}m
i=1
with βi ∈ {+1, −1}, ∀i, the kernel support vector machine [65] is
minimize
(38)

x,ξ,y

subject to

Pm
1
2
i=1 ξi ,
2 kxk2 + C
βi (x> φ(ai ) − y) ≥ 1

− ξi , ξi ≥ 0, ∀i ∈ [m],

where φ is a vector-to-vector map, mapping each data ai to a point in a
(possibly) higher-dimensional space. If φ(a) = a, then (38) reduces to the
linear support vector machine. The model (38) can be interpreted as finding
a hyperplane {w : x> w − y = 0} to separate two sets of points {φ(ai ) : βi =
1} and {φ(ai ) : βi = −1}.
The dual problem of (38) is
(39)

X
1
minimize s> Qs − e> s, subject to 0 ≤ si ≤ C, ∀i,
βi si = 0,
s
2
i

where Qij = βi βj k(ai , aj ), k(·, ·) is a so-called kernel function, and e =
(1, ..., 1)> . If φ(a) = a, then k(ai , aj ) = a>
i aj .
Unbiased case. If y = 0 is enforced in (38), then the solution hyperplane
{w : x> w = 0} passes through the origin and is called unbiased. Consequently, the dual problem (39) will no longer have the linear constraint
P
i βi si = 0, leaving it with the coordinate-wise separable box constraints
0 ≤ si ≤ C. To solve (39), we can apply the FBS operator T defined by (14).

32
Let d(s) := 12 s> Qs − e> s, A = prox[0,C] , and C = ∇d. The coordinate update based on FBS is
sk+1
= proj[0,C] (ski − γi ∇i d(sk )),
i
where we can take γi =

1
Qii .

Biased (general) case. In this case, the mode (38) has y ∈ R, so the
hyperplane {w : x> w − y = 0} may not pass the origin and
P is called biased.
Then, the dual problem (39) retains the linear constraint i βi si = 0. In this
case, we apply the primal-dual splitting scheme (28) or the three-operator
splitting scheme (13).
The coordinate update based on the full primal-dual splitting scheme (28)
is:
(40a)

tk+1 = tk + γ

(40b)

sk+1
i

m
X

βi ski ,

i=1



= proj[0,C] ski − η ∇i d(sk ) + βi (2tk+1 − tk ) ,

where t, sPare the primal and dual variables, respectively. Note that we can
let w := m
i=1 βi si and maintain it. With variable w and substituting (40a)
into (40b), we can equivalently write (40) into

(41)


if t is chosen (the index 0), then compute




tk+1 = tk + γwk ,

if si is chosen (an index i ∈ [m]), then compute


k+1
k − η q > sk − 1 + β (2γw k + tk )

s
s
=
proj

i
[0,C]
i
i
i


wk+1 = wk + βi (sk+1
− ski ).
i

We can also apply the
P three-operator splitting (13) as follows. Let D1 :=
[0, C]m and D2 := {s : m
i=1 βi si = 0}. Let A = projD2 , B = projD1 , and
C(x) = Qx − e, The full update corresponding to T = (I − ηk )I + ηk T3S is
(42a)
(42b)

sk+1 =projD2 (uk ),



uk+1 =uk + ηk projD1 2sk+1 − uk − γ(Qsk+1 − e) − sk+1 ,

where s is just an intermediate variable. Let β̃ :=

β
kβk2

and w := β̃ > u. Then

projD2 (u) = (I − β̃ β̃ > )u. Hence, sk+1 = uk − wk β̃. Plugging it into (42b)

Coordinate friendly structures, algorithms, and applications

33

yields the following coordinate update scheme:


if i ∈ [m] is chosen, then compute



 sk+1 = uk − wk β̃i ,
i
i




k+1
k+1
k+1
k+η
k − γ q > uk − w k (q > β̃) − 1
u
=
u
proj
2s
−
u
−
s

k
[0,C]
i
i
i
i
i
i
i



wk+1 = wk + β̃i (uk+1
− uki ),
i
where wk is the maintained variable and sk is the intermediate variable.
5.1.3. Group Lasso.
(43)

The group Lasso regression problem [84] is
minimize
f (x) +
n
x∈R

m
X

λi kxi k2 ,

i=1

where f is a differentiable convex function, often bearing the form 12 kAx −
bk22 , and xi ∈ Rni is a subvector of x ∈ Rn supported on Ii ⊂ [n], and
∪i Ii = [n]. If Ii ∩ Ij = ∅, ∀i 6= j, it is called non-overlapping group Lasso,
and if there are two different groups Ii and Ij with a non-empty intersection,
it is called overlapping group Lasso. The model finds a coefficient vector x
that minimizes the fitting (or loss) function f (x) and that is group sparse:
all but a few xi ’s are zero.
Let Ui be formed by the columns of the identity matrix I corresponding
> ] ∈ R(Σi ni )×n . Then, x = U > x.
to the indices in Ii , and let U = [U1> ; . . . ; Um
i
Pm i
Let hi (yi ) = λi kyi k2 , yi ∈ Rni for i ∈ [m], and h(y) =
h
(y
)
for
i
i
i=1
y = [y1 ; . . . ; ym ] ∈ RΣi ni . In this way, (43) becomes
(44)

minimize f (x) + h(U x).
x

Non-overlapping case [84]. In this case, we have Ii ∩ Ij = ∅, ∀i 6= j,
and can apply the FBS scheme (14) to (44). Specifically, let T1 = ∂(h ◦ U )
and T2 = ∇f . The FBS full update is
xk+1 = JγT1 ◦ (I − γT2 )(xk ).
The corresponding coordinate update is the following


(45)

if i ∈ [m] is chosen, then compute
xk+1
= arg minxi 21 kxi − xki + γi ∇i f (xk )k22 + γi hi (xi ),
i

34
where ∇i f (xk ) is the partial derivative of f with respect to xi and the step
1
size can be taken to be γi = kA:,i
k2 . When ∇f is either cheap or easy-tomaintain, the coordinate update in (45) is inexpensive.
Overlapping case [38]. This case allows Ii ∩ Ij 6= ∅ for some i 6= j,
causing the evaluation of JγT1 to be generally difficult. However, we can
apply the primal-dual update (28) to this problem as
(46a)

sk+1 = proxγh∗ (sk + γU xk ),

(46b)

xk+1 = xk − η(∇f (xk ) + U > (2sk+1 − sk )),

where s is the dual variable. Note that

0,
if ksi k2 ≤ λi , ∀i,
∗
h (s) =
+∞, otherwise,
is cheap. Hence, the corresponding coordinate update of (46) is
(47)

if si is chosen for some i ∈ [m], then compute



 sk+1 = projB (sk + γxk )
i
i
i
λi
if
x
is
chosen
for
some
i
∈
[m], then compute
i





 xk+1 = xk − η U T ∇f (xk ) + U T P T
k + γxk ) − sk ) ,
U
(2proj
(s
Bλj j
i
i
i
j
j
j,Ui Uj 6=0 j
i
where Bλ is the Euclidean ball of radius λ. When ∇f is easy-to-maintain,
the coordinate update in (47) is inexpensive. To the best of our knowledge,
the coordinate update method (47) is new.
5.2. Imaging
5.2.1. DRS for Image Processing in the Primal-dual Form [50].
Many convex image processing problems have the general form
minimize f (x) + g(Ax),
x

where A is a matrix such as a dictionary, sampling operator, or finite difference operator. We can reduce the problem to the system: 0 ∈ A(z) + B(z),
where z = [x; s],



 
∂f (x)
0
A> x
A(z) :=
, and B(z) :=
.
∂g ∗ (s)
s
−A
0

Coordinate friendly structures, algorithms, and applications

35

(see Appendix C for the reduction.) The work [50] gives their resolvents
JγA
JγB


proxγf
=
0

0
proxγg∗

0
−1
= (I + γB) =
0


,

  

>
0
I
I
2 >
−1
+
(I + γ A A)
,
I
γA
−γA

where JγA is often cheap or separable and we can explicitly form JγB as
a matrix or implement it based on a fast transform. With the defined JγA
and JγB , we can apply the RPRS method as z k+1 = TRPRS z k . The resulting RPRS operator is CF when JγB is CF. Hence, we can derive a new
RPRS coordinate update algorithm. We leave the derivation to the readers. Derivations of coordinate update algorithms for more specific image
processing problems are shown in the following subsections.
5.2.2. Total Variation Image Processing.
Total Variation (TV) image processing model
(48)

We consider the following

1
minimize λkxkTV + kA x − bk2 ,
x
2

where x ∈ Rn is the vector representation of the unknown image, A is an
m × n matrix describing the transformation from the image to the measurements b. Common A includes sampling matrices in MRI, CT, denoising, deblurring, etc. Let (∇hi , ∇vi ) be the discrete gradient at pixel i and
∇x = (∇h1 x, ∇v1 x, . . . , ∇hn x, ∇vn x)> . Then the TV semi-norm k · kTV in the
isotropic and anisotropic fashions are, respectively,
(49a)

kxkTV =

Xq

(∇hi x)2 + (∇vi x)2 ,

i

(49b)

kxkTV = k∇xk1 =

X


|∇hi x| + |∇vi x| .

i

For simplicity, we use the anisotropic TV for analysis and in the numerical experiment in § 6.2. It is slightly more complicated for the isotropic TV.
Introducing the following notation

 
∇
B: =
,
A

1
h(p, q): = λkpk1 + kq − bk2 ,
2

36
we can reformulate (48) as
minimize h(B x) = h(∇x, A x),
x

which reduces to the form of (24) with f = g = 0. Based on its definition,
the convex conjugate of h(p, q) and its proximal operator are, respectively,
(50)
(51)

1
1
h∗ (s, t) = ιk·k∞ ≤λ (s) + kt + bk2 − kbk2 ,
2
2
1
proxγh∗ (s, t) = projk·k∞ ≤λ (s) +
(t − γb).
1+γ

Let s, t be the dual variables corresponding to ∇x and Ax respectively, then
using (51) and applying (29) give the following full update:
(52a)
(52b)
(52c)

xk+1 = xk − η(∇> sk + A> tk ),


sk+1 = projk·k∞ ≤λ sk + γ∇(xk − 2η(∇> sk + A> tk )) ,

1 k
tk+1 =
t + γA(xk − 2η(∇> sk + A> tk )) − γb .
1+γ

To perform the coordinate updates as described in §4, we can maintain ∇> sk
and A> tk . Whenever a coordinate of (s, t) is updated, the corresponding
∇> sk (or A> tk ) should also be updated. Specifically, we have the following
coordinate update algorithm

if xi is chosen for some i ∈ [n], then compute





= xki − η(∇> sk + A> tk )i ;
xk+1

i



 if si is chosen for some i ∈ [2n], then compute



= projk·k∞ ≤λ ski + γ∇i (xk − 2η(∇> sk + A> tk ))
sk+1
i
(53)
and update ∇> sk to ∇> sk+1 ;




if ti is chosen for some i ∈ [m], then compute




1


tki + γAi,: (xk − 2η(∇> sk + A> tk )) − γbi
tk+1
= 1+γ

i


and update A> tk to A> tk+1 .
5.2.3. 3D Mesh Denoising. Following an example in [60], we consider
Y
Z
a 3D mesh described by their nodes x̄i = (x̄X
i , x̄i , x̄i ), i ∈ [n], and the
adjacency matrix A ∈ Rn×n , where Aij = 1 if nodes i and j are adjacent,
otherwise Aij = 0. We let Vi be the set of neighbours of node i. Noisy mesh
nodes zi , i ∈ [n], are observed. We try to recover the original mesh nodes by

Coordinate friendly structures, algorithms, and applications

37

solving the following optimization problem [60]:
(54)

minimize
x

n
X

fi (xi ) +

i=1

n
X

gi (xi ) +

i=1

XX
i

hi,j (xi − xj ),

j∈Vi

where fi ’s are differentiable data fidelity terms, gi ’s are the indicator funcP P
tions of box constraints, and i j∈Vi hi,j (xi − xj ) is the total variation on
the mesh.
We introduce a dual variable s with coordinates si,j , for all ordered
pairs of adjacent nodes (i, j), and, based on the overlapping-block coordinate
updating scheme (33), perform coordinate update:


select i from [n], then compute



k
k
k

s̃k+1

i,j = proxγh∗i,j (si,j + γxi − γxj ), ∀j ∈ Vi ,



k
k
k

s̃k+1

j,i = proxγh∗j,i (sj,i + γxj − γxi ), ∀j ∈ Vi ,
and update
P

k+1
k
k


xi k+1 = proxηgi (xki − η(∇fi (xki ) + j∈Vi (2s̃k+1

i,j − 2s̃j,i − si,j + sj,i ))),


1 k+1
k
k

sk+1

i,j = si,j + 2 (s̃i,j − si,j ), ∀j ∈ Vi ,


1 k+1
k
k
sk+1
j,i = sj,i + 2 (s̃j,i − sj,i ), ∀j ∈ Vi .
5.3. Finance
5.3.1. Portfolio Optimization. Assume that we have one unit of capital
and m assets to invest on. The ith asset has an expected return rate ξi ≥ 0.
Our goal is to find a portfolio with the minimal risk such that the expected
return is no less than c. This problem can be formulated as
1 >
x Qx,
2
m
m
X
X
subject to x ≥ 0,
xi ≤ 1,
ξi xi ≥ c,

minimize
x

i=1

i=1

where the objective function is a measure of risk, and the last constraint
√
√
imposes that the expected return is at least c. Let a1 = e/ m, b1 = 1/ m,
a2 = ξ/kξk2 , and b2 = c/kξk2 , where e = (1, . . . , 1)> , ξ = (ξ1 , . . . , ξm )> . The
above problem is rewritten as
(55)

1
>
minimize x> Qx, subject to x ≥ 0, a>
1 x ≤ b1 , a2 x ≥ b2 .
x
2

38
We apply the three-operator splitting scheme (13) to (55). Let f (x) =
1 >
2 x Qx,
{x : a>
1x

>
D1 = {x : x ≥ 0}, D2 = {x : a>
1 x ≤ b1 , a2 x ≥ b2 }, D21 =

= b1 }, and D22 = {x : a>
2 x = b2 }. Based on (13), the full update is

(56a)

y k+1 =projD2 (xk ),

(56b)


xk+1 =xk + ηk projD1 (2y k+1 − xk − γ∇f (y k+1 )) − y k+1 ,

where y is an intermediate variable. As the projection to D1 is simple, we
discuss how to evaluate the projection to D2 . Assume that a1 and a2 are
neither perpendicular nor co-linear, i.e., a>
1 a2 6= 0 and a1 6= λa2 for any
scalar λ. In addition, assume a>
1 a2 > 0 for simplicity. Let a3 = a2 −
b3 = b2 −

1
b1 ,
a>
1 a2

a4 = a1 −

1
a2 ,
a>
1 a2

and b4 = b1 −

1
b2 .
a>
1 a2

1
a1 ,
a>
1 a2

Then we can

partition the whole space into four areas by the four hyperplanes a>
i x = bi ,
>
i = 1, . . . , 4. Let Pi = {x : a>
i x ≤ bi , ai+1 x ≥ bi+1 }, i = 1, 2, 3 and P4 = {x :
>
a>
4 x ≤ b4 , a1 x ≥ b1 }. Then


x,



projD22 (x),
projD2 (x) =
projD21 ∩D22 (x),



projD21 (x),

if
if
if
if

x ∈ P1 ,
x ∈ P2 ,
x ∈ P3 ,
x ∈ P4 .

Let wi = a>
i x − bi , i = 1, 2, and maintain w1 , w2 . Let ã2 =
ã1 =

a1 −a2 (a>
1 a2 )
2 .
1−(a>
1 a2 )

Then

projD21 (x) = x − w1 a1 ,
projD22 (x) = x − w2 a2 ,
projD21 ∩D22 (x) = x − w1 ã1 − w2 ã2 ,

a2 −a1 (a>
1 a2 )
2 ,
1−(a>
1 a2 )

Coordinate friendly structures, algorithms, and applications

39

Hence, the coordinate update of (56) is
(57a)
xk ∈ P1 : xk+1
=(1 − ηk )xki + ηk max(0, xki − γqi> xk ),
i
xk ∈ P2 : xk+1
=(1 − ηk )xki + ηk w2k (a2 )i + ηk max (0,
i

xki − γqi> xk − w2k (2(a2 )i − γqi> a2 ) ,
(57b)

xk ∈ P3 : xk+1
=(1 − ηk )xki + ηk w1k (ã1 )i + w2k (ã2 )i + ηk max (0,
i


xki − γqi> xk − w1k (2(ã1 )i − γqi> ã1 ) − w2k (2(ã2 )i − γqi> ã2 ) ,

(57c)

xk ∈ P4 : xk+1
=(1 − ηk )xki + ηk w1k (a1 )i + ηk max (0,
i

(57d)
xki − γqi> xk − w1k (2(a1 )i − γqi> a1 ) ,
where qi is the ith column of Q. At each iteration, we select i ∈ [m], and
perform an update to xi according to (57) based on where xk is. We then
renew wjk+1 = wjk +aij (xk+1
−xki ), j = 1, 2. Note that checking xk in some Pj
i
requires only O(1) operations by using w1 and w2 , so the coordinate update
in (57) is inexpensive.
5.4. Distributed Computing
5.4.1. Network. Consider that m worker agents and one master agent
form a star-shaped network, where the master agent at the center connects
to each of the worker agents. The m + 1 agents collaboratively solve the
consensus problem:
m
X
fi (x),
minimize
x

i=1

Rd

where x ∈
is the common variable and each proximable function fi is
held privately by agent i. The problem can be reformulated as
(58)

minimize d F (x) :=

x1 ,...,xm ,y∈R

m
X

fi (xi ),

subject to xi = y, ∀i ∈ [m],

i=1

which has the KKT condition
 

  
0 0
I
x
x
∂F 0 0
>







y + 0 0 −e
y ,
(59)
0∈ 0 0 0
s
s
0 0 0
I −e
0
{z
}
|
|
{z
}
operator A
operator C

40
where s is the dual variable.
Applying the FBFS scheme (18) to (59) yields the following full update:

(60a)
(60b)

xk+1
= proxγfi (xki − γski ) + γ 2 xki − γ 2 y k − 2γski ,
i
m
m
X
X
k+1
2 k
k
2
y
= (1 + mγ )y + 3γ
sj − γ
xkj ,
j=1

(60c)

j=1

sk+1
= ski − 2γxki − γproxγfi (xki − γski ) + 3γy k + γ 2
i

m
X

skj ,

j=1

where (60a) and (60c) are applied to all i ∈ [m]. Hence, for each i, we
group xi and si together and assign them on agent i. We let the master
P
P
agent maintain j sj and j xj . Therefore, in the FBFS coordinate update,
P
updating any (xi , si ) needs only y and j sj from the master agent, and
updating y is done on the master agent. In synchronous parallel setting, at
, then the master
, xk+1
each iteration, each worker agent i computes sk+1
i
i
agent collects the updates from all of the worker agents and then updates
P
y and j sj . The above update can be relaxed to be asynchronous. In this
case, the master and worker agents work concurrently, the master agent
P
updates y and j sj as soon as it receives the updated si and xi from any
of the worker agents. It also periodically broadcasts y back to the worker
agents.
5.5. Dimension Reduction
5.5.1. Nonnegative Matrix Factorization. Nonnegative Matrix Factorization (NMF) is an important dimension reduction method for nonnegative data. It was proposed by Paatero and his coworkers in [51]. Given a
nonnegative matrix A ∈ Rp×n
+ , NMF aims at finding two nonnegative matrip×r
n×r
ces W ∈ R+ and H ∈ R+ such that W H > ≈ A, where r is user-specified
depending on the applications, and usually r  min(p, n). A widely used
model is

(61)

1
minimize F (W, H) := kW H > − Ak2F ,
W,H
2
n×r
subject to W ∈ Rp×r
+ , H ∈ R+ .

Coordinate friendly structures, algorithms, and applications

41

Applying the projected gradient method (21) to (61), we have

(62a)
W k+1 = max 0, W k − ηk ∇W F (W k , H k ) ,

(62b)
H k+1 = max 0, H k − ηk ∇H F (W k , H k ) .
In general, we do not know the Lipschitz constant of ∇F , so we have to
choose ηk by line search such that the Armijo condition is satisfied.
Partitioning the variables into 2r block coordinates: (w1 , . . . , wr , h1 , . . . , hr )
where wi and hi are the ith columns of W and H, respectively, we can apply
the coordinate update based on the projected-gradient method:

(63)


if wik is chosen for some ik ∈ [r], then compute




= max 0, wikk − ηk ∇wik F (W k , H k ) ;
wik+1
k
if hik −r is chosen for some ik ∈ {r + 1, ..., 2r}, then


 compute

k
k, Hk) .
hk+1
=
max
0,
h
−
η
∇
F
(W
k hik −r
ik −r
ik −r

It is easy to see that ∇wi F (W k , H k ) and ∇hi F (W k , H k ) are both Lipschitz
continuous with constants khki k22 and kwik k22 respectively. Hence, we can set

(
ηk =

1
,
khkik k22
1
,
kwikk −r k22

if 1 ≤ ik ≤ r,
if r + 1 ≤ ik ≤ 2r.

However, it is possible to have wik = 0 or hki = 0 for some i and k, and
thus the setting in the above formula may have trouble of being divided
by zero. To overcome this problem, one can first modify the problem (61)
by restricting W to have unit-norm columns and then apply the coordinate
update method in (63). Note that the modification does not change the
optimal value since W H > = (W D)(HD−1 )> for any r×r invertible diagonal
matrix D. We refer the readers to [82] for more details.
Note that ∇W F (W, H) = (W H > −A)H, ∇H F (W, H) = (W H > −A)> W
and ∇wi F (W, H) = (W H > − A)hi , ∇hi F (W, H) = (W H > − A)> wi , ∀i.
Therefore, the coordinate updates given in (63) are computationally worthy
(by maintaining the residual W k (H k )> − A).
5.6. Stylized Optimization
5.6.1. Second-Order Cone Programming (SOCP). SOCP extends
LP by incorporating second-order cones. A second-order cone in Rn is

Q = (x1 , x2 , . . . , xn ) ∈ Rn : k(x2 , . . . , xn )k2 ≤ x1 .

42
Given a point v ∈ Rn , let ρv1 := k(v2 , . . . , vn )k2 and ρv2 := 12 (v1 + ρv1 ). Then,
the projectionv of v to Q returns 0 if v1 < −ρv1 , returns v if v1 ≥ ρv1 , and
returns (ρv2 , ρρ2v · (v2 , . . . , vn )) otherwise. Therefore, if we define the scalar
1
couple:

v1 < −ρv1 ,

(0, 0),
(ξ1v , ξ2v ) = (1, 1),
v1 ≥ ρv1 ,

 v ρv2 
otherwise,
ρ2 , ρv ,
1

then we have u = projQ (v) = ξ1v v1 , ξ2v · (v2 , . . . , vn ) . Based on this, we
have
Proposition 2.
1. Let v ∈ Rn and v + := v + νei for any ν ∈ R. Then,
given ρv1 , ρv2 , ξ1v , ξ2v defined above, it takes O(1) operations to obtain
+
+
+
+
ρv1 , ρv2 , ξ1v , ξ2v .
2. Let v ∈ Rn and A = [a1 A2 ] ∈ Rm×n , where a1 ∈ Rm , A2 ∈ Rm×(n−1) .
Given ρv1 , ρv2 , ξ1v , ξ2v , we have
A(2 · projQ (v) − v) = ((2ξ1v − 1)v1 ) · a1 + (2ξ2v − 1) · A2 (v2 , . . . , vn )> .
By the proposition, if T1 is an affine operator, then in the composition
T1 ◦ projQ , the computation of projQ is cheap as long as we maintain
ρv1 , ρv2 , ξ1v , ξ2v .
Given x, c ∈ Rn , b ∈ Rm , and A ∈ Rm×n , the standard form of SOCP is
(64a)

minimize c> x, subject to Ax = b,
x

x ∈ X = Q1 × · · · × Qn̄ ,

(64b)

where each Qi is a second-order cone, and n̄ 6= n in general. The problem (64)
is equivalent to

minimize c> x + ιA·=b (x) + ιX (x),
x

to which we can apply the DRS iteration z k+1 = TDRS (z k ) (see (16)), in
which JγA = projX and TγB is a linear operator given by
JγB (x) = arg min c> y +
y

1
ky − xk2
2γ

subject to Ay = b.

Assume that the matrix A has full row-rank (otherwise, Ax = b has either
redundant rows or no solution). Then, in (16), we have RγB (x) = Bx + d,
where B := I − 2A> (AA> )−1 A and d := 2A> (AA> )−1 (b + γAc) − 2γc.

Coordinate friendly structures, algorithms, and applications

43

It is easy to apply coordinate updates to z k+1 = TDRS (z k ) following
Proposition 2. Specifically, by maintaining the scalars ρv1 , ρv2 , ξ1v , ξ2v for each
v = xi ∈ Qi during coordinate updates, the computation of the projection
can be completely avoided. We pre-compute (AA> )−1 and cache the matrix
B and vector d. Then, TDRS is CF, and we have the following coordinate
update method

(65)


 select i ∈ [n̄], then compute
yik+1 = Bi xk + di

xk+1
= projQi (yik+1 ) + 12 (xki − yik+1 ),
i

where Bi ∈ Rni ×n is the ith row block submatrix of B, and yik+1 is the
intermediate variable.
It is trivial to extend this method for SOCPs with a quadratic objective:
1
minimize c> x + x> Cx, subject to Ax = b, x ∈ X = Q1 × · · · × Qn̄ ,
x
2
because J2 is still linear. Clearly, this method applies to linear programs as
they are special SOCPs.
Note that many LPs and SOCPs have sparse matrices A, which deserve
further investigation. In particular, we may prefer not to form (AA> )−1 and
use the results in §4.2 instead.

6. Numerical Experiments
We illustrate the behavior of coordinate update algorithms for solving portfolio optimization, image processing, and sparse logistic regression problems.
Our primary goal is to show the efficiency of coordinate update algorithms
compared to the corresponding full update algorithms. We will also illustrate
that asynchronous parallel coordinate update algorithms are more scalable
than their synchronous parallel counterparts.
Our first two experiments run on Mac OSX 10.9 with 2.4 GHz Intel Core
i5 and 8 Gigabytes of RAM. The experiments were coded in Matlab. The
sparse logistic regression experiment runs on 1 to 16 threads on a machine
with two 2.5 Ghz 10-core Intel Xeon E5-2670v2 (20 cores in total) and
64 Gigabytes of RAM. The experiment was coded in C++ with OpenMP
enabled. We use the Eigen library7 for sparse matrix operations.
7

http://eigen.tuxfamily.org

44
6.1. Portfolio Optimization
In this subsection, we compare the performance of the 3S splitting scheme (56)
with the corresponding coordinate update algorithm (57) for solving the
portfolio optimization problem (55). In this problem, our goal is to distribute our investment resources to all the assets so that the investment
risk is minimized and the expected return is greater than c. This test uses
two datasets, which are summarized in Table 3. The NASDAQ dataset is
collected through Yahoo! Finance. We collected one year (from 10/31/2014
to 10/31/2015) of historical closing prices for 2730 stocks.

Number of assets (N)
Expected return rate
Asset return rate
Risk

Synthetic data

NASDAQ data

1000
0.02
3 * rand(N, 1) - 1
covariance matrix + 0.01 · I

2730
0.02
mean of 30 days return rate
positive definite matrix

Table 3: Two datasets for portfolio optimization
In our numerical experiments, for comparison purposes, we first obtain a
high accurate solution by solving (55) with an interior point solver. For both
full update and coordinate update, ηk is set to 0.8. However, we use different
2
, and for 3S
γ. For 3S full update, we used the step size parameter γ1 = kQk
2
2
coordinate update, γ2 = max{Q11 ,...,QN N } . In general, coordinate update can
benefit from more relaxed parameters. The results are reported in Figure 3.
We can observe that the coordinate update method converges much faster
than the 3S method for the synthetic data. This is due to the fact that γ2
is much larger than γ1 . However, for the NASDAQ dataset, γ1 ≈ γ2 , so 3S
coordinate update is only moderately faster than 3S full update.
6.2. Computed Tomography Image Reconstruction
We compare the performance of algorithm (52) and its corresponding coordinate version on Computed Tomography (CT) image reconstruction. We
generate a thorax phantom of size 284 × 284 to simulate spectral CT measurements. We then apply the Siddon’s algorithm [66] to form the sinogram
data. There are 90 parallel beam projections and, for each projection, there
are 362 measurements. Then the sinogram data is corrupted with Gaussian
noise. We formulate the image reconstruction problem in the form of (48).
The primal-dual full update corresponds to (52). For coordinate update, the

Coordinate friendly structures, algorithms, and applications
3S full up date vs 3S coordinate up date

−2

3S full up date vs 3S coordinate up date

−2

10

45

10

−4

10

−3

10

−6

10

fk − f∗

fk − f∗

−8

10

−10

10

−4

10

−12

10

−5

10

−14

10

3S full update
3S coordinate update

−16

10

0

20

40

60
Ep ochs

3S full update
3S coordinate update

−6

10
80

100

(a) Synthesis dataset

0

20

40

60
Ep ochs

80

100

(b) NASDAQ dataset

Figure 3: Compare the convergence of 3S full update with 3S coordinate
update algorithms.

block size for x is set to 284, which corresponds to a column of the image.
The dual variables s, t are also partitioned into 284 blocks accordingly. A
block of x and the corresponding blocks of s and t are bundled together as
a single block. In each iteration, a bundled block is randomly chosen and
updated. The reconstruction results are shown in Figure 4. After 100 epochs,
the image recovered by the coordinate version is better than that by (52).
As shown in Figure 4d, the coordinate version converges faster than (52).
6.3. `1 Regularized Logistic Regression
In this subsection, we compare the performance of sync-parallel coordinate
update and async-parallel coordinate update for solving the sparse logistic
regression problem
(66)

minimize
λkxk1 +
n
x∈R

N

1 X
log 1 + exp(−bj · a>
j x) ,
N
j=1

where {(aj , bj )}N
j=1 is the set of sample-label pairs with bj ∈ {1, −1}, λ =
0.0001, and n and N represent the numbers of features and samples, respec-

46

(a) Phantom image

(b) Recovered by PDS
10

6
PDS full update
PDS coordinate update

10 5

10 4

10 3

10 2
0

(c) Recovered by PDS coord

20

40

60

80

100

(d) Objective function value

Figure 4: CT image reconstruction.

tively. This test uses the datasets8 : real-sim and news20, which are summarized in Table 4.
We let each coordinate hold roughly 50 features. Since the total number
of features is not divisible by 50, some coordinates have 51 features. We let
each thread draw a coordinate uniformly at random at each iteration. We
stop all the tests after 10 epochs since they have nearly identical progress
per epoch. The step size is set to ηk = 0.9, ∀k. Let A = [a1 , . . . , aN ]> and
b = [b1 , ..., bN ]> . In global memory, we store A, b and x. We also store the
product Ax in global memory so that the forward step can be efficiently
computed. Whenever a coordinate of x gets updated, Ax is immediately
8

http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/

Coordinate friendly structures, algorithms, and applications
Name
real-sim
news20

# samples
72, 309
19,996

47

# features
20, 958
1,355,191

Table 4: Two datasets for sparse logistic regression.
updated at a low cost. Note that if Ax is not stored in global memory, every
coordinate update will have to compute Ax from scratch, which involves the
entire x and will be very expensive.
Table 5 gives the running times of the sync-parallel and async-parallel
implementations on the two datasets. We can observe that async-parallel
achieves almost-linear speedup, but sync-parallel scales very poorly as we
explain below.
In the sync-parallel implementation, all the running threads have to
wait for the last thread to finish an iteration, and therefore if a thread has a
large load, it slows down the iteration. Although every thread is (randomly)
assigned to roughly the same number of features (either 50 or 51 components
of x) at each iteration, their ai ’s have very different numbers of nonzeros,
and the thread with the largest number of nonzeros is the slowest. (Sparse
matrix computation is used for both datasets, which are very large.) As
more threads are used, despite that they altogether do more work at each
iteration, the per-iteration time may increase as the slowest thread tends
to be slower. On the other hand, async-parallel coordinate update does not
suffer from the load imbalance. Its performance grows nearly linear with the
number of threads.
Finally, we have observed that the progress toward solving (66) is mainly
a function of the number of epochs and does not change appreciably when
the number of threads increases or between sync-parallel and async-parallel.
Therefore, we always stop at 10 epochs.

7. Conclusions
We have presented a coordinate update method for fixed-point iterations,
which updates one coordinate (or a few variables) at every iteration and can
be applied to solve linear systems, optimization problems, saddle point problems, variational inequalities, and so on. We proposed a new concept called
CF operator. When an operator is CF, its coordinate update is computationally worthy and often preferable over the full update method, in particular in
a parallel computing setting. We gave examples of CF operators and also discussed how the properties can be preserved by composing two or more such
operators such as in operator splitting and primal-dual splitting schemes. In

48

# threads
1
2
4
8
16

time
async
81.6
45.9
21.6
16.1
7.1

real-sim
(s)
speedup
sync async sync
82.1
1.0
1.0
80.6
1.8
1.0
63.0
3.8
1.3
61.4
5.1
1.3
46.4
11.5
1.8

time
async
591.1
304.2
150.4
78.3
41.6

news20
(s)
speedup
sync async sync
591.3
1.0
1.0
590.1
1.9
1.0
557.0
3.9
1.1
525.1
7.5
1.1
493.2
14.2
1.2

Table 5: Running times of async-parallel and sync-parallel FBS implementations for `1 regularized logistic regression on two datasets. Sync-parallel
has very poor speedup due to the large distribution of coordinate sparsity
and thus the large load imbalance across threads.

addition, we have developed CF algorithms for problems arising in several
different areas including machine learning, imaging, finance, and distributed
computing. Numerical experiments on portfolio optimization, CT imaging,
and logistic regression have been provided to demonstrate the superiority
of CF methods over their counterparts that update all coordinates at every
iteration.

References
[1] Attouch, H., Bolte, J., Redont, P., Soubeyran, A.: Proximal alternating minimization and projection methods for nonconvex problems: An
approach based on the Kurdyka-Lojasiewicz inequality. Mathematics
of Operations Research 35(2), 438–457 (2010)
[2] Bahi, J., Miellou, J.C., Rhofir, K.: Asynchronous multisplitting methods for nonlinear fixed point problems. Numerical Algorithms 15(3-4),
315–345 (1997)
[3] Baudet, G.M.: Asynchronous iterative methods for multiprocessors. J.
ACM 25(2), 226–244 (1978).
[4] Bauschke, H.H., Borwein, J.M.: On the convergence of von Neumann’s
alternating projection algorithm for two sets. Set-Valued Analysis 1(2),
185–212 (1993)
[5] Bauschke, H.H., Combettes, P.L.: Convex analysis and monotone operator theory in Hilbert spaces. Springer Science & Business Media
(2011)

Coordinate friendly structures, algorithms, and applications

49

[6] Baz, D.E., Frommer, A., Spiteri, P.: Asynchronous iterations with flexible communication: contracting operators. Journal of Computational
and Applied Mathematics 176(1), 91 – 103 (2005)
[7] Baz, D.E., Gazen, D., Jarraya, M., Spiteri, P., Miellou, J.: Flexible
communication for parallel asynchronous methods with application to
a nonlinear optimization problem. In: E. D’Hollander, F. Peters, G. Joubert, U. Trottenberg, R. Volpel (eds.) Parallel Computing Fundamentals, Applications and New Directions, Advances in Parallel Computing,
vol. 12, pp. 429 – 436. North-Holland (1998)
[8] Beck, A., Tetruashvili, L.: On the convergence of block coordinate descent type methods. SIAM Journal on Optimization 23(4), 2037–2060
(2013)
[9] Bengio, Y., Delalleau, O., Le Roux, N.: Label propagation and
quadratic criterion. In: Semi-Supervised Learning, pp. 193–216. MIT
Press (2006)
[10] Bertsekas, D.P.: Distributed asynchronous computation of fixed points.
Mathematical Programming 27(1), 107–120 (1983)
[11] Bertsekas, D.P.: Nonlinear programming. Athena Scientific (1999)
[12] Bertsekas, D.P., Tsitsiklis, J.N.: Parallel and distributed computation:
numerical methods. Prentice hall Englewood Cliffs, NJ (1989)
[13] Bolte, J., Sabach, S., Teboulle, M.: Proximal alternating linearized minimization for nonconvex and nonsmooth problems. Mathematical Programming 146(1-2), 459–494 (2014)
[14] Bradley, J.K., Kyrola, A., Bickson, D., Guestrin, C.: Parallel coordinate
descent for l1-regularized loss minimization. In: Proceedings of the 28th
International Conference on Machine Learning (ICML-11), pp. 321–328
(2011)
[15] Briceño-Arias, L.M.: Forward-Douglas–Rachford splitting and forwardpartial inverse method for solving monotone inclusions. Optimization
64(5), 1239–1261 (2015)
[16] Briceno-Arias, L.M., Combettes, P.L.: Monotone operator methods for
Nash equilibria in non-potential games. In: Computational and Analytical Mathematics, pp. 143–159. Springer (2013)
[17] Chazan, D., Miranker, W.: Chaotic relaxation. Linear Algebra and its
Applications 2(2), 199–222 (1969)

50
[18] Combettes, P.L., Condat, L., Pesquet, J.C., Vu, B.C.: A forwardbackward view of some primal-dual optimization methods in image recovery. In: Proceedings of the 2014 IEEE International Conference on
Image Processing (ICIP), pp. 4141–4145 (2014)
[19] Combettes, P.L., Pesquet, J.C.: Stochastic quasi-Fejér block-coordinate
fixed point iterations with random sweeping. SIAM Journal on Optimization 25(2), 1221–1248 (2015).
[20] Condat, L.: A primal–dual splitting method for convex optimization
involving Lipschitzian, proximable and linear composite terms. Journal
of Optimization Theory and Applications 158(2), 460–479 (2013)
[21] Dang, C.D., Lan, G.: Stochastic block mirror descent methods for nonsmooth and stochastic optimization. SIAM Journal on Optimization
25(2), 856–881 (2015).
[22] Davis, D.: An o(n log(n)) algorithm for projecting onto the ordered
weighted `1 norm ball. arXiv preprint arXiv:1505.00870 (2015)
[23] Davis, D.: Convergence rate analysis of primal-dual splitting schemes.
SIAM Journal on Optimization 25(3), 1912–1943 (2015)
[24] Davis, D., Yin, W.: A three-operator splitting scheme and its optimization applications. arXiv preprint arXiv:1504.01032 (2015)
[25] Dhillon, I.S., Ravikumar, P.K., Tewari, A.: Nearest neighbor based
greedy coordinate descent. In: Advances in Neural Information Processing Systems, pp. 2160–2168 (2011)
[26] Douglas, J., Rachford, H.H.: On the numerical solution of heat conduction problems in two and three space variables. Transactions of the
American Mathematical Society 82(2), 421–439 (1956)
[27] El Baz, D., Gazen, D., Jarraya, M., Spiteri, P., Miellou, J.C.: Flexible
communication for parallel asynchronous methods with application to
a nonlinear optimization problem. Advances in Parallel Computing 12,
429–436 (1998)
[28] Fercoq, O., Bianchi, P.: A coordinate descent primal-dual algorithm
with large step size and possibly non separable functions. arXiv preprint
arXiv:1508.04625 (2015)
[29] Frommer, A., Szyld, D.B.: On asynchronous iterations. Journal of Computational and Applied Mathematics 123(1-2), 201–216 (2000)

Coordinate friendly structures, algorithms, and applications

51

[30] Gabay, D., Mercier, B.: A dual algorithm for the solution of nonlinear
variational problems via finite element approximation. Computers &
Mathematics with Applications 2(1), 17–40 (1976)
[31] Glowinski, R., Marroco, A.: Sur l’approximation, par éléments finis
d’ordre un, et la résolution, par pénalisation-dualité d’une classe de
problèmes de dirichlet non linéaires. Revue française d’automatique,
informatique, recherche opérationnelle. Analyse numérique 9(2), 41–76
(1975)
[32] Grippo, L., Sciandrone, M.: On the convergence of the block nonlinear
Gauss-Seidel method under convex constraints. Operations Research
Letters 26(3), 127–136 (2000)
[33] Han, S.: A successive projection method. Mathematical Programming
40(1), 1–14 (1988)
[34] Hastie, T., Tibshirani, R., Friedman, J., Franklin, J.: The elements of
statistical learning: data mining, inference and prediction. The Mathematical Intelligencer 27(2), 83–85 (2005)
[35] Hildreth, C.: A quadratic programming procedure. Naval Research Logistics Quarterly 4(1), 79–85 (1957)
[36] Hong, M., Wang, X., Razaviyayn, M., Luo, Z.Q.: Iteration complexity analysis of block coordinate descent methods. arXiv preprint
arXiv:1310.6957v2 (2015)
[37] Hsieh, C.j., Yu, H.f., Dhillon, I.: PASSCoDe: Parallel asynchronous
stochastic dual co-ordinate descent. In: Proceedings of the 32nd International Conference on Machine Learning (ICML-15), pp. 2370–2379
(2015)
[38] Jacob, L., Obozinski, G., Vert, J.P.: Group lasso with overlap and graph
lasso. In: Proceedings of the 26th International Conference on Machine
Learning (ICML-09), pp. 433–440. ACM (2009)
[39] Krasnosel’skii, M.A.: Two remarks on the method of successive approximations. Uspekhi Matematicheskikh Nauk 10(1), 123–127 (1955)
[40] Lebedev, V., Tynjanskiı, N.: Duality theory of concave-convex games.
In: Soviet Math. Dokl, vol. 8, pp. 752–756 (1967)
[41] Li, Y., Osher, S.: Coordinate descent optimization for `1 minimization
with application to compressed sensing; a greedy algorithm. Inverse
Problems and Imaging 3(3), 487–503 (2009)

52
[42] Liu, J., Wright, S.J.: Asynchronous stochastic coordinate descent: Parallelism and convergence properties. SIAM Journal on Optimization
25(1), 351–376 (2015)
[43] Liu, J., Wright, S.J., Ré, C., Bittorf, V., Sridhar, S.: An asynchronous
parallel stochastic coordinate descent algorithm. Journal of Machine
Learning Research 16, 285–322 (2015)
[44] Lu, Z., Xiao, L.: On the complexity analysis of randomized blockcoordinate descent methods. Mathematical Programming 152(1-2),
615–642 (2015).
[45] Luo, Z.Q., Tseng, P.: On the convergence of the coordinate descent
method for convex differentiable minimization. Journal of Optimization
Theory and Applications 72(1), 7–35 (1992)
[46] McLinden, L.: An extension of Fenchel’s duality theorem to saddle
functions and dual minimax problems. Pacific Journal of Mathematics
50(1), 135–158 (1974)
[47] Nedić, A., Bertsekas, D.P., Borkar, V.S.: Distributed asynchronous incremental subgradient methods. Studies in Computational Mathematics 8, 381–407 (2001)
[48] Nesterov, Y.: Efficiency of coordinate descent methods on huge-scale
optimization problems. SIAM Journal on Optimization 22(2), 341–362
(2012)
[49] Nutini, J., Schmidt, M., Laradji, I., Friedlander, M., Koepke, H.: Coordinate descent converges faster with the Gauss-Southwell rule than
random selection. In: Proceedings of the 32nd International Conference
on Machine Learning (ICML-15), pp. 1632–1641 (2015)
[50] O’Connor, D., Vandenberghe, L.: Primal-dual decomposition by operator splitting and applications to image deblurring. SIAM Journal on
Imaging Sciences 7(3), 1724–1754 (2014)
[51] Paatero, P., Tapper, U.: Positive matrix factorization: A non-negative
factor model with optimal utilization of error estimates of data values.
Environmetrics 5(2), 111–126 (1994)
[52] Passty, G.B.: Ergodic convergence to a zero of the sum of monotone
operators in Hilbert space. Journal of Mathematical Analysis and Applications 72(2), 383–390 (1979)

Coordinate friendly structures, algorithms, and applications

53

[53] Peaceman, D.W., Rachford Jr, H.H.: The numerical solution of
parabolic and elliptic differential equations. Journal of the Society for
Industrial and Applied Mathematics 3(1), 28–41 (1955)
[54] Peng, Z., Xu, Y., Yan, M., Yin, W.: ARock: an algorithmic framework for asynchronous parallel coordinate updates. ArXiv e-prints
arXiv:1506.02396 (2015)
[55] Peng, Z., Yan, M., Yin, W.: Parallel and distributed sparse optimization. In: Proceedings of the 2013 Asilomar Conference on Signals, Systems and Computers, pp. 659–646 (2013)
[56] Pesquet, J.C., Repetti, A.: A class of randomized primal-dual algorithms for distributed optimization. Journal of Nonlinear and Convex
Analysis 16(12), 2453–2490 (2015)
[57] Polak, E., Sargent, R., Sebastian, D.: On the convergence of sequential
minimization algorithms. Journal of Optimization Theory and Applications 12(6), 567–575 (1973)
[58] Razaviyayn, M., Hong, M., Luo, Z.Q.: A unified convergence analysis
of block successive minimization methods for nonsmooth optimization.
SIAM Journal on Optimization 23(2), 1126–1153 (2013)
[59] Recht, B., Re, C., Wright, S., Niu, F.: Hogwild: A lock-free approach
to parallelizing stochastic gradient descent. In: Advances in Neural
Information Processing Systems, pp. 693–701 (2011)
[60] Repetti, A., Chouzenoux, E., Pesquet, J.C.: A random block-coordinate
primal-dual proximal algorithm with application to 3d mesh denoising.
In: Proceedings of the 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 3561–3565 (2015)
[61] Richtárik, P., Takáč, M.: Iteration complexity of randomized blockcoordinate descent methods for minimizing a composite function. Mathematical Programming 144(1-2), 1–38 (2014)
[62] Richtárik, P., Takáč, M.: Parallel coordinate descent methods for
big data optimization. Mathematical Programming 156(1), 433–484
(2016).
[63] Rockafellar, R.T.: Convex analysis. Princeton University Press (1997)
[64] Rue, H., Held, L.: Gaussian Markov random fields: theory and applications. CRC Press (2005)

54
[65] Scholkopf, B., Smola, A.J.: Learning with kernels: support vector machines, regularization, optimization, and beyond. MIT press (2001)
[66] Siddon, R.L.: Fast calculation of the exact radiological path for a threedimensional CT array. Medical Physics 12(2), 252–255 (1985)
[67] Strikwerda, J.C.: A probabilistic analysis of asynchronous iteration.
Linear Algebra and its Applications 349(1), 125–154 (2002)
[68] Tseng, P.: Applications of a splitting algorithm to decomposition in
convex programming and variational inequalities. SIAM Journal on
Control and Optimization 29(1), 119–138 (1991)
[69] Tseng, P.: On the rate of convergence of a partially asynchronous gradient projection algorithm. SIAM Journal on Optimization 1(4), 603–619
(1991)
[70] Tseng, P.: Dual coordinate ascent methods for non-strictly convex minimization. Mathematical Programming 59(1), 231–247 (1993)
[71] Tseng, P.: A modified forward-backward splitting method for maximal
monotone mappings. SIAM J. Control and Optimization 38(2), 431–
446 (2000).
[72] Tseng, P.: Convergence of a block coordinate descent method for nondifferentiable minimization. Journal of Optimization Theory and Applications 109(3), 475–494 (2001)
[73] Tseng, P., Yun, S.: Block-coordinate gradient descent method for linearly constrained nonsmooth separable optimization. Journal of Optimization Theory and Applications 140(3), 513–535 (2009)
[74] Tseng, P., Yun, S.: A coordinate gradient descent method for nonsmooth separable minimization. Mathematical Programming 117(1-2),
387–423 (2009)
[75] Von Neumann, J.: On rings of operators. reduction theory. Annals of
Mathematics 50(2), 401–485 (1949)
[76] Vũ, B.C.: A splitting algorithm for dual monotone inclusions involving
cocoercive operators. Advances in Computational Mathematics 38(3),
667–681 (2013)
[77] Warga, J.: Minimizing certain convex functions. Journal of the Society
for Industrial and Applied Mathematics 11(3), 588–593 (1963)
[78] Wright, S.J.: Coordinate descent algorithms. Mathematical Programming 151(1), 3–34 (2015)

Coordinate friendly structures, algorithms, and applications

55

[79] Wu, T.T., Lange, K.: Coordinate descent algorithms for lasso penalized
regression. The Annals of Applied Statistics 2(1), 224–244 (2008)
[80] Xu, Y.: Alternating proximal gradient method for sparse nonnegative
Tucker decomposition. Mathematical Programming Computation 7(1),
39–70 (2015).
[81] Xu, Y., Yin, W.: A block coordinate descent method for regularized
multiconvex optimization with applications to nonnegative tensor factorization and completion. SIAM Journal on Imaging Sciences 6(3),
1758–1789 (2013)
[82] Xu, Y., Yin, W.: A globally convergent algorithm for nonconvex
optimization based on block coordinate update.
arXiv preprint
arXiv:1410.1386 (2014)
[83] Xu, Y., Yin, W.: Block stochastic gradient iteration for convex and
nonconvex optimization. SIAM Journal on Optimization 25(3), 1686–
1716 (2015).
[84] Yuan, M., Lin, Y.: Model selection and estimation in regression with
grouped variables. Journal of the Royal Statistical Society: Series B
(Statistical Methodology) 68(1), 49–67 (2006)
[85] Zadeh, N.: A note on the cyclic coordinate ascent method. Management
Science 16(9), 642–644 (1970)

Appendix A. Some Key Concepts of Operators
In this section, we go over a few key concepts in monotone operator theory
and operator splitting theory.
Definition 8 (monotone operator). A set-valued operator T : H ⇒ H is
monotone if hx − y, u − vi ≥ 0, ∀x, y ∈ H, u ∈ T x, v ∈ T y. Furthermore, T
is maximally monotone if its graph Grph(T ) = {(x, u) ∈ H × H : u ∈ T x}
is not strictly contained in the graph of any other monotone operator.
Example 18. An important maximally monotone operator is the subdifferential ∂f of a closed proper convex function f .
Definition 9 (nonexpansive operator). An operator T : H → H is nonexpansive if kT x − T yk ≤ kx − yk, ∀x, y ∈ H. We say T is averaged,
or α-averaged, if there is one nonexpansive operator R such that T =
(1 − α)I + αR for some 0 < α < 1. A 21 -averaged operator T is also called
firmly-nonexpansive.

56
By definition, a nonexpansive operator is single-valued. Let T be averaged. If T has a fixed point, the iteration (2) converges to a fixed point;
otherwise, the iteration diverges unboundedly. Now let T be nonexpansive.
The damped update of T : xk+1 = xk −η(xk −T xk ), is equivalent to applying
the averaged operator (1 − η)I + ηT .
Example 19. A common firmly-nonexpansive operator is the resolvent of
a maximally monotone map T , written as
(67)

JA := (I + A)−1 .

Given x ∈ H, JA (x) = {y : x ∈ y + Ay}. (By monotonicity of A, JA is a
singleton, and by maximality of A, JA (x) is well defined for all x ∈ H. ) A
reflective resolvent is
(68)

RA := 2JA − I.

Definition 10 (proximal map). The proximal map for a function f is a
special resolvent defined as:
(69)


1
proxγf (y) = arg min f (x) +
kx − yk2 ,
2γ
x

where γ > 0. The first-order variational condition of the minimization yields
proxγf (y) = (I + γ∂f )−1 ; hence, proxγf is firmly-nonexpansive. When
x ∈ Rm and proxγf can be computed in O(m) or O(m log m) operations,
we call f proximable.
Examples of proximable functions include `1 , `2 , `∞ -norms, several matrix norms, the owl-norm [22], (piece-wise) linear functions, certain quadratic
functions, and many more.
Example 20. A special proximal map is the projection map. Let X be a
nonempty closed convex set, and ιS be its indicator function. Minimizing
ιS (x) enforces x ∈ S, so proxγιS reduces to the projection map projS for
any γ > 0. Therefore, projS is also firmly nonexpansive.
Definition 11 (β-cocoercive operator). An operator T : H → H is βcocoercive if hx − y, T x − T yi ≥ βkT x − T yk2 , ∀x, y ∈ H.
Example 21. A special example of cocoercive operator is the gradient of a
smooth function. Let f be a differentiable function. Then ∇f is β-Lipschitz
continuous if and only if ∇f is β1 -cocoercive [5, Corollary 18.16].

Coordinate friendly structures, algorithms, and applications

57

Appendix B. Derivation of ADMM from the DRS Update
We derive the ADMM update in (23) from the DRS update
(70a)
(70b)

sk = JηB (tk ),


1
1
k+1
t
=
(2JηA − I) ◦ (2JηB − I) + I (tk ),
2
2

where A = −∂f ∗ (−·) and B = ∂g ∗ .
Note (70a) is equivalent to tk ∈ sk +η∂g ∗ (sk ), i.e., there is a y k ∈ ∂g ∗ (sk )
such that tk = sk + ηy k , so
tk − ηy k = sk ∈ ∂g(y k ).

(71)

In addition, (70b) can be written as
tk+1 = JηA (2sk − tk ) + tk − sk
= sk + (JηA − I)(2sk − tk )
= sk + (I − (I + η∂f ∗ )−1 )(tk − 2sk )
= sk + η(ηI + ∂f )−1 (tk − 2sk )
= sk + η(ηI + ∂f )−1 (ηy k − sk ),

(72)

where in the fourth equality, we have used the Moreau’s Identity [63]: (I +
∂h)−1 + (I + ∂h∗ )−1 = I for any closed convex function h. Let
(73)

1
1
xk+1 = (ηI + ∂f )−1 (ηy k − sk ) = (I + ∂f )−1 (y k − sk ).
η
η

Then (72) becomes
tk+1 = sk + ηxk+1 ,
and
(74)

(71)

sk+1 = tk+1 − ηy k+1 = sk + ηxk+1 − ηy k+1 ,

which together with sk+1 ∈ ∂g(y k+1 ) gives
(75)

1
1
y k+1 = (ηI + ∂g)−1 (sk + ηxk+1 ) = (I + ∂g)−1 (xk+1 + sk ).
η
η

Hence, from (73), (74), and (75), the ADMM update in (23) is equivalent
to the DRS update in (70) with η = γ1 .

58

Appendix C. Representing the Condat-Vũ Algorithm as a
Nonexpansive Operator
We show how to derive the Condat-Vũ algorithm (28) by applying a forwardbackward operator to the optimality condition (27):

 
 
  
∇f 0
∂g 0
x
0 A>
,
+
+
0 0
0 ∂h∗
s
−A 0
| {z } |
{z
} |{z}
z
operator A
operator B

 
(76)

0∈

 
x
It can be written as 0 ∈ Az + Bz after we define z =
. Let M be a
s
symmetric positive definite matrix, we have
0 ∈ Az + Bz
⇔M z − Az ∈ M z + Bz
⇔z − M −1 Az ∈ z + M −1 Bz
⇔z = (I + M −1 B)−1 ◦ (I − M −1 A)z.
Convergence and other results can be found in [23]. The last equivalent
relation is due to M −1 B being a maximally monotone operator under the
norm induced by M . We let
"
#
1
>
ηI A
M=
0
A γ1 I
and iterate
z k+1 = T z k = (I + M −1 B)−1 ◦ (I − M −1 A)z k .
We have M z k+1 + Bz k+1 = M z k − Az k :
(
1 k
1 k+1
> k
k
+ A> sk+1 + A> sk+1 + ∂g(xk+1 ),
η x + A s − ∇f (x ) ∈ η x
1 k
1 k+1
k
∈ γs
+ A xk+1 − A xk+1 + ∂h∗ (sk+1 ),
γs + A x
which is equivalent to
 k+1
s
= proxγh∗ (sk + γAxk ),
xk+1 = proxηg (xk − η(∇f (xk ) + A> (2sk+1 − sk ))).

Coordinate friendly structures, algorithms, and applications

59

Now we derived the Condat-Vũ algorithm. With proper choices of η and γ,
the forward-backward operator T = (I + M −1 B)−1 ◦ (I − M −1 A) can be
shown to be α-averaged if we use the inner product
hz1 , z2 iM = z1> M z2 and
 
√
x
norm kzkM = z > M z on the space of z =
. More details can be found
s
in [23].
"
#
1
>
I
−A
If we change the matrix M to η
, the other algorithm (29)
1
−A
γI
can be derived similarly.

Appendix D. Proof of Convergence for Async-parallel
Primal-dual Coordinate Update Algorithms
Algorithms 1 and 2 differ from that in [54] in the following aspects:
1. the operator TCV is nonexpansive under a norm induced by a symmetric positive definite matrix M (see Appendix C), instead of the
standard Euclidean norm;
2. the coordinate updates are no longer orthogonal to each other under
the norm induced by M ;
3. the block coordinates may overlap each other.
Because of these differences, we make two major modifications to the proof
in [54, Section 3]: (i) adjusting parameters in [54, Lemma 2] and modify
its proof to accommodate for the new norm; (2) modify the inner product
and induced norm used in [54, Theorem 2] and adjust the constants in [54,
Theorems 2 and 3].
We assume the same inconsistent case as in [54], i.e., the relationship
between ẑ k and z k is
X
(77)
ẑ k = z k +
(z d − z d+1 ),
d∈J(k)

where J(k) ⊆ {k − 1, ..., k − τ } and τ is the maximum number of other
updates to z during the computation of the update. Let S = I − TCV .
ηk
Then the coordinate update can be rewritten as z k+1 = z k − (m+p)q
Sik ẑ k ,
i
k

k , (S ẑ k ) , ẑ k , . . . , ẑ k
where Si ẑ k = (ẑ1k , . . . , ẑi−1
i i+1
m+p ) for Algorithm 1. For Algorithm 2, the update is

(78)

z k+1 = z k −

ηk
Si ẑ k ,
mqik k

60
where


0
 ..

.


0


IHi



Si ẑ k = 


















 k
 S ẑ .









0
..

.
0
ρi,1 IG1
..

.
ρi,p IGp

Let λmax and λmin be the maximal and minimal eigenvalues of the matrix
be the condition number. Then we have the
M , respectively, and κ = λλmax
min
following lemma.
Lemma 1. For both Algorithms 1 and 2,
X
(79)
Si ẑ k = S ẑ k ,
i

X

(80)

kSi ẑ k k2M ≤ κkS ẑ k k2M ,

i

where i runs from 1 to m + p for Algorithm 1 and 1 to m for Algorithm 2.
Proof. The first part comes immediately from the definition of S for both
algorithms. For the second part, we have
(81)

X

kSi ẑ k k2M ≤

i

X

λmax kSi ẑ k k2 = λmax kS ẑ k k2 ≤

i

λmax
kS ẑ k k2M ,
λmin

for Algorithm 1. For Algorithm 2, the equality is replaced by “≤”.
At last we define
(82)

z̄ k+1 := z k − ηk S ẑ k ,

qmin = mini qi > 0, and |J(k)| be the number of elements in J(k). It is shown
in [23] that with proper choices of η and γ, TCV is nonexpansive under the
norm induced by M . Then Lemma 2 shows that S is 1/2-cocoercive under
the same norm.

Coordinate friendly structures, algorithms, and applications

61

Lemma 2. An operator T : F → F is nonexpansive under the induced norm
by M if and only if S = I − T is 1/2-cocoercive under the same norm, i.e.,
(83)

1
hz − z̃, Sz − S z̃iM ≥ kSz − S z̃k2M ,
2

∀ z, z̃ ∈ F.

The proof is the same as that of [5, Proposition 4.33].
We state the complete theorem for Algorithm 2. The theorem for Algorithm 1 is similar (we need to change m to m + p when necessary).
Theorem 2. Let Z ∗ be the set of optimal solutions of (24) and (z k )k≥0 ⊂ F
be the sequence generated by Algorithm 2 (with proper choices of η and γ
such that TCV is nonexpansive under the norm induced by M ), under the
following conditions:
(i) f, g, h∗ are closed proper convex functions. In addition, f is differentiable and ∇f is Lipschitz continuous with β;
min
(ii) ηk ∈ [ηmin , ηmax ] for certain 0 < ηmax < 2τ √mq
κqmin +κ and any 0 < ηmin ≤
ηmax .
Then (z k )k≥0 converges to a Z ∗ -valued random variable with probability 1.
The proof directly follows [54, Section 3]. Here we only present the key
modifications. Interested readers are referred to [54] for the complete procedure.
The next lemma shows that the conditional expectation of the distance
between z k+1 and any z ∗ ∈ FixTCV = Z ∗ for given Z k = {z 0 , z 1 , · · · , z k }
has an upper bound that depends on Z k and z ∗ only.
Lemma 3. Let (z k )k≥0 be the sequence generated by Algorithm 2. Then for
any z ∗ ∈ FixTCV , we have

(84)


σ X
E kz k+1 − z ∗ k2M Z k ≤kz k − z ∗ k2M +
kz d − z d+1 k2M
m
d∈J(k)


κ
1
1 |J(k)|
+
+
−
kz k − z̄ k+1 k2M
m
σ
mqmin ηk

where E(· | Z k ) denotes conditional expectation on Z k and σ > 0 (to be
optimized later).

62
Proof. We have

(85)



E kz k+1 − z ∗ k2M | Z k


(78)
ηk
k − z ∗ k2 | Z k
= E kz k − mq
ẑ
S
i
k
M
ik

2ηk
=kz k − z ∗ k2M + E mq
Sik ẑ k , z ∗ − z k
ik
Pm
k ∗
k
k
=kz k − z ∗ k2M + 2η
i=1 Si ẑ , z − z
m
=kz k − z ∗ k2M +

2ηk
m

S ẑ k , z ∗ − z k

M

+

M

+

ηk2
k 2
m2 qi2k kSik ẑ kM
2

M

ηk2
m2

ηk
+m
2
Pm

Zk



Pm

1
k 2
i=1 qi kSi ẑ kM

1
k 2
i=1 qi kSi ẑ kM ,

where the third equality holds because the probability of choosing i is qi .
Note that

Pm

1
k 2
i=1 qi kSi ẑ kM

(86)

≤

m
1 X

qmin

kSi ẑ k k2M

i=1

m
κ X

(80)

≤

qmin

kS ẑ k k2M

i=1

κ
kz k − z̄ k+1 k2M ,
= 2
ηk qmin

(82)

and
(87)
hS ẑ k , z ∗ − z k iM
=hS ẑ k , z ∗ − ẑ k +

P

d∈J(k) (z

(82)

d

− z d+1 )iM

= hS ẑ k , z ∗ − ẑ k iM +
≤hS ẑ k − Sz ∗ , z ∗ −

(83)

≤ − 21 kS ẑ k k2M +

(82)

= −

1
k
2ηk2 kz

1 P
k
k+1 , z d − z d+1 i
M
d∈J(k) hz − z̄
ηk
1 P
1
k
k+1
k
k2M
ẑ iM + 2ηk d∈J(k) σ kz − z̄

1
2ηk

+ σkz d − z d+1 k2M



1
k
d∈J(k) ( σ kz

− z̄ k+1 k2M + σkz d − z d+1 k2M )
P
k
k+1 k2 + σ
d
d+1 k2 ,
+ |J(k)|
M
M
d∈J(k) kz − z
2σηk kz − z̄
2ηk

P

− z̄ k+1 k2M

where the first inequality follows from the Young’s inequality. Plugging (86)
and (87) into (85) gives the desired result.
Let Fτ +1 =
product:

Qτ

i=0 F

be a product space and h· | ·i be the induced inner

h(z 0 , . . . , z τ ) | (z̃ 0 , . . . , z̃ τ )i =

τ
X
i=0

hz i , z̃ i iM ,

∀(z 0 , . . . , z τ ), (z̃ 0 , . . . , z̃ τ ) ∈ Fτ +1 .

Coordinate friendly structures, algorithms, and applications

63

Define a (τ + 1) × (τ + 1) matrix U 0 by



1 0 ···
0 0 · · ·

U 0 :=  . . .
..
 .. ..
0 0 ···


0
r
0
qmin

..  +
κ

.
0

τ
−τ








−τ
2τ − 1 1 − τ
1 − τ 2τ − 3 2 − τ
..
..
.
.
−2






,
..

.

3 −1
−1 1

and let U = U 0 ⊗ IF . Here ⊗ represents the Kronecker product. For a given
(y 0 , · · · , y τ ) ∈ Fτ +1 , (z 0 , · · · , z τ ) = U (y 0 , · · · , y τ ) is given by:

p
0
1
z 0 = y 0 + τ qmin
κ (y − y ),
p
i−1 + (2τ − 2i + 1)y i + (i − τ )y i+1 ), if 1 ≤ i ≤ τ − 1,
z i = qmin
κ ((i − τ − 1)y
p
τ
τ −1 ).
z τ = qmin
κ (y − y
Then U is a self-adjoint and positive definite linear operator since U 0 is
symmetric and positive definite, and we define h· | ·iU = h· | U ·i as the U weighted inner product and k · kU the induced norm.
Let
zk = (z k , z k−1 , . . . , z k−τ ) ∈ Fτ +1 , k ≥ 0, z∗ = (z ∗ , . . . , z ∗ ) ∈ Z∗ ⊆ Fτ +1 ,
where z k = z 0 for k < 0. With
(88)
p
Pk−1
i
i+1 k2 ,
ξk (z∗ ) := kzk −z∗ k2U = kz k −z ∗ k2M + qmin
M
i=k−τ (i−(k−τ )+1)kz −z
κ
we have the following fundamental inequality:
Theorem 3 (fundamental inequality). Let (z k )k≥0 be the sequence generated
by Algorithm 2. Then for any z∗ ∈ Z∗ , it holds that



∗

E ξk+1 (z ) Z

k



1
≤ ξk (z ) +
m
∗




√
2τ κ
κ
1
+
−
kz̄ k+1 − z k k2M .
√
m qmin mqmin ηk

64

p
Proof. Let σ = m qmin
κ . We have
E(ξk+1 (z∗ )|Z k )
(88)

= E(kz k+1 − z ∗ k2M |Z k ) + σ

Pk

i−(k−τ )
E(kz i
m

− z i+1 k2M |Z k )
Pk−1
(78)
ηk2
i−(k−τ )
k 2
k
kz i − z i+1 k2M
= E(kz k+1 − z ∗ k2M |Z k ) + στ
i=k+1−τ
m E( m2 qi2k kSik ẑ kM |Z ) + σ
m
P
i−(k−τ )
κ
kz k − z̄ k+1 k2M + σ k−1
≤E(kz k+1 − z ∗ k2M |Z k ) + mστ
kz i − z i+1 k2M
3q
i=k+1−τ
m
min


(84)
|J(k)|
κ
1
στ κ
1
k
k+1 k2
+
−
+
≤ kz k − z ∗ k2M + m
2
M
σ
m qmin
mqmin
ηk kz − z̄
Pk−1
i−(k−τ )
σ P
d
d+1
2
i
i+1
kz − z k2M
+ m d∈J(k) kz − z kM + σ i=k+1−τ
m


1
τ
στ κ
κ
1
k
k+1 k2
≤kz k − z ∗ k2M + m
+
+
−
2
M
σ
m qmin
mqmin
ηk kz − z̄
Pk−1
i−(k−τ )
σ Pk−1
i
i+1
2
i
+ m i=k−τ kz − z kM + σ i=k+1−τ
kz − z i+1 k2M
m

 √
(88)
2τ κ
κ
1
1
k
k+1 k2 .
√
+
−
= ξk (x∗ ) + m
M
m qmin
mqmin
ηk kz − z̄
i=k+1−τ

The first inequality follows from the computation of the conditional expectation on Z k and (86), the third inequality holds because J(k) ⊂ {k − 1, k −
p
2, · · · , k − τ }, and the last equality uses σ = m qmin
κ , which minimizes
στ κ
τ
σ + m2 qmin over σ > 0. Hence, the desired inequality holds.

Zhimin Peng
PO Box 951555
UCLA Math Department
Los Angeles, CA 90095
E-mail address: zhimin.peng@math.ucla.edu
Tianyu Wu
PO Box 951555
UCLA Math Department
Los Angeles, CA 90095
E-mail address: wuty11@math.ucla.edu
Yangyang Xu
207 Church St SE
University of Minnesota, Twin Cities
Minneapolis, MN 55455
E-mail address: yangyang@ima.umn.edu

Coordinate friendly structures, algorithms, and applications

65

Ming Yan
Department of Computational Mathematics, Science and Engineering
Department of Mathematics
Michigan State University
East Lansing, MI 48824
E-mail address: yanm@math.msu.edu
Wotao Yin
PO Box 951555
UCLA Math Department
Los Angeles, CA 90095
E-mail address: wotaoyin@math.ucla.edu
Received January 1, 2016

