A Distributed, Asynchronous and Incremental Algorithm for Nonconvex Optimization: An ADMM Based Approach

Mingyi Hong

I Introduction

Consider the following nonconvex and nonsmooth problem

where gkg_{k}’s are a set of smooth, possibly nonconvex functions; h(x)h(x) is a convex nonsmooth regularization term. In this paper we consider the scenario where the component functions gkg_{k}’s are located at different distributed computing nodes. We seek an algorithm that is capable of computing high quality solutions for problem (1) in a distributed, asynchronous and incremental manner.

Dealing with asynchrony is a central theme in designing distributed algorithms. Indeed, often in a completely decentralized setting, there is no clock synchronization, little coordination among the distributed nodes, and minimum mechanism to ensure reliable communication. Therefore an ideal distributed algorithm should be robust enough to handle different sources of asynchrony, while still producing high quality solutions in a reasonable amount of time. Since the seminal work of Bertsekas and Tsitsiklis , there has been a large body of literature focusing on asynchronous implementation of various distributed schemes; see, e.g., for the developments by the optimization and signal processing communities. In , an incremental and asynchronous gradient-based algorithm is proposed to solve a convex problem, where at each step certain outdated gradients can be used for update. In , the authors show that the well-known iterative water-filling algorithm can be implemented in a totally asynchronous manner, as long as the interference among the users are weak enough.

The recent interest in optimization and machine learning for problems with massive amounts of data introduces yet another compelling reason for dealing with asynchrony; see [10, Chapter 10]. When large amounts of data are distributedly located at computing nodes, local computations can be costly and time consuming. If synchronous algorithms are used, then the slowest nodes can drag the performance of the entire system. To make distributed learning algorithms scalable and efficient, the machine learning community has also started to deal with asynchrony; see recent results in . For example in , an asynchronous randomized block coordinate descent method is developed for solving convex block structured problem, where the per-block update can utilize delayed gradient information. In , the authors show that it is also possible to tolerate asynchrony in stochastic optimization. Further, they prove that the rate of the convergence is more or less independent of the maximum allowable delay, which is an improvement over earlier results in .

In this paper, we show that through the lens of the ADMM method, the nonconvex and nonsmooth problem (1) can be optimized in an asynchronous, distributed, and incremental manner. The ADMM, originally developed in early 1970s , has been extensively studied in the last two decades . It is known to be effective in solving large-scale linearly constrained convex optimization problems. Its application includes machine learning, computer vision, signal and image processing, networking, etc; see . However, despite various successful numerical attempts (see, e.g., ), little is known about whether ADMM is capable of handling nonconvex optimization problems, or whether it can be used in an asynchronous setting. There are a few recent results that start to fill these gaps. Reference shows that the ADMM converges when applied to certain nonconvex consensus and sharing problems, provided that the stepsize is chosen large enough. However it is not clear whether asynchrony will destroy the convergence. Reference proposes an asynchronous implementation for convex global consensus problem, where the distributed worker nodes can use outdated information for updates. Two conditions are imposed on the protocol, namely the partial barrier and bounded delay. The algorithm cannot deal with the asynchrony cause by loss/delay in the communication link, nor does it cover nonconvex problems. In randomized versions of ADMM are proposed for consensus problems, where the nodes are allowed to be randomly activated for updates. We note that the algorithms in still require the nodes to use up-to-date information whenever they update, therefore they are more in line with randomized algorithms than asynchronous algorithms. Further, it is not known whether the analysis carries over to the case when the problem is nonconvex.

The algorithm proposed in this work is a generalization of the flexible proximal ADMM algorithm proposed in [43, Section 2.3]. The key feature of the proposed algorithm is that it can deal with asynchrony arises from the heterogeneity of the computing nodes as well as the loss/delay caused by unreliable communication links. The basic requirement here is that the combined effects of these sources leads to a bounded delay on the component gradient evaluation, and that the stepsize of the algorithm is chosen appropriately. Further, we show that the framework studied here can be viewed as an (possibly asynchronous) incremental scheme for nonconvex problem, where at each iteration only a subset of (possibly delayed) component gradients are updated. To the best of our knowledge, asynchronous incremental schemes of this kind hasn’t been studied in the literature; see for recent works on synchronous incremental algorithm for nonconvex problems.

II The ADMM-based Framework

Consider the optimization problem (1). In many practical applications, gkg_{k}’s need to be handled by a single distributed node, such as a thread or a processor, which motivates the so-called global consensus formulation [49, Section 7]. Suppose there is a master node and KK distributed nodes available. Let us introduce a set of new variables {xk}k=1K\{x_{k}\}_{k=1}^{K}, and transform problem (1) to the following linearly constrained problem

The augmented Lagrangian function is given by

where ρk>0\rho_{k}>0 is some constant, and y:={y1,⋯ ,yK}y:=\{y_{1},\cdots,y_{K}\}. Applying the vanilla ADMM algorithm, listed below in (4), one obtains a distributed solution where each function gkg_{k} is only handled by a single node kk at any iteration t=0,1,2,…t=0,1,2,\ldots:

Under suitable conditions the algorithm converges to the set of stationary solutions of (1); see .

At this point, it is important to note that the algorithm described in (4) uses a synchronous protocol, that is

The set of agents that are selected to update at each iteration act in a coordinated way;

There is no communication delay and/or loss between the agents and the master node;

All local updates are performed assuming that the most up-to-date information is available.

However, in many practical large-scale networks, these assumptions are hardly true. Nodes may have different computational capacity, or they may be assigned jobs that have different computational requirements. Therefore the time consumed to complete local computation can vary significantly among the nodes. This makes them difficult to coordinate with each other in terms of when to update, which information to use for the update and so on. Further, the communication links between the distributed and the master nodes can have delays or may even be lossy.

Additionally, we want to mention that in certain machine learning and signal processing problems when there is a large number of component functions, it is desirable that the algorithm is incremental, meaning at each iteration only a subset of gk′sg_{k}^{\prime}s are used for update; see . Clearly the vanilla ADMM described in (4) does not belong to this type of algorithm.

II-B The Proposed Algorithm

There are two key features that we want to build into the ADMM-based algorithm. One is to allow the nodes to use staled information for local computation, as long as such information is not “too old” (this notion will be made precise shortly). This enables the nodes to have varying update frequency, therefore faster nodes do not need to wait for the slower ones. The other feature is to take into account scenarios where the communication links among the node are lossy or have delays. Below we give a high level description of the proposed scheme.

Suppose there is a master node and KK distributed nodes in the system. Let the index t=1,⋯t=1,\cdots denote the total number of updates that have been performed on the variable xx. The master node takes care of updating all the primal and dual variables, while the distributed nodes compute the gradients for each component function gkg_{k}. At each iteration t+1t+1, the master node first updates xx. Then it waits a fixed period of time, collects a few (possibly staled) gradients of component functions returned by a subset of local nodes {\mbox{\mathcal{C}}}^{t+1}\subseteq\{1,\cdots,K\}, then proceed to the next round of update. On the other hand, each node kk is in charge of a local component function gkg_{k}. Based on the copy of xx passed along by the master node, node kk computes and returns the gradient of gkg_{k} to the master node. Note that for data intensive applications, the computation of the gradient can be time consuming. Also there can be delays of communication between two different nodes in the network. Therefore there is no guarantee that during the period of computation and communication of the gradient of gkg_{k}, the xx variable at the master node will always remain the same.

To characterize the possible delay involved in the computation and communication, we define a new sequence {t(k)}\{t(k)\}, where each t(k)t(k) represents the index of the copy of xx that evaluates the ∇gk\nabla g_{k} used by the master node at iteration tt.

The proposed algorithm, named Asynchronous Proximal ADMM (Async-PADMM), is given in the following table.

Algorithm 1. The Async-PADMM for Problem (2) S1) At each iteration t+1t+1, compute: \displaystyle\begin{split}x^{t+1}&={\rm arg}\!\min_{x\in X}\;L(\{x_{k}^{t}\},x;y^{t})\\ &=\mbox{prox}_{\iota(X)+h}\left[\frac{\sum_{k=1}^{K}\rho_{k}x_{k}^{t}+\sum_{k=1}^{K}y_{k}^{t}}{\sum_{k=1}^{K}\rho_{k}}\right].\end{split} (5) S2) Pick a set {\mbox{\mathcal{C}}}^{t+1}\subseteq\{1,\cdots,K\}, for all k\in{\mbox{\mathcal{C}}}^{t+1}, update index [t+1](k)[t+1](k); for all k\notin{\mbox{\mathcal{C}}}^{t+1}, let [t+1](k)=[t](k)[t+1](k)=[t](k) S3) Update xkx_{k} by solving: xkt+1\displaystyle x^{t+1}_{k} =arg⁡ ⁣min⁡xk  ⟨∇gk(x[t+1](k)),xk−xt+1⟩+⟨ykt,xk−xt+1⟩\displaystyle=\arg\!\min_{x_{k}}\;\langle\nabla g_{k}(x^{[t+1](k)}),x_{k}-x^{t+1}\rangle+\langle y^{t}_{k},x_{k}-x^{t+1}\rangle +ρk2∥xk−xt+1∥2,  ∀ k=1,⋯ ,K.\displaystyle+\frac{\rho_{k}}{2}\|x_{k}-x^{t+1}\|^{2},\;\forall~{}k=1,\cdots,K. (6) S4) Update the dual variable: ykt+1=ykt+ρk(xkt+1−xt+1),  ∀ k=1,⋯ ,K.\displaystyle y^{t+1}_{k}=y_{k}^{t}+\rho_{k}\left(x^{t+1}_{k}-x^{t+1}\right),\;\forall~{}k=1,\cdots,K. (7)

We note that in Step S2, {\mbox{\mathcal{C}}}^{t+1} defines the subset of component functions whose gradients have arrived during iteration t+1t+1; again [t+1](k)[t+1](k) is the index of the copy of xx that evaluates the ∇gk\nabla g_{k} used by the master node at iteration t+1t+1. For those component functions without new gradient information available, the old gradients will continue to be used (indeed, note that we have for all k\notin{\mbox{\mathcal{C}}}^{t+1}, [t+1](k)=[t](k)[t+1](k)=[t](k)). In Step S3, all the variables {xk}\{x_{k}\}, regardless k\in{\mbox{\mathcal{C}}}^{t+1} or not, are updated according to the following gradient-type scheme:

Despite the fact that the gradients of all the component functions are used at each step t+1t+1, only a subset of them (i.e., thosed indexed by {\mbox{\mathcal{C}}}^{t+1}) differ from those at the previous iteration. Therefore the algorithm can be classified as incremental algorithm; see for related incremental algorithms for convex problems.

To highlight the asynchronous aspect of the algorithm, below we present an equivalent version of Algorithm 1, from the perspective of the distributed nodes and the master node, respectively. We use rkr_{k}, k=1,⋯ ,Kk=1,\cdots,K to denote the clock at node kk, and use r0r_{0} to denote the clock at the master node.

Algorithm 1(a). Async-PADMM at the Master Node S0) Set r0=1r_{0}=1, initialize {xk1,yk1}\{x^{1}_{k},y^{1}_{k}\}, x1x^{1}. S1) Update xx: xr0+1\displaystyle x^{r_{0}+1} =arg ⁣min⁡x∈X  L({xkr0},x;yr0)\displaystyle={\rm arg}\!\min_{x\in X}\;L(\{x_{k}^{r_{0}}\},x;y^{r_{0}}) =\mboxproxι(X)+h[∑k=1Kρkxkr0+∑k=1Kykr0∑k=1Kρk].\displaystyle=\mbox{prox}_{\iota(X)+h}\left[\frac{\sum_{k=1}^{K}\rho_{k}x_{k}^{r_{0}}+\sum_{k=1}^{K}y_{k}^{r_{0}}}{\sum_{k=1}^{K}\rho_{k}}\right]. (10) S2) Broadcast xr0+1x^{r_{0}+1} to all agents. S3) Wait for a fixed period of time. S4) Collect a set {\mbox{\mathcal{C}}}^{r_{0}+1}\subseteq\{1,\cdots,K\} of new local gradients, denoted as \{z^{r_{0}+1}_{k}\}_{k\in{\mbox{\mathcal{C}}}^{r_{0}+1}}, arrived during S3). If multiple gradients arrive from the same node, pick the one with the smallest local time stamp. S5) Let ∇gkr0+1=zkr0+1\nabla g^{r_{0}+1}_{k}=z^{r_{0}+1}_{k}, \forall~{}k\in{\mbox{\mathcal{C}}}^{r_{0}+1}. S6) Let ∇gkr0+1=∇gkr0\nabla g^{r_{0}+1}_{k}=\nabla g^{r_{0}}_{k}, \forall~{}k\notin{\mbox{\mathcal{C}}}^{r_{0}+1}. S7) Compute \displaystyle\begin{split}x_{k}^{r_{0}+1}&=x^{r_{0}+1}-\frac{1}{\rho_{k}}\left(\nabla g_{k}^{r_{0}+1}+y^{r_{0}}_{k}\right),\;\forall~{}k\\ y^{r_{0}+1}_{k}&=y_{k}^{r_{0}}+\rho_{k}\left(x^{r_{0}+1}_{k}-x^{r_{0}+1}\right),\;\forall~{}k.\end{split} (11) S8) Set r0=r0+1r_{0}=r_{0}+1, go to step S1).

Algorithm 1(b). The Async- PADMM at Node kk S0) Set rk=1r_{k}=1. S1) Wait until a new xx is arrived, mark it as xrkx^{r_{k}}. S2) Compute the gradient ∇gk(xrk)\nabla g_{k}(x^{r_{k}}). S3) Send ∇gk(xrk)\nabla g_{k}(x^{r_{k}}) and the local time stamp rkr_{k} to the master node. S4) Set rk=rk+1r_{k}=r_{k}+1, go to step S1).

It is not hard to see that the scheme described here is equivalent to Algorithm 1, except that in Algorithm 1 every step is measured using the clock at the master node. We have the following remarks regarding to the above algorithm descriptions.

(Blocking Events) There is a minimal number of blocking events for both the master node and the distributed agents. In Algorithm 1(a), the master node only needs to wait for a given period of time in step S3). After the waiting period, it collects the set of new gradients that has arrived during that period. Note that {\mbox{\mathcal{C}}}^{r_{0}+1} is allowed to be an empty set, meaning the master node is not blocking on the arrival of any local gradients. Similarly, each node kk does not need to wait for the rest of the agents to perform computation: once it obtains a new copy of xrk+1x^{r_{k}+1} the computation starts immediately. As soon as the computation is done node kk can send out the new gradient, without checking whether that gradient has arrived at the master node. Admittedly, in Step S1 of Algorithm 1(b), node kk needs to wait for a new xx, but this is reasonable because otherwise there is nothing it can do.

(Characterization on the Delays) The proposed algorithm allows communication delays and packet loss between the master and the distributed nodes. For example, the vector xt+1x^{t+1} broadcasted by the master node may arrive at the different distributed nodes at different time instances; it may even arrive at a given node out of order, i.e., xt+1x^{t+1} arrives before xtx^{t}. Further, xt+1x^{t+1} may get lost during the transmission and never reaches a given node. All these scenarios can happen in the reverse communication direction as well. Comparing Algorithm 1 and Algorithm 1(a)–(b), we see that if k\in{\mbox{\mathcal{C}}}^{t+1}, then the difference (t+1)−[t+1](k)(t+1)-[t+1](k) is the total computation time and the round-trip communication delay, starting from broadcasting x[t+1](k)x^{[t+1](k)} until the updated ∇gk(x[t+1](k))\nabla g_{k}(x^{[t+1](k)}) is received by the master node. If k\notin{\mbox{\mathcal{C}}}^{t+1}, then the difference (t+1)−[t+1](k)(t+1)-[t+1](k) is the number of times that the gradient ∇gk(x[t+1](k))\nabla g_{k}(x^{[t+1](k)}) has been used so far (or equivalently the number of iterations since the last gradient from node kk has arrived). Clearly, when there is no delay at all , then the system is synchronous and we have [t+1](k)=t+1[t+1](k)=t+1. In Fig. 1, we illustrate the relationship tt and t(k)t(k), and different types of asynchronous events covered by the algorithm.

(Connection to Existing Algorithms) To the best of our knowledge, the proposed algorithm can tolerate the highest degree of asynchrony, among all known asynchronous variants of ADMM. For example, the scheme proposed in corresponds to the case where there is no communication delay or loss (all messages sent are received instantaneously by the intended receiver). It is not clear whether the scheme in can be generalized to our caseIn fact, no proof is provided in . Therefore it becomes difficult to see whether it is possible to extend their analysis.. The schemes proposed in and require the nodes to use the most up-to-date information, hence hardly asynchronous. The second major difference with the existing literature is about the tasks performed by the distributed nodes: in each node directly optimizes the augmented Lagrangian, while here each node computes the gradient of their respective component functions. The third difference is on the assumptions made on problem (1): the schemes in handle convex problem but each component function gig_{i} can be nonsmooth, while we can handle nonconvex functions, but there can be only a single nonsmooth function hh (see Assumption A1 below). The fourth difference is on the assumed network topology: the schemes in deal with general topology, where nodes are interconnected according to certain graphs; our work and are restricted to the “star” network topology where all distributed nodes communicate directly with the master node.

(Incrementalism) Algorithm 1 can be viewed as an incremental algorithm, as long as each |{\mbox{\mathcal{C}}}^{t+1}| is a strict subset of {1,⋯K}\{1,\cdots K\}, in which case the gradients of only a subset of component functions are updated. This is in the same spirit of several recent incremental algorithms for convex problems , despite the fact that our algorithm has a different form, and we can further handle nonconvexity and asychrony.

It is worth noting that Algorithm 1 can be modified to resemble the more traditional incremental algorithm , where each iteration only those variables with “fresh” gradients are updated. That is, steps S3 and S4 are replaced with the following steps:

S3)’ Update xkx_{k} by solving: xkt+1\displaystyle x^{t+1}_{k} =arg⁡ ⁣min⁡xk  ⟨∇gk(x[t+1](k)),xk−xt+1⟩+⟨ykt,xk−xt+1⟩\displaystyle=\arg\!\min_{x_{k}}\;\langle\nabla g_{k}(x^{[t+1](k)}),x_{k}-x^{t+1}\rangle+\langle y^{t}_{k},x_{k}-x^{t+1}\rangle \displaystyle+\frac{\rho_{k}}{2}\|x_{k}-x^{t+1}\|^{2},\;\forall~{}k\in{\mbox{\mathcal{C}}}^{t+1}. S4)’ Update the dual variable: \displaystyle y^{t+1}_{k}=y_{k}^{t}+\rho_{k}\left(x^{t+1}_{k}-x^{t+1}\right),\;\forall~{}k\in{\mbox{\mathcal{C}}}^{t+1}.

However, we found that this variant leads to much more complicated analysis To analyze this version, we need to define a few additional sequences, one for each node kk, to characterize the iteration indices in which each component variable xkx_{k} is updated. We will also need to impose that the xkx_{k}’s are updated often enough; see [1, Chapter 7]., stringent requirement on the range of stepsizes ρk\rho_{k}’s, and most importantly, slow convergence. Therefore we choose not to discuss the related variants in the paper. We also note that recent works in incremental-type algorithms for solving (1) either do not deal with nonconvex problem , or they do not consider asynchrony .

III Convergence Analysis

In order to reduce the notational burden, our analysis will be based on Algorithm 1, which uses a global clock. We first make a few assumptions.

(On the Problem) There exists a positive constant Lk>0L_{k}>0 such that

Moreover, hh is convex (possibly nonsmooth); XX is a closed, convex and compact set. f(x)f(x) is bounded from below over XX.

(On the Asynchrony) The total delays are bounded, i.e., for each node kk there exists finite constants TkT_{k} such that t−t(k)≤Tkt-t(k)\leq T_{k} for all tt and kk.

(On the Algorithm) For all kk, the stepsize ρk\rho_{k} is chosen large enough such that:

By Assumption A2, we see that the only requirement on the asynchrony is that when each xkx_{k} is updated, the information used to compute the gradient should be one of xx generated within last TkT_{k} iterations. So it is perfectly legitimate if copies of xx or copies of the gradients get lost due to unsuccessful communication. Also there is nothing preventing copies of xx from arriving at the same node kk with reversed order (e.g., xtx^{t} arrives after xt+1x^{t+1}). Due to this assumption on the boundedness of the asynchrony, Algorithms 1 belongs to the family of “partially asynchronous algorithm”, as opposed to the “totally asynchronous algorithm” in which the delays can potentially be unbounded In short, the only requirement for the totally asynchronous algorithm is that no nodes quits forever. ; see the definitions and discussions in .

From Assumption A3, it is clear that when the system is synchronous, i.e., when Tk=0T_{k}=0, the bound for αk\alpha_{k} becomes

Suppose Assumption A is satisfied. Then for Algorithm 1, the following is true for all kk

Proof. From the update of xkx_{k} in (II-B), we observe that the following is true

Note that both xkx_{k} and yky_{k} are updated at each iteration, so we have the following equality for iteration tt as well

Suppose k\notin{\mbox{\mathcal{C}}}^{t+1}, which means that no new gradient information arrives for node kk. In this case, we have [t+1](k)=t(k)[t+1](k)=t(k), therefore

It follows that for k\notin{\mbox{\mathcal{C}}}^{t+1}, (15) is true.

Suppose that k\in{\mbox{\mathcal{C}}}^{t+1}, then we have

Therefore we have, for all k\in{\mbox{\mathcal{C}}}^{t+1}

Next, we upper bound the successive difference of the augmented Lagrangian. To this end, let us define a few new functions, given below

Using these short-handed definitions, we have

Suppose Assumption A is satisfied. Let {xkt,xt,yt}\{x^{t}_{k},x^{t},y^{t}\} be generated by Algorithm 1. Then we have the following

Observe that xkt+1x^{t+1}_{k} is generated according to (24). Combined with the strong convexity of uˉk(xk;xt+1,yt)\bar{u}_{k}(x_{k};x^{t+1},y^{t}) with respect to xkx_{k}, we have

Using the strong convexity of uku_{k}, we have the series of inequalities given below

Further, we have the following series of inequalities

where the first two inequalities follow from Assumption A1. Combining (26) – (III) we obtain

Next, we bound the difference of the augmented Lagrangian function values.

Assume the same set up as in Lemma III.2. Then we have

where αk\alpha_{k} is the constant defined in (12).

Proof. We first bound the successive difference L({xkt+1},xt+1;yt+1)−L({xkt},xt;yt)L(\{x^{t+1}_{k}\},x^{t+1};y^{t+1})-L(\{x^{t}_{k}\},x^{t};y^{t}). We first decompose the difference by

The first term in (33) can be expressed as

To bound the second term in (33), we use Lemma III.2. We have the series of inequalities in (34), where the last inequality follows from Lemma III.2 and the strong convexity of L({xkt},x;yt)L(\{x^{t}_{k}\},x;y^{t}) with respect to the variable xx (with modulus γ=∑k=1Kρk\gamma=\sum_{k=1}^{K}\rho_{k}) at x=xt+1x=x^{t+1}.

Combining the above two inequalities and use Lemma III.1, we obtain the inequality below:

Then for any given tt, the difference L({xkt+1},xt+1;yt+1)−L({xk1},x1;y1)L(\{x^{t+1}_{k}\},x^{t+1};y^{t+1})-L(\{x^{1}_{k}\},x^{1};y^{1}) is obtained by summing (III) over all iterations:

We conclude that to make the augmented Lagrangian decrease at each iteration, it is sufficient to require that αk>0\alpha_{k}>0 and ρk−7Lk>0\rho_{k}-7L_{k}>0 for all kk. Note that one can always find a set of ρk\rho_{k}’s large enough such that the above condition is satisfied.

Next we show that L({xkt},xt;yt)L(\{x^{t}_{k}\},x^{t};y^{t}) is convergent.

Suppose Assumption A is satisfied. Then Algorithm 1 generates a sequence of augmented Lagrangian that satisfies

where \mboxdiam(X):=sup⁡{∥x1−x2∥∣x1,x2∈X}\mbox{diam}(X):=\sup\{\|x_{1}-x_{2}\|\mid x_{1},x_{2}\in X\}, which is the diameter of the set XX.

Proof. Observe that the augmented Lagrangian can be expressed as

In the above series of inequalities, (a)\rm(a) is from (17); (b)\rm(b) is due to the Cauchy-Schwartz inequality, Assumption A1, and the following inequalities

The inequality in (c)(c) is due to the assumption that ρk≥4Lk\rho_{k}\geq 4L_{k}, and by the definition of the \mboxdiam(X)\mbox{diam}(X); (d)(d) is because of the assumption that f(x)f(x) is bounded over all XX, and that XX is a compact set. It follows from Lemma III.3 that whenever the stepsize ρk\rho_{k}’s are chosen sufficiently large (as per Assumption A), L({xkt+1},xt+1;yt+1)L(\{x^{t+1}_{k}\},x^{t+1};y^{t+1}) will monotonically decrease and is convergent. This completes the proof.

Using Lemmas III.1–III.4, we arrive at the following convergence result.

Suppose that Assumption A holds. Then the following is true for Algorithm 1.

We have lim⁡t→∞∥xt+1−xkt+1∥=0\lim_{t\to\infty}\|x^{t+1}-x^{t+1}_{k}\|=0, k=1,⋯ ,Kk=1,\cdots,K. That is, the primal feasibility is always satisfied in the limit.

The sequence {{xkt+1},xt+1,yt+1}\{\{x^{t+1}_{k}\},x^{t+1},y^{t+1}\} converges to the set of stationary solutions of problem (2). Moreover, the sequence {{xkt+1},xt+1}\{\{x^{t+1}_{k}\},x^{t+1}\} converges to the set of stationary solutions of problem (1).

Proof. Combining Lemma III.3 – III.4 we must have

Taking limit on both sides of (15) and use the above two results, we immediately obtain

Once we can show that the primal feasibility gap goes to zero, the proof for stationarity is straightforward. We refer the readers to for detailed arguments.

It turns out that for some special cases of gkg_{k}’s, the requirement on the stepsize can be further relaxed.

Suppose Assumption A1 and A3 are true. We have the following:

If gkg_{k} is a convex function, then the corresponding ρk\rho_{k} should satisfy:

If gkg_{k} is a concave function, then the corresponding ρk\rho_{k} should satisfy:

where the last inequality comes from the convexity of gkg_{k}. Similarly, we can replace the series of inequalities in (III) by

where in (a)(a) we have again used the Cauchy-Schwartz inequality and the convexity of gkg_{k}. Then by simple manipulation we arrive at the claimed result.

where (a)(a) comes from the concavity of gkg_{k}.

(On the Bounded Delays) Our convergence results critically dependent on the choice of the stepsizes {ρk}\{\rho_{k}\}, which in turn is a function of the bounds {Tk}\{T_{k}\}. Clearly all TkT_{k}’s must be finite, therefore the scheme proposed here is reminiscent to the family of partially asynchronous algorithm discussed in [1, Chapter 7]. Clearly those bounds on ρk\rho_{k}’s are developed for the worst case delay scenarios. If we model different delays {t−t(k)}\{t-t(k)\} as random variables following certain probability distributions with finite supports, we can slightly modify the analysis so that the final bounds on ρk\rho_{k}’s are dependent on the statistical properties of the random variables. Such modification is minor so we do not intend to go over it in this paper. The more interesting case would be when the delays {t−t(k)}\{t-t(k)\} follow distributions with finite means and variances but unbounded supports. However our current approach cannot be directly used.

(On the Relationship with ) The analysis presented above follows the general recipe first alluded in and later generalized in , for dealing with nonconvex ADMM-type algorithms. The same three-step approach is used here: 1) Bounding the size of the successive difference of the dual variable; 2) Bounding the successive difference of the augmented Lagrangian; 3) Bounding the sequence of the augmented Lagrangian. However, several important improvements have been made to both the algorithm and the analysis in order to better incorporate asynchrony. For example, compared with the flexible Proximal ADMM algorithm in , we have increased the stepsize for updating xkx_{k} from 1ρk+Lk\frac{1}{\rho_{k}+L_{k}} to 1ρk\frac{1}{\rho_{k}}. This change significantly simplifies the analysis and leads to a better bound for ρk\rho_{k}. Second, in the flexible Proximal ADMM, a given tuple (xk,yk)(x_{k},y_{k}) is only updated when the new gradient is available, while here (xk,yk)(x_{k},y_{k}) is updated at every iteration regardless of the availability of new gradients. This also leads to better bound for ρk\rho_{k} and faster algorithm. Third, different analysis techniques have been used throughout to take into consideration the changes in the algorithm as well as the presence of staled gradients.

(On the Necessity of Lemma III.4) As a technical remark, we emphasize that lower-bounding the sequence of the augmented Lagrangian, as we have done in Lemma III.4, is a key step in the entire analysis. The reason is that the compactness of the set XX only guarantees that the sequence of {xt}\{x^{t}\} is bounded, but not the sequences {xkt,ykt}\{x^{t}_{k},y^{t}_{k}\} (note that xktx^{t}_{k} is generated by solving an unconstrained problem). Without the boundedness of {xkt,ykt}\{x^{t}_{k},y^{t}_{k}\}, the augmented Lagrangian L({xkt},xt;yt)L(\{x^{t}_{k}\},x^{t};y^{t}) can go to −∞-\infty, therefore one cannot claim that ∥xt+1−xt∥→0\|x^{t+1}-x^{t}\|\to 0 and ∥xkt−xkt+1∥→0\|x^{t}_{k}-x^{t+1}_{k}\|\to 0.

IV Numerical Results

In this section we conduct numerical experiments to validate the performance of the proposed algorithm.

We consider the following nonconvex problem:

To formulate above sparse PCA problem in the form of (2), we introduce a set of new variable {xk}\{x_{k}\}:

It is straightforward to see that when applying ADMM or the Async-PADMM, each subproblem can be solved in closed form. It is also worth noting that each smooth term in the objective of (46), −12xkTBkTBkxk-\frac{1}{2}x_{k}^{T}B^{T}_{k}B_{k}x_{k}, is a concave function, so the refined the stepsize rule (41) can be used for Async-PADMM.

In our experiment, we compare the Async-PADMM with the following two algorithms

Synchronous ADMM Algorithm: This is the vanilla ADMM algorithm discussed in [43, Section 2.2]. The algorithm can handle nonconvex problems, but its protocol is synchronous. Therefore when the nodes have different computational time, the master node has to wait for all the distributed nodes to complete one iteration of computation before proceeding to the next step. The downside of this approach is that fast nodes have to wait for the slow nodes. The choice of the stepsize ρk\rho_{k} follows the condition given in [43, Assumption A].

Synchronous PADMM Algorithm: This is the period-1 proximal ADMM algorithm discussed in [43, Section 2.3]. Again the algorithm is synchronous. The choice of the stepsize ρk\rho_{k} follows the condition given in [43, Assumption B].

It is worth noting that by simply waiting for the slowest nodes at each iteration, the ADMM and PADMM are capable of handling the asynchrony caused by the computational delay, albeit in a rather inefficient way. However neither of them can deal with the asynchrony caused by imperfect communication link (i.e., loss of messages, out-of-sequence messages, etc). Therefore for fair comparison, in our experiments we only consider scenarios where the communication links are perfect. That is, we require all messages sent by the nodes are perfectly received, in sequence.

where ak(i,j)a_{k}(i,j) and ck(i,j)c_{k}(i,j) follow uniform distribution: ak(i,j)∼\mboxUniform(0,1)a_{k}(i,j)\sim\mbox{Uniform}(0,1) and ck(i,j)∼\mboxUniform(0,1)c_{k}(i,j)\sim\mbox{Uniform}(0,1). In words, there is a probability pk>0p_{k}>0 such that bk(i,j)b_{k}(i,j) is nonzero and follows a Gaussian distribution with random mean and variance. The asynchrony in the system is simulated as follows. For each node kk we assign a distinct TkT_{k} which characterizes the maximum delay for that node. Each time node kk starts to perform computation (i.e., Step S2 in Algorithm 1(b)), the computational delay is randomly drawn from the distribution \mboxUniform(0,Tk)\mbox{Uniform}(0,T_{k}).

To measure the progress of different algorithms, we need the following definitions. For a given iterate xtx^{t}, it is known that the size of the proximal gradient, expressed below, can be used to measure the optimality:

All the algorithms we tested will be terminated when e(xt,{xkt})e\left(x^{t},\{x^{t}_{k}\}\right) reaches below 10−310^{-3}.

IV-B The Results

We first graphically illustrate the convergence behavior of different algorithms. We set N=500N=500, K=10K=10, λ=0\lambda=0, Tk=5T_{k}=5, Mk=100M_{k}=100, pk=0.1p_{k}=0.1 for all kk. In Fig. 2 – 3, the progress of the algorithm is shown by the sequences of the augmented Lagrangian L(xt,{xkt};yt)L(x^{t},\{x^{t}_{k}\};y^{t}) as well as the optimality measure e(xt,{xkt})e\left(x^{t},\{x^{t}_{k}\}\right). First we observe that as predicted in Lemma III.3 and Lemma III.4, the augmented Lagrangian generated by the Async-PADMM algorithm is a decreasing and lower bounded sequence. Second we see that the augmented Lagrangian generated by ADMM (or the PADMM) resembles a stair function. The reason is that between two successive updates, the master node has to wait for a few iterations for the slowest nodes to finish the computation. Nothing is done during such waiting period, leading to constant augmented Lagrangian. Third, we see that the ADMM seems to be able to reduce the augmented Lagrangian quickly, but in terms of the overall optimality measure it takes longer to converge compared with Async-PADMM.

Next we show the averaged convergence behavior of different algorithms, under various different scenarios. Note that each number in the following table is the average of 5050 independent runs of the respective algorithm.

In Table I, we compare the behavior of different algorithms with varying number of distributed nodes. We observe that for all three algorithms, the iteration required for convergence is increasing with the number of nodes KK. It is also clear that the proposed Async-PADMM performs the best (we use underlines to highlight the best result for each scenario).

In Table II, we compare different algorithms under varying degrees of asynchrony. More specifically, in the first four scenarios, we change TkT_{k} from to 99, for all kk. In the last two scenarios, we set all the TkT_{k}’s to be zero, except for Node 10, whose TkT_{k} is either 55 or 1010. This is to test how the algorithms react to a system with a single slow node. We observe that all three algorithms perform well when there is no computational delay (i.e., when Tk=0T_{k}=0). However once delay starts to increase, ADMM and PADMM become slow, and the transition is quite abrupt (for example ADMM triples its convergence time when TkT_{k}’s change from to 33). This is reasonable as ADMM and PADMM are not designed to deal with asynchrony.

In the last set of experiments, we increase the dimensions of the unknown variables and the value of penalization parameter λ\lambda. The results are in Tables III – IV. Again we see that the proposed algorithm works well in both cases.

V Conclusion

In this paper, we propose an ADMM-based algorithm that is capable of solving the nonsmooth and nonconvex problem (1) in a distributed, asynchronous and incremental manner. We show that as long as the stepsize of the primal and dual updates are chosen sufficiently large, the algorithm converges to the set of stationary solutions of the problem. Numerically we show that the proposed algorithm can efficiently deal with the asynchrony arises from distributedly solving certain sparse PCA problem. In the future, we are interested in analyzing the iteration complexity of the algorithm proposed in this paper. That is, we want to bound the maximum number of iterations needed to reach an ϵ\epsilon-stationary solution for problem (1). To the best of our knowledge, such iteration complexity analysis for nonconvex ADMM-type algorithm is not available yet. We are also interested in extending the analysis to asynchronous ADMM algorithm without using the proximal step. Our current analysis is critically dependent on the availability of the gradient information for each component function, therefore cannot be directly applied to the aforementioned case.

VI Acknowledgement

The author wish thank Zhi-Quan Luo from University of Minnesota, and Tsung-Hui Chang from National Taiwan University of Science and Technology, and Xiangfeng Wang from East China Normal University, for helpful discussions.

References