ARock: an Algorithmic Framework for Asynchronous Parallel Coordinate Updates

Introduction

Technological advances in data gathering and storage have led to a rapid proliferation of big data in diverse areas such as climate studies, cosmology, medicine, the Internet, and engineering house2014big . The data involved in many of these modern applications are large and grow quickly. Therefore, parallel computational approaches are needed. This paper introduces a new approach to asynchronous parallel computing with convergence guarantees.

In a synchronous(sync) parallel iterative algorithm, the agents must wait for the slowest agent to finish an iteration before they can all proceed to the next one (Figure 1(a)). Hence, the slowest agent may cripple the system. In contract, the agents in an asynchronous(async) parallel iterative algorithm run continuously with little idling (Figure 1(b)). However, the iterations are disordered, and an agent may carry out an iteration without the newest information from other agents.

Asynchrony has other advantages bertsekas1991some : the system is more tolerant to computing faults and communication glitches; it is also easy to incorporate new agents.

On the other hand, it is more difficult to analyze asynchronous algorithms and ensure their convergence. It becomes impossible to find a sequence of iterates that one completely determines the next. Nonetheless, we let any update be a new iteration and propose an async-parallel algorithm (ARock) for the generic fixed-point iteration. It converges if the fixed-point operator is nonexpansive (Def. 1) and has a fixed point.

Let ${\mathcal{H}}_{1},\ldots,{\mathcal{H}}_{m}$ be Hilbert spaces and ${\mathcal{H}}:={\mathcal{H}}_{1}\times\cdots\times{\mathcal{H}}_{m}$ be their Cartesian product. For a nonexpansive operator $T:{\mathcal{H}}\to{\mathcal{H}}$, our problem is to

Finding a fixed point to $T$ is equivalent to finding a zero of $S\equiv I-T,$ denoted by $x^{*}$ such that $0=Sx^{*}$. Hereafter, we will use both $S$ and $T$ for convenience.

Problem (1) is widely applicable in linear and nonlinear equations, statistical regression, machine learning, convex optimization, and optimal control. A generic framework for problem (1) is the Krasnosel’skiĭ–Mann (KM) iteration krasnosel1955two :

where $\alpha\in(0,1)$ is the step size. If $\operatorname*{Fix}T$ — the set of fixed points of $T$ (zeros of $S$) — is nonempty, then the sequence $(x^{k})_{k\geq 0}$ converges weakly to a point in $\operatorname*{Fix}T$ and $(Tx^{k}-x^{k})_{k\geq 0}$ converges strongly to 0. The KM iteration generalizes algorithms in convex optimization, linear algebra, differential equations, and monotone inclusions. Its special cases include the following iterations: alternating projection, gradient descent, projected gradient descent, proximal-point algorithm, Forward-Backward Splitting (FBS) passty1979ergodic , Douglas-Rachford Splitting (DRS) lions1979splitting , a three-operator splitting davis2015three , and the Alternating Direction Method of Multipliers (ADMM) lions1979splitting ; glowinski1975approximation .

In ARock, a set of $p$ agents, $p\geq 1$, solve problem (1) by updating the coordinates $x_{i}\in{\mathcal{H}}_{i}$, $i=1,\ldots,m$, in a random and asynchronous fashion. Algorithm 1 describes the framework. Its special forms for several applications are given in Section 2 below.

Whenever an agent updates a coordinate, the global iteration counter $k$ increases by one. The $k$th update is applied to $x_{i_{k}}\in{\mathcal{H}}_{i_{k}}$, where $i_{k}\in\{1,\ldots,m\}$ is an independent random variable. Each coordinate update has the form:

where $\eta_{k}>0$ is a scalar whose range will be set later, $S_{i_{k}}{x}:=(0,...,0,(Sx)_{i_{k}},0,...,0),$ and $mp_{i_{k}}$ is used to normalize nonuniform selection probabilities. In the uniform case, namely, $p_{i}\equiv\frac{1}{m}$ for all $i$, we have $mp_{i_{k}}\equiv 1$, which simplifies the update (3) to

Here, the point $\hat{x}^{k}$ is what an agent reads from global memory to its local cache and to which $S_{i_{k}}$ is applied, and $x^{k}$ denotes the state of $x$ in global memory just before the update (3) is applied. In a sync-parallel algorithm, we have $\hat{x}^{k}=x^{k}$, but in ARock, due to possible updates to $x$ by other agents, $\hat{x}^{k}$ can be different from $x^{k}$. This is a key difference between sync-parallel and async-parallel algorithms. In Subsection 1.2 below, we will establish the relationship between $\hat{x}^{k}$ and $x^{k}$ as

The update (3) is only computationally worthy if $S_{i}x$ is much cheaper to compute than $Sx$. Otherwise, it is more preferable to apply the full KM update (2). In Section 2, we will present several applications that have the favorable structures for ARock. The recent work peng2016coordinate studies coordinate friendly structures more thoroughly.

The convergence of ARock (Algorithm 1) is stated in Theorems 3.2 and 3.3. Here we include a shortened version, leaving detailed bounds to the full theorems:

Let $T:{\mathcal{H}}\to{\mathcal{H}}$ be a nonexpansive operator that has a fixed point. Let $(x^{k})_{k\geq 0}$ be the sequence generated by Algorithm 1 with properly bounded step sizes $\eta_{k}$. Then, with probability one, $(x^{k})_{k\geq 0}$ converges weakly to a fixed point of $T$. This convergence becomes strong if ${\mathcal{H}}$ has a finite dimension.

In addition, if $T$ is demicompact (see Definition 2 below), then with probability one, $(x^{k})_{k\geq 0}$ converges strongly to a fixed point of $T$ .

In the theorem, the weak convergence result only requires $T$ to be nonexpansive and has a fixed point. In addition, the computation requires: (a) bounded step sizes; (b) random coordinate selection; and (c) a finite maximal delay $\tau$. Assumption (a) is standard, and we will see the bound can be $O(1)$. Assumption (b) is essential to both the analysis and the numerical performance of our algorithms. Assumption (c) is not essential; an infinite delay with a light tail is allowed (but we leave it to future work). The strong convergence result applies to all the examples in Section 2, and the linear convergence result applies to Examples 2.2 and 2.4 when the corresponding operator $S$ is quasi-strongly monotone. Step sizes $\eta_{k}$ are discussed in Remarks 2 and 4.

ARock employs random coordinate selection. This subsection discusses its advantages and disadvantages.

Its main disadvantage is that an agent cannot caching the data associated with a coordinate. The variable $x$ and its related data must be either stored in global memory or passed through communication. A secondary disadvantage is that pseudo-random number generation takes time, which becomes relatively significant if each coordinate update is cheap. (The network optimization examples in Subsections 2.3 and 2.6.2 are exceptions, where data are naturally stored in a distributed fashion and random coordinate assignments are the results of Poisson processes.)

There are several advantages of random coordinate selection. It realizes the user-specified update frequency $p_{i}$ for every component $x_{i}$, $i=1,\ldots,m$, even when different agents have different computing powers and different coordinate updates cost different amounts of computation. Therefore, random assignment ensures load balance. The algorithm is also fault tolerant in the sense that if one or more agents fail, it will still converge to a fixed-point of $T$. In addition, it has been observed numerically on certain problems chang2008coordinate that random coordinate selection accelerates convergence.

2 Uncoordinated memory access

In ARock, since multiple agents simultaneously read and update $x$ in global memory, $\hat{x}^{k}$ — the result of $x$ that is read from global memory by an agent to its local cache for computation — may not equal $x^{j}$ for any $j\leq k$, that is, $\hat{x}^{k}$ may never be consistent with a state of $x$ in global memory. This is known as inconsistent read. In contrast, consistent read means that $\hat{x}^{k}=x^{j}$ for some $j\leq k$, i.e., $\hat{x}^{k}$ is consistent with a state of $x$ that existed in global memory.

Even with inconsistent read, each component is consistent under the atomic coordinate update assumption, which will be defined below. Therefore, we can express what has been read in terms of the changes of individual coordinates. In the above example, the first change is $x^{1}_{1}-x^{0}_{1}=1$, which is added to $x_{1}$ just before time $t_{1}$ by agent 2, and the second change is $x^{2}_{4}-x^{1}_{4}=2$, added to $x_{4}$ just before time $t_{2}$ by agent 3. The inconsistent read by agent 1, which gives the result $^{T}$, equals $x^{0}+0\times(x^{1}-x^{0})+1\times(x^{2}-x^{1})$.

We have demonstrated that $\hat{x}^{k}$ can be inconsistent, but each of its coordinates is consistent, that is, for each $i$, $\hat{x}^{k}_{i}$ is an ever-existed state of $x_{i}$ among $x_{i}^{k},\ldots,x_{i}^{k-\tau}$. Suppose that $\hat{x}^{k}_{i}=x_{i}^{\underline{d}}$, where $\underline{d}\in\{k,k-1,\ldots,k-\tau\}$. Therefore, $\hat{x}^{k}_{i}$ can be related to $x^{k}_{i}$ through the interim changes applied to $x_{i}$. Let $J_{i}(k)\subset\{k-1,\ldots,k-\tau\}$ be the index set of these interim changes. If $J_{i}(k)\not=\emptyset$, then $\underline{d}=\min\{d\in J_{i}(k)\}$; otherwise, $\underline{d}=k$. In addition, we have $\hat{x}^{k}_{i}=x_{i}^{\underline{d}}=x^{k}_{i}+\sum_{d\in J_{i}(k)}(x_{i}^{d}-x_{i}^{d+1})$. Since the global counter $k$ is increased after each coordinate update, updates to $x_{i}$ and $x_{j}$, $i\not=j$, must occur at different $k$’s and thus $J_{i}(k)\cap J_{j}(k)=\emptyset,\,\forall i\neq j$. Therefore, by letting $J(k):=\cup_{i}J_{i}(k)\subset\{k-1,\ldots,k-\tau\}$ and noticing $(x_{i}^{d}-x_{i}^{d+1})=0$ for $d\in J_{j}(k)$ where $i\not=j$, we have $\textstyle\hat{x}^{k}_{i}=x^{k}_{i}+\sum_{d\in J(k)}(x_{i}^{d}-x_{i}^{d+1}),\forall i=1,\ldots,m$, which is equivalent to (5). Here, we have made two assumptions:

atomic coordinate update: a coordinate is not further broken to smaller components during an update; they are all updated at once.

bounded maximal delay $\tau$: during any update cycle of an agent, $x$ in global memory is updated at most $\tau$ times by other agents.

When each coordinate is a single scalar, updating the scalar is a single atomic instruction on most modern hardware, so the first assumption naturally holds, and our algorithm is lock-free. The case where a coordinate is a block that includes multiple scalars is discussed in the next subsection.

In the “block coordinate” case (updating a block of several coordinates each time), the atomic coordinate update assumption can be met by either employing a per-coordinate memory lock or taking the following dual-memory approach: Store two copies of each coordinate $x_{i}\in{\mathcal{H}}_{i}$ in global memory, denoting them as $x_{i}^{(0)}$ and $x_{i}^{(1)}$; let a bit $\alpha_{i}\in\{0,1\}$ point to the active copy; an agent will only read $x_{i}$ from the active copy $x_{i}^{(\alpha_{i})}$; before an agent updates the components of $x_{i}$, it obtains a memory lock to the inactive copy $x_{i}^{(1-\alpha_{i})}$ to prevent other agents from simultaneously updating it; then after it finishes updating $x_{i}^{(1-\alpha_{i})}$, flip the bit $\alpha_{i}$ so that other agents will begin reading from the updated copy. This approach never blocks any read of $x_{i}$, yet it eliminates inconsistency.

3 Straightforward generalization

Our async-parallel coordinate update scheme (3) can be generalized to (overlapping) block coordinate updates after a change to the step size. Specifically, the scheme (3) can be generalized to

where $U_{i_{k}}$ is randomly drawn from a set of operators $\{U_{1},\ldots,U_{n}\}$ ($n\leq m$), $U_{i}:{\mathcal{H}}\to{\mathcal{H}}$, following the probability $P(i_{k}=i)=p_{i}$, $i=1,\ldots,n$ ($p_{i}>0$, and $\sum_{i=1}^{n}p_{i}=1$). The operators must satisfy $\sum_{i=1}^{n}U_{i}=I_{\mathcal{H}}$ and $\sum_{i=1}^{n}\|U_{i}x\|^{2}\leq C\|x\|^{2}$ for some $C>0$.

Let $U_{i}:x\mapsto(0,\ldots,0,x_{i},0,\ldots,0),~{}i=1,\ldots,m$, which has $C=1$; then (6) reduces to (3). If ${\mathcal{H}}$ is endowed with a metric $M$ such that $\rho_{1}\|x\|^{2}\leq\|x\|_{M}^{2}\leq\rho_{2}\|x\|^{2}$ (e.g., the metric in the Condat-Vũ primal-dual splitting condat2013primal ; vu2013splitting ), then we have

In general, multiple coordinates can be updated in (6). Consider linear $U_{i}:x\mapsto(a_{i1}x_{1},\cdots,a_{im}x_{m})$, $i=1,\ldots,m$, where $\sum_{i=1}^{n}a_{ij}=1$ for each $j$. Then, for $C:=\max\left\{\sum_{i=1}^{n}a_{i1}^{2},\cdots,\sum_{i=1}^{n}a_{im}^{2}\right\}$, we have

4 Special cases

If there is only one agent ($p=1$), ARock (Algorithm 1) reduces to randomized coordinate update, which includes the special case of randomized coordinate descent nesterov2012rcd for convex optimization. Sync-parallel coordinate update is another special case of ARock corresponding to $\hat{x}^{k}\equiv x^{k}$. In both cases, there is no delay, i.e., $\tau=0$ and $J(k)=\emptyset$. In addition, the step size $\eta_{k}$ can be more relaxed. In particular, if $p_{i}=\frac{1}{m}$, $\forall i$, then we can let $\eta_{k}=\eta$, $\forall k$, for any $\eta<1$, or $\eta<1/\alpha$ when $T$ is $\alpha$-averaged (see Definition 2 for the definition of an $\alpha$-averaged operator).

5 Related work

Chazan and Miranker chazan1969chaotic proposed the first async-parallel method in 1969. The method was designed for solving linear systems. Later, async-parallel methods have been successful applied in many fields, e.g., linear systems avron2014revisiting ; bethune2014performance ; FSS1997asyn-addSch ; rosenfeld1969case , nonlinear problems BMR1997asyn-multisplit ; baudet1978asynchronous , differential equations aharoni2000parallel ; AAI1998implicit ; Chau20081126 ; donzis2014asynchronous , consensus problems LMS1986asynchronous ; leifang_information_2005 , and optimization hsieh2015passcode ; liu2014asynchronous ; liu2013asynchronous ; tai2002convergence ; zhang2014asynchronous . We review the theory for async-parallel fixed-point iteration and its applications.

The above works assign coordinates in a deterministic manner. Different from them, ARock is stochastic, works for nonexpansive operators, and is more applicable.

Linear, nonlinear, and differential equations. The first async-parallel method for solving linear equations was introduced by Chazan and Miranker in chazan1969chaotic . They proved that on solving linear systems, P-contraction was necessary and sufficient for convergence. The performance of the algorithm was studied by Iain et al. bethune2014performance ; rosenfeld1969case on different High Performance Computing (HPC) architectures. Recently, Avron et al. avron2014revisiting revisited the async-parallel coordinate update and showed its linear convergence for solving positive-definite linear systems. Tarazi and Nabih el1982some extended the poineering work chazan1969chaotic to solving nonlinear equations, and the async-parallel methods have also been applied for solving differential equations, e.g., in aharoni2000parallel ; AAI1998implicit ; Chau20081126 ; donzis2014asynchronous . Except for avron2014revisiting , all these methods are totally async-parallel with the P-contraction condition or its variants. On solving a positive-definite linear system, avron2014revisiting made assumptions similar to ours, and it obtained better linear convergence rate on that special problem.

Our framework differs from the recent surge of the aforementioned sync-parallel and async-parallel coordinate descent algorithms (e.g., peng2013parallel ; kyrola2011parallel ; liu2013asynchronous ; liu2014asynchronous ; hsieh2015passcode ; richtarik2015parallel ). While they apply to convex function minimization, ARock covers more cases (such as ADMM, primal-dual, and decentralized methods) and also provides sequence convergence. In Section 2, we will show that some of the existing async-parallel coordinate descent algorithms are special cases of ARock, through relating their optimality conditions to nonexpansive operators. Another difference is that the convergence of ARock only requires a nonexpansive operator with a fixed point, whereas properties such as strong convexity, bounded feasible set, and bounded sequence, which are seen in some of the recent literature for async-parallel convex minimization, are unnecessary.

Others. Besides solving equations and optimization problems, there are also applications of async-parallel algorithms to optimal control problems LMS1986asynchronous , network flow problems ESMG1996asyn-flex , and consensus problems of multi-agent systems leifang_information_2005 .

6 Contributions

Our contributions and techniques are summarized below:

ARock is the first async-parallel coordinate update framework for finding a fixed point to a nonexpansive operator.

By introducing a new metric and establishing stochastic Fejér monotonicity, we show that, with probability one, ARock converges to a point in the solution set; linear convergence is obtained for quasi-strongly monotone operators.

Based on ARock, we introduce an async-parallel algorithm for linear systems, async-parallel ADMM algorithms for distributed or decentralized computing problems, as well as async-parallel operator-splitting algorithms for nonsmooth minimization problems. Some problems are treated in they async-parallel fashion for the first time in history. The developed algorithms are not straightforward modifications to their serial versions because their underlying nonexpansive operators must be identified before applying ARock.

7 Notation, definitions, background of monotone operators

Throughout this paper, ${\mathcal{H}}$ denotes a separable Hilbert space equipped with the inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$, and $(\Omega,{\mathcal{F}},P)$ denotes the underlying probability space, where $\Omega$, ${\mathcal{F}}$, and $P$ are the sample space, $\sigma$-algebra, and probability measure, respectively. The map $x:(\Omega,{\mathcal{F}})\rightarrow({\mathcal{H}},{\mathcal{B}})$, where ${\mathcal{B}}$ is the Borel $\sigma$-algebra, is an ${\mathcal{H}}$-valued random variable. Let $(x^{k})_{k\geq 0}$ denote either a sequence of deterministic points in ${\mathcal{H}}$ or a sequence of ${\mathcal{H}}$-valued random variables, which will be clear from the context, and let $x_{i}\in{\mathcal{H}}_{i}$ denote the $i$th coordinate of $x$. In addition, we let ${\mathcal{X}}^{k}:=\sigma(x^{0},\hat{x}^{1},x^{1},...,\hat{x}^{k},x^{k})$ denote the smallest $\sigma$-algebra generated by $x^{0},\hat{x}^{1},x^{1},...,\hat{x}^{k},x^{k}$. “Almost surely” is abbreviated as “a.s.”, and the $n$ product space of ${\mathcal{H}}$ is denoted by ${\mathcal{H}}^{n}$. We use $\to$ and $\rightharpoonup$ for strong convergence and weak convergence, respectively.

We define $\operatorname*{Fix}T:=\{x\in{\mathcal{H}}~{}|~{}Tx=x\}$ as the set of fixed points of operator $T$, and, in the product space, we let ${\mathbf{X}}^{*}:=\{(x^{*},x^{*},...,x^{*})~{}|~{}x^{*}\in\operatorname*{Fix}T\}\subseteq{\mathcal{H}}^{\tau+1}$.

An operator $T:{\mathcal{H}}\rightarrow{\mathcal{H}}$ is $c$-Lipschitz, where $c\geq 0$, if it satisfies $\|Tx-Ty\|\leq c\|x-y\|$, $\forall x,y\in{\mathcal{H}}$. In particular, $T$ is nonexpansive if $c\leq 1$, and contractive if $c<1$.

Consider an operator $T:{\mathcal{H}}\rightarrow{\mathcal{H}}$.

$T$ is $\alpha$-averaged with $\alpha\in(0,1)$, if there is a nonexpansive operator $R:{\mathcal{H}}\rightarrow{\mathcal{H}}$ such that $T=(1-\alpha)I_{\mathcal{H}}+\alpha R$, where $I_{\mathcal{H}}:{\mathcal{H}}\rightarrow{\mathcal{H}}$ is the identity operator.

$T$ is $\beta$-cocoercive with $\beta>0$, if $\langle x-y,Tx-Ty\rangle\geq\beta\|Tx-Ty\|^{2},~{}\forall x,y\in{\mathcal{H}}.$

$T$ is $\mu$-strongly monotone, where $\mu>0$, if it satisfies $\langle x-y,Tx-Ty\rangle\geq\mu\|x-y\|^{2},~{}\forall x,y\in{\mathcal{H}}.$ When the inequality holds for $\mu=0$, $T$ is monotone.

$T$ is quasi-$\mu$-strongly monotone, where $\mu>0$, if it satisfies $\langle x-y,Tx\rangle\geq\mu\|x-y\|^{2},~{}\forall x\in{\mathcal{H}},y\in\operatorname*{zer}T:=\{y\in{\mathcal{H}}\mid Ty=0\}$. When the inequality holds for $\mu=0$, $T$ is quasi-monotone.

$T$ is demicompact petryshyn1966construction at $x\in{\mathcal{H}}$ if for every bounded sequence $(x^{k})_{k\geq 0}$ in ${\mathcal{H}}$ such that $Tx^{k}-x^{k}\to x$, there exists a strongly convergent subsequence.

Averaged operators are nonexpansive. By the Cauchy-Schwarz inequality, a $\beta$-cocoercive operator is $\frac{1}{\beta}$-Lipschitz; the converse is generally untrue, but true for the gradients of convex differentiable functions. Examples are given in the next section.

Applications

In this section, we provide some applications that are special cases of the fixed-point problem (1). For each application, we identify its nonexpansive operator $T$ (or the corresponding operator $S$) and implement the conditions in Theorem 1.1. For simplicity, we use the uniform distribution, $p_{1}=\cdots=p_{m}=1/m$, and apply the simpler update (4) instead of (3).

(bauschke2011convex, , Example 22.5) Suppose that $T$ is $c$-Lipschitz continuous with $c\in[0,1)$. Then, $I-T$ is $(1-c)$-strongly monotone.

Suppose $\|M\|_{2}<1$. Since $T$ is $\|M\|_{2}$-Lipschitz continuous, by Proposition 1, $S$ is $(1-\|M\|_{2})$-strongly monotone. By Theorem 3.3, Algorithm 2 converges linearly.

2 Minimize convex smooth function

where $f$ is a closed proper convex differentiable function and $\nabla f$ is $L$-Lipschitz continuous, $L>0$. Let $S:=\frac{2}{L}\,\nabla f$. As $f$ is convex and differentiable, $x$ is a minimizer of $f$ if and only if $x$ is a zero of $S$. Note that $S$ is $\frac{1}{2}$-cocoercive. By Lemma 1, $T\equiv I-S$ is nonexpansive. Applying ARock, we have the following iteration:

where $S_{i_{k}}x=\frac{2}{L}(0,...,0,\nabla_{i_{k}}f(x),0,...,0)^{T}$. Note that $\nabla f$ needs a structure that makes it cheap to compute $\nabla_{i_{k}}f(\hat{x}^{k})$. Let us give two such examples: (i) quadratic programming: $f(x)=\frac{1}{2}x^{T}Ax-b^{T}x$, where $\nabla f(x)=Ax-b$ and $\nabla_{i_{k}}f(\hat{x}^{k})$ only depends on a part of $A$ and $b$; (ii) sum of sparsely supported functions: $f=\sum_{j=1}^{N}f_{j}$ and $\nabla f=\sum_{j=1}^{N}\nabla f_{j}$, where each $f_{j}$ depends on just a few variables.

Theorem 3.2 below guarantees the convergence of $(x^{k})_{k\geq 0}$ if $\eta_{k}\in[\eta_{\min},\frac{1}{2\tau/\sqrt{m}+1})$. In addition, If $f(x)$ is restricted strongly convex, namely, for any $x\in{\mathcal{H}}$ and $x^{*}\in X^{*}$, where $X^{*}$ is the solution set to (7), we have $\langle x-x^{*},\nabla f(x)\rangle\geq\mu\|x-x^{*}\|^{2}$ for some $\mu>0$, then $S$ is quasi-strongly monotone with modulus $\mu$. According to Theorem 3.3, iteration (8) converges at a linear rate if the step size meets the condition therein.

3 Decentralized consensus optimization

If the agents are $p$ independent Poisson processes and that each agent $i$ has activation rate $\lambda_{i}$, then the probability that agent $i$ activates before other agents is equal to $\frac{\lambda_{i}}{\sum_{i=1}^{p}\lambda_{i}}$ larson1981urban and therefore our random sample scheme holds and ARock applies naturally. The algorithm is summarized as follows:

4 Minimize smooth ++ nonsmooth functions

+ nonsmooth functions Consider the problem

where $f$ is closed proper convex and $g$ is convex and $L$-Lipschitz differentiable with $L>0$. Problems in the form of (10) arise in statistical regression, machine learning, and signal processing and include well-known problems such as the support vector machine, regularized least-squares, and regularized logistic regression. For any $x\in{\mathcal{H}}$ and scalar $\gamma\in(0,\frac{2}{L})$, define the proximal operator $\mathbf{prox}_{f}:{\mathcal{H}}\rightarrow{\mathcal{H}}$ and the reflective-proximal operator $\mathbf{refl}_{f}:{\mathcal{H}}\rightarrow{\mathcal{H}}$ as

Assume that $f$ is a closed proper convex function, and $g$ is $L$-Lipschitz differentiable and strongly convex with modulus $\mu>0$. Let $\gamma\in(0,\frac{2}{L})$. Then, both $I-\gamma\nabla g$ and $\mathbf{prox}_{\gamma f}\circ(I-\gamma\nabla g)$ are quasi-contractive operators.

We first show that $I-\gamma\nabla g$ is a quasi-contractive operator. Note

where the first inequality follows from the Baillon-Haddad theoremLet $g$ be a convex differentiable function. Then, $\nabla g$ is $L$-Lipschitz if and only if it is $\frac{1}{L}$-cocoercive. and the second one from the strong convexity of $g$. Hence, $I-\gamma\nabla g$ is quasi-contractive if $0<\gamma<2/L$. Since $f$ is convex, $\mathbf{prox}_{\gamma f}$ is firmly nonexpansive, and thus we immediately have the quasi-contractiveness of $\mathbf{prox}_{\gamma f}\circ(I-\gamma\nabla g)$ from that of $I-\gamma\nabla g$.

5 Minimize nonsmooth ++ nonsmooth functions

where both $f(x)$ and $g(x)$ are closed proper convex and their $\mathbf{prox}$ maps are easy to compute. Define the Peaceman-Rachford lions1979splitting operator:

where $\hat{x}^{k}$ and $\hat{y}^{k}$ are intermediate variables. Note that the order in which the proximal operators are applied to $f$ and $g$ affects both $z^{k}$ YanYin2014 and whether coordinate-wise updates can be efficiently computed. Next, we present two special cases of (13) in Subsections 2.5.1 and 2.6 and discuss how to efficiently implement the update (15).

Suppose that $C_{1},...,C_{m}$ are closed convex subsets of ${\mathcal{H}}$ with a nonempty intersection. The problem is to find a point in the intersection. Let ${\mathcal{I}}_{C_{i}}$ be the indicator function of the set $C_{i}$, that is, ${\mathcal{I}}_{C_{i}}(x)=0$ if $x\in C_{i}$ and $\infty$ otherwise. The feasibility problem can be formulated as the following

Let $z^{k}=(z^{k}_{1},\ldots,z^{k}_{m})\in{\mathcal{H}}^{m}$, $\hat{z}^{k}=(\hat{z}^{k}_{1},\ldots,\hat{z}^{k}_{m})\in{\mathcal{H}}^{m}$, and $\hat{\bar{z}}^{k}\in{\mathcal{H}}$. We can implement (15) as follows (see Appendix A for the step-by-step derivation):

The update (16) can be implemented as follows. Let global memory hold $z_{1},\ldots,z_{m}$, as well as $\bar{z}=\frac{1}{m}\sum_{i=1}^{m}z_{i}$. At the $k$th update, an agent independently generates a random number $i_{k}\in\{1,\ldots,m\}$, then reads $z_{i_{k}}$ as $\hat{z}_{i_{k}}^{k}$ and $\bar{z}$ as ${\hat{\bar{z}}}^{k}$, and finally computes $\hat{y}_{i_{k}}$ and updates $z_{i_{k}}$ in global memory according to (16). Since $\bar{z}$ is maintained in global memory, the agent updates $\bar{z}$ according to $\bar{z}^{k+1}=\bar{z}^{k}+\frac{1}{m}(z_{i_{k}}^{k+1}-z_{i_{k}}^{k})$. This implementation saves each agent from computing (16a) or reading all $z_{1},\ldots,z_{m}$. Each agent only reads $z_{i_{k}}$ and $\bar{z}$, executes (16b), and updates $z_{i_{k}}$ (16c) and $\bar{z}$.

6 Async-parallel ADMM

This is another application of (15). Consider

where ${\mathcal{H}}_{1}$ and ${\mathcal{H}}_{2}$ are Hilbert spaces, $A$ and $B$ are bounded linear operators. We apply the update (15) to the Lagrange dual of (17) (see gabay1983chapter for the derivation):

where $d_{f}(w):=f^{*}(A^{*}w)$, $d_{g}(w):=g^{*}(B^{*}w)-\langle w,b\rangle$, and $f^{*}$ and $g^{*}$ denote the convex conjugates of $f$ and $g$, respectively. The proximal maps induced by $d_{f}$ and $d_{g}$ can be computed via solving subproblems that involve only the original terms in (17): $z^{+}=\mathbf{prox}_{\gamma d_{f}}(z)$ can be computed by (see Appendix A for the derivation)

and $z^{+}=\mathbf{prox}_{\gamma d_{g}}(z)$ by

Plugging (19) and (20) into (15) yields the following naive implementation

Note that $2\eta_{k}$ in (15c) becomes $\eta_{k}$ in (21e) because ADMM is equivalent to the Douglas-Rachford operator, which is the average of the Peaceman-Rachford operator and the identity operator lions1979splitting . Under favorable structures, (21) can be implemented efficiently. For instance, when $A$ and $B$ are block diagonal matrices and $f,g$ are corresponding block separable functions, steps (21a)–(21d) reduce to independent computation for each $i$. Since only $\hat{w}_{f,i_{k}}^{k}$ and $\hat{w}_{g,i_{k}}^{k}$ are needed to update the main variable $z^{k}$, we only need to compute (21a)–(21d) for the $i_{k}$th block. This is exploited in distributed and decentralized ADMM in the next two subsections.

Consider the consensus optimization problem:

where $f_{i}(x_{i})$ are proper close convex functions. Rewrite (22) to the ADMM form:

where $g=0$. Now apply the async-parallel ADMM (21) to (2.6.1) with dual variables $z_{1},...,z_{m}\in{\mathcal{H}}$. In particular, the update (21a), (21b), (21c), (21d) reduce to

Therefore, we obtain the following async-parallel ADMM algorithm for the problem (22). This algorithm applies to all the distributed applications in boyd2011distributed .

6.2 Async-parallel ADMM for decentralized optimization

Let $V=\{1,...,m\}$ be a set of agents and E=\{(i,j)~{}|~{}\text{if agenti$connects to agent$j},i<j\} be the set of undirected links between the agents. Consider the following decentralized consensus optimization problem on the graph $G=(V,E)$:

for proper matrices $A$ and $B$. Applying the async-parallel ADMM (21) to (29) gives rise to the following simplified update: Let $E(i)$ be the set of edges connected with agent $i$ and $|E(i)|$ be its cardinality. Let $L(i)=\{j~{}|~{}(j,i)\in E(i),j<i\}$ and $R(i)=\{j~{}|~{}(i,j)\in E(i),j>i\}$. To every pair of constraints $x_{i}=y_{ij}$ and $x_{j}=y_{ij}$, $(i,j)\in E$, we associate the dual variables $z_{ij,i}$ and $z_{ij,j}$, respectively. Whenever some agent $i$ is activated, it calculates

We present the algorithm based on (30) for problem (27) in Algorithm 5.

Algorithm 5 activates one agent at each iteration and updates all the dual variables associated with the agent. In this case, only one-sided communication is needed, for sending the updated dual variables in the last step. We allow this communication to be delayed in the sense that agent $i$’s neighbors may be activated and start their computation before receiving the latest dual variables from agent $i$.

Our algorithm is different from the asynchronous ADMM algorithm by Wei and Ozdaglar wei2013on . Their algorithm activates an edge and its two associated agents at each iteration and thus requires two-sided communication at each activation. We can recover their algorithm as a special case by activating an edge $(i,j)\in E$ and its associated agents $i$ and $j$ at each iteration, updating the dual variables $z_{ij,i}$ and $z_{ij,j}$ associated with the edge, as well as computing the intermediate variables $x_{i}$, $x_{j}$, and $y_{ij}$. The updates are derived from (29) with the orders of $x$ and $y$ swapped. Note that wei2013on does not consider the situation that adjacent edges are activated in a short period of time, which may cause overlapped computation and delay communication. Indeed, their algorithm corresponds to $\tau=0$ and the corresponding stepsize $\eta_{k}\equiv 1$. Appendix B presents the steps to derive the algorithms in this subsection.

Convergence

We establish weak and strong convergence in Subsection 3.1 and linear convergence in Subsection 3.2. Step size selection is also discussed.

Throughout the our analysis, we assume $p_{\min}:=\min_{i}p_{i}>0$ and

We let $|J(k)|$ be the number of elements in $J(k)$ (see Subsection 1.2). Only for the purpose of analysis, we define the (never computed) full update at $k$th iteration:

Lemma 1 below shows that $T$ is nonexpansive if and only if $S$ is ${1/2}$-cocoercive.

Operator $T:{\mathcal{H}}\to{\mathcal{H}}$ is nonexpansive if and only if $S=I-T$ is ${1}/{2}$-cocoercive, i.e., $\langle x-y,Sx-Sy\rangle\geq\frac{1}{2}\|Sx-Sy\|^{2},{\forall~{}x,y\in{\mathcal{H}}}$.

See textbook (bauschke2011convex, , Proposition 4.33) for the proof of the “if” part, and the “only if” part, though missing there, follows by just reversing the proof.

The lemma below develops an an upper bound for the expected distance between $x^{k+1}$ and any $x^{*}\in\operatorname*{Fix}T$.

Let $(x^{k})_{k\geq 0}$ be the sequence generated by Algorithm 1. Then for any $x^{*}\in\operatorname*{Fix}T$ and $\gamma>0$ (to be optimized later), we have

where the first inequality follows from the Young’s inequality. Plugging (35) and (36) into (34) gives the desired result.

We need the following lemma on nonnegative almost supermartingales robbins1985convergence .

Let ${\mathcal{H}}^{\tau+1}=\prod_{i=0}^{\tau}{\mathcal{H}}$ be a product space and $\langle\cdot\,|\,\cdot\rangle$ be the induced inner product:

Let $M^{\prime}$ be a symmetric $(\tau+1)\times(\tau+1)$ tri-diagonal matrix with its main diagonal as $\sqrt{p_{\min}}[\frac{1}{\sqrt{p_{\min}}}+\tau,2\tau-1,2\tau-3,\ldots,1]$ and first off-diagonal as $-\sqrt{p_{\min}}[\tau,\tau-1,\ldots,1]$, and let $M=M^{\prime}\otimes I_{\mathcal{H}}$. Here $\otimes$ represents the Kronecker product. For a given $(y^{0},\cdots,y^{\tau})\in{\mathcal{H}}^{\tau+1}$, $(z^{0},\cdots,z^{\tau})=M(y^{0},\cdots,y^{\tau})$ is given by:

Then $M$ is a self-adjoint and positive definite linear operator since $M^{\prime}$ is symmetric and positive definite, and we define $\langle\cdot\,|\,\cdot\rangle_{M}=\langle\cdot\,|\,M\cdot\rangle$ as the $M$-weighted inner product and $\|\cdot\|_{M}$ the induced norm. Let

where we set $x^{k}=x^{0}$ for $k<0$. With

we have the following fundamental inequality:

Let $(x^{k})_{k\geq 0}$ be the sequence generated by ARock. Then for any ${\mathbf{x}}^{*}\in{\mathbf{X}}^{*}$, it holds that

Let $\gamma=m\sqrt{p_{\min}}$. Since $J(k)\subset\{k-1,\cdots,k-\tau\}$, then (33) indicates

Let us check our step size bound $\frac{mp_{\min}}{2\tau\sqrt{p_{\min}}+1}$. Consider the uniform case: $p_{\min}\equiv p_{i}\equiv\frac{1}{m}$. Then, the bound simplifies to $\frac{1}{1+2\tau/\sqrt{m}}$. If the max delay is no more than the square root of coordinates, i.e., $\tau=O(\sqrt{m})$, then the bound is $O(1)$. In general, $\tau$ depends on several factors such as problem structure, system architecture, load balance, etc. If all updates and agents are identical, then $\tau$ is proportional to $p$, the number of agents. Hence, ARock takes an $O(1)$ step size for solving a problem with $m$ coordinates by $p=\sqrt{m}$ agents under balanced loads.

The next lemma is a direct consequence of the invertibility of the metric $M$.

A sequence $({\mathbf{z}}^{k})_{k\geq 0}\subset{\mathcal{H}}^{\tau+1}$ (weakly) converges to ${\mathbf{z}}\in{\mathcal{H}}^{\tau+1}$ under the metric $\langle\cdot\,|\,\cdot\rangle$ if and only if it does so under the metric $\langle\cdot\,|\,\cdot\rangle_{M}$.

In light of Lemma 4, the metric of the inner product for weak convergence in the next lemma is not specified. The lemma and its proof are adapted from combettes2014stochastic .

Let $(x^{k})_{k\geq 0}\subset{\mathcal{H}}$ be the sequence generated by ARock with $\eta_{k}\in[\eta_{\min},\frac{cmp_{\min}}{2\tau\sqrt{p_{\min}}+1}]$ for any $\eta_{\min}>0$ and $0<c<1$. Then we have:

$\sum_{k=0}^{\infty}\|x^{k}-\bar{x}^{k+1}\|^{2}<\infty$ a.s..

$x^{k}-x^{k+1}\rightarrow 0$ a.s. and $\hat{x}^{k}-x^{k+1}\rightarrow 0$ a.s..

The sequence $({\mathbf{x}}^{k})_{k\geq 0}\subset{\mathcal{H}}^{\tau+1}$ is bounded a.s..

Let $\mathscr{Z}({\mathbf{x}}^{k})$ be the set of weakly convergent cluster points of $({\mathbf{x}}^{k})_{k\geq 0}$. Then, $\mathscr{Z}({\mathbf{x}}^{k})\subseteq{\mathbf{X}}^{*}$ a.s..

(i): Note that $\inf_{k}\left(\frac{1}{\eta_{k}}-\frac{2\tau}{m\sqrt{p_{\min}}}-\frac{1}{mp_{\min}}\right)>0$. Also note that, in (42), $\|\bar{x}^{k+1}-x^{k}\|^{2}=\|\eta_{k}S\hat{x}^{k}\|^{2}$ is ${\mathcal{X}}^{k}$-measurable. Hence, applying Lemma 3 with $\xi_{k}=\eta_{k}=0$ and $\alpha_{k}=\xi_{k}({\mathbf{x}}^{*}),\,\forall k,$ to (42) gives this result directly.

(ii) From (i), we have $x^{k}-\bar{x}^{k+1}\rightarrow 0$ a.s.. Since $\|x^{k}-x^{k+1}\|\leq\frac{1}{mp_{\min}}\|x^{k}-\bar{x}^{k+1}\|$, we have $x^{k}-x^{k+1}\rightarrow 0$ a.s.. Then from (5), we have $\hat{x}^{k}-x^{k}\rightarrow 0$ a.s..

(iii): From Lemma 3, we have that $(\|{\mathbf{x}}^{k}-{{\mathbf{x}}}^{*}\|_{M}^{2})_{k\geq 0}$ converges a.s. and so does $(\|{\mathbf{x}}^{k}-{{\mathbf{x}}}^{*}\|_{M})_{k\geq 0}$, i.e., $\lim_{k\rightarrow\infty}\|{\mathbf{x}}^{k}-{{\mathbf{x}}}^{*}\|_{M}=\gamma$ a.s., where $\gamma$ is a $[0,+\infty)$-valued random variable. Hence, $(\|{\mathbf{x}}^{k}-{{\mathbf{x}}}^{*}\|_{M})_{k\geq 0}$ must be bounded a.s. and so is $({\mathbf{x}}^{k})_{k\geq 0}$.

(v): By (ii), there exists $\hat{\Omega}\in{\mathcal{F}}$ such that $P(\hat{\Omega})=1$ and

For any $\omega\in\hat{\Omega}$, let $({\mathbf{x}}^{k_{n}}(\omega))_{n\geq 0}$ be a weakly convergent subsequence of $({\mathbf{x}}^{k}(\omega))_{k\geq 0}$, i.e., ${\mathbf{x}}^{k_{n}}(\omega)\rightharpoonup{\mathbf{x}}$, where ${\mathbf{x}}^{k_{n}}(\omega)=(x^{k_{n}}(\omega),x^{k_{n}-1}(\omega)...,x^{k_{n}-\tau}(\omega))$ and ${\mathbf{x}}=(u^{0},...,u^{\tau})$. Note that ${\mathbf{x}}^{k_{n}}(\omega)\rightharpoonup{\mathbf{x}}$ implies $x^{k_{n}-j}(\omega)\rightharpoonup u^{j},\,\forall j.$ Therefore, $u^{i}=u^{j}$, for any $i,j\in\{0,\cdots,\tau\}$ because $x^{k_{n}-i}(\omega)-x^{k_{n}-j}(\omega)\rightarrow 0$.

Furthermore, observing $\eta_{k}\geq\eta_{\min}>0$, we have

From the triangle inequality and the nonexpansiveness of $T$, it follows that

From (44), (45), and the above inequality, it follows $\lim_{n\rightarrow\infty}x^{k_{n}}(\omega)-Tx^{k_{n}}(\omega)=0.$ Finally, the demiclosedness principle (bauschke2011convex, , Theorem 4.17) implies $u^{0}\in\operatorname*{Fix}T$.

Under the assumptions of Lemma 5, the sequence $({\mathbf{x}}^{k})_{k\geq 0}$ weakly converges to an ${\mathbf{X}}^{*}$-valued random variable a.s.. In addition, if $T$ is demicompact at 0, $({\mathbf{x}}^{k})_{k\geq 0}$ strongly converges to an ${\mathbf{X}}^{*}$-valued random variable a.s..

For the generalization in Section 1.3, we need to replace (35) by

and update the step size condition to $\eta_{k}\in[\eta_{\min},\frac{cmp_{\min}}{2\tau\sqrt{p_{\min}}+C}]$. Then the proofs of Theorem 3.2 and Lemma 5 will go through and yield the same convergence result.

2 Linear convergence

In this section, we establish linear convergence under the assumption that $S$ is quasi-strongly monotone. We first present a key lemma.

Assume that the step size is fixed, i.e., $\eta_{k}=\eta$, and satisfies

for some $\rho>1$. Then we have, for all $k\geq 1$,

We prove (47) by induction. First, based on the inequality $\|a\|^{2}-\|b\|^{2}\leq 2\|a\|\|b-a\|$ we observe that, for any $k\geq 1$,

Applying the triangle inequality and (5) yields

For the basic case, we have $\hat{x}^{0}=x^{0}$, $\hat{x}^{1}\in\{x^{0},x^{1}\}$. Letting $k=1$ in (48) gets us

For the induction step, applying Young’s inequality gives us

Taking the expectation on (51) and combining it with (48) yield

With this lemma, we are ready to derive the linear convergence rate of ARock.

Assume that $S$ is quasi-$\mu$-strongly monotone with $\mu>0$. Let $\beta\in(0,1)$ and $(x^{k})_{k\geq 0}$ be the sequence generated by ARock with a constant stepsize $\eta\in(0,\min\{\underline{\eta}_{1},\underline{\eta}_{2}\}]$, where $\underline{\eta}_{1}$ is given in (46) and

Following the proof of Lemma 2 and starting from (36), we have

where the second inequality holds because $S$ is $\frac{1}{2}$-cocoercive and also quasi-$\mu$-strongly monotone, and the last one comes from the Cauchy-Schwartz inequality. Plugging the above inequality and (35) into (34) and noting $|J(k)|\subset\{k-\tau,\ldots,k-1\}$ gives

where we have let $\gamma=m\sqrt{\frac{\tau(\rho-1)p_{\min}}{\rho(\rho^{\tau}-1)}}$ in the second equality, and the last inequality holds because of the choice of $\eta$. Therefore, (55) holds.

Assume $i_{k}$ is chosen uniformly at random, so $p_{\min}=\frac{1}{m}$. We consider the case when $m$ and $\tau$ are large. Let $\sqrt{\rho}=1+\frac{1}{\tau}$. Then from the fact that $(1+\frac{1}{k})^{k}$ increasingly converges to the natural number $e$, we have from (46) that $\underline{\eta}_{1}=O(\frac{\sqrt{m}}{\tau^{2}})$. In addition, note from (54) that $a=O(b^{2})=O(\frac{\tau^{2}}{m})$, and thus $\underline{\eta}_{2}=O(\frac{\sqrt{m}}{\tau})$. Therefore, if $\tau=O(m^{\frac{1}{4}})$, then the stepsize in Theorem 3.3 can be $\eta=O(1)$. Hence, linear speedup can be achieved.

Experiments

Our experiments run on 1 to 32 threads on a machine with eight Quad-Core AMD OpteronTM Processors (32 cores in total) and $64$ Gigabytes of RAM. All of the experiments were coded in C++ and OpenMP. We use the Eigen libraryhttp://eigen.tuxfamily.org for sparse matrix operations. Our codes as well as numerical results for other applications will be publicly available on the authors’ website.

The running times and speedup ratios of both sync-parallel and async-parallel algorithms are sensitive to a number of factors, such as the size of each coordinate update (granularity), sparsity of the problem data, compiler optimization flags, and operations that affect cache performance and memory access contention. In addition, since all agents in the sync-parallel implementation must wait for the last agent to finish an iteration, a large load imbalance will significantly degrade the performance. We do not have the space in this paper to present numerical results under all variations of these cases.

where $\{(a_{i},b_{i})\}_{i=1}^{N}$ is the set of sample-label pairs with $b_{i}\in\{1,-1\}$, $\lambda=0.0001$, and $n$ and $N$ represent the numbers of features and samples, respectively. This test uses the datasetshttp://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/: rcv1 and news20, which are summarized in Table 1.

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 agent draw a coordinate uniformly at random at each iteration. We stop all the tests after 100 epochs since they have nearly identical progress per iteration. The step size is set to $\eta_{k}=0.9,\,\forall k$. Let $A=[a_{1},\ldots,a_{N}]^{T}$ and $b=[b_{1},...,b_{N}]^{T}$. 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 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 2 gives the running times of the sync-parallel and ARock (async-parallel) implementations on the two datasets. We can observe that ARock achieves almost-linear speedup, but sync-parallel scales very poorly as we explain below.

In the sync-parallel implementation, all the running cores have to wait for the last core to finish an iteration, and therefore if a core has a large load, it slows down the iteration. Although every core is (randomly) assigned to roughly the same number of features (either 50 or 51 components of $x$) at each iteration, their $a_{i}$’s have very different numbers of nonzeros (see Figure 3 for the distribution), and the core with the largest number of nonzeros is the slowest (Sparse matrix computation is used for both datasets, which are very large.) As more cores are used, despite that they altogether do more work at each iteration, the per-iteration time reduces as the slowest core tends to be slower. The very large imbalance of load explains why the 32 cores only give speedup ratios of 4.0 and 1.3 in Table 2.

On the other hand, being asynchronous, ARock does not suffer from the load imbalance. Its performance grows nearly linear with the number of cores. In theory, a large load imbalance may cause a large $\tau$, and thus a small $\eta_{k}$. However, the uniform $\eta_{k}=0.9$ works well in all the tests, possibly because the $a_{i}$’s are sparse.

Finally, we have observed that the progress toward solving (56) is mainly a function of the number of epochs and does not change appreciably when the number of cores increases or between sync-parallel and async-parallel. Therefore, we always stop at 100 epochs.

Conclusion

We have proposed an async-parallel framework, ARock, for finding a fixed-point of a nonexpansive operator by coordinate updates. We establish the almost sure weak and strong convergence, linear convergence rate and almost-linear speedup of ARock under certain assumptions. Preliminary numerical results on real data illustrate the high efficiency of the proposed framework compared to the traditional parallel (sync-parallel) algorithms.

Acknowledgements

We would like to thank Brent Edmunds for offering invaluable suggestions on the organization and writing of this paper. We would also like to thank Robert Hannah for coming up with the dual-memory approach. The authors are grateful to Kun Yuan for helpful discussions on decentralized optimization.

References

Appendix A Derivation of certain updates

We show in details how to obtain the updates in (16) and (21).

Let $x=(x_{1},\ldots,x_{m})\in{\mathcal{H}}^{m}$,

where $g(x)$ equals if $x_{1}=\cdots=x_{m}$ and $\infty$ otherwise. Then (15a) reduces to

where the last equality is obtained by noting that $z_{1}=\frac{1}{m}\sum_{i=1}^{m}\hat{z}_{i}^{k}$ is the unique minimizer of $\sum_{i=1}^{m}\|z_{1}-\hat{z}_{i}^{k}\|^{2}$. Next, (15b) reduces to

Since (15c) only updates the $i_{k}$th coordinate of $z$, we only need $\hat{x}_{i_{k}}^{k}$ and $\hat{y}_{i_{k}}^{k}$, and thus in (16a) and (16b), we only compute $\hat{x}_{i_{k}}^{k}$ and $\hat{y}_{i_{k}}^{k}$. Plugging the above $\hat{x}^{k}$ and $\hat{y}^{k}$ into (15c) gives (16c) directly.

A.2 Derivation of (21)

We first show how to get (18). The Lagrangian of (17) is $L(x,y,w)=f(x)+g(y)-\langle w,Ax+By-b\rangle,$ and the Lagrange dual function is

where the last equality is from the definition of convex conjugate: $f^{*}(z)=\max_{x}\langle z,x\rangle-f(x).$ Hence, the dual problem is $\max_{w}d(w)$, which is equivalent to (18).

Secondly, we show why $z^{+}=\mathbf{prox}_{\gamma\cdot d_{g}}(z)$ is given by (20). Note

where the fifth equality holds because $s^{*}=z-\gamma(By-b)=\arg\min_{s}\langle s,By-b\rangle+\frac{1}{2\gamma}\|s-z\|^{2}.$ Hence, by the definition of the proximal operator and the above arguments, we have that $z^{+}=\mathbf{prox}_{\gamma\cdot d_{g}}(z)$ can be obtained from (20). Then (19) is from (20) through replacing $g$ to $f$, $B$ to $A$, and $b$ to .

Finally, it is straightforward to have (21) by plugging (19) and (20) into (15).

Appendix B Derivation of async-parallel ADMM for decentralized optimization

This section describes how to implement the updates (21) for the model (29).

In (29), $g(y)$ and $b$ vanish and, corresponding to the two constraints $x_{i}=y_{ij}$ and $x_{j}=y_{ij}$, the two rows of matrices $A$ and $B$ are $\begin{bmatrix}\cdots&1&\cdots&\cdots&\cdots\\ \cdots&\cdots&\cdots&1&\cdots\end{bmatrix}\quad\begin{bmatrix}\cdots&-1&\cdots\\ \cdots&-1&\cdots\end{bmatrix},$ where $\cdots$ are zeros, the two coefficients 1 correspond to $x_{i}$ and $x_{j}$, and the two coefficients $-1$ correspond to $y_{ij}$. Then, (21a) and (21b) can be calculated as