Asynchronous Optimization Over Heterogeneous Networks via Consensus ADMM

Sandeep Kumar, Rahul Jain, Ketan Rajawat

I Introduction

Multi-agent networked systems arise in a number of engineering disciplines such as tactical ad hoc networks, environmental monitoring networks , multi-robot control and tracking , internet-scale monitoring, and large-scale learning . The estimation, resource allocation, and network control tasks required in these applications are often formulated as distributed optimization problems, where each node is associated with a local cost function, determined from set of local and possibly private measurements . The nodes must nevertheless cooperate in order to minimize the network objective function, which is the sum of local costs. The algorithm design becomes challenging in the absence of a centralized controller or a fusion center, where nodal interactions are limited to their neighborhoods, and global network state information is largely unavailable . In high-dimensional problems, even local message passing may be prohibitive, since each node may only be interested in a subset of the optimization variables.

In general, distributed optimization algorithms are designed either in the primal , [13, Chap. 10] or dual domain . A popular dual approach is the distributed alternating direction method of multipliers (ADMM), where the nodal variables are decoupled through the introduction of the consensus constraints, and the updates are carried out in the dual domain . High-dimensional problems are handled through the so-called general-form consensus formulation, where local updates depend only on a subset of optimization variables[8, Chap. 7]. The ADMM algorithm is also applicable to a class of non-convex problems, where it has been shown to converge to a local minimum .

Practical networks are also heterogeneous with respect to their processing powers, energy availability, and communication capabilities, giving rise to asynchrony . Indeed, a key feature required of the distributed algorithms is their tolerance to processing and communication delays arising due to slow or energy-starved nodes . In general, distributed optimization algorithms such as ADMM must be appropriately modified to allow updates to be skipped or delayed, and the convergence of the asynchronous variants is neither obvious, nor guaranteed [18, 16, Chap. 6]. Asynchronous or randomized variants of the distributed ADMM are well-known for the case when the cost functions are convex; see e.g., and references therein. Likewise, if a master node or a fusion center is available, the asynchronous distributed ADMM variant proposed in , is applicable to non-convex cost functions and has been shown to converge to a local optimum.

This work considers the non-convex general-form consensus optimization problem arising in a multi-agent networked system. The first contribution is the development of an asynchronous and distributed ADMM framework that runs without a fusion center. Two variants, namely, the proximal and the majorized ADMM, are proposed, each handling the non-convex objective function in a different way. While the ADMM iterations at all nodes still occur according to a common schedule, both algorithms allow the nodes to, at times, skip the computationally intensive steps and/or the transmission of updates. As the second contribution, it is shown that both variants converge to a local minimum, under certain regularity conditions. The convergence analysis reveals that with appropriately chosen parameters, the algorithms can tolerate any bounded level of asynchrony. Finally, the third contribution is the application of the proposed ADMM algorithm to the problem of cooperative localization for distributed networks . Detailed simulations and comparisons with existing distributed and asynchronous localization algorithms are carried out, establishing the superior performance of the ADMM algorithm.

This paper is organized as follows. The problem formulation and related examples are provided in Sec. II. The proposed proximal asynchronous ADMM and the associated convergence results are presented in Sec. III, while the proposed majorized asynchronous ADMM is detailed in Sec. IV. Finally, the simulation results for the proposed algorithm applied to the cooperative localization problem are provided in Sec. V and sec. VI concludes the paper.

II Problem Formulation

This section details the partially separable non-convex problem formulation considered here, and motivates the need for a distributed optimization algorithm via two examples. Before describing the problem at hand, some background regarding the general multi-agent optimization problem is first presented.

Consider a network represented by the undirected graph G=(K,E){\mathcal{G}}=({\mathcal{K}},{\mathcal{E}}), where K:={1,2,…,K}{\mathcal{K}}:=\{1,2,\ldots,K\} denotes the set of agents or nodes and E{\mathcal{E}}, the set of edges that represent communication links. A node k∈Kk\in{\mathcal{K}} may only communicate with its neighbors Nk′:={j∣(j,k)∈E}{\mathcal{N}}^{\prime}_{k}:=\{j|(j,k)\in{\mathcal{E}}\}. Consider first the general multi-agent problem where the nodes want to cooperatively solve the following optimization problem:

In general, non-convex problems such as (1) are solved in a distributed manner using the first order gradient or subgradient descent, dual methods such as ADMM, convex relaxation (such as semidefinite relaxation), successive convex optimization, or leveraging the problem structure; see and references therein. These approaches result in algorithms that are parallelizable to various extents, with different computational and message passing requirements. Of particular interest here is the high-dimensional regime, where large NN prohibits nodes from operating over and exchanging the full vector x\mathbf{x}. To this end, the next section considers a partially separable form of (1) which is amenable to a distributed optimization algorithm.

II-B Partially separable form

Consider a special case of (1) where the optimization variables {xn}n=1N\{x_{n}\}_{n=1}^{N} are also partitioned among nodes, and each variable is of interest to exactly one node. To this end, let {Sk}k=1K\{\mathcal{S}_{k}\}_{k=1}^{K} denote disjoint subsets such that the variables {xn∣xn∈Sk}\{x_{n}|x_{n}\in S_{k}\} are local to node kk. Further, the component function gk(⋅)g_{k}(\cdot) at node kk depends only on variables that are local either to node kk or its neighbors. The overall problem considered here takes the following form:

where the set Sk′:=⋃j∈{k}∪Nk′Sj\mathcal{S}_{k}^{\prime}:=\bigcup_{j\in\{k\}\cup{\mathcal{N}}^{\prime}_{k}}\mathcal{S}_{j}. Observe here that for each node kk, the function hkh_{k} and the constraint set Xk{\mathcal{X}}_{k} depend only on the variables in Sk\mathcal{S}_{k}. For this partially separable form, it is now possible to express the objective function as a bipartite factor graph, with check nodes representing the summands gkg_{k}, and the variable nodes representing the sets Sk\mathcal{S}_{k}. For the four-node example network shown in Fig. 1, the factor graph is shown in Fig. 2. From the perspective of algorithm design, the dependence structure imposed by (2) can be exploited to eliminate message passing between non-neighboring nodes. The partially separable form considered here occurs commonly in the context of distributed estimation, where the parameter of interest is a collection of node-specific quantities such as temperatures, node locations, harmful algal blooms, pH, and temperature . Cooperation between nodes is still required here, since parameters at neighboring nodes are often coupled or correlated. The subsequent examples detail two specific applications where the partially separable form problem structure in (2) arises.

II-C Examples

It can be observed that (3) is a special case of (2) with gk(⋅)=ln⁡pk(⋅)g_{k}(\cdot)=\ln p_{k}(\cdot) and hk(⋅)=ln⁡p(yk∣⋅)h_{k}(\cdot)=\ln p({\mathbf{y}}_{k}|\cdot).

where, gk({xj}j∈Nk)=∑j∈Nk′wkj(δkj−dkj(xk,xj))2g_{k}(\{\mathbf{x}_{j}\}_{j\in{\mathcal{N}}_{k}})=\sum_{j\in{\mathcal{N}}^{\prime}_{k}}w_{kj}\left(\delta_{kj}-d_{kj}(\mathbf{x}_{k},\mathbf{x}_{j})\right)^{2}. When the location xkp\mathbf{x}_{k}^{p} of node kk is known a priori, a regularization term of the form ri∥xk−xkp∥r_{i}\left\|\mathbf{x}_{k}-\mathbf{x}^{p}_{k}\right\| may also be added to gk(⋅)g_{k}(\cdot). Since (4) is non-convex, it is often solved via majorization or via SDP relaxation . A distributed and incremental algorithm for solving (4) via majorization was first detailed in . A distributed and asynchronous algorithm using SDP relaxation was described in . The subsequent sections describe the proximal and majorized ADMM algorithms that are also applicable to (4).

III Distributed Asynchronous ADMM

This section details the proposed distributed asynchronous algorithm for solving (2) via ADMM, and provides the relevant convergence results. To begin with, the next subsection describes a distributed synchronous implementation, which is motivated from the so called general-form consensus algorithm for solving convex optimization problems . The synchronous version serves as a starting point for the asynchronous algorithm described in Sec. III-B.

where xk\mathbf{x}_{k} collects the variables required at node kk, i.e., [xk]j=xkj[\mathbf{x}_{k}]_{j}=x_{kj} for j∈Nkj\in{\mathcal{N}}_{k}, and zero otherwise.

The idea of introducing consensus variables, in order to make the updates separable in the optimization variables, is well known [10, Chap. 5]. It is now possible to apply the ADMM method by associating dual variables ykjy_{kj} for each constraint in (6) and writing the augmented Lagrangian as

From (III-A), it is clear that LL is separable in {xk}\{\mathbf{x}_{k}\}. The Lagrangian is also separable in {zj}\{z_{j}\}, since it holds that

Together, (III-A) and (III-A) allow the ADMM updates to be carried out in a distributed fashion. In particular, starting with arbitrary {xk1}\{\mathbf{x}_{k}^{1}\} and {ykj1}\{{\mathbf{y}}^{1}_{kj}\}, the update for {zjt+1}\{z_{j}^{t+1}\} are evaluated as

Similarly updates for xkt+1\mathbf{x}_{k}^{t+1} can be obtained by minimizing (III-A) with respect to xk\mathbf{x}_{k}, which yields

where the optimization is with respect to {xkj}j∈Nk\{x_{kj}\}_{j\in{\mathcal{N}}_{k}}. Observe that since the component functions gk(⋅)g_{k}(\cdot) are not necessarily convex, the update in (11) is difficult to carry out. As suggested in , the update xkt+1\mathbf{x}_{k}^{t+1} can however be calculated approximately as follows

where the vector [zk]j:=zj[{\mathbf{z}}_{k}]_{j}:=z_{j} for all j∈Nkj\in{\mathcal{N}}_{k} and zero otherwise. Since nodal functions gk(⋅)g_{k}(\cdot) depend only on {xn}n∈Nk\{x_{n}\}_{n\in{\mathcal{N}}_{k}}, the gradient vector is defined as

The approximate update of xkjt+1x_{kj}^{t+1} thus becomes

Algorithm 1 summarizes the implementation of the distributed ADMM described here. The main feature of Algorithm 1 is that it does not require a master node, and all message passing is limited to the neighboring nodes only. The distributed implementation also requires that each node must be capable of carrying out the updates in (9) and (III-A). Algorithm 1 may be viewed as the proximal variant of the distributed ADMM algorithm , applied to the non-convex problem (2). The stopping criterion in Algorithm 1 is simply ∣zkt+1−zkt∣≤δ\lvert z_{k}^{t+1}-z_{k}^{t}\rvert\leq\delta for all k∈Kk\in{\mathcal{K}} and a small δ>0\delta>0.

Algorithm 1 is a synchronous protocol since all updates must necessarily be carried out at every iteration by every node. Its applicability to heterogeneous networks is therefore limited, since the progress of the algorithm is determined by the slowest node in the network. For instance, the following issues may arise when Algorithm 1 is implemented on a wireless network with energy-constrained, low-cost devices.

Each node is required to transmit two messages to each of its neighbors per-iteration. This might be excessive for nodes operating on a power-budget.

Delays in function computation or communication may also arise because of the heterogeneity of nodes in a wireless network. For instance, mission critical nodes may attempt to extend their battery lives by operating under a power-saving mode, thus deliberately reducing their computational capabilities. On the other hand, the available power at some energy-harvesting nodes may vary throughput the day, depending on, for instance, the received solar energy.

Sec. III-B describes an asynchronous version of Algorithm 1 that overcomes (S1)-(S2) by allowing nodes to carry out the resource-intensive operations in Steps 5 and/or 7 intermittently. Each time an update is skipped, the node saves on both computational and communication costs. It is shown that the asynchronous algorithm converges, as long as the updates are performed “often enough,” a notion that is made precise through some constraints on the algorithm parameters.

III-B Distributed Asynchronous Algorithm with Optional and Delayed Updates

Whenever j∉Stj\notin\mathcal{S}^{t}, the subsequent transmission may not be carried out either, and the non-updating node may simply stay silent. The neighboring nodes will then wait for a fixed amount of time to receive an update, and assume that zjt+1=zjtz_{j}^{t+1}=z_{j}^{t} holds for all nodes j∈Nkj\in{\mathcal{N}}_{k} that do not transmit anything.

The transmission of updates in Step 3 is again optional and in the event that no update is received at a neighbor j∈Nkj\in{\mathcal{N}}_{k}, (16) is used instead.

The proposed asynchronous algorithm for node k∈Kk\in{\mathcal{K}} is summarized in Algorithm 2. The salient features of the proposed asynchronous algorithm are as follows.

All resource-intensive steps, such as gradient calculation, proximal function calculation, and transmission of updates are now optional.

Nodes operating on a power budget may only carry out the updates in Steps 11 and 12 at every iteration. When using the old gradient, these updates amount to simple addition/subtraction operations. It is also possible to defer Steps 11 and 12 to a later time slot when the transmission in Step 3 occurs.

The nodes may implement a timeout mechanism when listening for updates [cf. Steps 3 and 5]. The appropriate default actions must be triggered if nothing is heard from one or more neighbors.

Fig. 3 shows an illustration where nodes 2 and 3 are neighbors of a node 1, that go to sleep at iteration t=2t=2. During the sleep state, a node is only able to receive updates, but not transmit or carry out any computations. Observe that in the second sub-slot of the 2-nd iteration, z23z^{3}_{2} and z33z^{3}_{3} are also not updated at nodes 22 and 33 since no update for x122x_{12}^{2} and y122y_{12}^{2} was received. Consequently, it will hold that zkt=zkt−1z_{k}^{t}=z_{k}^{t-1} for k=1k=1, 22, and 33. More generally, each sleeping node will force all its neighbors to use the update in Step 8 of Algorithm 2. Since the update frequency at each node cannot be too small, it is required that several nodes should be awake at any given time.

The next section provides the convergence analysis of Algorithm 2, ensuring that for any δ>0\delta>0, the stopping criteria in Step 13 of Algorithm 2 is eventually met. Note that the convergence results also apply to Algorithm 1, which is simply a special case of Algorithm 2. Before concluding this subsection, the following remark about related work in the context of asynchronous algorithms is due.

The key feature of the proposed algorithms is that some of the updates and transmissions are entirely optional. From the perspective of algorithm design, this results in two sources of asynchrony, namely, delayed gradients, and skipped updates. The asynchronous algorithm in can handle bounded delays in the gradient calculation, but does not consider skipping updates. This is because the DA-ADMM algorithm in utilizes a fusion center that is not energy constrained. Note that the absence of such a fusion center significantly complicates the algorithm design, also restricting the application of Algorithms 2 to (2), which is a special case of (1).

An asynchronous ADMM algorithm was also proposed and applied to non-convex problem in . However, no proof of convergence was provided, and the performance of the algorithm was only tested via simulations. A number of asynchronous variants exist for convex problems. The optional-update idea used in (16) is in fact inspired from the asynchronous ADMM proposed in . Different from Algorithm 2 however, the asynchronous algorithm in also uses random network states, and is applicable only to convex problems. For partially separable problems such as (2), it may also be possible to develop an asynchronous variant of the coordinate descent algorithm. The stochastic asynchronous coordinate descent algorithm proposed in also allows erroneous gradients but is applicable only to convex problems. Finally, asymptotic results on the effect of asynchrony on the stochastic gradient descent algorithm were provided in .

III-C Convergence Analysis for Algorithm 2

In order to establish the convergence of the asynchronous algorithm, some assumptions regarding the problem structure and algorithm parameters must be made. In particular, it is shown that the algorithm converges as long as the optional updates happen “often enough.” Specifically, recall that t+1−[t+1]≤Tkt+1-[t+1]\leq T_{k}, which implies that the gradients may be calculated using updates that are at most TkT_{k}-old. Similarly, let the frequency of update in (9) being carried out at node kk be denoted by 0<fk≤10<f_{k}\leq 1. For instance, if the update occurs once every KK time slots, fk=1/Kf_{k}=1/K. Then the following assumptions are required.

The set X\mathcal{X} is a closed, convex, and compact. The functions gk(x)g_{k}(x) is bounded from below over X\mathcal{X}.

For node kk, the step size ρk\rho_{k} is chosen large enough such that, it holds that αk>0\alpha_{k}>0 and βk>0\beta_{k}>0, where

Of these, Assumptions (A1) and (A2) are standard in the context of non-convex optimization and are satisfied for most problems of interest. Intuitively, the use of old gradients is permitted only because they are Lipschitz continuous, and therefore change slowly over iterations. The boundedness assumption (A2) is required to ensure that the updates {xt,zt,yt}\{\mathbf{x}^{t},\mathbf{z}^{t},{\mathbf{y}}^{t}\} generated by Algorithm 2 stay bounded. This assumption may be dropped for certain problems where boundedness of these iterates may arise naturally. Finally, Assumption (A3) specifies the exact relationship that the algorithm and problem parameters {Lk,ρk,fk,Tk}k=1K\{L_{k},\rho_{k},f_{k},T_{k}\}_{k=1}^{K} must satisfy in the worst case. Interestingly, the choice of ρk\rho_{k} is node-specific since nodes may differ with respect to their computational capabilities resulting in different levels of asynchrony. While not stated explicitly, the convergence also requires that ρk<∞\rho_{k}<\infty. Thus, for (19) to hold, it is necessary that fk>0f_{k}>0 and Tk<∞T_{k}<\infty. In other words, the delay cannot be unbounded in the worst case.

The convergence of the asynchronous algorithm is established through the following intermediate lemma that holds under Assumptions (A1)-(A3).

Starting from any time t=t0t=t_{0}, there exists T<∞T<\infty such that

The augmented Lagrangian values in (III-A) are bounded from below, i.e., for any time t≥1t\geq 1, it holds that Lagrangian satisfies

The proof of Lemma 1 is provided in Appendix A. Lemma 1(a) establishes that there exists some finite TT such that the augmented Lagrangian values are non-increasing after TT iterations. In practice, the value of TT depends on the update frequencies {fk}k=1K\{f_{k}\}_{k=1}^{K}. For instance, TT could be the minimum number of iterations in which each node kk updates zktz_{k}^{t} at least fkTf_{k}T times. Lemma 1(b) establishes that the Lagrangian is bounded from below. One way to interpret Lemma 1 is to define the sequence Ln:=L({xknT+1};znT+1,{yknT+1})L^{n}:=L(\{\mathbf{x}_{k}^{nT+1}\};{\mathbf{z}}^{nT+1},\{{\mathbf{y}}_{k}^{nT+1}\}) for all n≥0n\geq 0, and observe that {Ln}\{L^{n}\} is non-increasing and bounded from below, and therefore convergent. The subsequent theorem establishes the final convergence result and related properties.

The iterates generated by Algorithm 2 converges in the following sense

For each k∈Kk\in{\mathcal{K}} and j∈Nkj\in{\mathcal{N}}_{k}, denote limit points of the sequences {zkt}\{z_{k}^{t}\}, {xkjt}\{x_{kj}^{t}\}, and {ykjt}\{y_{kj}^{t}\} by zk⋆z_{k}^{\star}, xkj⋆x_{kj}^{\star}, and ykj⋆y_{kj}^{\star}, respectively. Then {{zk⋆},{xkj⋆},{ykj⋆}}\{\{z_{k}^{\star}\},\{x_{kj}^{\star}\},\{y_{kj}^{\star}\}\} is a stationary point of (5) and satisfies

The proof of Theorem 1 is provided in Appendix B. Note that it suffices to show that Algorithm 2 converges to a stationary solution of (5) which is equivalent to (2). In other words, {zn⋆}n=1N\{z_{n}^{\star}\}_{n=1}^{N} can be used as a solution to (2). It is emphasized that Algorithm 2 may not necessarily converge to a globally optimum solution to (2).

Further, using assumptions (A1)-(A3) alone, it is also difficult to quantify the convergence rate of the ADMM algorithm for the non-convex case. As will be shown in Sec. V however, rate of convergence for the localization problem is low whenever ρk\rho_{k} is large. Assuming that this observation applies to (2) generally, it makes sense to choose ρk\rho_{k} as small as possible, while respecting (19). Interestingly, this result also matches with the intuition that the Algorithm 2 converges slowly if the asynchrony is high, i.e., when Tk≫1T_{k}\gg 1 and fk≪1f_{k}\ll 1, since it would require choosing a larger ρk\rho_{k} for each k∈Kk\in{\mathcal{K}}.

IV Majorized Asynchronous ADMM

In many problems, it is possible to upper bound the non-convex component functions gk(xk)g_{k}(\mathbf{x}_{k}) with an appropriate convex surrogate function. Given zk∈X\mathbf{z}_{k}\in{\mathcal{X}}, the surrogate fk(xk,zk)f_{k}(\mathbf{x}_{k},\mathbf{z}_{k}), also referred to as the majorizing function, is such that for all xk\mathbf{x}_{k}, gk(xk)≤fk(xk,zk)g_{k}(\mathbf{x}_{k})\leq f_{k}(\mathbf{x}_{k},\mathbf{z}_{k}), and satisfies

When such a majorizer exists and can be found easily, the update of xkt+1\mathbf{x}_{k}^{t+1} in (14) can be carried out more accurately, without appreciable increase in the per-iteration complexity.

This section develops a provably convergent variant of Algorithm 2 that utilizes majorization while updating xkt+1\mathbf{x}_{k}^{t+1}. Specifically, while (15) and (16) stay the same, the update for xkt+1\mathbf{x}_{k}^{t+1} becomes

where, as in (17), zk[t+1]\mathbf{z}_{k}^{[t+1]} is used for calculating the surrogate function. Interestingly, the majorized ADMM proposed here enjoys the flexibility afforded by its proximal counterpart, and can be implemented as in Algorithms 2 or LABEL:async_algowsn.

In order to show that the majorized ADMM converges, the following assumptions are required in addition to (A2).

For each node kk, there exists a constant Lk≥0L_{k}\geq 0, such that for all x,x′,z,z′\mathbf{x},\mathbf{x}^{\prime},\mathbf{z},\mathbf{z}^{\prime}, the following inequalities are satisfied:

For node kk, the step size ρk\rho_{k} is chosen large enough such that αk>0\alpha_{k}>0 and βk>0\beta_{k}>0, where

Note that Assumption (A4) utilizes the same Lipschitz constant in (26) for the sake of simplicity. In general, if there exist constants LkgL^{g}_{k}, LkxL^{x}_{k}, and LkzL^{z}_{k} in (26), it is always possible to define Lk:=max⁡{Lkg,Lkx,Lkz}L_{k}:=\max\{L^{g}_{k},L^{x}_{k},L^{z}_{k}\} for all kk. Likewise, the convergence proof provided here can also be developed for the case when the three constants are different, resulting in slightly tighter bounds.

The conditions required in (A4) are not very restrictive. Consider for instance a component function gk(x)g_{k}(\mathbf{x}) that is expressible as a sum of a convex function gk1(x)g^{1}_{k}(\mathbf{x}) and a concave function gk2(x)g^{2}_{k}(\mathbf{x}). The concavity of gk2(x)g^{2}_{k}(\mathbf{x}) allows it to be majorized by its supporting hyperplane, i.e., given any z∈dom gk2(x)\mathbf{z}\in\textrm{dom}~{}g^{2}_{k}(\mathbf{x}),

where it can be verified that fk(x,z)f_{k}(\mathbf{x},\mathbf{z}) satisfies (23) and (24). Then, observe that ∇fk(x,z)=∇gk1(x)+∇gk2(x)\nabla f_{k}(\mathbf{x},\mathbf{z})=\nabla g_{k}^{1}(\mathbf{x})+\nabla g_{k}^{2}(\mathbf{x}) for all x\mathbf{x} and z\mathbf{z}. If Lk1L^{1}_{k} and Lk2L^{2}_{k} are Lipschitz constants of ∇gk1\nabla g^{1}_{k} and ∇gk2\nabla g^{2}_{k} respectively, it may be seen that (26a) and (26b) hold with Lipschitz constant Lk1+Lk2L^{1}_{k}+L^{2}_{k}, while the left-hand side of (26c) is identically zero. The following Theorem, whose proof is provided in Appendix C, summarizes the main result of this section.

The iterates generated by the majorized ADMM ((16), (IV), and (15)) converge to a stationary point of (5).

V Distributed Localization in Networks

This section builds upon the localization examples introduced in Sec. II-C and provides simulation results for comparing two different related algorithms. Monte-Carlo simulations are performed over randomly generated networks with N=25N=25 nodes. To this end, nodes are uniformly distributed over a unit two-dimensional area R=2\mathcal{R}=^{2}. Of these, m=5m=5 nodes are located at (0.25,0.25)(0.25,0.25), (0.75,0.25)(0.75,0.25), (0.25.0.75)(0.25.0.75), (0.5,0.5)(0.5,0.5), at (0.75,0.75)(0.75,0.75), and serve as anchor nodes for others. Distance measurement and communication between nodes kk and jj is possible only if they are within a distance of R=0.5R=0.5 of each other. For neighboring nodes, the weights wkjw_{kj} are all set to unity.

It is remarked that the classical SMACOF algorithm used for solving in is not distributed, and is therefore not applicable in the present context. The distributed weighted MDS (DwMDS) framework introduced in is an incremental algorithm that utilizes component wise majorization in order to circumvent the non-convexity of the objective function. Different from the proposed algorithms, incremental algorithms utilize a cyclic message passing routine and are synchronous in nature.

In the present case, the following modified version of the objective function in (4) is considered dkj(xk,xj)=∥xk−xj∥+ϵd_{kj}(\mathbf{x}_{k},\mathbf{x}_{j})=\sqrt{\left\|\mathbf{x}_{k}-\mathbf{x}_{j}\right\|+\epsilon}, where ϵ>0\epsilon>0 is a small number introduced to make the objective function differentiable everywhere. This modification also makes ∇gk\nabla g_{k} Lipschitz continuous, as required in (A1), and verified in Appendix D. Recall that Algorithm 2 allows two modes of asynchrony, namely, old gradient and skipped updates. In the implementation described here, it is assumed that at each time instant, some nodes go to sleep and do not carry out any update. Further, at each iteration, each node randomly chooses an older available gradient or calculates a more recent one, while ensuring that the gradient used is at most TkT_{k}-old.

Fig. 4 compares the NRMSE performance of DwMDS algorithm , the synchronous (SyncE-ML) and asynchronous (AsyncE-ML) versions of the E-ML algorithm , and the Algorithms 1 (SyncDADMM), 2 (AsyncDADMM) . For the asynchronous E-ML algorithm, 125 of the 130 edges are activated per iteration. For Algorithm 2, we set Tk=8T_{k}=8 and fk=0.75f_{k}=0.75, while ρk\rho_{k} is chosen so as to satisfy (A3).

From Fig. 4, it is clear that both, Algorithms 1, 2 outperform all the other state-of-the-art algorithms. As expected, the asynchronous versions of the E-ML and ADMM algorithms have poorer performance than their synchronous counterparts. Interestingly, asynchrony has a greater effect on the asymptotic performance of the E-ML algorithm than the convergence rate of the ADMM algorithm. It is also worthwhile to compare the implementation complexity of the two asynchronous algorithms. Since the E-ML algorithm uses SDR, its computational complexity is approximately O(∣Nk∣3)\mathcal{O}(\lvert{\mathcal{N}}_{k}\rvert^{3}) for node kk per iteration. This is because each iteration requires solving an ∣Nk∣\lvert{\mathcal{N}}_{k}\rvert-sized convex semidefinite program. On the other hand, gradient calculation is the most computationally demanding step in Algorithms 1 and 2, which translates to a total computational complexity of O(∣Nk∣)\mathcal{O}(\lvert{\mathcal{N}}_{k}\rvert) for node kk per iteration. Note further that unlike Algorithm 2, the updates in AsyncE-ML algorithm are also not optional, irrespective of the number of edges activated per iteration. The communication complexity for the two algorithms is of the same order, i.e., O(∣Nk∣)\mathcal{O}(\lvert{\mathcal{N}}_{k}\rvert) for node kk per iteration. However, the message transmitted to each neighbor in AsyncE-ML consists of nine real numbers, as opposed to four real numbers in Algorithm 2.

V-B Convergence rates

As stated earlier, while Theorem 1 establishes convergence, it does not provide any indication regarding the convergence rate exhibited by Algorithm 2. This subsection compares the convergence rates of the synchronous proximal ADMM (SyncADMM), syncrhonous majorized ADMM (SyncMDADMM), and the consensus-based distributed gradient descent (C-DGD) method proposed in [13, Chap. 10]. While C-DGD has only been proposed for convex problems, it has been implemented in a distributed and asynchronous manner , and is therefore an interesting candidate for the non-convex localization problem. The goal is to compare the three algorithms in terms of how fast they approach a local minimum, while ignoring their distances to the global minimum. To this end, Fig. 5 shows the evolution of the convergence criterion ϕ(t):=∥1N∑i=1N(xit+1−xit)∥F\phi(t):=\left\|\frac{1}{N}\sum_{i=1}^{N}({\mathbf{x}_{i}^{t+1}-\mathbf{x}_{i}^{t}})\right\|_{F} with the iteration index tt. It can be observed that the majorized ADMM performs slightly better than the proximal ADMM. On the other hand, the C-DGD algorithm convergence relatively slowly. In particular, unlike the other simulations, the plots shown in Fig. 5 are generated for a network with communication range R=0.8R=0.8, since the C-DGD algorithm does not converge for R=0.5R=0.5. It is remarked that as the range decreases, fewer measurements are available, making network localization more difficult.

At this point, it may be useful to also compare the complexity and run-times of the three algorithms. In general, the majorized ADMM is the most complex, since it involves solving a convex optimization problem at every iteration. On the other hand, both the C-DGD and the proximal ADMM algorithms require simple calculations with the gradients of gkg_{k}. For the localization problem considered here, the per-iteration majorization update at node kk requires solving a system of ∣Nk∣×∣Nk∣)|\mathcal{N}_{k}|\times|\mathcal{N}_{k}|) linear equations, incurring almost four times as much CPU-time as the proximal ADMM and the C-DGD algorithms. This extra complexity more than compensates for the reduced number of iterations afforded by the majorized ADMM. In summary, the proximal ADMM algorithm takes the least wall-time to reach a certain accuracy, say ϕ(t)=10−15\phi(t)=10^{-15}.

V-C Choice of parameters

In order to further study the convergence rates, performance of Algorithm 2 is studied for different parameter values. The convergence rate is analyzed by plotting the stopping critereon given by ψ(t):=∥zt+1−zt∥\psi(t):=\left\|\mathbf{z}^{t+1}-\mathbf{z}^{t}\right\| against the iteration index tt.

Theorem 1 guarantees the convergence of Algorithm 2 whenever ρk\rho_{k} is chosen in accordance with (A3). In practice however, it may be possible to improve the convergence rate of Algorithm 2 by choosing smaller values of ρk\rho_{k}. Fig. 6 shows the evolution of ψ(t)\psi(t) and NRMSE with iterations, for different values of ρk\rho_{k}, while Tk=4T_{k}=4 and fk=0.75f_{k}=0.75 are kept fixed. As expected, the algorithm takes longer to converge for larger values of ρk\rho_{k}, although the asymptotic NRMSE performance for all cases remains the same.

V-C2 Effect of Asynchrony

Fig. 7 shows the effect of choosing the parameters TkT_{k} and fkf_{k} on the rate of convergence. To this end, we set ρk=10\rho_{k}=10 for all k∈Nk\in{\mathcal{N}}, and vary TkT_{k} and fkf_{k} separately. Both figures confirm the intuition that introduction of asynchrony results in slower convergence, even when ρk\rho_{k} stays the same. Note however that in order to guarantee convergence, it is necessary to choose increasingly larger values of ρk\rho_{k} for larger values of TkT_{k} or smaller values of fkf_{k}. Such a choice may therefore result in even slower convergence but allow higher asynchrony.

VI Conclusion

This paper develops an asynchronous distributed ADMM algorithm that is applicable to a class of non-convex optimization problems. The non-convexity of the cost functions is handled either by making a first order approximation or via majorization, resulting in two variants of the proposed algorithm. Both variants converge to a stationary point of the optimization problem as long as the ADMM updates are applied “often enough.” The proposed algorithms find applications in distributed in-network estimation and localization. Comparisons with state-of-the-art distributed algorithms for the problem of cooperative localization in ad hoc networks demonstrates the superior performance of the proposed algorithm.

Appendix A Proof of Lemma 1

This appendix provides the convergence analysis for the Lagrangian in Algorithm 2. Before proceeding with the proof, some notation is introduced. Recall the definitions of xk\mathbf{x}_{k}, yk{\mathbf{y}}_{k}, and zk\mathbf{z}_{k}, and similarly define [xˇj]k:=xkj[\check{\mathbf{x}}_{j}]_{k}:=x_{kj} and [yˇj]k:=xkj[\check{{\mathbf{y}}}_{j}]_{k}:=x_{kj} for all k∈Njk\in{\mathcal{N}}_{j}. Since the Lagrangian is separable in both {xk}\{\mathbf{x}_{k}\} and {zj}\{z_{j}\}, the following notation is introduced for the summands in (III-A) and (III-A):

With this definition, observe that the update for xkt+1\mathbf{x}_{k}^{t+1} in (III-A) is given by

where the gradient with respect to xk\mathbf{x}_{k} is defined similar to that in (13). In order to show that the Lagrangian decreases over several iterations, express the difference between consecutive Lagrangian values as

The subsequent lemma establishes bounds on the different terms in (A), and will be utilized to prove Lemma 1(a).

Define St\mathcal{S}^{t} as the set of nodes for which the update in (16) is carried out at time tt. Then it holds that

The proof begins by establishing a bound on difference between successive dual values. Rearranging (17),

which, together with the update for ykjt+1y_{kj}^{t+1} [cf. (15) ] yields

Similarly, it holds that ykjt=−[∇gk(zk[t])]jy_{kj}^{t}=-[\nabla g_{k}({\mathbf{z}}_{k}^{[t]})]_{j} for all j∈Nkj\in{\mathcal{N}}_{k}. Therefore, for each k∈Kk\in{\mathcal{K}}, the following bound applies:

where (38a) follows from Assumption (A1) while (38b) and (38c) follow from the use of triangle inequality. Next, the different bounds in (34), (35), and (36) are proved.

(a) The bound in (34) follows from the following equalities utilizing (15).

Finally, the bound in (34) follows by substituting (38c) into the right-hand side of (39).

where (43b) follows from (32), (43c) follows from (A1), and the rest follow from the use of triangle inequality similar to that in (38).

Similar to (40), the following quadratic upper bound is also implied by (A1),

where (46a) follows from (A1) and (46b) from the use of triangle inequality.

Having derived inequalities (A)-(46), the difference between the summands of (35) can finally be bounded. Towards this end, it holds that

where (48) follows from using the update for ykjt+1y_{kj}^{t+1} and (49) from (38). Finally, summing (49) over k=1,2,…,Kk=1,2,\ldots,K, the result in (35) follows.

(c) Define the indicator function 11X(x)=111_{{\mathcal{X}}}(x)=1 if x∈Xx\in{\mathcal{X}} and zero otherwise, the observe that update for zjz_{j} can be written as [cf. (9)],

for all j∈Stj\in\mathcal{S}^{t}. The first order optimality condition for (51) is that

Since 11Xj(zjt+1)=11Xj(zjt)=011_{{\mathcal{X}}_{j}}(z_{j}^{t+1})=11_{{\mathcal{X}}_{j}}(z_{j}^{t})=0 for j∈Stj\in\mathcal{S}^{t}, from (52) and (53) it holds that

For all j∉Stj\notin\mathcal{S}^{t}, since zjt+1=zjtz_{j}^{t+1}=z_{j}^{t}, the right-hand side is clearly zero. Finally, the result in (36) is obtained by summing both sides in (54) over j=1,2,…,Kj=1,2,\ldots,K. ∎

From Lemma 2, the decrease in the Lagrangian over consecutive time slots is given by

Given t0≥1t_{0}\geq 1 and T≥1T\geq 1, define the inverse mapping ST−1(k):={t0≤t≤t0+T∣k∈St}\mathcal{S}^{-1}_{T}(k):=\{t_{0}\leq t\leq t_{0}+T\mid k\in\mathcal{S}^{t}\} as the set of iterations for which update (9) is applied at node kk and let fk≥∣ST−1(k)∣/Tf_{k}\geq\lvert\mathcal{S}^{-1}_{T}(k)\rvert/T. Then, summing both sides of (A) over t=t0t=t_{0}, t0+1t_{0}+1, …\ldots, T+t0T+t_{0},

where αk\alpha_{k} and βk\beta_{k} are defined as in (19), yielding the desired result. Note that from Assumption (A4), the term on the right-hand side of (A) is negative.

: Using the Lipschitz continuity of ∇gk(⋅)\nabla g_{k}(\cdot) and applying triangle inequality, it follows that

Next, using the relationship from (37), observe that the

Next, using Cauchy-Schwarz inequality on the last term, it follows that

where the last inequality follows from the fact that ρk≥7Lk\rho_{k}\geq 7L_{k} [cf. (A4)], and that Xj{\mathcal{X}}_{j} is compact [cf. (A2)]. ∎

Appendix B Proof of Theorem 1

From Lemma 1, it follows that L({xkt},zt;{ykt})L(\{\mathbf{x}_{k}^{t}\},{\mathbf{z}}^{t};\{{\mathbf{y}}_{k}^{t}\}) converges as t→∞t\rightarrow\infty. Therefore, it holds from ((a)) that,

From (65), it follows that the limit points {{xkj⋆},{zj⋆},{ykj⋆}}\{\{x_{kj}^{\star}\},\{z_{j}^{\star}\},\{y_{kj}^{\star}\}\} exist, and satisfy

which is the primal feasibility condition for (5). Since zjt+1∈Xjz_{j}^{t+1}\in{\mathcal{X}}_{j} for all tt, it should hold that zj⋆∈Xjz_{j}^{\star}\in{\mathcal{X}}_{j} and xkj⋆=xlj⋆x_{kj}^{\star}=x_{lj}^{\star} for all j,l,k∈Nkj,l,k\in{\mathcal{N}}_{k}. The first order optimality condition for xk\mathbf{x}_{k} can be obtained from (37), which implies (22a). Finally, the first order optimality condition for zjz_{j} can be obtained from (52), which implies that

The desired result follows by using (66) and noting that 0∈∂11z∈Xj∣z=zj⋆0\in\partial 11_{z\in{\mathcal{X}}_{j}}\mid_{z=z_{j}^{\star}}. ∎

Appendix C Convergence of the Majorized ADMM

From the update in (IV), it holds that ∇xkuk(xkt+1,zk[t+1],zkt+1,ykt)=0\nabla_{\mathbf{x}_{k}}u_{k}(\mathbf{x}_{k}^{t+1},\mathbf{z}_{k}^{[t+1]},\mathbf{z}_{k}^{t+1},{\mathbf{y}}_{k}^{t})=0, and

Further, from the update of ykjt+1y_{kj}^{t+1} and from (68),

Similarly, it holds that ykjt=−[∇xkfk(xkt,zk[t])]jy_{kj}^{t}=-[\nabla_{\mathbf{x}_{k}}f_{k}(\mathbf{x}_{k}^{t},\mathbf{z}_{k}^{[t]})]_{j} for all j∈Nkj\in{\mathcal{N}}_{k}. Therefore, for each k∈Kk\in{\mathcal{K}}, the following bound applies:

where (70a) follows from Assumption, and (70b) follows similarly as in (38c).

As in Appendix A, the Lagriangian is split into three summands [cf. Lemma 2], each of which must be separately bounded. The bound on the first summand follows directly from (70).

Utilizing (68) and following steps (43a)-(43d), (74b) is obtained, which follow similarly as in (43).

The Lipschitz continuity of ∇gk(⋅)\nabla g_{k}(\cdot) and fk(⋅,zkt+1)f_{k}(\cdot,\mathbf{z}_{k}^{t+1}) also imply the following upper bounds

where (78c) follows from (A4) and (78d) from the use of triangle inequality. Finally, it holds that

where (81) follows from using the update for ykjt+1y_{kj}^{t+1} and (82) from (70). Finally, summing (82) over k=1,2,…,Kk=1,2,\ldots,K, the required bound is obtained. Since the expression for the third summand in (36) remains the same, the decrease in the Lagrangian over consecutive time slots t=t0+1,…,t0+Tt=t_{0}+1,\ldots,t_{0}+T, is given by

where, αk\alpha_{k} and βk\beta_{k} are given in (27) yielding the desired result. Note that from Assumption (A5), the term on the right-hand side of (A) is negative. Next, the boundedness of the Lagrangian follows as shown in Appendix A. Finally, it is possible to apply Theorem 1 to this case, yielding the desired result.

Appendix D Lipschitz continuity of objective function in (4)

Recall that gk({xj}j∈Nk)=∑jwkj(δkj−dkj(xk,xj))2g_{k}(\{\mathbf{x}_{j}\}_{j\in{\mathcal{N}}_{k}})=\sum_{j}w_{kj}(\delta_{kj}-d_{kj}(\mathbf{x}_{k},\mathbf{x}_{j}))^{2}, where the modified definition of dkj(xk,xj)=∥xk−xj∥+ϵd_{kj}(\mathbf{x}_{k},\mathbf{x}_{j})=\sqrt{\left\|\mathbf{x}_{k}-\mathbf{x}_{j}\right\|+\epsilon} is utilized and for a pair of node (k,j)(k,j) weight wkj=1w_{kj}=1 for j∈Nkj\in{\mathcal{N}}_{k} otherwise 0. The gradient of gkg_{k} is given by ∇xlgk({xj}j∈Nk)=\nabla_{\mathbf{x}_{l}}g_{k}(\{x_{j}\}_{j\in{\mathcal{N}}_{k}})=

For the present localization example, we have max⁡k,j{δkj}≤1\max_{k,j}\{\delta_{kj}\}\leq 1, max⁡k∣Nk∣≤N\max_{k}|{\mathcal{N}}_{k}|\leq N, wkj≤1w_{kj}\leq 1, ∥xk∥≤1∀k=1,…,K\|\mathbf{x}_{k}\|\leq 1\quad\forall\quad k=1,\ldots,K, and ∥1dkj(xk,xj)∥≤1ϵ\|\frac{1}{d_{kj}(\mathbf{x}_{k},\mathbf{x}_{j})}\|\leq\frac{1}{\sqrt{\epsilon}}. The application of triangle inequality therefore implies that ∥C(xk,xm)∥≤8(2/ϵ+1)\left\|C(\mathbf{x}_{k},\mathbf{x}_{m})\right\|\leq 8(2/\sqrt{\epsilon}+1) for all 1≤k,m≤N1\leq k,m\leq N. Therefore it follows that ∇xl,xm2gk({xj}j∈Nk)\nabla^{2}_{\mathbf{x}_{l},\mathbf{x}_{m}}g_{k}(\{x_{j}\}_{j\in{\mathcal{N}}_{k}}) is bounded, as claimed.

References