---
abstract: 'We propose a distributed algorithm, named D-ADMM, for solving separable optimization problems in networks of interconnected nodes or agents. In a separable optimization problem, the cost function is the sum of all the agents’ private cost functions, and the constraint set is the intersection of all the agents’ private constraint sets. We require the private cost function and constraint set of a node to be known by that node only, during and before the execution of the algorithm. The application of our algorithm is illustrated with problems from signal processing and control, namely average consensus, compressed sensing, and support vector machines. It is well known that communicating in distributed environments is the most energy/time-demanding operation. Thus, algorithms using less communications are more prone to make networks live longer, e.g., sensor networks, or to execute faster, e.g., in supercomputing platforms. Through simulations for several network types and problems, we show that our algorithm requires less communications than the state-of-the-art algorithms.'
author:
- 'João F. C. Mota, João M. F. Xavier, Pedro M. Q. Aguiar, and Markus Püschel [^1]'
bibliography:
- 'paper.bib'
title: 'D-ADMM: A Communication-Efficient Distributed Algorithm For Separable Optimization'
---

Distributed algorithms, alternating direction method of multipliers, consensus, compressed sensing, machine learning, sensor networks.

Introduction {#Sec:intro}
============

Recently, there has been a growing interest in distributed methods for data processing. Such interest was triggered both by the need of processing large quantities of data in distributed platforms, such as supercomputers, and by applications where data comes from spatially different locations and central processing is costly or impractical, as in sensor networks. On the other hand, over the last decades, convex optimization theory has proven to be an invaluable tool for deriving data processing algorithms [@Boyd:ConvexOpti]. Standard algorithms like interior-point methods are, however, unable to handle either large-scale or distributed problems. To solve these problems we thus need new algorithms.

In this paper, we aim to solve separable optimization problems in a distributed, decentralized way. In a separable optimization problem, the cost function can be decomposed as a sum of $P$ functions $f_p$ and the constraint set can be decomposed as the intersection of $P$ sets $X_p$: $$\label{Eq:IntroSeparableProb}
        \begin{array}{ll}
            \underset{x}{\text{minimize}} & f_1(x) + f_2(x) + \cdots + f_P(x) \\
            \text{subject to} & x \in X_1 \cap X_2 \cap \cdots \cap X_P\,,
        \end{array}$$ where $x \in \mathbb{R}^n$ is the global optimization variable. We will denote any solution of  with $x^\star$. We associate with problem  a network of $P$ nodes, where only node $p$ has access to its private cost function $f_p$ and to its private constraint set $X_p$. This situation is illustrated in Figure \[Fig:IntroFig7Nodes\]. Each node can communicate with its neighbors using the network infrastructure. A method that does not use any kind of special or central node will be called *distributed algorithm*. Thus, in a distributed algorithm, aggregating data is not allowed; instead, each node has to exchange information only with its neighbors in order to achieve the common goal: arriving at a solution of . We propose a novel distributed algorithm based on the alternating direction method of multipliers (ADMM) for solving .

In applications such as sensor networks, nodes are battery-operated devices and thus energy is a scarce resource. It is well known that communication is the most energy-consuming operation in such a device [@Akyildiz02-WirelessNetworksASurvey]; therefore, algorithms using less communications can increase the lifespan of the network. Similarly, communication is known to be the algorithmic bottleneck in applications with more controlled environments, for example, in network protocols [@Kurose05-ComputerNetworking] or supercomputing platforms [@Graham05-FutureOfSuperComputing]. Algorithms using few communications are thus required for such applications.

We use our proposed distributed algorithm to solve several instances of  from several areas in engineering, namely average consensus, compressed sensing problems, and support vector machines. Our simulations show that, for each of these problems and for almost all network types, our proposed algorithm requires significantly less communication than prior state-of-the-art algorithms, including those that are specifically designed for a particular problem and are not applicable to the entire problem class .

Next we formally state the problem we solve; then we give an overview of related work.

[[**Formal problem statement.**]{}]{} Given a network with $P$ nodes, we associate $f_p$ and $X_p$ in  with the $p$th node of the network, $p=1,\ldots,P$. We make the following assumptions:

\[Ass:FunctionsAndSets\] Each $f_p: \mathbb{R}^n \xrightarrow{} \mathbb{R}$ is a convex function over $\mathbb{R}^n$, and each set $X_p$ is closed and convex.

\[Ass:ProblemSolvable\] Problem  is solvable.

\[Ass:Network\] The network is connected and its topology does not vary with time.

Assumption \[Ass:FunctionsAndSets\] involves the concepts of convexity of a function and convexity of a set. We assume the reader is familiar with these concepts, as well as with Lagrangian duality [@Boyd:ConvexOpti; @Bertsekas:Convex]. The second assumption states that there is at least one vector $x^\star$ that solves . This implies that the optimal cost $f^\star := \sum_{p=1}^P f_p(x^\star)$ is finite (and that $\cap_{p=1}^P X_p \neq \emptyset$). In Assumption \[Ass:Network\], a network is connected if there is a path between every pair of nodes.

Under the previous assumptions, we solve the following problem: *given a network, design a distributed algorithm that solves *. By distributed algorithm, we mean that there is no notion of a central or special node; also, no node in the network, except node $p$, has access to $f_p$ or $X_p$ at any time before or during the algorithm.

[[**Related work.**]{}]{} The first method addressing the generic problem  was based on an average consensus algorithm coupled with a subgradient method [@Lobel10-DistributedMultiAgentOptimization; @CXVOptimInSPandComm]. Such method has been extensively studied and is known to be robust to noise and link failures; however, it generally requires too many iterations to converge.

Another approach for solving  in a distributed way is via augmented Lagrangian methods, which consist of two loops: an outer loop updating the dual variables, and an inner loop updating the primal variables. There are some options for the algorithms in each loop. For the outer loop, the most common is the gradient method [@Ruszczynski; @Spars09; @Jakovetic11-CooperativeConvexOptimization], which degenerates in the so-called method of multipliers [@Bertsekas97-ParallelDistributed]. For the inner loop, common choices are the nonlinear Gauss-Seidel method [@Spars09; @Jakovetic11-CooperativeConvexOptimization] and Jacobi type algorithms, such as the diagonal quadratic approximation [@Ruszczynski]. In [@Mota11ICASSP] we applied fast gradient algorithms [@BeckTeboulleFISTA] in both loops. We mention that, although [@Ruszczynski; @Spars09; @Mota11ICASSP] do not address  specifically, they can be easily adapted to solve it, as explained in section \[Sec:RelatedAlgs\].

A special augmented Lagrangian-based algorithm, which uses one loop only, is the alternating direction method of multipliers (ADMM) [@Glowinski75-ADMMFirst; @Gabay76-ADMMFirst; @Bertsekas97-ParallelDistributed]. In fact, ADMM is the method of multipliers concatenated with one iteration of the nonlinear Gauss-Seidel, and therefore, it consists of just one loop. ADMM is not directly applicable to : one has to reformulate that problem first. There are several options for such reformulation. For example, [@Zhu09-DistributedInNetworkChannelCoding] proposes a distributed algorithm for solving specific instance of , based on a particular reformulation (the works [@Bazerque10-DistributedSpectrumSensing; @Kim11-CooperativeSpectrumSensingKrigedKalmanFiltering; @Forero10-ConsensusBasedDistributedSVMs; @Erseghe11-FastConsensusByADMM] apply the same algorithm to other instances of ). Another possible reformulation of  is explored in [@Schizas08-ConsensusAdHocWSNsPartI; @Schizas08-ConsensusAdHocWSNsPartII], which also proposes an ADMM-based algorithm, but in contrast with [@Zhu09-DistributedInNetworkChannelCoding], it uses two communication steps per iteration, or in other words, each node uses information from its second-order neighborhood. A comparison between [@Zhu09-DistributedInNetworkChannelCoding] and [@Schizas08-ConsensusAdHocWSNsPartI] for the average consensus problem is provided in [@Erseghe11-FastConsensusByADMM]; it is reported that both methods have similar convergence rates (and thus [@Schizas08-ConsensusAdHocWSNsPartI] uses twice as much communications as [@Zhu09-DistributedInNetworkChannelCoding]), while [@Zhu09-DistributedInNetworkChannelCoding] shows more resilience to noise. The algorithm we propose in this paper is also based on ADMM, but applied to a new reformulation of . We show through extensive simulations that the proposed algorithm requires less communications than any of the previous approaches.

All the above algorithms can solve  in a distributed way and will be briefly described in section \[Sec:RelatedAlgs\]. There are, however, other algorithms that can solve , but are not distributed. For example, [@Boyd11-ADMM] applies ADMM to distributed scenarios, but the resulting algorithms are not distributed in the sense that they require a central node or a special network topology. Our algorithm and the ones described before, in contrast, are decentralized and can run on any connected network topology.

[[**Contributions.**]{}]{} We propose a distributed algorithm for solving . This algorithm is an extension of our prior work [@Mota12-DistributedBP], where we solved a particular instance of , namely Basis Pursuit [@AtomicDecompBP]. Here, we first show that the algorithm in [@Mota12-DistributedBP] can be generalized to the problem class . Then we apply it to following problems: average consensus, two compressed sensing problems (different from the Basis Pursuit), and support vector machines. We also provide extensive simulations of several distributed algorithms solving . Our simulations show that the proposed algorithm outperforms, in terms of the number of used communications, any of the previous algorithms. One of the above problems, the average consensus, has been extensively studied due to its importance and simplicity. For example, many distributed algorithms, which cannot be trivially generalized to solve , have been proposed to solve it. Our algorithm, designed to solve the generic problem , can also solve the average consensus problem and, as our simulations show, it outperforms a state-of-the-art consensus algorithm [@Olshevsky11-ConvergenceSpeedDistributedConsensusAveraging], for many networks of interest.

[[**Organization.**]{}]{} The paper is organized as follows. In section \[Sec:AlgDerivation\] we derive the algorithm, by manipulating the problem and then applying ADMM. Section \[Sec:RelatedAlgs\] is dedicated to related algorithms; in particular, we describe the competing algorithms and show how we compare them with ours. In section \[Sec:Applications\] we take problems from several areas and see how they can be recast as . Additionally, we analyze the resulting algorithm for each one of these problems. Finally, in section \[Sec:Experiments\], we describe and show the results of our simulations. In Appendix \[App:ADMMProof\] we provide a proof of convergence of the particular version of ADMM we use, making this paper self-contained.

Algorithm Derivation {#Sec:AlgDerivation}
====================

In this section, we propose a distributed algorithm for solving . We begin by introducing some graph theory concepts and notation. Next, we address the scenario of a bipartite network, where a proof of convergence is given and the exposition is simpler; then, we generalize the algorithm for arbitrary networks. In this case, although we do not provide a proof of convergence, the algorithm always converges in practice, as shown in section \[Sec:Experiments\].

[[**Network notation.**]{}]{} We represent a network as an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V} = \{1,2,\ldots,P\}$ is the set of nodes and $\mathcal{E} \subset \mathcal{V} \times \mathcal{V}$ is the set of edges. The cardinality of these sets is represented respectively by $P$, the number of nodes, and $E$, the number of edges. An edge is represented by $\{i,j\} = \{j,i\}$, and $\{i,j\} \in \mathcal{E}$ means that nodes $i$ and $j$ can exchange data with each other. We define the neighborhood $\mathcal{N}_p$ of a node $p$ as the set of nodes connected to node $p$, but excluding it; the cardinality of this set, $D_p := |\mathcal{N}_p|$, is the degree of node $p$.

[[**Node coloring.**]{}]{} Given a network $\mathcal{G}$, a (node) coloring scheme is an assignment of numbers, which we call colors, to the nodes of the network such that no adjacent nodes have the same color. The colors will be represented as a set $\mathcal{C} = \{1,2,\ldots,C\}$, where $C$ is the total number of colors. For each color $c$, the set $\mathcal{C}_c \subset \mathcal{V}$ represents the nodes with color $c$, and $C_c$ is the number of nodes with that color, i.e., $|\mathcal{C}_c| = C_c$. A network $\mathcal{G}$ can have several coloring schemes, and $\chi(\mathcal{G})$ denotes the minimum number of colors required to color it; it is called its chromatic number. Finding a coloring scheme that uses only $\chi(\mathcal{G})$ colors is NP-hard for $\chi(\mathcal{G}) > 2$. Node-coloring is very important in wireless networks, since it is used in some medium access (MAC) protocols; consequently, many distributed node-coloring algorithms have been proposed, e.g., [@Kuhn06ComplexityDistributedGraphColoring; @Leith06DistributedLearningAlgorithms; @Duffy09ComplexityAnalysisDecentralizedGraphColoring; @Linial92LocalityDistributedGraphAlgorithms]. Henceforth, we will assume that the network $\mathcal{G}$ is given together with a coloring scheme $\mathcal{C}$.

[[**Problem manipulations.**]{}]{} We start by assuming the network $\mathcal{G}$ is bipartite or, equivalently, that $\chi(\mathcal{G}) = 2$, i.e., it can be colored with two colors. An example of a bipartite network is a grid. We start with this assumption for two reasons: the exposition is simpler, and we prove that our algorithm converges for this case. Without loss of generality, assume the nodes are ordered such that $\mathcal{C}_1 = \{1,2,\ldots,C_1\}$ and $\mathcal{C}_2 = \{C_1+1,C_1+2,\ldots,C\}$.

A common technique to decouple problem  is to assign copies of the global variable $x$ to each node and then constrain all copies to be equal. Denoting the copy held by node $p$ with $x_p \in \mathbb{R}^n$, problem  is written equivalently as $$\label{Eq:BipClonedProb}
    \begin{array}{cl}
      \underset{\bar{x} = (x_1,\ldots,x_P)}{\text{minimize}} & f_1(x_1) + f_2(x_2) +\cdots + f_P(x_P) \\
      \text{subject to} & x_p \in X_p\,,\quad p = 1,\ldots,P \\
                        & x_i = x_j \,,\quad \{i,j\} \in \mathcal{E}\,,
    \end{array}$$ where $\bar{x}= (x_1,\ldots,x_P) \in (\mathbb{R}^n)^P$ is the optimization variable. Problem  is no longer coupled by the intersection of the sets $X_p$ or by a common variable in all $f_p$’s, as , but instead by the new equations $x_i=x_j$, for all pairs of edges in the network $\{i,j\} \in \mathcal{E}$. These equations enforce all copies to be equal since the network is connected, due to Assumption \[Ass:Network\]. Note that they can be written more compactly as $(B^\top \otimes I_n) \bar{x} = 0$, where $B \in \mathbb{R}^{P\times E}$ is the node arc-incidence matrix of the graph, $I_n$ is the identity matrix in $\mathbb{R}^n$, and $\otimes$ is the Kronecker product. Each column of $B$ is associated with an edge $\{i,j\} \in \mathcal{E}$ and has $1$ and $-1$ in the $i$th and $j$th entry, respectively; the remaining entries are zeros. We can partition $B$ into two matrices $B_1$ and $B_2$, where $B_1$ contains the first $C_1$ rows and $B_2$ contains the remaining rows. Therefore, $(B^\top \otimes I_n)\bar{x} = (B_1^\top \otimes I_n)\bar{x}_1 + (B_2^\top \otimes I_n)\bar{x}_2$, where $\bar{x}_1 = (x_1,\ldots,x_{C_1}) \in (\mathbb{R}^{n})^{C_1}$ and $\bar{x}_2 = (x_{C_1 + 1},\ldots,x_{C}) \in (\mathbb{R}^{n})^{C_2}$. This enables rewriting  as $$\label{Eq:BipClonedProb2}
    \begin{array}{ll}
      \underset{\bar{x}_1,\bar{x}_2}{\text{minimize}} & \sum_{p \in \mathcal{C}_1}f_p(x_p) + \sum_{p \in \mathcal{C}_2}f_p(x_p) \\
      \text{subject to} & \bar{x}_1 \in \bar{X}_1\,\,,\,\, \bar{x}_2  \in \bar{X}_2\\
                        & (B_1^\top \otimes I_n) \bar{x_1} + (B_2^\top \otimes I_n) \bar{x_2} = 0\,,
    \end{array}$$ where $\bar{X}_i = \bigcap_{p \in \mathcal{C}_i} X_p$, for $i=1,2$. Problem  can be solved with ADMM, explained next.

[[**ADMM.**]{}]{} The alternating direction method of multipliers (ADMM) was proposed in the seventies by [@Glowinski75-ADMMFirst; @Gabay76-ADMMFirst]. Given two functions $g_1$ and $g_2$, two sets $X_1$ and $X_2$, and two matrices $A_1$ and $A_2$, ADMM solves $$\label{Eq:ProblemSolvedByADMM}
    \begin{array}{cl}
      \underset{x_1,x_2}{\text{minimize}} & g_1(x_1) + g_2(x_2) \\
      \text{subject to} & x_1 \in X_1\,\,,\,\, x_2\in X_2 \\
                        & A_1x_1 + A_2x_2 = 0\,.
    \end{array}$$ It uses the method of multipliers concatenated with an iteration of the Gauss-Seidel algorithm [@Bertsekas97-ParallelDistributed], i.e., it iterates on $k$: $$\begin{aligned}
        x_1^{k+1} &\in \arg\min_{x_1 \in X_1} L_{\rho}(x_1, x_2^{k};\lambda^k) 
                \label{Eq:ADMMAlg1}
                \\
        x_2^{k+1} &\in \arg\min_{x_2 \in X_2} L_{\rho}(x_1^{k+1}, x_2; \lambda^k)
                \label{Eq:ADMMAlg2}
                \\
        \lambda^{k+1} &= \lambda^k + \rho(A_1x_1^{k+1} + A_2x_2^{k+1})\,,
        \label{Eq:ADMMAlg3}
  \end{aligned}$$ where is the augmented Lagrangian of , $\lambda$ is the dual variable, and $\rho>0$ is a predefined parameter. More information about ADMM, including variations, applications, and related algorithms can be found, for example, in [@Bertsekas97-ParallelDistributed; @Boyd11-ADMM; @Combettes10-ProximalSplittingInSP]. The following theorem establishes its convergence.

\[Teo:ConvergenceADMM\] Assume $g_i:\mathbb{R}^{n_i}\xrightarrow{} \mathbb{R}$ is a convex function over $\mathbb{R}^{n_i}$, $X_i \subset \mathbb{R}^{n_i}$ is a closed convex set, and $A_i$ is a full column-rank matrix, for $i=1,2$. Also, assume problem  is solvable. Then, the sequence $\{(x_1^k,x_2^k,\lambda^k)\}$ generated by - converges to $(x_1^\star, x_2^\star, \lambda^\star)$, where

1.  $(x_1^\star,x_2^\star)$ solves 

2.  $\lambda^\star$ solves the dual problem of : $$\label{Eq:ADMMDualProb}
                        \underset{\lambda}{\text{minimize}}\,\,\, G_1(\lambda) + G_2(\lambda)\,,$$ where $G_i(\lambda) = \inf_{x_i \in X_i} (g_i(x_i) + \lambda^\top A_ix_i)$, $i=1,2$.

This theorem is more general than the versions in [@Boyd11-ADMM] and [@Bertsekas97-ParallelDistributed Prop.4.2]. Namely, [@Boyd11-ADMM] only proves objective convergence, but since it does not assume $A_1$ and $A_2$ to be full column-rank, the sequence $\{(x_1^k,x_2^k,\lambda^k)\}$ might not even converge. On the other hand, [@Bertsekas97-ParallelDistributed] proves claims $1)$ and $2)$ only for the special case when one of the $A_i$’s is the identity. Since we could not find a proof of this version of the theorem in literature, we provide it in Appendix \[App:ADMMProof\]. The proof is, however, very similar to the ones in [@Bertsekas97-ParallelDistributed; @Boyd11-ADMM].

[[**Applying ADMM.**]{}]{} Clearly, has the same format as : just make the following associations $$\label{Eq:BipAssociations}
        g_i(x_i) = \sum_{p \in \mathcal{C}_i} f_p(x_p)\,,\quad X_i = \bar{X}_i\,,\quad A_i = B_i^\top \otimes I_n\,,$$ for $i=1,2$. The following theorem guarantees that, under Assumptions \[Ass:FunctionsAndSets\]-\[Ass:Network\], problem  satisfies all the conditions of Theorem \[Teo:ConvergenceADMM\]; thus, we can use ADMM to solve .

\[Teo:BipApplyADMM\] Let Assumptions \[Ass:FunctionsAndSets\]-\[Ass:Network\] hold and consider problem . Then, for $i=1,2$, the function $\sum_{p \in \mathcal{C}_i} f_p(x_p)$ is convex over $\mathbb{R}^n$, the set $\bar{X}_i$ is closed and convex, and the matrix $B_i^\top \otimes I_n$ has full column-rank. Furthermore, problem  is solvable.

All conclusions of Theorem \[Teo:BipApplyADMM\], except that $B_i^\top \otimes I_n$ has full column-rank, are derived straightforwardly Assumptions \[Ass:FunctionsAndSets\] and \[Ass:ProblemSolvable\] and the equivalence between  and . To see how Assumption \[Ass:Network\] implies that each $B_i^\top \otimes I_n$ has full column-rank for $i=1,2$, note that it is sufficient to prove that $B_i^\top$ has full column-rank. If we prove that $B_iB_i^\top$ has full rank, then the results follows because $\text{rank}(B_iB_i^\top) = \text{rank}(B_i^\top)$. Now note that $B_iB_i^\top$ is a diagonal matrix, where in the diagonal are the degrees of the nodes belonging to the subnetwork composed by the nodes in $\mathcal{C}_i$. Assumption \[Ass:Network\] states that the network is connected and thus no node can have degree $0$; consequently, $B_iB_i^\top$ has full rank.

[[**Algorithm for bipartite networks.**]{}]{} We now apply ADMM, i.e., - directly to . By analyzing the augmented Lagrangian , we will see that  yields $C_1$ optimization problems that can be solved in parallel. Similarly, yields $C_2$ optimization problems that can be solved in parallel. In fact, developing the squared term in  we obtain for the $x_1$-minimization : $$\label{Eq:BipMinx1}
        x_1^{k+1} \in \arg\min_{x_1 \in X_1} g_1(x_1) + \bar{\eta}_1^\top x_1 + \frac{\rho}{2} x_1^\top (A_1^\top A_1) x_1\,,$$ where $\bar{\eta}_1 :=A_1^\top \lambda^k + \rho (A_1^\top A_2)x_2^k$. Using , the quadratic term in  becomes $(\rho/2)\bar{x}_1^\top(B_1B_1^\top \otimes I_n) \bar{x}_1$. As we had seen before, $B_1B_1^\top$ is a diagonal matrix with the degree of node $p$, $D_p$, in the $(pp)$th entry. This quadratic term can then be rewritten as $(\rho/2)\sum_{p \in \mathcal{C}_1}D_p\|x_p\|^2$. Regarding the linear term, $$\begin{aligned}
          \bar{\eta}_1^\top \bar{x}_1 
        &=
          \Bigl((B_1\otimes I_n)\lambda^k  + \rho(B_1B_2^\top \otimes I_n)\bar{x}_2^k\Bigr)^\top \bar{x}_1
        \notag
        \\
        &=
          \sum_{p \in \mathcal{C}_1} 
                \Bigl(
                    \sum_{j \in \mathcal{N}_p}\text{sign}(j-p) \lambda_{\{p,j\}}^k - \rho\sum_{j \in \mathcal{N}_p}x_j^k
                \Bigr)^\top x_p\,,
        \label{Eq:BipGamma}
    \end{aligned}$$ where $\text{sign}(w) = 1$ if $w \geq 0$ and $\text{sign}(w) = -1$ otherwise. We decomposed the dual variable $\lambda$ as $(\ldots,\lambda_{\{i,j\}},\ldots)$, where $\lambda_{\{i,j\}}= \lambda_{\{j,i\}}$ is associated with the edge $\{i,j\} \in \mathcal{E}$. To see why  holds, note that the $(ij)$th entry of $B_1B_2^\top$ is $-1$ if $\{i,j\} \in \mathcal{E}$ and $0$ otherwise. In conclusion, using  on , we get $C_1$ optimization problems that can be solved in parallel: $$x_p^{k+1} \in \arg\min_{x_p \in X_p} f_p(x_p) + ({\gamma_p^k}^\top -\rho\sum_{j \in \mathcal{N}_p}x_j^k)^\top x_p + \frac{\rho D_p}{2}\|x_p\|^2\,,$$ where $\gamma_p^k := \sum_{j \in \mathcal{N}_p}\text{sign}(j-p) \lambda_{\{p,j\}}^k$, for $p \in \mathcal{C}_1$, i.e., all nodes with color $1$. A similar analysis for  yields the same problem for each node in $\mathcal{C}_2$, but using the neighbors’ estimates at $k+1$.

Note that the optimization problem each node solves depends on the dual variable $\lambda$ through $\gamma_p$. In fact, we can derive a formula for updating $\gamma_p$ from : write  edge-wise, $\lambda_{\{i,j\}}^{(k+1)} = \lambda_{\{i,j\}}^{(k)} + \rho \,\text{sign}(j-p)(x_i^{k+1} - x_j^{k+1})$, and insert this formula in the definition of $\gamma_p^k$; the result is $
        \gamma_p^{k+1} = \gamma_p^k + \rho \sum_{j \in \mathcal{N}_p}(x_p^{k+1} - x_j^k)\,.
    $ In sum, we arrive at the following algorithm, named D-ADMM, after *Distributed*-ADMM.

\[Alg:ADMMBipartite\]

for all $p \in \mathcal{V}$, set $\gamma_{p}^{1} = x_p^{1} = 0$ and $k=1$ \[SubAlg:ADMMBipLoop\] Set $v_p^{k} =  \gamma_p^{k} -\rho \sum_{j \in \mathcal{N}_p} x_j^{k}$ and find \[SubAlg:ADMMBipOptimProb\] $$x_p^{k+1}   = \begin{array}[t]{cl}
                           \textrm{argmin} & f_p(x_p) + {v_p^{k}}^\top x_p + \frac{D_p \rho}{2}\|x_p\|^2\\
                            \textrm{s.t.} & x_p \in X_p
                          \end{array}$$ Send $x_p^{k+1}$ to $\mathcal{N}_p$ \[SubAlg:ADMMBipCommunication\]

\[SubAlg:ADMMBipEndLoop\]

Repeat \[SubAlg:ADMMBipLoop\]-\[SubAlg:ADMMBipEndLoop\] for all $p \in \mathcal{C}_2$, replacing $x_j^{k}$ by $x_j^{k+1}$ in $v_p^k$

$
            \gamma_p^{k+1} = \gamma_p^{k} + \rho \sum_{j \in \mathcal{N}_p} (x_p^{k+1} -  x_j^{k+1})
        $ \[SubAlg:ADMMBipDualVarUp\] $k \gets k+1$

In Algorithm \[Alg:ADMMBipartite\], there are two groups of nodes: $\mathcal{C}_1$ and $\mathcal{C}_2$. In each group, the nodes are not neighbors between themselves and operate at the same time. In particular, one node receives $x_j^k$ or $x_j^{k+1}$ from all its neighbors and then solves the optimization problem of step \[SubAlg:ADMMBipOptimProb\], finding its estimate $x_p^{k+1}$; it then broadcasts $x_p^{k+1}$ to its neighbors. After having received all the estimates from its neighbors, node $p$ can update $\gamma_p$ as in step \[SubAlg:ADMMBipDualVarUp\].

Algorithm \[Alg:ADMMBipartite\] only applies to bipartite networks and its convergence is guaranteed by Theorems \[Teo:ConvergenceADMM\] and \[Teo:BipApplyADMM\]. Next, we generalize it for any type of network.

[[**Algorithm for arbitrary networks.**]{}]{} So far we assumed the network was bipartite, i.e., it could be colored with just two colors: $\chi(\mathcal{G})=2$. Now, suppose the network is colored with $C\geq 2$ colors. We derived problem  from problem  by partitioning the node-arc incidence matrix $B$ into $\begin{bmatrix}B_1^\top & B_2^\top\end{bmatrix}^\top$, where $B_1$ and $B_2$ contain the rows corresponding to nodes in $\mathcal{C}_1$ and $\mathcal{C}_2$, respectively. For convenience, we had assumed the nodes were ordered such that $\mathcal{C}_1 = \{1,2,\ldots,C_1\}$ and $\mathcal{C}_2 = \{C_1+1,C_1+2,\ldots,C\}$. Without loss of generality, assume a similar ordering for the nodes in the current case, i.e., $\mathcal{C}_1 = \{1,2,\ldots,C_1\}$, $\mathcal{C}_2 = \{C_1+1,\ldots,C_1+C_2\}$, …, $\mathcal{C}_C = \{\sum_{c=1}^{C-1}C_c+1,\ldots,C\}$. This induces a natural partition of $B$ as $\begin{bmatrix}B_1^\top & B_2^\top & \cdots & B_C^\top\end{bmatrix}^\top$. The counterpart of  for this case is then $$\label{Eq:GenNetCloned}
            \begin{array}{ll}
                \underset{\bar{x}_1,\ldots,\bar{x}_C}{\text{minimize}} & \sum_{c=1}^C \sum_{p \in \mathcal{C}_p} f_p(x_p) \\
                \text{subject to} & \bar{x}_c \in \bar{X}_c\,,\quad c = 1,\ldots,C \\
                                  & \sum_{c=1}^C (B_c^\top \otimes I_n)\bar{x}_c= 0\,,
            \end{array}$$ where the variable is $(\bar{x}_1,\ldots,\bar{x}_C) \in (\mathbb{R}^n)^P$. To solve , we cannot apply ADMM - directly, but its generalization to solve $$\begin{array}{cl}
        \underset{x_1,\ldots,x_C}{\text{minimize}} & g_1(x_1) + \cdots + g_C(x_C) \\
        \text{subject to} & x_c \in X_c\,,\quad c=1,\ldots,C \\
                          & A_1x_1 + \cdots + A_Cx_C = 0\,.
      \end{array}$$ Such a generalization is straightforward, but there are no convergence guarantees in the sense that there is no known proof of Theorem \[Teo:ConvergenceADMM\] for the resulting algorithm. Yet, experimental results, as the ones presented here, suggest that Theorem \[Teo:ConvergenceADMM\] might still hold.

The following algorithm is a straightforward generalization of Algorithm \[Alg:ADMMBipartite\] to solve .

\[Alg:DADMMGeneral\]

for all $p \in \mathcal{V}$, set $\gamma_{p}^{1} = x_p^{1} = 0$ and $k=1$ $$v_p^{k} = \gamma_p^{k}-
                \rho \sum_{\begin{subarray}{c}
                             j \in \mathcal{N}_p \\
                             j < p
                           \end{subarray}
                }x_j^{k+1} - \rho \sum_{\begin{subarray}{c}
                             j \in \mathcal{N}_p \\
                             j > p
                           \end{subarray}
                }x_j^{k}$$ \[SubAlg:GeneralDADMMv\] and find $$x_p^{k+1} = \begin{array}[t]{cl}
                            \textrm{argmin} & f_p(x_p) + {v_p^{k}}^\top x_p + \frac{D_p \rho}{2}\|x_p\|^2\\
                            \textrm{s.t.} & x_p \in X_p
                          \end{array}$$ \[SubAlg:ADMMGenProb\] Send $x_p^{k+1}$ to $\mathcal{N}_p$ \[SubAlg:ADMMGenComm\]

$
            \gamma_p^{k+1} = \gamma_p^{k} + \rho \sum_{j \in \mathcal{N}_p} (x_p^{k+1} -  x_j^{k+1})
        $ \[SubAlg:ADMMGenDualVarUp\] $k \gets k+1$

In Algorithm \[Alg:DADMMGeneral\] there are $C$ groups of nodes, arranged by colors. As in the bipartite case, in each group, the nodes are not neighbors between themselves and operate at the same time. In practice, they cannot operate at the same time because there is no central coordination. As seen in the expression for $v_p^k$ (step \[SubAlg:GeneralDADMMv\]), one node, knowing the color of its neighbors, can operate only after having received the estimates $x_k^{k+1}$ from the nodes with lower color. Excepting that, the algorithm is very similar to Algorithm \[Alg:DADMMGeneral\].

As said before, although there is no convergence guarantee for Algorithm \[Alg:DADMMGeneral\], it converges in practice for several kinds of problems, and we have never seen a nonconvergent example. We point out that, using a reasoning similar to that of the discussion after Theorem \[Teo:BipApplyADMM\], each matrix $B_i^\top$ can be proved to be full column-rank.

Related Algorithms {#Sec:RelatedAlgs}
==================

We now review some algorithms for solving . We start with describing the performance measure we will use and its practical relevance.

[[**Performance measure: communication steps.**]{}]{} We divide the algorithms that solve  in two kinds: single-looped and double-looped. The single-looped algorithms have just one (nested) loop, and in each iteration, each node receives the neighbors’ estimates, computes its on estimate, and broadcasts it to its neighbors. The double-looped algorithms have two nested loops and they are based on augmented Lagrangian methods (see ). The outer loop serves to update the dual variable of the augmented Lagrangian; the inner loop serves for each node to receive the neighbors’ estimates, compute its own estimate, and broadcast it to the neighbors, as one iteration of a single-looped algorithm. The natural way of comparing all algorithms is then through a *communication step*, which will denote one iteration (for the single-looped algorithms) or one inner loop iteration (for the double-looped algorithms). We say that an algorithm is more efficient than another one if it uses less communication steps for the same solution accuracy.

Communication steps is a fair measure because it compares, not only similar computational operations at each node (as we will soon see), but also the amount of communications each node requires to achieve a solution with a given precision. Furthermore, one of the main concerns in designing algorithms and protocols for wireless networks is reducing the number of communications; the reason is because communicating is the slowest and most energy-consuming task [@Akyildiz02-WirelessNetworksASurvey; @Krishnamachari05-NetworkingWirelessSensors]. The same concern is shared in supercomputing environments [@Graham05-FutureOfSuperComputing]. In contrast with execution time, communication steps are application and implementation independent. For example, given a distributed algorithm, we get different execution times if we change any of the following: medium-access protocol, transmission medium, scheduling, platform for the nodes, or even the way the algorithm is coded. Communication steps are, on the other hand, invariant to any of these. Note that, given a network with $E$ edges, if we multiply $2E$ by the number of communication steps, we obtain the total number of communications in the network.

We point out, however, that our algorithm, D-ADMM, is asynchronous, i.e., not all nodes perform the same tasks at the same time, in contrast with most of the algorithms presented next. More concretely, D-ADMM assumes a coloring for the nodes, and each node operates according to its color. Other algorithms do not rely on any coloring scheme, enabling them to be synchronous, i.e., all nodes can perform the same tasks at the same time.[^2] Therefore, if one communication step in D-ADMM takes $T$ time units, it would take $T/C$ units in a synchronous algorithm, where $C$ is the number of colors, in case the nodes were able to transmit messages at the same time. [[**Gradient/Subgradient method.**]{}]{} The classical approach for solving  is to use a subgradient method combined with an averaging consensus algorithm. This approach was taken in [@Lobel10-DistributedMultiAgentOptimization], and consists of iterating in $k$ $$\label{Eq:OASubgradAlg}
        x_p^{k+1} = \Bigl[a_{pp}^{k} x_p^{k} + \sum_{j \in \mathcal{N}_p} a_{pj}^{k} x_j^{k} - \alpha^{k} g_p^{k}\Bigr]_{X_p}^+\,,$$ for each node $p$. In , $[s]_S$ represents the projection of the point $s$ onto the (convex) set $S$, $x_p^k$ the estimate of node $p$ at iteration $k$, $g_p^k$ the gradient (or a subgradient, in case the objective is not differentiable) of $f_p$ at the point $a_{pp}^{k} x_p^{k} + \sum_{j \in \mathcal{N}_p} a_{pj}^{k} x_j^{k}$, and $a_{ij}^k$ represents the influence node $j$ exerts on node $i$ at iteration $k$. There are some conditions the weights $a_{ij}^k$ have to satisfy [@Lobel10-DistributedMultiAgentOptimization], for example, for each $k$, the matrix with entries $a_{ij}^k$ should be doubly stochastic. In each iteration, node $p$ receives the estimates $x_{j}^k$ from its neighbors, computes , and then broadcasts $x_p^{k+1}$ to its neighbors. It is thus a (synchronous) single-looped algorithm. Although this algorithm is provably robust to noise, it is known to perform very slowly, requiring lots of iterations to converge.

[[**ADMM-based algorithms.**]{}]{} We now describe two algorithms that are based on ADMM, like ours. Both of them are based on different reformulations of , yielding two different algorithms. The first algorithm is [@Zhu09-DistributedInNetworkChannelCoding] (it also appears in [@Bazerque10-DistributedSpectrumSensing; @Kim11-CooperativeSpectrumSensingKrigedKalmanFiltering; @Forero10-ConsensusBasedDistributedSVMs; @Erseghe11-FastConsensusByADMM]), described in Algorithm \[Alg:Giannakis\].

\[Alg:Giannakis\]

for all $p \in \mathcal{V}$, set $\gamma_{p}^{1} = x_p^{1} = 0$ and $k=1$

set $
                v_p^{k} = \gamma_p^{k}-
                \rho \sum_{j \in \mathcal{N}_p \cup \{p\}} x_j^{k}
            $ and find $$x_p^{k+1} = \begin{array}[t]{cl}
                            \textrm{argmin} & f_p(x_p) + {v_p^{k}}^\top x_p + \rho D_p\|x_p\|^2\\
                            \textrm{s.t.} & x_p \in X_p
                          \end{array}$$ \[SubAlg:Giannakis-Optm\]

Send $x_p^{k+1}$ to $\mathcal{N}_p$, and receive $x_j^{k+1}$, $j \in \mathcal{N}_p$

$
            \gamma_p^{k+1} = \gamma_p^{k} + \rho \sum_{j \in \mathcal{N}_p} (x_p^{k+1} -  x_j^{k+1})
        $ \[SubAlg:ADMMGenDualVarUp\] $k \gets k+1$

Although derived from a different reformulation of , Algorithm \[Alg:Giannakis\] algorithm is very similar to D-ADMM. The three main differences are:

1.  Algorithm \[Alg:Giannakis\] is synchronous, i.e., all the nodes can perform the same tasks at the same time;

2.  while in Algorithm \[Alg:Giannakis\] node $p$ uses $x_p^k$ to construct the vector $v_p^k$ (see step \[SubAlg:Giannakis-Optm\]), D-ADMM does not;

3.  Algorithm \[Alg:Giannakis\] is proven to converge for any connected network, and D-ADMM is only proven to converge in bipartite networks (although in practice it converges in any connected network).

We will see that D-ADMM requires almost always less communication steps to achieve a given accuracy than Algorithm \[Alg:Giannakis\], for several instances of .

Another ADMM-based algorithm is [@Schizas08-ConsensusAdHocWSNsPartI] (also appearing in [@Schizas08-ConsensusAdHocWSNsPartII; @Erseghe11-FastConsensusByADMM]). That algorithm, in contrast with D-ADMM and Algorithm \[Alg:Giannakis\], uses two communication steps per iteration. Since the number of iterations it takes to achieve a given accuracy in the average consensus problem is similar to Algorithm \[Alg:Giannakis\] [@Erseghe11-FastConsensusByADMM], it takes the double of the communication steps.

[[**Double-looped algorithms.**]{}]{} Several double-looped algorithms were proposed to solve particular optimization problems, which as D-ADMM, can be generalized to solve . In particular, they can solve , by solving the following dual problem of : $$\label{Eq:DualProblemCloned}
    \underset{\lambda}{\text{maximize}} \,\,\, L_{\rho}(\lambda)\,,$$ where $L_\rho(\lambda)$ is the (augmented) dual function $$\label{Eq:LagrangianProblemCloned}
    L_\rho(\lambda) = 
    \begin{array}[t]{cl}
      \underset{\bar{x}}{\inf} & \sum_{p=1}^P f_p(x_p) + \sum_{\{i,j\}\in \mathcal{E}} \phi_{\lambda_{\{i,j\}}}^\rho (x_i - x_j) \\
      \text{s.t.} & x_p \in X_p\,,\quad p=1,\ldots,P\,,
    \end{array}$$ and $\phi_{v}^\rho(\eta) := v^\top \eta + \frac{\rho}{2}\|\eta\|^2$. It can be proved that $L_\rho(\lambda)$ is differentiable, and thus  can be solved with a gradient method, or even with a fast (Nesterov) gradient method. The outer loop consists of applying such algorithm: it updates the dual variables $\lambda_{\{i,j\}}$. However, to compute the gradient of $L_\rho(\lambda)$, one needs to solve the optimization problem in ; and this requires another iterative algorithm, consisting of the second nested loop. Several approaches have been used to solve that problem, e.g., nonlinear Gauss-Seidel methods [@Spars09; @Jakovetic11-CooperativeConvexOptimization], Jacobi-type methods [@Ruszczynski], and Nesterov methods [@Mota11ICASSP].

Of all the algorithms described in this section, only gradient/subgradient method does not solve an optimization problem at each iteration. All the remaining algorithms solve the same problem as D-ADMM (step \[SubAlg:ADMMGenProb\] of Algorithm \[Alg:DADMMGeneral\]), but possibly with different parameters. As illustrated in section \[Sec:Experiments\], in practice, our D-ADMM performs generally better than any of algorithms described in this subsection. The only true competitor is Algorithm \[Alg:Giannakis\].

Applications {#Sec:Applications}
============

We now take some optimization problems from signal processing and see how they can be recast as . We also review some of the conventional methods to solve them in a distributed way. Then, in section \[Sec:Experiments\], we present simulation results of D-ADMM and other algorithms applied to these problems.

[[**Consensus.**]{}]{} Consensus is one of the most fundamental problems in networks. Given a network with $P$ nodes, node $p$ generates a number or measurement, say $\theta_p$, and the goal is to compute the average $\theta^\star = (1/P)\sum_{p=1}^P \theta_p$ in every node. The classical approach to consensus [@DeGroot74ReachingConsensus; @Olshevsky11-ConvergenceSpeedDistributedConsensusAveraging] is $$\label{Eq:ApplicationConsensusLinearSys}
    x^{k+1} = A x^k\,,\qquad x^0 = \theta\,,$$ where the $p$th entry of $x^k = (x_1^k,\ldots,x_P^k) \in \mathbb{R}^P$ contains the estimate of node $p$ at time $k$ and $\theta = (\theta_1,\ldots,\theta_P)\in \mathbb{R}^P$ is the vector collecting the measurements. The matrix $A$ reflects the topology of the network by having $a_{ij} = 0$ for $\{i,j\}\not\in \mathcal{E}$. Thus, is equivalent to $x_p^{k+1} = \sum_{j \in \mathcal{N}_p \cup \{p\}}a_{pj} x_j^k$, i.e., the estimate of node $p$ at $k+1$ will be a weighted average of its estimate and the neighbors’ estimates at $k$. There are some constraints on choosing the weights, for example, $\sum_{j \in \mathcal{N}_p \cup \{p\}}a_{pj} = 1$. See [@Olshevsky11-ConvergenceSpeedDistributedConsensusAveraging] for a recent survey on consensus algorithms using the approach .

Another approach to consensus, recently adopted in [@Erseghe11-FastConsensusByADMM], is to solve the optimization problem $$\label{Eq:ApplicationConsensusOptimProb}
      \underset{x}{\text{minimize}} \,\,\, \sum_{p=1}^P (x - \theta_p)^2$$ in a distributed way. The solution of  is $\theta^\star = (1/P)\sum_{p=1}^P\theta_p$. Clearly, has the format of  by assigning $f_p(x) = (x-\theta_p)^2$, $X_p = \mathbb{R}$; thus, it can be solved by D-ADMM or any of the methods described in the last section. The work [@Erseghe11-FastConsensusByADMM] applies Algorithm \[Alg:Giannakis\] for solving .

In all algorithms, node $p$ has to solve at each iteration $$\underset{x}{\text{minimize}} \,\,\, (x - \theta_p)^2 + v^\top x + c x^2\,,$$ which has the closed-form solution $x^\star = (2\theta_p - v)/(2(1+c))$.

[[**Compressed sensing problems.**]{}]{} Compressed sensing is a new paradigm for the acquisition of signals. For an introductory survey see, for example, [@IntroCompressiveSamplingCandesWakin]. In compressed sensing, a signal is compressed at the same time as it is acquired, in contrast to the classical acquire-then-compress paradigm. Although compressed sensing makes the acquisition process simpler, the reconstruction of a compressed signal is more complicated, since an optimization problem has to be solved. The goal of this optimization problem is to find the sparsest solution of a linear system $Ax = b$, where $A \in \mathbb{R}^{m\times n}$ is a matrix related with the acquisition process and $b \in \mathbb{R}^m$ is the acquired/compressed signal. If $A$ satisfies some technical conditions [@IntroCompressiveSamplingCandesWakin], then the solution to the *Basis Pursuit* (BP) [@AtomicDecompBP], $$\label{Eq:BP}
            \begin{array}{ll}
                \underset{x}{\text{minimize}} & \|x\|_1 \\
                \text{subject to} & Ax = b\,,
            \end{array}$$ recovers the sparsest solution of $Ax = b$. In , $\|x\|_1 = |x_1| + \cdots + |x_n|$ is the $\ell_1$-norm of $x$. A variation of  is the LASSO [@IntroCompressiveSamplingCandesWakin]: $$\label{Eq:Lasso}
        \begin{array}{ll}
            \underset{x}{\text{minimize}} & \|x\|_1 \\
            \text{subject to} & \|Ax - b\| \leq \sigma\,,
        \end{array}$$ where $\sigma > 0$ is a fixed parameter and $\|\cdot\|$ the $\ell_2$-norm. Problem  is more appropriate when there is noise in the acquisition process and we have an upper bound $\sigma$ on the norm of the noise. Closely related with , is $$\label{Eq:BPDN}
        \underset{x}{\text{minimize}} \,\,\, \|Ax - b\|^2 + \beta\|x\|_1\,,$$ called *Basis Pursuit Denoising* (BPDN) [@AtomicDecompBP]. The parameter $\beta > 0$ regulates the tradeoff between the sparsity of the signal and the weight of the noise.

Here, we are interested in solving distributed versions of -. Namely, we assume the matrix $A$ is partitioned into $P$ blocks, either by rows or by columns, as depicted in Figure \[Fig:PartitionOfA\]. We assign the $p$th block to the $p$th node of a network with $P$ nodes. For applications of this problem see, for example, our previous work [@Mota12-DistributedBP]. In [@Mota12-DistributedBP], we solved BP with D-ADMM for both the row and the column partition. In this paper, we solve  for the column partition, and  for the row partition. The reverse cases, i.e., for the row partition or  for the column partition cannot be trivially recast as .

[[**LASSO: column partition.**]{}]{} Assume the matrix $A$ is partitioned by columns, $A = \begin{bmatrix}A_1 & \cdots & A_P\end{bmatrix}$, where the $p$th block of $A$ is only known at node $p$. Assume vector $b$, parameter $\sigma$, and the number of nodes $P$ are available at all nodes. This problem cannot be directly recast as ; we will have to do it through duality. However, if we solve the “simple dual” of , we will not be able to recover solution to the primal problem afterwards. The reason is that the objective of  is not strictly convex. We thus start by regularizing  by adding a small term to its objective, making it strictly convex: $$\label{Eq:LassoRegularized}
        \begin{array}{ll}
            \underset{x}{\text{minimize}} & \|x\|_1 + \frac{\delta}{2} \|x\|^2 \\
            \text{subject to} & \|Ax - b\| \leq \sigma\,,
        \end{array}$$ where $\delta >0$ is sufficiently small. This regularization is inspired by [@Friedlander07ExactRegularizationConvexPrograms], which establishes some conditions for regularizations of this type to be exact. By exact, we mean that there exists $\bar{\delta} > 0$ such that the solution of  is always a solution of , for $\delta \leq \bar{\delta}$. One of the conditions for exact regularization is that the objective is linear and the constraint set is the intersection of a linear system with a closed polyhedral cone. LASSO  can be recast as $$\label{Eq:LassoManipulated}
    \begin{array}{ll}
      \underset{x,t,u,v}{\text{minimize}} & 1_n^\top t \\
      \text{subject to} & \|u\| \leq v \\
                        & u = Ax - b\,,\,\,\, v = \sigma \\
                        & x \leq t\,,\,\,\, -x\leq t\,,
    \end{array}$$ where $1_n \in \mathbb{R}^n$ is the vector of ones, and $(x,t,u,v)$ the variable. Although the objective of  is linear and its constraint set is the intersection of a linear system with a closed convex cone $\mathcal{K} := \{(u,v)\,:\, \|u\|\leq v\}$, there is not a proof of exact regularization for , because $\mathcal{K}$ is not polyhedral. However, experimental results in [@Friedlander07ExactRegularizationConvexPrograms] suggest that exact regularization might occur for non-polyhedral cones. We will also see in the results of our simulations that approximating  with  can still yield very accurate solutions.

Before computing the dual of , we first introduce a new variable $y \in \mathbb{R}^{m}$ and make the column partition explicit: $$\label{Eq:LassoWithY}
    \begin{array}{ll}
      \underset{x,y}{\text{minimize}} & \sum_{p=1}^P (\|x_p\|_1 + \frac{\delta}{2}\|x_p\|^2) \\
      \text{subject to} & \|y\| \leq \sigma \\
                        & y = \sum_{p=1}^P A_p x_p - b\,.
    \end{array}$$ The dual problem of , by dualizing the last constraint only, can be shown to be equivalent to $$\underset{\lambda}{\text{minimize}} \,\,\, \sum_{p=1}^P g_p(\lambda)\,,$$ where $g_p(\lambda) := \frac{1}{P}(b^\top \lambda + \sigma\|\lambda\|) - \inf_{x_p}(\|x_p\|_1 + (A_p^\top \lambda)^\top x_p + \frac{\delta}{2}\|x_p\|^2)$ is the function associated to node $p$. Each function $g_p$ is convex because it is the pointwise supremum of convex functions [@Bertsekas:Convex Prop.1.2.4]. The problem each node has to solve in each iteration is addressed in Appendix \[App:ProbLasso\].

We point out that problem  can be ill-conditioned due to the small value that has to be used for $\delta$, if we want to get a good approximation of the original problem. In particular, the problem solved at each node can be very time-consuming, since it may require lots of internal iterations. We used $\delta = 10^{-3}$ in our simulations. While such value for $\delta$ required a considerable amount of iterations for solving the internal problem of each node (we set $500$ as the number of maximum iterations), it also allowed to get solutions of the original problem  with relative errors of $0.1\%$, using few communications. BPDN, addressed next, handles the same type of problems as LASSO, but its distributed implementation does not require solving ill-conditioned problems.

[[**BPDN: row partition.**]{}]{} In BPDN , we assume that both the matrix $A$ and the vector $b$ are partitioned by rows: $A = \begin{bmatrix}A_1^\top & \cdots & A_P^\top\end{bmatrix}^\top$, $b = \begin{bmatrix}b_1^\top & \cdots & b_P^\top\end{bmatrix}^\top$. Therefore, can be readily rewritten as $$\label{Eq:BPDN2}
        \underset{x}{\text{minimize}} \,\,\, \sum_{p=1}^P \Bigl( \|A_px - b_p\|^2 + \frac{\beta}{P}\|x\|_1\Bigr)\,.$$ The identification between  and  is immediate: set $f_p(x) := \|A_px - b_p\|^2 + \frac{\beta}{P}\|x\|_1$. The problem that has to be solved at node $p$, $$\label{Eq:BPDNAtEachNode}
        \underset{x}{\text{minimize}}\,\,\, \|A_px - b_p\|^2 + \frac{\beta}{P}\|x\|_1 + v^\top x + c\|x\|^2\,,$$ can be solved, for example, with the GPSR solver [@GPSR].

[[**Distributed support vector machines.**]{}]{} In a *Support Vector Machine* (SVM) [@Bishop06-PatternRecognitionMachineLearning Ch.7], we are given data in a matrix $A \in \mathbb{R}^{m \times n}$, where each row of $A$ represents a data point in $\mathbb{R}^n$. Each data point can be classified as belonging to group $\mathcal{X}$ or to group $\mathcal{Y}$. If data point $a_i$ belongs to $\mathcal{X}$ (resp. $\mathcal{Y}$), we set the $i$th entry of a vector $d \in \mathbb{R}^m$ to $+1$ (resp. $-1$). The goal is to find the parameters $(s,r) \in \mathbb{R}^n\times \mathbb{R}$ of a linear separator $z(a) = s^\top a + r$ such that $z(a) > 0$ if $a \in \mathcal{X}$, and $z(a) < 0$ if $a \in \mathcal{Y}$. Once the pair $(s,r)$ is found, we can attempt to determine if a new point $a$ belongs to $\mathcal{X}$ or $\mathcal{Y}$ just by analyzing the sign of $z(a)$. There are several ways of finding $(s,r)$; one of them is solving the quadratic program $$\label{Eq:SVMProb}
    \begin{array}{ll}
      \underset{s,r}{\text{minimize}} & \|s\|^2 \\
      \text{subject to} & D(As - 1_m r) \geq 1_m\,,
    \end{array}$$ where $D = \text{Diag}(d_1,\ldots,d_m)$, and $1_m$ is the vector of ones in $\mathbb{R}^m$. Many times, the data is distributed among several entities, and processing it in a distributed way is prohibitive for complexity or privacy reasons. Here, we assume the model of [@Mangasarian07-PrivacyPreservingClassificationHP; @Forero10-ConsensusBasedDistributedSVMs]: there are $P$ entities, and each entity $p$ has $m_p$ data points, or equivalently, a block $A_p \in \mathbb{R}^{m_p \times n}$ of rows of $A$, and their corresponding classification, i.e., the $p$th diagonal block $D_p\in \mathbb{R}^{m_p \times m_p}$ of $D$, such that $m_1 + \cdots + m_P = m$. This enables us to rewrite  as $$\label{Eq:SVMProb2}
    \begin{array}{ll}
      \underset{s,r}{\text{minimize}} & \sum_{p=1}^P\|s\|^2 \\
      \text{subject to} & D_p(A_ps - 1_{m_p}r) \geq 1_{m_p}\,,\quad p = 1,\ldots,P\,.
    \end{array}$$ Problem  has the format of : just set $f_p(s,r) = \|s\|^2$ and $X_p = \{(s,r) \in \mathbb{R}^n\times \mathbb{R}\,:\,  D_p(A_ps - 1_{m_p}r) \geq 1_{m_p}\}$. The problem each node has to solve at each iteration is a quadratic program, which we solve with Matlab’s `quadprog` function. Algorithm \[Alg:Giannakis\] was applied in [@Forero10-ConsensusBasedDistributedSVMs] to solve a version of  that takes into account possible violations of the linear separability of the data. Here, for simplicity, we just consider  and assume the data can be separated linearly.

Experimental Results {#Sec:Experiments}
====================

In this section, we show some simulation results of D-ADMM, and the algorithms from section \[Sec:RelatedAlgs\], applied to the problems just described. Of these applications, consensus is the one that assesses the algorithms the best, since neither the problem solution, nor the problem each node solves at each iteration, requires iterative algorithms or any kind of solver: they have closed-form solutions. We will, therefore, give more emphasis to consensus than to the other applications. We start with describing the network models.

  Name                                                        Parameters     Description
  ------------------------------------------------------- -- ------------ -- --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  Erdős-Rényi [@Erdos59RandomGraphs]                             $p$         Every pair of nodes $\{i,j\}$ is connected or not with probability $p$
  Watts-Strogatz [@Watts98CollectiveDynamicsSmallWorld]        $(n,p)$       First, it creates a lattice where every node is connected to $n$ nodes; then, it rewires every link with probability $p$. If link $\{i,j\}$ is to be rewired, it removes the link, and connects node $i$ or node $j$ (chosen with equal probability) to another node in the network, chosen uniformly.
  Barabasi-Albert [@Barabasi99EmergenceOfScaling]                 ——         It starts with one node. At each step, one node is added to the network by connecting it to $2$ existing nodes: the probability to connect it to node $p$ is proportional to $D_p$.
  Geometric [@Penrose]                                           $d$         It drops $P$ points, corresponding to the nodes of the network, randomly in a $[0,1]^2$ square; then, it connects nodes whose (Euclidean) distance is less than $d$.
  Lattice                                                         ——         Creates a lattice of dimensions $m\times n$; $m$ and $n$ are chosen to make the lattice as square as possible.

[@clc&gt;r&gt;rrrrrrr@]{} Number & Model & Parameters &\
& & &\
(r)[4-11]{} & & & 10 & 50 & 100 & 200 & 500 & 700 & 1000 & 2000\
1 & Erdős-Rényi & $0.25$ & & & & & & & &\

2 & Erdős-Rényi & $0.75$ & & & & & & & &\

3 & Watts-Strogatz & $(2,0.8)$ & & & & & & & &\

4 & Watts-Strogatz & $(4,0.6)$ & & & & & & & &\

5 & Barabasi-Albert & —— & & & & & & & &\

6 & Geometric & $0.2$& & & & & & & &\

7 & Lattice & ——& & & & & & & &\

[[**Network models.**]{}]{} We generated networks according to the models described in Table \[Tab:NetworkModels\]. All networks were generated using the igraph package in R [@igraph; @R]. Table \[Tab:NetworkParametersADNC\] shows the $56$ different networks we generated: given a model and fixed its parameters, we generated $8$ networks with different number of nodes, ranging from $10$ to $2000$; we used $7$ different combinations of models and parameters. All networks, except Lattice networks, were generated randomly, and we had to guarantee that all networks were connected (cf. Assumption \[Ass:Network\]). So, every time we would get a non-connected network, we would generate another network with slightly different parameters, changed in the direction to make the next network connected with greater probability. Therefore, there were some exceptions to the parameters described in Table \[Tab:NetworkParametersADNC\]. For example, network number $1$ with $10$ nodes was generated with $p=0.27$, instead of with $p=0.25$, or network number $6$ with $10$ nodes was generated with $d=0.36$, instead of with $d=0.2$.

Table \[Tab:NetworkParametersADNC\] also shows the average node degree and the number of colors used for each network. The average node degree can be as low as $3$, for example in the Lattice network with $50$ nodes, or as high as $1499$, in the Erdős-Rényi network with $2000$ nodes. The number of colors each network was colored with is also varied. Note that only the Lattice networks were colored with $2$ colors, and thus these are the only networks for which D-ADMM is proven to converge.

We use all networks of Table \[Tab:NetworkParametersADNC\] for the consensus problem, since we study that problem in more detail. To illustrate the algorithms for the remaining applications, we will only use the networks with $10$ and $50$ nodes.

[[**Choosing .**]{}]{} In section \[Sec:RelatedAlgs\] we reviewed the distributed algorithms in literature that can solve . All these algorithms, including D-ADMM, but excluding the gradient/subgradient method, are based on augmented Lagrangian duality. This means that, for each algorithm, we have to choose the parameter $\rho$ of the augmented Lagrangian (see, e.g., ). We are unaware of any scheme that chooses this parameter in an optimal way before the execution of the algorithm. There are, however, schemes for updating it during the execution of the algorithm [@Boyd11-ADMM; @He00ADMWithSlefAdaptativePenaltyForVIs], but it is not straightforward to implement them in a distributed scenario. We will thus adopt the following scheme for choosing $\rho$: given a fixed network, we execute several times all the algorithms that depend on $\rho$, each for a different $\rho$, chosen from the set $\{10^{-4},10^{-3},10^{-2},10^{-1},10^{0},10^{1},10^{2}\}$; then, we pick the $\rho$ that leads to the best result, i.e., the least number of communication steps. This scheme was adopted in all experiments.

[[**Consensus.**]{}]{} We now present simulation results for the consensus problem. The data was generated the following way: node $p$ generates, independently of the other nodes, one number $\theta_p$ from a Gaussian distribution with mean $10$ and standard deviation $100$. We chose such a large standard deviation to make sure that the nodes’ numbers differed significantly.

Figure \[Fig:ConsAllAlgs50nodes\] shows the results of the simulations of D-ADMM, the algorithms described in section \[Sec:RelatedAlgs\], and the additional algorithm from [@Olshevsky11-ConvergenceSpeedDistributedConsensusAveraging], which uses the approach . Note that [@Olshevsky11-ConvergenceSpeedDistributedConsensusAveraging] only solves the consensus problem and cannot be trivially generalized to solve . Figure \[Fig:ConsAllAlgs50nodes\] depicts the number of communication steps, which is equal to the number of iterations for the single-looped algorithms and to the number of inner iterations for the double-looped algorithms, for all the networks with $50$ nodes. All algorithms stopped after reaching a relative error of $10^{-2} \%$, i.e., $\|x^k - 1_P\theta^\star\|/(\sqrt{P}|\theta^\star|) \leq 10^{-4}$, or after reaching the maximum number of $1000$ communication steps. The optimal solution $\theta^\star$ had been computed before. It can been seen from the plot of Figure \[Fig:ConsAllAlgs50nodes\] that D-ADMM is the algorithm requiring the least communication steps to converge, i.e., to achieve a precision of $10^{-2}\%$, for all networks except the first two, where [@Olshevsky11-ConvergenceSpeedDistributedConsensusAveraging] was better. Figure \[Fig:ConsAllAlgs50nodes\] shows the error evolution along the communication steps for network $3$ (cf. with point $3$ of Figure \[Fig:ConsAllAlgs50nodes\]). It can be observed that D-ADMM required always the least number of communication steps to reach any accuracy between $10\%$ and $10^{-2}\%$.

We only consider the three best algorithms (D-ADMM, [@Zhu09-DistributedInNetworkChannelCoding], and [@Olshevsky11-ConvergenceSpeedDistributedConsensusAveraging]) for a more thorough analysis, in Figure \[Fig:ConsSize\]; here, we simulate these algorithms in all the networks of Table \[Tab:NetworkParametersADNC\]. In each of the plots of Figure \[Fig:ConsSize\] it is shown the number of communication steps required to achieve $10^{-2}\%$ of accuracy as a function of the number of nodes. It can be seen that D-ADMM always required less communication steps than [@Zhu09-DistributedInNetworkChannelCoding], except for the smallest network in almost all scenarios. D-ADMM also required less communication steps than [@Olshevsky11-ConvergenceSpeedDistributedConsensusAveraging], except for the Erdős-Rényi networks (Figures \[Fig:ConsSize\] and \[Fig:ConsSize\]). [@Olshevsky11-ConvergenceSpeedDistributedConsensusAveraging] was the only algorithm achieving the maximum number of $1000$ iterations, which occurred for the Geometric and Lattice networks (Figures \[Fig:ConsSize\] and \[Fig:ConsSize\]). We also note a curious fact from these plots: the number of communication steps is almost invariant to the size of the networks (for the considered range).

We thus conclude that D-ADMM is a competing algorithm for solving the consensus problem, since it requires, in general, less communication steps to converge than other state-of-the-art algorithm, and it has approximately the same behavior for very different networks.

[[**BPDN.**]{}]{} We also executed some experiments comparing D-ADMM and [@Zhu09-DistributedInNetworkChannelCoding] (Algorithm \[Alg:Giannakis\]), the two best algorithms, solving BPDN . Figure \[Fig:BPDN\] shows the results of such experiments, where we only used the networks with $50$ nodes. For generating data, we used the Sparco toolbox [@Sparco], namely problems with ID’s $7$ and $902$. Problem with ID $7$ contains a $600\times 2560$ matrix $A$ with Gaussian entries, and problem with ID $902$ contains a $200\times 1000$ partial DCT matrix. In both problems, we added a small amount of noise to the provided vector $b$, in order to make the application of  more realistic. The trade-off parameter $\beta$ was set to $1$ for problem ID $7$ and to $0.3$ for problem ID $902$. We ran the algorithms until they reached a relative error of $10^{-2}\%$ or the maximum number of iterations $1000$, only reached by [@Zhu09-DistributedInNetworkChannelCoding] in network $6$. At each iteration, each node solved problem  with GPSR [@GPSR].

The results for problems with ID’s $7$ and $902$ are shown respectively in Figures \[Fig:BPDN\] and \[Fig:BPDN\]. These show the number of communication steps for each network model. Again, each point in the plots is the best result for $\rho \in \{10^{-4},10^{-3},\ldots,10^2\}$. It can be observed that D-ADMM always required less communication steps than [@Zhu09-DistributedInNetworkChannelCoding]; however, there were some cases where both algorithms required almost the same communication steps, namely, for networks $2$ and $3$ of Figure \[Fig:BPDN\], and for networks $2$ and $4$ of Figure \[Fig:BPDN\]. The similarity of the shape of the curves of Figures \[Fig:BPDN\] and \[Fig:BPDN\] indicates that changing the network models causes more variability in the difficulty of the problem than the problem data.

[[**LASSO.**]{}]{} The results of our experiments for LASSO  are shown in Figure \[Fig:LASSO\]. Again, we just compared the best two algorithms: D-ADMM and [@Zhu09-DistributedInNetworkChannelCoding]. For the matrix $A$ and vector $b$, we used the same data as in BPDN: problems $7$ and $902$ from the Sparco toolbox. The bound for the noise was set to $\sigma=0.5$ for problem ID $7$, and to $\sigma=0.1$ for problem ID $902$. In all experiments, which were executed for the networks with $P=10$ nodes, the regularizing parameter in  was always $\delta = 10^{-3}$. With such a choice, both algorithms were able to achieve the relative error of $10^{-1}\%$ in all instances (the maximum number of iterations was never achieved). Note that the error was measured with respect to the solution of the original problem , which was obtained beforehand using the algorithm SPGL1 [@spgl1].

In Figure \[Fig:LASSO\], we can see that, in spite of solving a slightly ill-conditioned problem, D-ADMM and [@Zhu09-DistributedInNetworkChannelCoding] required considerably few communication steps to converge. Also note that the behavior of the algorithms is approximately uniform over all the types of networks. In contrast with BPDN, it was changing the data what produced different curves: the curves for ID $902$ are on top of the curves for ID $7$, which reveals that the former was harder to solve. For each problem, D-ADMM always took less communications steps than [@Zhu09-DistributedInNetworkChannelCoding] to converge.

[[**SVM.**]{}]{} Quite surprisingly, the problem that required the largest number of communication steps was SVM. For this problem, we used data from the UCI machine learning repository [@UCIMachineLearningRepository], namely the Iris dataset, which contains two clusters of points that are linearly separable. Each point has $n = 4$ features and there are $m = 100$ total points. We first solved  with the Matlab function `quadprog`. Next, we executed D-ADMM and [@Zhu09-DistributedInNetworkChannelCoding], which stopped after achieving a relative error of $10^{-2}\%$ or $10^4$ communication steps.

Figure \[Fig:SVM\] shows the results for the networks with $50$ nodes. Again, D-ADMM required always less communication steps than [@Zhu09-DistributedInNetworkChannelCoding] to converge; however, both algorithms required an amount of communications one order of magnitude larger than all the previous applications. And this happened with a relatively small problem. This might indicate that problem  is inherently difficult to solve in a distributed way.

Conclusions
===========

We proposed an algorithm for solving separable problems in a distributed way. This means that, given a connected network, each node has a cost function and a constraint set associated. While keeping their costs and constraints sets private all the time, all nodes cooperate in order to solve a global optimization problem, which minimizes the sum of all the nodes’ costs, and constrains the solution to be in the intersection of all sets. We showed how several problems from signal processing and control can be solved with the proposed algorithm, named D-ADMM, and presented extensive simulation results. These results show that D-ADMM requires almost always less communications to converge than any other state-of-the-art distributed optimization algorithm.

Proof of Theorem \[Teo:ConvergenceADMM\] {#App:ADMMProof}
========================================

First, note that the assumptions imply that strong duality holds for . This follows directly from [@Bertsekas:Convex Prop.6.4.2] and it means that the optimal cost of  and of  have the same (finite) value, which will be denoted by $p^\star$. It also means that  is solvable. Let then $(x_1^\star,x_2^\star,\lambda^\star)$ denote any primal-dual solution. We will use the notation $p^k = g_1(x_1^k) + g_2(x_2^k)$, $r^k = A_1x_1^k + A_2x_2^k$. Using only strong duality and our assumptions on $X_1$ and $X_2$, we can prove that $$\begin{aligned}
            &p^\star - p^{k+1} \leq {\lambda^{\star}}^\top r^{k+1}
            \label{Eq:ADMMProofIneq1}
            \\
          &p^{k+1} - p^{\star} \leq -(\lambda^{k+1})^{\top} r^{k+1}
                {\ifthenelse{\boolean{draft}}{}{\notag
                \\ 
                &\phantom{aaaaa}}}-\rho (A_2(x_2^{k+1} - x_2^k))^{\top}(A_2(x_2^{k+1} - x_2^\star) - r^{k+1})
            \label{Eq:ADMMProofIneq2}
            \\
          &V^{k+1} \leq V^k - \rho \|r^{k+1}\|^2 - \rho \|A_2(x_2^{k+1} - x_2^k)\|^2\,,
            \label{Eq:ADMMProofIneq3}
        \end{aligned}$$ where $V^k$ is the Lyapunov function $$V^k := \frac{1}{\rho}\|\lambda^k - \lambda^\star\|^2 + \rho\|A_2(x_2^k - x_2^\star)\|^2\,.$$ This is done in [@Boyd11-ADMM] and we omit the proof here. In the course of that proof, the following fact is used: if $\phi$ and $\Psi$ are convex functions and $\Psi$ is continuously differentiable, then $x^\star \in \arg\min_x\{\phi(x) + \Psi(x)\,:\, x \in X\}$ implies $x^\star \in \arg\min_x\{\phi(x) + \nabla \Psi(x^\star)^\top (x - x^\star)\,:\, x \in X\}$, for any closed convex set $X$. Using $\lambda^{k} = \lambda^{k-1} + \rho r^{k}$ and applying the previous fact to  and , we get, respectively, $$\begin{aligned}
            x_1^{k} &\in \arg\min_{x_1 \in X_1} g_1(x_1) + (\lambda^{k} - \rho A_2(x_2^{k} - x_2^{k-1}))^\top A_1 x_1
            \label{Eq:ADMMProofIneq4}
            \\
            x_2^{k} &\in \arg\min_{x_2 \in X_2} g_2(x_2) + {\lambda^{k}}^\top A_2 x_2\,.
            \label{Eq:ADMMProofIneq5}
        \end{aligned}$$

[[**Limit points of are primal optimal.**]{}]{} From , we see that $0 \leq V^k \leq V^0$ for any $k$, and thus $\{\lambda^k\}$ and $\{A_2x_2^k\}$ are bounded, showing that they have limit points. Since $A_2$ has full column-rank, $\{x_2^k\}$ is also bounded and thus has limit points. To show existence of limit points of $\{x_1^k\}$, we first have to show that $r^k \xrightarrow{} 0$. Iterating  we get $$\rho \sum_{k=0}^\infty (\|r^{k+1}\|^2 + \|A_2(x_2^{k+1} - x_2^k)\|^2) \leq V^0\,,$$ from which we conclude $r^k \xrightarrow{} 0$ and $A_2(x_2^{k+1} - x_2^k) \xrightarrow{} 0$ as $k \xrightarrow{} 0$. Now, $r^k = A_1x_1^k + A_2x_2^k \xrightarrow{} 0$, together with the boundedness of $\{A_2x_2^k\}$, shows that $\{A_1x_1^k\}$ is bounded and thus also $\{x_1^k\}$, because $A_1$ has full column-rank. Therefore, $\{x_1^k\}$ has limit points.

We showed the existence of limit points of $\{(x_1^k,x_2^k)\}$; now, we show their optimality. In fact, taking $k \xrightarrow{} \infty$ in  and , using the boundedness of $\{A_2x_2^k\}$ and $\{\lambda^k\}$, and the facts that $r^k \xrightarrow{} 0$ and $A_2(x_2^{k+1} - x_2^k) \xrightarrow{} 0$, we conclude that $p^k = g_1(x_1^k) + g_2(x_2^k) \xrightarrow{} p^\star$.

[[**Limit points of are dual optimal.**]{}]{} We had seen that $\{\lambda^k\}$ has limit points because it is bounded. Let $\bar{\lambda}$ be such limit point, and let $\mathcal{K}$ be the set of indices such that $\{\lambda^k\}_{k \in \mathcal{K}} \xrightarrow{} \bar{\lambda}$ and also that $\{(x_1^k,x_2^k)\}_{k \in \mathcal{K}}\xrightarrow{}(\bar{x}_1,\bar{x}_2)$, where $(\bar{x}_1,\bar{x}_2)$ is a limit point of $\{(x_1^k,x_2^k)\}$. Now, define $\hat{\lambda}^k = \lambda^k - \rho A_2(x_2^k - x_2^{k-1})$. Since $A_2(x_2^k - x_2^{k-1}) \xrightarrow{} 0$, we have $\lim_{k\xrightarrow{} \infty}\{\hat{\lambda}^k\}_{k \in \mathcal{K}} = \lim_{k\xrightarrow{} \infty}\{\lambda^k\}_{k \in \mathcal{K}} = \bar{\lambda}$. Now, using the definition of $G_1$ and $G_2$, and -, we see that $$\begin{aligned}
              G_1(\hat{\lambda}^{k}) &= \inf_{x_1 \in X_1} g_1(x_1) + (\hat{\lambda}^{k})^\top A_1x_1
              \notag
            \\
            &= 
              g_1(x_1^{k}) + (\hat{\lambda}^{k})^\top A_1x_1^{k}
            \label{Eq:ADMMProofDualProb1}
            \\
            &\leq
              g_1(x_1) + (\hat{\lambda}^{k})^\top A_1x_1\,, \qquad \forall_{x_1 \in X_1}\,,
            \label{Eq:ADMMProofDualProb2}
            \\              G_2(\lambda^{k}) &= \inf_{x_2 \in X_2} g_2(x_2) + (\lambda^{k})^\top A_2x_2
              \notag
            \\
            &=
              g_2(x_2^{k}) + (\lambda^{k})^\top A_2x_2^{k}
            \label{Eq:ADMMProofDualProb3}
            \\
            &\leq
              g_2(x_2) + (\lambda^{k})^\top A_2x_2\,, \qquad \forall_{x_2 \in X_2}\,.
            \label{Eq:ADMMProofDualProb4}
        \end{aligned}$$ Adding  and  and taking the limit $k\xrightarrow{} +\infty$ ($k \in \mathcal{K}$), where the second-to-last equality follows from primal optimality (and feasibility) of $(\bar{x}_1,\bar{x}_2)$, and the last equality follows from strong duality. Adding  and  and taking the limit $k\xrightarrow{} +\infty$ ($k \in \mathcal{K}$), $$\lim_{\mathcal{K}\ni k\xrightarrow{}\infty} (G_1(\hat{\lambda}^{k}) + G_2(\lambda^{k}))
              \leq
              \sum_{i=1}^2 (g_i(x_i) + \bar{\lambda}^\top A_ix_i) \,,$$ for all $x_i \in X_i$, $i=1,2$. In particular, taking the infimum on the right-hand side: $$\label{Eq:ADMMProofDualProb6}
             \lim_{\mathcal{K}\ni k\xrightarrow{}\infty} (G_1(\hat{\lambda}^{k}) + G_2(\lambda^{k}))
              \leq
           G_1(\bar{\lambda}) + G_2(\bar{\lambda})\,.$$ From  and  we conclude that $G_1(\bar{\lambda}) + G_2(\bar{\lambda}) \geq L(\lambda^\star)$, showing that $\bar{\lambda}$ is dual optimal.

[[**The sequence converges.**]{}]{} We have seen that $\{(x_1^k,x_2^k,\lambda^k)\}$ has limit points and they are primal-dual optimal. Now, we will see that this sequence has only one limit point, and thus converges. Let $(\bar{x}_1,\bar{x}_2,\bar{\lambda})$ be a limit point of $\{(x_1^k,x_2^k,\lambda^k)\}$. Since it is primal-dual optimal, - hold with $(x_1^\star,x_2^\star,\lambda^\star)$ replaced by $(\bar{x}_1,\bar{x}_2,\bar{\lambda})$. We had seen that $V^k$ was convergent because it was bounded and nonincreasing. With the previous replacement, we now see that its limit is $0$. Therefore, $\lambda^k \xrightarrow{} \bar{\lambda}$ and $x_2^k \xrightarrow{} \bar{x}_2$ (because $A_2$ has full column-rank). Since $r^k = A_1(x_1 - \bar{x}_1) + A_2(x_2 - \bar{x}_2) \xrightarrow{} 0$, we also have $x_1^k \xrightarrow{} \bar{x}_1$ (because $A_1$ has full-column rank), i.e., $\{(x_1^k,x_2^k,\lambda^k)\}$ converges.

Problem For Each Node in LASSO {#App:ProbLasso}
==============================

When LASSO  is solved in the column partition case, node $p$ has to solve, in each iteration, $$\label{Eq:AppLassoProb}
    \underset{\lambda}{\text{minimize}} \,\,\, g_p(\lambda) + v^\top \lambda + c\|\lambda\|^2\,,$$ for some vector $v$ and scalar $c$, and with $g_p(\lambda) =  \frac{1}{P}(b^\top \lambda + \sigma\|\lambda\|) - \inf_{x_p}(\|x_p\|_1 + (A_p^\top \lambda)^\top x_p + \frac{\delta}{2}\|x_p\|^2)$. Define $$\label{Eq:AppLasso1}
    \phi_p(\lambda) = -\inf_{x_p} \,\, \|x_p\|_1 + (A_p^\top \lambda)^\top x_p + \frac{\delta}{2}\|x_p\|^2\,,$$ and let $x_p(\lambda)$ denote the solution of the optimization problem in  for a given $\lambda$. Since the objective in that problem is strictly convex, $\phi_p$ is differentiable and its gradient is given by $\nabla \phi_p(\lambda) = -A_px_p(\lambda)$. Furthermore, each component of $x_p(\lambda)$ can be found in closed-form: $-(r_i + 1)/\delta$ if $r_i < -1$, $-(r_i - 1)/\delta$ if $r_i > 1$, and $0$ otherwise, where $r_i$ denotes the $i$th row of $A_p^\top \lambda$. Problem  then becomes equivalent to $$\underset{\lambda}{\text{minimize}} \,\,\, \phi_p(\lambda) + \frac{1}{P}b^\top\lambda + \frac{\sigma}{P}\|\lambda\| + v^\top \lambda + c\|\lambda\|^2\,,$$ or introducing an epigraph variable $t$, $$\label{Eq:AppLasso2}
        \begin{array}{ll}
            \underset{x,t}{\text{minimize}} & \phi_p(\lambda) + \frac{1}{P}b^\top\lambda + \frac{\sigma}{P}t + v^\top \lambda + c\|\lambda\|^2 \\
            \text{subject to} & \|\lambda\| \leq t\,.
        \end{array}$$ Problem  can be solved using a projected gradient method. We will use FISTA [@BeckTeboulleFISTA §4], whose complexity is $O(1/k^2)$. FISTA solves $$\label{Eq:FistaProb}
            \underset{\lambda}{\text{minimize}} \,\,\, f(\lambda) + g(\lambda)\,,$$ where $f$ is a convex continuously differentiable function with a Lipschitz continuous gradient $\nabla f$ with constant $L$ (i.e., $\|\nabla f(\eta) - \nabla f(\lambda)\| \leq L\|\eta - \lambda\|$ for any $\lambda,\eta$), and $g$ is a convex function. FISTA consists of iterating $
        \lambda^{k+1} = \bigl[\lambda^k - \frac{1}{L} \nabla f(\lambda^k))\bigr]^+\,,
    $ where $$\label{Eq:FistaProjection}
        [\eta]^+ = \arg\min_z \,\, g(z) + \frac{L}{2} \|z - \eta\|^2\,.$$ We set $f(\lambda)$ to be the objective of , and $h(\lambda)$ to be the indicator function of the set $\mathcal{S} := \{(\lambda,t)\,:\, \|\lambda\|\leq t\}$. In this case, becomes the projection onto $\mathcal{S}$, which can be shown to be $$[(\lambda,t)]^+_{\mathcal{S}} = 
        \left\{
            \begin{array}{lcl}
                (\lambda,t) &,& t\geq \|\lambda\| \\
                (0,0) &,& t \leq -\|\lambda\| \\
                \frac{t+\|\lambda\|}{2} (\frac{\lambda}{\|\lambda\|},1) &,& -\|x\|<t< \|x\|\,.
            \end{array}
        \right.$$ The Lipschitz constant of $\nabla f$ can also be shown to be $\sigma_{\max}^2(A_p)/\delta + 2c$, where $\sigma_{\max}(A_p)$ is the largest singular value of $A_p$.

[^1]: This work was supported by the following grants from Fundação para a Ciência e Tecnologia (FCT): CMU-PT/SIA/0026/2009, PTDC/EEA-ACR/73749/2006, PEst-OE/EEI/LA0009/2011, and SFRH/BD/33520/2008 (through the Carnegie Mellon/Portugal Program managed by ICTI).

[^2]: At our abstraction level, we are overlooking packet collision and other medium-access problems that may potentially destroy synchronism.
