Coordinate Friendly Structures, Algorithms and Applications

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.

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, 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

such that the limit of the sequence $\{x^{k}\}$ exists and is a fixed point of ${\mathcal{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 $x_{i}$ 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 ${\mathcal{T}}$ that make the following coordinate update algorithm computationally worthy

where $\eta_{k}$ is a step size and $i\in[m]:=\{1,\ldots,m\}$ is arbitrary. Specifically, the cost of performing (3) is roughly $\frac{1}{m}$, or lower, of that of performing (2). We call such ${\mathcal{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 second-order cone programming, variational image processing, support vector machine, empirical risk minimization, portfolio optimization, distributed computing, and nonnegative matrix factorization. For each application, we present 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 $\infty$, and all sets $X,Y,Z$ are nonempty closed convex. The indicator function $\iota_{X}(x)$ returns if $x\in X$, and $\infty$ elsewhere. For a positive integer $m$, we let $[m]:=\{1,\ldots,m\}$.

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

update $x^{k+1}_{i}$ for $i={i_{k}}$ while keeping $x_{i}^{k+1}=x_{i}^{k}$, $\forall\,i\not={i_{k}}$;

Next we review the index rules and the methods to update $x_{i}$.

In this framework, there is a sequence of coordinate indices $i_{1},i_{2},\ldots$ chosen according to one of the following rules: cyclic, cyclic permutation, random, and greedy rules. At iteration $k$, only the $i_{k}$th coordinate is updated:

Sequential updates have been applied to many problems such as the Gauss-Seidel iteration for solving a linear system of equations, alternating projection for finding a point in the intersection of two sets, ADMM for solving monotropic programs, and Douglas-Rachford Splitting (DRS) for finding a zero to the sum of two operators.

In optimization, coordinate descent algorithms, at each iteration, minimize the function $f(x_{1},\ldots,x_{m})$ by fixing all but one variable $x_{i}$. Let

collect all but the $i$th coordinate of $x$. Coordinate descent solves one of the following subproblems:

which are called direct update, proximal update, gradient update, and prox-gradient update, respectively. The last update applies to the function

Sequential-update literature. Coordinate descent algorithms date back to the 1950s , 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 . Proximal updates are studied in and developed into prox-gradient updates in and mixed updates in .

The random index rule first appeared in and then . Recently, 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 .

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 and the save in the total number of iterations significantly outweighs the extra calculation .

A simple example. We present the coordinate update algorithms under different index rules for solving a simple least squares problem:

The four tested index rules are: cyclic, cyclic permutation, random, and greedy under the Gauss-Southwellit selects $i_{k}=\operatorname*{arg\,max}_{i}\|\nabla_{i}f(x^{k})\|$. 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 $\eta_{k}$ is set to the theoretical upper bound $\frac{2}{\|A\|_{2}^{2}}$, where $\|A\|_{2}$ denotes the matrix operator norm and equals the largest singular value of $A$. For each coordinate update to $x_{i}$, the step size $\eta_{k}$ is set to $\frac{1}{(A^{\top}A)_{ii}}$. All of the full and coordinate updates have the same per-epoch complexity, so we plot the objective errors in Figure 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.

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 $x_{i}$ back to the shared memory (or through communicating with other agents without shared memory):

Unlike before, $k$ increases whenever any agent completes an update.

The lack of synchronization often results in computation with out-of-date information. During the computation of the $k$th update, other agents make $d_{k}$ updates to $x$ in the shared memory; when the $k$th update is written, its input is already $d_{k}$ iterations out of date. This number is referred to as the asynchronous delay. In (5), the agent reads $x^{k-d_{k}}$ and commits the update to $x_{i_{k}}^{k}$. Here we have assumed consistent reading, i.e., $x^{k-d_{k}}$ lying in the set $\{x^{j}\}_{j=1}^{k}$. 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 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 for systems of linear equations. For function minimization, introduced an async-parallel gradient projection method. Convergence rates are obtained in . Recently, developed parallel randomized methods.

For fixed-point problems, async-parallel methods date back to in 1978. In the pre-2010 methods and the review , each agent updates its own subset of coordinates. Convergence is established under the $P$-contraction condition and its variants . Papers show convergence for async-parallel iterations with simultaneous reading and writing to the same set of components. Unbounded but stochastic delays are considered in .

Recently, random coordinate selection appeared in for fixed-point problems. The works introduced async-parallel stochastic methods for function minimization. For fixed-point problems, introduced async-parallel stochastic methods, as well as several applications.

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 coordinate-update 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 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.

Coordinate Friendly Operators

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:

Note that $x^{+}_{j}=x_{j}$ for all $j\not=i$. Hence, $x^{+}=(x_{1},\ldots,x_{i}+\delta_{i},\ldots,x_{m})$.

We let $\mathfrak{M}\left[{a}\mapsto{b}\right]$ denote the number of basic operations that it takes to compute the quantity $b$ from the input $a$.

For example, $\mathfrak{M}\left[{x}\mapsto{({\mathcal{T}}x)_{i}}\right]$ denotes the number of operations to compute the $i$th component of ${\mathcal{T}}x$ given $x$. We explore the possibility to compute $({\mathcal{T}}x)_{i}$ with much fewer operations than what is needed to first compute ${\mathcal{T}}x$ and then take its $i$th component.

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 $A_{i,:}$ and $A_{:,j}$ be the $i$th row and $j$th column of a matrix $A$, respectively. Let $A^{\top}$ be the transpose of $A$ and $A^{\top}_{i,:}$ be $(A^{\top})_{i,:}$, i.e., the $i$th row of the transpose of $A$.

Assuming that $A^{\top}A$ and $A^{\top}b$ are already computed, we have $\mathfrak{M}\left[{x}\mapsto{{\mathcal{T}}x}\right]=O(m^{2})$. The coordinate update at the $k$th iteration performs

and $x_{j}^{k+1}=x_{j}^{k},\forall j\neq i_{k}$, where $i_{k}$ is some selected coordinate.

Since for all $i$, $\nabla_{i}f(x^{k})=\left(A^{\top}(Ax-b)\right)_{i}=(A^{\top}A)_{i,:}\cdot x-(A^{\top}b)_{i}$, we have $\mathfrak{M}\left[{x}\mapsto{({\mathcal{T}}x)_{i}}\right]=O(m)$ and thus $\mathfrak{M}\left[{x}\mapsto{({\mathcal{T}}x)_{i}}\right]=O(\frac{1}{m}\mathfrak{M}\left[{x}\mapsto{{\mathcal{T}}x}\right])$. Therefore, the coordinate gradient descent is computationally worthy.

The operator ${\mathcal{T}}$ in the above example is a special Type-I CF operator.

We can implement the coordinate update in Example 1 in a different manner by maintaining the result ${\mathcal{T}}x^{k}$ in the memory. This approach works when $m=\Theta(p)$ or $p\gg m$. The full update (8) is unchanged. At each coordinate update, from the maintained quantity ${\mathcal{T}}x^{k}$, we immediately obtain $x_{i_{k}}^{k+1}=({\mathcal{T}}x^{k})_{i_{k}}$. But we need to update ${\mathcal{T}}x^{k}$ to ${\mathcal{T}}x^{k+1}$. Since $x^{k+1}$ and $x^{k}$ differ only over the coordinate $i_{k}$, this update can be computed as

which is a scalar-vector multiplication followed by vector addition, taking only $O(m)$ operations. Computing ${\mathcal{T}}x^{k+1}$ from scratch involves a matrix-vector multiplication, taking $O(\mathfrak{M}\left[{x}\mapsto{{\mathcal{T}}(x)}\right])=O(m^{2})$ operations. Therefore,

The operator ${\mathcal{T}}$ in the above example is a special Type-II CF operator.

An operator ${\mathcal{T}}$ is called Type-II CF (denoted as ${\mathcal{F}}_{2}$) if, for any $i,x$ and x^{+}:=\big{(}x_{1},\ldots,({\mathcal{T}}x)_{i},\ldots,x_{m}\big{)}, the following holds

The next example illustrates an efficient coordinate update by maintaining certain quantity other than ${\mathcal{T}}x$.

For the case $p\ll m$, we should avoid pre-computing the relative large matrix $A^{\top}A$, and it is cheaper to compute $A^{\top}(Ax)$ than $(A^{\top}A)x$. Therefore, we change the implementations of both the full and coordinate updates in Example 1. In particular, the full update

pre-multiplies $x^{k}$ by $A$ and then $A^{\top}$. Hence, $\mathfrak{M}\left[{x^{k}}\mapsto{{\mathcal{T}}(x^{k})}\right]=O(mp)$.

We change the coordinate update to maintain the intermediate quantity $Ax^{k}$. In the first step, the coordinate update computes

by pre-multiplying $Ax^{k}$ by $A^{\top}_{i_{k},:}$. Then, the second step updates $Ax^{k}$ to $Ax^{k+1}$ by adding $(x^{k+1}_{i_{k}}-x^{k}_{i_{k}})A_{:,i_{k}}$ to $Ax^{k}$. Both steps take $O(p)$ operations, so

Combining Type-I and Type-II CF operators with the last example, we arrive at the following CF definition.

where ${\mathcal{M}}(x)$ is some quantity maintained in the memory to facilitate each coordinate update and refreshed to ${\mathcal{M}}(x^{+})$. ${\mathcal{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 ${\mathcal{M}}(x)$ to ${\mathcal{M}}(x^{+})$) while the right-hand side measures the average per-coordinate cost of updating all the coordinates together. When (10) holds, ${\mathcal{T}}$ is amenable to coordinate updates.

By definition, a Type-I CF operator ${\mathcal{T}}$ is CF without maintaining any quantity, i.e., ${\mathcal{M}}(x)=\emptyset$.

A Type-II CF operator ${\mathcal{T}}$ satisfies (10) with ${\mathcal{M}}(x)={\mathcal{T}}x$, so it is also CF. Indeed, given any $x$ and $i$, we can compute $x^{+}$ by immediately letting $x^{+}_{i}=({\mathcal{T}}x)_{i}$ (at $O(1)$ cost) and keeping $x^{+}_{j}=x_{j},\,\forall j\neq i$; then, by (9), we update ${\mathcal{T}}x$ to ${\mathcal{T}}x^{+}$ at a low cost. Formally, letting ${\mathcal{M}}(x)={\mathcal{T}}x$,

In general, the set of CF operators is much larger than the union of Type-I and Type-II CF operators.

non-separable operator: ${\mathcal{C}}_{3}:=\mathfrak{T}\setminus({\mathcal{C}}_{1}\cup{\mathcal{C}}_{2})$. If ${\mathcal{T}}\in{\mathcal{C}}_{3}$, there exists some $i$ such that $({\mathcal{T}}x)_{i}$ depends on many coordinates of $x$.

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.

Clearly ${\mathcal{T}}:x\mapsto Ax$ is separable.

Then, both $\nabla f$ and $\mathbf{prox}_{\gamma f}$ are separable, in particular,

Here, $\mathbf{prox}_{\gamma f}(x)$ ($\gamma>0$) is the proximal operator that we define in Definition 10 in Appendix A.

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.

Let $E$ be a class of index sets and every $e\in E$ be a small subset of $[m]$, $|e|\ll m$. In addition $\#\{e:i\in e\}\ll\#\{e\}$ for all $i\in[m]$. Let $x_{e}:=(x_{i})_{i\in e}$, and

The gradient map $\nabla f$ is nearly-separable.

An application of this example arises in wireless communication over a graph of $m$ nodes. Let each $x_{i}$ be the spectrum assignment to node $i$, each $e$ be a neighborhood of nodes, and each $f_{e}$ be a utility function. The input of $f_{e}$ is $x_{e}$ 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 .

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

Let us maintain ${\mathcal{M}}(x)=a^{\top}x$. For arbitrary $x$ and $i$, let

and $x^{+}_{j}:=x_{j},\,\forall j\neq i$. Then, computing $x^{+}_{i}$ from $x$ and $a^{\top}x$ takes $O(1)$ (as $a^{\top}x$ is maintained), and computing $a^{\top}x^{+}$ from $x^{+}_{i}-x_{i}$ and $a^{\top}x$ costs $O(1)$. Formally, we have

On the other hand, $\mathfrak{M}\left[{x}\mapsto{{\mathcal{T}}x}\right]=O(m)$. Therefore, (10) holds, and ${\mathcal{T}}$ defined in (11) is CF.

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 $x^{k}$ to $x^{k+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.

We start by an example with numerous applications. It is a generalization of Example 9.

Assume that evaluating $\phi^{\prime}_{j}$ costs $O(1)$ for each $j$. Then, $\nabla f$ is CF. Indeed, let

Since $\mathfrak{M}\left[{x}\mapsto{\nabla f(x)}\right]=O(pm)$, therefore $\nabla f={\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}\circ{\mathcal{T}}_{3}$ is CF.

There are general rules to preserve Type-I and Type-II CF. For example, ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$ is still Type-I CF, and ${\mathcal{T}}_{2}\circ{\mathcal{T}}_{3}$ is still CF, but there are counter examples where ${\mathcal{T}}_{2}\circ{\mathcal{T}}_{3}$ 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 splitting schemes in §3.2 below will be based on ${\mathcal{T}}_{1}+{\mathcal{T}}_{2}$ or ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$, as well as a sequence of such combinations. If ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{2}$ are both CF, ${\mathcal{T}}_{1}+{\mathcal{T}}_{2}$ remains CF, but ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$ is not necessarily so. This subsection discusses how ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$ inherits the properties from ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{2}$. Our results are summarized in Tables 1 and 2 and explained in detail below.

The combination ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$ generally inherits the weaker property from ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{2}$.

The separability (${\mathcal{C}}_{1}$) property is preserved by composition. If ${\mathcal{T}}_{1},\ldots,{\mathcal{T}}_{n}$ are separable, then ${\mathcal{T}}_{1}\circ\cdots\circ{\mathcal{T}}_{n}$ is separable. However, combining nearly-separable (${\mathcal{C}}_{2}$) 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 ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$ inherits the CF properties from ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{2}$. For simplicity, we only use matrix-vector multiplication as examples to illustrate the ideas; more interesting examples will be given later.

If ${\mathcal{T}}_{1}$ is separable or nearly-separable (${\mathcal{C}}_{1}\cup{\mathcal{C}}_{2}$), then as long as ${\mathcal{T}}_{2}$ is CF (${\mathcal{F}}$), ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$ remains CF. In addition, if ${\mathcal{T}}_{2}$ is Type-I CF (${\mathcal{F}}_{1}$), so is ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$.

Assume that ${\mathcal{T}}_{2}$ is separable (${\mathcal{C}}_{1}$). It is easy to see that if ${\mathcal{T}}_{1}$ is CF (${\mathcal{F}}$), then ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$ remains CF. In addition if ${\mathcal{T}}_{1}$ is Type-II CF (${\mathcal{F}}_{2}$), so is ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$; see Example 10.

Note that, if ${\mathcal{T}}_{2}$ is nearly-separable, we do not always have CF properties for ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$. This is because ${\mathcal{T}}_{2}x$ and ${\mathcal{T}}_{2}x^{+}$ can be totally different (so updating ${\mathcal{T}}_{2}x$ is expensive) even if $x$ and $x^{+}$ only differ over one coordinate; see the footnote 3 on Page 3.

Assume that ${\mathcal{T}}_{1}$ is Type-I CF (${\mathcal{F}}_{1}$). If ${\mathcal{T}}_{2}$ is Type-II CF (${\mathcal{F}}_{2}$), then ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$ is CF (${\mathcal{F}}$).

Assume that one of ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{2}$ is cheap. If ${\mathcal{T}}_{2}$ is cheap, then as long as ${\mathcal{T}}_{1}$ is Type-I CF (${\mathcal{F}}_{1}$), ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$ is Type-I CF. If ${\mathcal{T}}_{1}$ is cheap, then as long as ${\mathcal{T}}_{2}$ is Type-II CF (${\mathcal{F}}_{2}$), ${\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$ is CF (${\mathcal{F}}$); see Example 13.

We will see more examples of the above cases in the rest of the paper.

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 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 ${\mathcal{A}}$, the resolvent operator ${\mathcal{J}}_{{\mathcal{A}}}$ and the reflective-resolvent operator ${\mathcal{R}}_{{\mathcal{A}}}$ are also defined there, in (67) and (68), respectively.

Consider the following problem: given three operators ${\mathcal{A}},{\mathcal{B}},{\mathcal{C}}$, possibly set-valued,

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).

Indeed, by setting $\gamma\in(0,2\beta)$, ${\mathcal{T}}_{3S}$ is $(\frac{2\beta}{4\beta-\gamma})$-averaged (think it as a property weaker than the Picard contraction; in particular, ${\mathcal{T}}$ may not have a fixed point). Following the standard convergence result (cf. textbook ), provided that ${\mathcal{T}}$ has a fixed point, the sequence from (2) converges to a fixed-point $x^{*}$ of ${\mathcal{T}}$. Note that, instead of $x^{*}$, ${\mathcal{J}}_{\gamma{\mathcal{B}}}(x^{*})$ is a solution to (12).

for solving the problem $0\in{\mathcal{A}}x+{\mathcal{C}}x$.

introduced in for solving the problem $0\in{\mathcal{A}}x+{\mathcal{B}}x$. A more general splitting is the Relaxed Peaceman-Rachford Splitting (RPRS) with $\lambda\in$:

Forward-Douglas-Rachford Splitting (FDRS): Let $V$ be a linear subspace, and ${\mathcal{N}}_{V}$ and ${\mathcal{P}}_{V}$ be its normal cone and projection operator, respectively. The FDRS

where $X$ is the feasible set and $f$ and $g$ are objective functions. We present examples of operator splitting methods discussed above.

A special case of (20) with $g=\iota_{X}$ is the projected gradient iteration:

If $\nabla f$ is CF and $\mathbf{prox}_{\gamma g}$ is (nearly-)separable (e.g., $g(x)=\|x\|_{1}$ or the indicator function of a box constraint) or if $\nabla f$ is Type-II CF and $\mathbf{prox}_{\gamma g}$ is cheap (e.g., $\nabla f(x)=Ax-b$ and $g=\|x\|_{2}$), then the FBS iteration (20) is CF. In the latter case, we can also apply the BFS iteration (15) (i.e, compute $\mathbf{prox}_{\gamma g}$ and then perform the gradient update), which is also CF.

(The iteration can be generalized to handle the constraint $Ax-By=b$.) The dual problem of (22) is $\min_{s}f^{*}(-s)+g^{*}(s)$, where $f^{*}$ is the convex conjugate of $f$. Letting ${\mathcal{A}}=-\partial f^{*}(-\cdot)$ and ${\mathcal{B}}=\partial g^{*}$ in (16) recovers the iteration (23) through (see the derivation in Appendix B)

From the results in §3.1, a sufficient condition for the above iteration to be CF is that ${\mathcal{J}}_{\gamma{\mathcal{A}}}$ is (nearly-)separable and ${\mathcal{J}}_{\gamma{\mathcal{B}}}$ being CF.

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

Primal-dual Coordinate Friendly Operators

Let $u^{0}$ be an image, where $u_{i}^{0}\in$, and $B$ be the blurring linear operator. Let $\|\nabla u\|_{1}$ be the anisotropicGeneralization to the isotropic case is straightforward by grouping variables properly. total variation of $u$ (see (49) for definition). Suppose that $b$ is a noisy observation of $Bu^{0}$. Then, we can try to recover $u^{0}$ by solving

which can be written in the form of $\eqref{pdproblem}$ with $f=\frac{1}{2}\|B\cdot-b\|^{2}$, $g=\iota_{}$, $A=\nabla$, and $h=\lambda\|\cdot\|_{1}$.

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 $\|\nabla u\|_{1}$.

In many applications, although $h$ is proximable, $h\circ A$ is generally non-proximable and non-differentiable. To avoid using slow subgradient methods, we can consider the primal-dual splitting approaches to separate $h$ and $A$ so that $\mathbf{prox}_{h}$ can be applied. We derive that the equivalent form (for convex cases) of $\eqref{pdproblem}$ is to find $x$ such that

Problem $\eqref{pdkkt}$ can be solved by the Condat-Vũ algorithm :

Switching the orders of $x$ and $s$ yields the following algorithm:

It is known from that both $\eqref{vucondat}$ and $\eqref{vucondat2}$ reduce to iterations of nonexpansive operators (under a special metric), i.e., ${{\mathcal{T}}_{\textnormal{CV}}}$ is nonexpansive; see Appendix C for the reasoning.

Similar primal-dual algorithms can be used to solve other problems such as saddle point problems and variational inequalities . Our coordinate update algorithms below apply to these problems as well.

In this subsection, we make the following assumption.

Functions $g$ and $h^{*}$ in the problem (24) are separable and proximable. Specifically,

when $p=O(m)$, the Condat-Vu operator ${{\mathcal{T}}_{\textnormal{CV}}}$ in (28) is CF, more specifically,

when $m\ll p$ and $\mathfrak{M}\left[{x}\mapsto{\nabla f(x)}\right]=O(m)$, the Condat-Vu operator ${\mathcal{T}}^{\prime}_{\textnormal{CV}}$ in (29) is CF, more specifically,

Computing $z^{k+1}={{\mathcal{T}}_{\textnormal{CV}}}z^{k}$ involves evaluating $\nabla f$, $\mathbf{prox}_{g}$, and $\mathbf{prox}_{h^{*}}$, applying $A$ and $A^{\top}$, and adding vectors. It is easy to see $\mathfrak{M}\left[{z^{k}}\mapsto{{{\mathcal{T}}_{\textnormal{CV}}}z^{k}}\right]=O(mp+m+p)+\mathfrak{M}[x\to\nabla f(x)]$, and $\mathfrak{M}\left[{z^{k}}\mapsto{{\mathcal{T}}^{\prime}_{\textnormal{CV}}z^{k}}\right]$ is the same. (a) We assume $\nabla f\in{\mathcal{F}}_{1}$ for simplicity, and other cases are similar.

If $({{\mathcal{T}}_{\textnormal{CV}}}z^{k})_{j}=s^{k+1}_{i}$, computing it involves: adding $s^{k}_{i}$ and $\gamma(Ax^{k})_{i}$, and evaluating $\mathbf{prox}_{\gamma h^{*}_{i}}$. In this case $\mathfrak{M}\left[{\{z^{k},Ax\}}\mapsto{\{z^{+},Ax^{+}\}}\right]=O(1)$.

If $({{\mathcal{T}}_{\textnormal{CV}}}z^{k})_{j}=x^{k+1}_{i}$, computing it involves evaluating: the entire $s^{k+1}$ for $O(p)$ operations, $(A^{\top}(2s^{k+1}-s^{k}))_{i}$ for $O(p)$ operations, $\mathbf{prox}_{\eta g_{i}}$ for $O(1)$ operations, $\nabla_{i}f({x}^{k})$ for $O(\frac{1}{m}\mathfrak{M}\left[{x}\mapsto{\nabla f(x)}\right])$ operations, as well as updating $Ax^{+}$ for $O(p)$ operations. In this case $\mathfrak{M}\left[{\{z^{k},Ax\}}\mapsto{\{z^{+},Ax^{+}\}}\right]=O(p+\frac{1}{m}\mathfrak{M}\left[{x}\mapsto{\nabla f(x)}\right])$.

Therefore, \mathfrak{M}\left[{\{z^{k},Ax\}}\mapsto{\{z^{+},Ax^{+}\}}\right]=O\big{(}\frac{1}{m+p}\mathfrak{M}\left[{z^{k}}\mapsto{{{\mathcal{T}}_{\textnormal{CV}}}z^{k}}\right]\big{)}. (b) When $m\ll p$ and $\mathfrak{M}\left[{x}\mapsto{\nabla f(x)}\right]=O(m)$, following arguments similar to the above, we have $\mathfrak{M}\left[{\{z^{k},A^{\top}s\}}\mapsto{\{z^{+},A^{\top}s^{+}\}}\right]=O(1)+\mathfrak{M}\left[{x}\mapsto{\nabla_{i}f(x)}\right]$ if $({\mathcal{T}}^{\prime}_{\textnormal{CV}}z^{k})_{j}=x_{i}^{k+1}$; and $\mathfrak{M}\left[{\{z^{k},A^{\top}s\}}\mapsto{\{z^{+},A^{\top}s^{+}\}}\right]=O(m)+\mathfrak{M}\left[{x}\mapsto{\nabla f(x)}\right]$ if $({\mathcal{T}}^{\prime}_{\textnormal{CV}}z^{k})_{j}=s_{i}^{k+1}$. In both cases $\mathfrak{M}\left[{\{z^{k},A^{\top}s\}}\mapsto{\{z^{+},A^{\top}s^{+}\}}\right]=O(\frac{1}{m+p}\mathfrak{M}\left[{z^{k}}\mapsto{{\mathcal{T}}^{\prime}_{\textnormal{CV}}z^{k}}\right])$. ∎

2 Extended Monotropic Programming

We develop a primal-dual coordinate update algorithm for the extended monotropic program:

where $U$ is a symmetric positive semidefinite matrix and $X=\{x:x_{i}\geq 0~{}\forall i\}$. Then, (31) is a special case of (30) with $g_{i}(x_{i})=\iota_{\cdot\geq 0}(x_{i})$, $f(x)=\frac{1}{2}x^{\top}Ux+c^{\top}x$ and $h=\iota_{\{b\}}$.

Applying iteration $\eqref{vucondat}$ to problem $\eqref{emp}$ and eliminating $s^{k+1}$ from the second row yield the Jacobi-style update (denoted as ${\mathcal{T}}_{\textnormal{emp}}$):

To the best of our knowledge, this update is never found in the literature. Note that $x^{k+1}$ no longer depends on $s^{k+1}$, making it more convenient to perform coordinate updates.

In general, when the $s$ update is affine, we can decouple $s^{k+1}$ and $x^{k+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 ${\mathcal{T}}_{\textnormal{emp}}$ to be CF is $\mathbf{prox}_{g}\in{\mathcal{C}}_{1}$ i.e., separable. Indeed, we have ${\mathcal{T}}_{\textnormal{emp}}={\mathcal{T}}_{1}\circ{\mathcal{T}}_{2}$, where

Following Case 5 of Table 2, ${\mathcal{T}}_{\textnormal{emp}}$ is CF. When $m=\Theta(p)$, the separability condition on $\mathbf{prox}_{g}$ can be relaxed to $\mathbf{prox}_{g}\in{\mathcal{F}}_{1}$ since in this case ${\mathcal{T}}_{2}\in{\mathcal{F}}_{2}$, and we can apply Case 7 of Table 2 (by maintaining $\nabla f(x)$, $A^{\top}s$, $Ax$ and $A^{\top}Ax$.)

3 Overlapping-Block Coordinate Updates

In the coordinate update scheme based on (28), if we select $x_{i}$ to update then we must first compute $s^{k+1}$, because the variables $x_{i}$’s and $s_{j}$’s are coupled through the matrix $A$. However, once $x_{i}^{k+1}$ is obtained, $s^{k+1}$ is discarded. It is not used to update $s$ or cached for further use. This subsection introduces ways to utilize the otherwise wasted computation.

The use of relaxation parameters $\rho_{i,j}$ makes our scheme different from that in .

Following the assumptions and arguments in §4.1, if we maintain $Ax$, the cost for each block coordinate update is $O(p)+\mathfrak{M}\left[{x}\mapsto{\nabla_{i}f(x)}\right]$, which is $O(\frac{1}{m}\mathfrak{M}\left[{z}\mapsto{{{\mathcal{T}}_{\textnormal{CV}}}z}\right])$. Therefore the coordinate update scheme (33) is computationally worthy.

The recent paper proposes a different primal-dual coordinate update algorithm. The authors produce a new matrix $\bar{A}$ based on $A$, with only one nonzero entry in each row, i.e. $m_{j}=1$ for each $j$. They also modify $h$ to $\bar{h}$ so that the problem

has the same solution as $\eqref{pdproblem}$. Then they solve $\eqref{fbpdproblem}$ by the scheme $\eqref{pdoverlap}$. Because they have $m_{j}=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 $s_{j}$ multiple times. The computation cost for each block coordinate update of their algorithm is the same as $\eqref{pdoverlap}$, but more memory is needed for the duplicated copies of each $s_{j}$.

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 and state their convergence results. When there is only one agent, all algorithms proposed in this section reduce to stochastic coordinate update algorithms , and their convergence is a direct consequence of Theorem 1. Moreover, our convergence analysis also applies to sync-parallel algorithms.

The two algorithms are based on §4.1 and §4.3, respectively.

Whenever an agent updates a coordinate, the global iteration number $k$ increases by one. The $k$th update is applied to $z_{i_{k}}$, with $i_{k}$ being independent random variables: $z_{i}=x_{i}$ when $i\leq m$ and $z_{i}=s_{i-m}$ when $i>m$. Each coordinate update has the form:

where $\eta_{k}$ is the step size, $z^{k}$ denotes the state of $z$ in global memory just before the update (35) is applied, and $\hat{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 $(\hat{z}^{k}_{i_{k}}-({{\mathcal{T}}_{\textnormal{CV}}}\hat{z}^{k})_{i_{k}})$ is being computed, asynchronous parallel computing allows other agents to make updates to $z$, introducing so-called asynchronous delays. Therefore, $\hat{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 update is introduced to save computation).

If shared memory is used, it is recommended to set all but one $\rho_{i,j}$’s to for each $i$.

$f,g,h^{*}$ are closed proper convex functions, $f$ is differentiable, and $\nabla f$ is Lipschitz continuous with constant $\beta$;

the delay for every coordinate is bounded by a positive number $\tau$, i.e. for every $1\leq i\leq m+p$, $\hat{z}^{k}_{i}=z_{i}^{k-d_{i}^{k}}$ for some $0\leq d_{i}^{k}\leq\tau$;

$\eta_{k}\in[\eta_{\min},\eta_{\max}]$ for certain $\eta_{\min},\eta_{\max}>0$.

Then $(z^{k})_{k\geq 0}$ converges to a $Z^{*}$-valued random variable with probability 1.

The formulas for $\eta_{\min}$ and $\eta_{\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.

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 ${\mathcal{T}}$ and how to efficiently calculate $({\mathcal{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.

We consider the following regularized empirical risk minimization problem

where $a_{j}$’s are sample vectors, $\phi_{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 . 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.2 Support Vector Machine

Given the training data $\{(a_{i},\beta_{i})\}_{i=1}^{m}$ with $\beta_{i}\in\{+1,-1\},\,\forall i$, the kernel support vector machine is

where $\phi$ is a vector-to-vector map, mapping each data $a_{i}$ to a point in a (possibly) higher-dimensional space. If $\phi(a)=a$, then (38) reduces to the linear support vector machine. The model (38) can be interpreted as finding a hyperplane $\{w:x^{\top}w-y=0\}$ to separate two sets of points $\{\phi(a_{i}):\beta_{i}=1\}$ and $\{\phi(a_{i}):\beta_{i}=-1\}$.

where $Q_{ij}=\beta_{i}\beta_{j}k(a_{i},a_{j})$, $k(\cdot,\cdot)$ is a so-called kernel function, and $e=(1,...,1)^{\top}$. If $\phi(a)=a$, then $k(a_{i},a_{j})=a_{i}^{\top}a_{j}$.

Unbiased case

If $y=0$ is enforced in (38), then the solution hyperplane $\{w:x^{\top}w=0\}$ passes through the origin and is called unbiased. Consequently, the dual problem (39) will no longer have the linear constraint $\sum_{i}\beta_{i}s_{i}=0$, leaving it with the coordinate-wise separable box constraints $0\leq s_{i}\leq C$. To solve (39), we can apply the FBS operator ${\mathcal{T}}$ defined by (14). Let $d(s):=\frac{1}{2}s^{\top}Qs-e^{\top}s$, ${\mathcal{A}}=\mathbf{prox}_{[0,C]}$, and ${\mathcal{C}}=\nabla d$. The coordinate update based on FBS is

where we can take $\gamma_{i}=\frac{1}{Q_{ii}}$.

Biased (general) case

The coordinate update based on the full primal-dual splitting scheme (28) is:

where $t,s$ are the primal and dual variables, respectively. Note that we can let $w:=\sum_{i=1}^{m}\beta_{i}s_{i}$ and maintain it. With variable $w$ and substituting (40a) into (40b), we can equivalently write (40) into

where $w^{k}$ is the maintained variable and $s^{k}$ is the intermediate variable.

1.3 Group Lasso

Non-overlapping case [84]

The corresponding coordinate update is the following

where $\nabla_{i}f(x^{k})$ is the partial derivative of $f$ with respect to $x_{i}$ and the step size can be taken to be $\gamma_{i}=\frac{1}{\|A_{:,i}\|^{2}}$. When $\nabla f$ is either cheap or easy-to-maintain, the coordinate update in (45) is inexpensive.

Overlapping case [38]

is cheap. Hence, the corresponding coordinate update of (46) is

where $B_{\lambda}$ is the Euclidean ball of radius $\lambda$. When $\nabla 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.

2 Imaging

Many convex image processing problems have the general form

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

(see Appendix C for the reduction.) The work gives their resolvents

where ${\mathcal{J}}_{\gamma{\mathcal{A}}}$ is often cheap or separable and we can explicitly form ${\mathcal{J}}_{\gamma{\mathcal{B}}}$ as a matrix or implement it based on a fast transform. With the defined ${\mathcal{J}}_{\gamma{\mathcal{A}}}$ and ${\mathcal{J}}_{\gamma{\mathcal{B}}}$, we can apply the RPRS method as $z^{k+1}={\mathcal{T}}_{\text{RPRS}}z^{k}$. The resulting RPRS operator is CF when ${\mathcal{J}}_{\gamma{\mathcal{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.

2.2 Total Variation Image Processing

We consider the following Total Variation (TV) image processing model

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

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,

Let $s,t$ be the dual variables corresponding to $\nabla x$ and $Ax$ respectively, then using (51) and applying (29) give the following full update:

To perform the coordinate updates as described in §4, we can maintain $\nabla^{\top}s^{k}$ and $A^{\top}t^{k}$. Whenever a coordinate of $(s,t)$ is updated, the corresponding $\nabla^{\top}s^{k}$ (or $A^{\top}t^{k})$ should also be updated. Specifically, we have the following coordinate update algorithm

2.3 3D Mesh Denoising

where $f_{i}$’s are differentiable data fidelity terms, $g_{i}$’s are the indicator functions of box constraints, and $\sum_{i}\sum_{j\in{\mathcal{V}}_{i}}h_{i,j}(x_{i}-x_{j})$ is the total variation on the mesh.

We introduce a dual variable $s$ with coordinates $s_{i,j}$, for all ordered pairs of adjacent nodes $(i,j)$, and, based on the overlapping-block coordinate updating scheme $\eqref{pdoverlap}$, perform coordinate update:

3 Finance

Assume that we have one unit of capital and $m$ assets to invest on. The $i$th asset has an expected return rate $\xi_{i}\geq 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

where the objective function is a measure of risk, and the last constraint imposes that the expected return is at least $c$. Let $a_{1}=e/\sqrt{m}$, $b_{1}=1/\sqrt{m}$, $a_{2}=\xi/\|\xi\|_{2}$, and $b_{2}=c/\|\xi\|_{2}$, where $e=(1,\dots,1)^{\top},\xi=(\xi_{1},\dots,\xi_{m})^{\top}$. The above problem is rewritten as

We apply the three-operator splitting scheme (13) to (55). Let $f(x)=\frac{1}{2}x^{\top}Qx$, $D_{1}=\{x:x\geq 0\}$, $D_{2}=\{x:a_{1}^{\top}x\leq b_{1},\,a_{2}^{\top}x\geq b_{2}\}$, $D_{21}=\{x:a_{1}^{\top}x=b_{1}\}$, and $D_{22}=\{x:a_{2}^{\top}x=b_{2}\}$. Based on (13), the full update is

where $y$ is an intermediate variable. As the projection to $D_{1}$ is simple, we discuss how to evaluate the projection to $D_{2}$. Assume that $a_{1}$ and $a_{2}$ are neither perpendicular nor co-linear, i.e., $a_{1}^{\top}a_{2}\neq 0$ and $a_{1}\neq\lambda a_{2}$ for any scalar $\lambda$. In addition, assume $a_{1}^{\top}a_{2}>0$ for simplicity. Let $a_{3}=a_{2}-\frac{1}{a_{1}^{\top}a_{2}}a_{1}$, $b_{3}=b_{2}-\frac{1}{a_{1}^{\top}a_{2}}b_{1}$, $a_{4}=a_{1}-\frac{1}{a_{1}^{\top}a_{2}}a_{2}$, and $b_{4}=b_{1}-\frac{1}{a_{1}^{\top}a_{2}}b_{2}$. Then we can partition the whole space into four areas by the four hyperplanes $a_{i}^{\top}x=b_{i}$, $i=1,\ldots,4$. Let $P_{i}=\{x:a_{i}^{\top}x\leq b_{i},a_{i+1}^{\top}x\geq b_{i+1}\},\,i=1,2,3$ and $P_{4}=\{x:a_{4}^{\top}x\leq b_{4},a_{1}^{\top}x\geq b_{1}\}$. Then

where $q_{i}$ is the $i$th column of $Q$. At each iteration, we select $i\in[m]$, and perform an update to $x_{i}$ according to (57) based on where $x^{k}$ is. We then renew $w_{j}^{k+1}=w_{j}^{k}+a_{ij}(x_{i}^{k+1}-x_{i}^{k}),j=1,2$. Note that checking $x^{k}$ in some $P_{j}$ requires only $O(1)$ operations by using $w_{1}$ and $w_{2}$, so the coordinate update in (57) is inexpensive.

4 Distributed Computing

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:

Applying the FBFS scheme (18) to (59) yields the following full update:

where (60a) and (60c) are applied to all $i\in[m]$. Hence, for each $i$, we group $x_{i}$ and $s_{i}$ together and assign them on agent $i$. We let the master agent maintain $\sum_{j}s_{j}$ and $\sum_{j}x_{j}$. Therefore, in the FBFS coordinate update, updating any $(x_{i},s_{i})$ needs only $y$ and $\sum_{j}s_{j}$ from the master agent, and updating $y$ is done on the master agent. In synchronous parallel setting, at each iteration, each worker agent $i$ computes $s_{i}^{k+1},x_{i}^{k+1}$, then the master agent collects the updates from all of the worker agents and then updates $y$ and $\sum_{j}s_{j}$. The above update can be relaxed to be asynchronous. In this case, the master and worker agents work concurrently, the master agent updates $y$ and $\sum_{j}s_{j}$ as soon as it receives the updated $s_{i}$ and $x_{i}$ from any of the worker agents. It also periodically broadcasts $y$ back to the worker agents.

5 Dimension Reduction

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

In general, we do not know the Lipschitz constant of $\nabla F$, so we have to choose $\eta_{k}$ by line search such that the Armijo condition is satisfied.

Partitioning the variables into $2r$ block coordinates: $(w_{1},\ldots,w_{r},h_{1},\ldots,h_{r})$ where $w_{i}$ and $h_{i}$ are the $i$th columns of $W$ and $H$, respectively, we can apply the coordinate update based on the projected-gradient method:

It is easy to see that $\nabla_{w_{i}}F(W^{k},H^{k})$ and $\nabla_{h_{i}}F(W^{k},H^{k})$ are both Lipschitz continuous with constants $\|h_{i}^{k}\|_{2}^{2}$ and $\|w_{i}^{k}\|_{2}^{2}$ respectively. Hence, we can set

However, it is possible to have $w_{i}^{k}=0$ or $h_{i}^{k}=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 $WH^{\top}=(WD)(HD^{-1})^{\top}$ for any $r\times r$ invertible diagonal matrix $D$. We refer the readers to for more details.

Note that $\nabla_{W}F(W,H)=(WH^{\top}-A)H,\nabla_{H}F(W,H)=(WH^{\top}-A)^{\top}W$ and $\nabla_{w_{i}}F(W,H)=(WH^{\top}-A)h_{i},\nabla_{h_{i}}F(W,H)=(WH^{\top}-A)^{\top}w_{i},\,\forall i.$ Therefore, the coordinate updates given in (63) are computationally worthy (by maintaining the residual $W^{k}(H^{k})^{\top}-A$).

6 Stylized Optimization

then we have u={\mathbf{proj}}_{Q}(v)=\big{(}\xi_{1}^{v}v_{1},\,\xi_{2}^{v}\cdot(v_{2},\ldots,v_{n})\big{)}. Based on this, we have

By the proposition, if ${\mathcal{T}}_{1}$ is an affine operator, then in the composition ${\mathcal{T}}_{1}\circ{\mathbf{proj}}_{Q}$, the computation of ${\mathbf{proj}}_{Q}$ is cheap as long as we maintain $\rho_{1}^{v},\rho_{2}^{v},\xi_{1}^{v},\xi_{2}^{v}$.

where each $Q_{i}$ is a second-order cone, and $\bar{n}\not=n$ in general. The problem (64) is equivalent to

Assume that the matrix $A$ has full row-rank (otherwise, $Ax=b$ has either redundant rows or no solution). Then, in (16), we have ${\mathcal{R}}_{\gamma{\mathcal{B}}}(x)=Bx+d$, where $B:=I-2A^{\top}(AA^{\top})^{-1}A$ and $d:=2A^{\top}(AA^{\top})^{-1}(b+\gamma Ac)-2\gamma c$.

It is trivial to extend this method for SOCPs with a quadratic objective:

because ${\mathcal{J}}_{2}$ is still linear. Clearly, this method applies to linear programs as they are special SOCPs.