Stochastic Proximal Gradient Consensus Over Random Networks

I Introduction

This problem has found applications in various domains such as distributed consensus , distributed and parallel machine learning and distributed signal processing ; see for a recent survey. The key research question is: How to compute an optimal solution of (1), through a distributed process where each agent only utilizes local gradient information about the objective.

Let each agent $i$ keep a local copy of $y$, say $y_{i}$. The well-known distributed subgradient (DSG) method is given by

where $r$ denotes the iteration counter; $d^{r}_{i}\in\partial f_{i}(y^{r}_{i})$ denotes a subgradient of the local function $f_{i}$ evaluated at $y^{r}_{i}$; $w^{r}_{ij}\geq 0$ denotes the weight for the link $e_{ij}\in\mathcal{E}$ at iteration $r$; and $\gamma^{r}>0$ denotes some stepsize parameter. Let $\bar{y}^{r}_{i}:=\frac{1}{r}\sum_{t=1}^{r}y_{i}^{t}$.

The convergence of the DSG iteration (2) was first analyzed in by Nedić and Ozdaglar. It was shown that if the subdifferential is bounded, and that the weights $\{w_{ij}\}$ and the graph $\mathcal{G}$ satisfy certain regularity assumptions, then each $\bar{y}^{r}_{i}$ converges to a neighborhood of the optimal solution (resp. the exact optimal solution) if $\gamma^{r}$ is a constant (resp. a diminishing sequence). As a special case, when $f(x)\equiv 0$ (only the consensus among the agents is sought for), then the convergence of the iteration (2) was first studied by Tsitsiklis . The DSG iteration has been extended to scenarios where there is a local constraint for each agent , or the messages exchanged among the agents are quantized , or the communication among the agents is noisy . Also see for other related methods for solving (1).

The rate of convergence analysis of the DSG-type method has been a central research issue. In its most general form, it is known that when appropriate diminishing stepsizes are chosen, DSG converges with a rate of $\mathcal{O}(\ln(r)/\sqrt{r})$ in terms of the differences between the local objective functions and the optimal objective function , for both static and time-varying networks. Duchi et al. propose a distributed dual-averaging algorithm and show that it converges with a rate of $\mathcal{O}(\ln(r)/\sqrt{r})$. Jakovetic et al. show that when the objective has Lipschitz continuous and bounded gradient, and when the graph is static, it is possible to accelerate the DSG to achieve an $\mathcal{O}(1/r^{2})$ rate, but at the expense of solving more complicated subproblems, each of which involves multiple rounds of communication and computation. If only simple computation/communication steps are performed, the rate becomes $\mathcal{O}(\ln(r)/r)$. A related acceleration scheme has also been proposed in , which further works for time-varying $B$-connected graphs A $B$-connected graph is a time-varying graph in which at each iteration the graph is not necessarily connected, but the union of the graphs across every $B>0$ consecutive iterations is connected.. Under the smoothness assumption on $f$, Shi et al. propose an interesting algorithm called EXTRA, which adds certain error-correction terms to the DSG (2). By adding such correction, EXTRA uses constant stepsize and achieves an $\mathcal{O}(1/r)$ rate for smooth convex problem and linear convergence for certain smooth strongly convex problems. This method has also been generalized to solve nonsmooth problems , but both algorithms in can only work for static networks. Other recent developments can be found in and the references therein.

Another popular approach for distributed optimization is to use the alternating direction method of multipliers (ADMM) . Applying the ADMM to distributed optimization has been first suggested in , and subsequently popularized in . The $\mathcal{O}(1/r)$ sublinear rate of convergence for decentralized consensus ADMM (C-ADMM) has been shown by Wei and Ozdaglar , where it is assumed that the underlying graph is generated according to certain stochastic mechanism. When the problem is smooth, the linear convergence of C-ADMM is shown in . Recently a broadcast based C-ADMM has been proposed in . However the C-ADMM usually requires solving local optimization problems exactly (cf. ), which can be expensive in certain applications. This requirement has been relaxed by two recent works and . In particular, Chang et al. develop an inexact C-ADMM (IC-ADMM) algorithm which uses a simple (proximal) gradient step at each ADMM iteration. Ling et al. also propose to replace the exact minimization by certain proximal gradient steps. While we are finalizing the paper, we were made aware of an independent work that also proposes linearlized ADMM method for consensus composite optimization. In particular, for a static network, convergence rates of $\mathcal{O}(1/{r})$ and $\mathcal{O}(1/\sqrt{r})$ are shown for certain deterministic and stochastic linearlized distributed ADMM. Recently, Hong et al. show that the ADMM-based method (with exact or inexact update) can be used to solve certain nonconvex global consensus problem, with a convergence rate of $\mathcal{O}(1/\sqrt{r})$.

There has been a few works that design distributed optimization algorithms in the primal-dual perspective. For example, propose random coordinate primal-dual algorithms, with possible applications in distributed and asynchronous optimization. However no convergence rate has been provided. In , the authors propose an augmented Lagrangian method based algorithm for distributed optimization and analyzed its linear convergence, but the algorithm and analysis only works for smooth and strongly convex problems. Further, the algorithm has double loops, and requires some global knowledge about the objective function and the underlying graph. These requirements can be restrictive in practical applications. Reference develops a primal-dual algorithm in which each agent is updated by performing local stochastic averaging gradients.

Below we provide a high level comparison of the DSG-based and ADMM-based algorithms.

(Problem types) The DSG can solve convex problems with only subgradient information about the objective, while to our best knowledge the ADMM does not directly work for this case.

(Gradient Information) The DSG only needs (stochastic) subgradients of the objective , while the ADMM usually requires subproblems to have some nice structures so that they can be solved in closed-form .

(Convergence rates) When the objective function $f$ has certain additional structures (e.g., smooth or a smooth plus a simple nonsmooth function), the distributed ADMM generally converges faster in practice (which has a convergence rate of $\mathcal{O}(1/r)$) than its DSG counterpart (which has a convergence rate of $\mathcal{O}(\ln(r)/\sqrt{r})$). Theoretically, it is possible to modify the iteration of the DSG algorithm to improve its rate to $\mathcal{O}(1/r)$; see recent developments in .

(Network structures) The DSG generally works when the underlying network is time-varying and follows the so-called $B$-connected structure . However the ADMM-based method only works for static network, except for the recent variants proposed in , both of which work for certain randomized networks.

I-B Contribution of This Work

In this work, we consider the following structured version of the global consensus problem (1)

We propose an ADMM based method, named dynamic stochastic proximal-gradient consensus (DySPGC), that has the following key features:

• When only an unbiased estimate of each $\nabla g_{i}$ is known, the algorithm converges with a rate $\mathcal{O}(1/\sqrt{r})$;

• When the exact $\nabla g_{i}$ is known, the rate becomes $\mathcal{O}(1/r)$;

• The algorithm works for both the static and certain random time-varying networks.

What is more interesting is our insight on the connection between the C-ADMM-type methods and a few DSG-type methods. In particular, we show that the EXTRA/PG-EXTRA , despite being posed as error-corrected DSGs, can be viewed as special cases of the proposed DySPGC (for static network with symmetric weights and exact gradients). This observation explains the relative fast practical convergence performance of these two algorithms compared with the DSG [for structured problems (3)]. Further, we also establish a close connection between the DSG (2) and the proposed DySPGC. Additionally our method generalizes other distributed ADMM-type methods such as the DLM and the IC-ADMM .

II System Model

We make the following blanket assumptions for (3).

Let $g(y):=\sum_{i=1}^{N}g_{i}(y)$, and $h(y):=\sum_{i=1}^{N}h_{i}(y)$. The $h_{i}$’s prox operators, defined below

Each $\nabla g_{i}$ is Lipschitz continuous (with constant $P_{i}>0$)

As have been mentioned in the introduction, we consider a collection of $N$ agents defined over a connected undirected graph {\mbox{\mathcal{G}}}=\{\mathcal{V},\mathcal{E}\}, with $|\mathcal{V}|=N$ vertices and $|\mathcal{E}|=E$ edges. Define a companion symmetric directed graph given by {\mbox{\mathcal{G}}}_{d}=\{\mathcal{V},\mathcal{A},W\}, where $\mathcal{A}$ is a set of directed arcs with $|\mathcal{A}|=2E$, and for every edge in $\mathcal{E}$ which connects nodes $i,j$, we have both $e_{ij},e_{ji}\in\mathcal{A}$. Note that using a companion graph to represent the original graph $\mathcal{G}$ is conventional in the consensus ADMM literature; see e.g., and . It helps to simplify the definition of the consensus constraint (to be provided shortly).

Let us use $\mathcal{N}_{i}$ to denote the neighborhood of node $i$, i.e.,

$W$ is a row stochastic matrix, i.e., $\{W[i,j]\geq 0\}$, $\sum_{j}W[i,j]=1,\;\forall~{}i$;

The diagonal elements of $W$ are all positive, and its off-diagonal elements all satisfy

Later we will provide explicit expressions for $W$.

Consider an equivalent reformulation of problem (3) (equivalent when $\mathcal{G}$ is connected)

II-B Randomly Time-Varying Graph Structure

We assume that the edges of the graph $\mathcal{G}$ are activated according to certain randomly time-varying patterns. To describe such random pattern, at a given time $r$, define a new graph {\mbox{\mathcal{G}}}^{r}=\{{\mbox{\mathcal{V}}}^{r},\mathcal{E}^{r}\}, and its companion graph {\mbox{\mathcal{G}}}_{d}^{r}=\{{\mbox{\mathcal{V}}}^{r},{\mbox{\mathcal{A}}}^{r},W^{r}\} where {\mbox{\mathcal{V}}}^{r}\subseteq{\mbox{\mathcal{V}}}, $\mathcal{E}^{r}\subseteq\mathcal{E}$ and {\mbox{\mathcal{A}}}^{r}\subseteq{\mbox{\mathcal{A}}}, and each weight matrix $W^{r}$ is a stochastic matrix satisfying (7). Again {\mbox{\mathcal{G}}}_{d}^{r} is symmetric, meaning if $e$ connects nodes $i$ and $j$ with $e\in\mathcal{E}^{r}$, then e_{ij},e_{ji}\in{\mbox{\mathcal{A}}}^{r}. The precise specification of the random graphs \{{\mbox{\mathcal{G}}}_{d}^{r}\} and \{{\mbox{\mathcal{G}}}^{r}\} is given below .

(Randomly Activated Graph) At each time $r$, each link pair (i,j),(j,i)\in{\mbox{\mathcal{A}}} has a probability $p_{ij}=p_{ji}\in(0,1]$ of being active. The set of active nodes {\mbox{\mathcal{V}}}^{r} is given by:

Effectively at each time $r$ a node i\in{\mbox{\mathcal{V}}} has a probability $\alpha_{i}>0$ of being active, while such $\alpha_{i}$ is a function of $\{p_{ij}\mid j\in\mathcal{N}_{i}\}$. Let us collect these probabilities and define

where $\mbox{diag}\{\alpha_{i}\}$ represents a diagonal matrix whose diagonal entries are the elements in the set $\{\alpha_{i}\}$. Further, assume that {\mbox{\mathcal{G}}}_{d} is strongly connected, and realizations of the graphs {\mbox{\mathcal{G}}}_{d}^{r} and {\mbox{\mathcal{G}}}_{d}^{t} are independent and identically distributed across all $r\neq t$. $\blacksquare$

In practice, the random network pattern can be used to model communication and/or node failures . It is the stochastic variant of the so-called $B$-strongly connected network which has been widely considered in the literature, under very different context . The connection between such randomly generated graph and popular communication protocols such as the gossip protocol and asynchronous protocols has been explored in . Note the graph $\mathcal{G}$ is required to be connected, but {\mbox{\mathcal{G}}}^{r}’s are not necessarily so. At a given iteration $r$, we can define the neighborhood $\mathcal{N}_{i}^{r}$ for each node $i$ similarly as in (6), and define the matrices $A^{r}$ and $B^{r}$ similarly as in (13), making all quantities conforming to the instantaneous graph structure.

II-C The Gradient Information

Define the gradient of the smooth part of the objective as $G(x):=[\nabla g_{1}(x_{1});\cdots;\nabla g_{N}(x_{N})]$. In this work, we will consider situations in which only an estimate of $\nabla g_{i}(x_{i})$, denoted by $\widetilde{g}_{i}(x_{i},\xi_{i})$, is available for each agent $i$. In this case, the estimate $\widetilde{g}_{i}(x_{i},\xi_{i})$ will satisfy the following

where each $\xi_{i}$ is a random variable following an unknown distribution, and $\xi_{i}$, $\xi_{j}$ are not necessarily independent for any $i\neq j$. Further, when time is involved (cf. Section II-B), we will assume $\xi_{i}$ to be independent over time. Each $\widetilde{g}_{i}(x_{i},\xi_{i})$ is assumed to be a measurable function; the constant $\sigma^{2}$ represents the maximum expected deviation of the gradient estimate.

III The Proposed Algorithms

Our proposed algorithms are based on the ADMM. To describe the algorithm in its general form, let us first define a vector of positive penalty constants \rho:=\{\rho_{ij}>0\mid e_{ij}\in{\mbox{\mathcal{A}}}\}, i.e., each $\rho_{ij}$ corresponds to a link variable $z_{ij}$. For a given graph {\mbox{\mathcal{G}}}_{d}, we can construct a diagonal matrix $\Gamma\succeq 0$ by

Using the above definition, let us write the augmented Lagrangian of (P):

To proceed, we need the following definitions. For each i\in{\mbox{\mathcal{V}}} and some $\omega_{i}\geq 0$, define

where the latter matrix is a block diagonal matrix with diagonal blocks being $\Omega_{1},\cdots,\Omega_{N}$. Define the following matrices

It can be verified that $\frac{1}{2}M_{-}M^{T}_{-}$ and $\frac{1}{2}M_{+}M^{T}_{+}$ represent the signed and signless graph Laplacian matrices, respectively (see, e.g., [31, Section II] for detailed discussion on these matrices).

To illustrate various quantities related to the graphs, let us consider a simple graph with $3$ nodes and two edges connecting nodes $\{1,2\}$ and nodes $\{2,3\}$. Suppose that $M=1$. In this case, $\mathcal{A}=\{(1,2),(2,3),(2,1),(3,2)\}$. Let us order the links as $(1,2),(2,3),(2,1),(3,2)$, then the matrices $A_{1}$ and $A_{2}$ are given below

The matrices $M_{+}$ and $M_{-}$ are given by

Define h^{r+1}(x):=\sum_{i\in{\mbox{\mathcal{V}}}^{r+1}}h_{i}(x_{i}). Let $\{\eta^{r}\geq 0\}$ denote a sequence of iteration-dependent parameters, whose values will be given shortly.

Using these definitions, we present in the table below the proposed algorithm in its general form, named the dynamic stochastic proximal-gradient consensus (DySPGC) algorithm.

Algorithm 1. DySPGC Over Random Graphs At iteration , select $\lambda^{0},z^{0},x^{0}$ such that $B^{T}\lambda^{0}=0,\;z^{0}=\frac{1}{2}M^{T}_{+}x^{0}.$ At each iteration $r+1$, update the variable blocks by: $\displaystyle x^{r+1}$ \displaystyle=\arg\min_{x}\;\bigg{\{}\left\langle\widetilde{G}^{r+1}(x^{r},\xi^{r+1}),x-x^{r}\right\rangle+h^{r+1}(x) $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{1}{2}\left\|A^{r+1}x+B^{r+1}z^{r}+\Gamma^{-1}{\lambda^{r}}\right\|_{\Gamma}^{2}$ \displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{1}{2}\|x-x^{r}\|^{2}_{\Omega+{\eta^{r+1}}I_{MN}}\bigg{\}} (20a) $\displaystyle x^{r+1}_{i}$ \displaystyle=x^{r}_{i},\quad\mbox{if}~{}i\notin{\mbox{\mathcal{V}}}^{r+1} (20b) $\displaystyle z^{r+1}$ $\displaystyle=\arg\min_{z}\;\frac{1}{2}\left\|A^{r+1}x^{r+1}+B^{r+1}z+\Gamma^{-1}{\lambda^{r}}\right\|_{\Gamma}^{2}$ (20c) $\displaystyle z^{r+1}_{ij}$ \displaystyle=z_{ij}^{r},\quad\mbox{if}~{}e_{ij}\notin{\mbox{\mathcal{A}}}^{r+1} (20d) $\displaystyle\lambda^{r+1}$ $\displaystyle=\lambda^{r}+\Gamma\left(A^{r+1}x^{r+1}+B^{r+1}z^{r+1}\right)$ (20e)

Let us make a few comments about DySPGC. First, the penalty parameter used for the $x$-update for the proximal term $\|x-x^{r}\|^{2}$ is given by $\Omega+\eta^{r+1}I_{MN}$. Here $\Omega$ is a fixed constant matrix defined in (18); the iteration-dependent parameter $\eta^{r+1}$, when being chosen as an appropriate increasing sequence (to be specified in Theorem IV.2), is used to deal with the stochasticity in the gradient. Second, when the gradients are precisely known, we can set $\eta^{r+1}=0$ for all $r$, in which case the $x$-update rule (20a) becomes

where $G^{r+1}(x^{r})$ is defined similarly as $\widetilde{G}^{r+1}(x^{r},\xi^{r+1})$ (with inexact gradients replaced by the exact gradients).

When we assume that the graph is static and the exact gradients are known, i.e., {\mbox{\mathcal{G}}}_{d}^{r}={\mbox{\mathcal{G}}}_{d} and $\widetilde{G}^{r+1}(x^{r},\xi^{r+1})=G(x^{r})$ for all $r$, then the DySPGC reduces to a simplified version named the proximal gradient consensus (PGC) algorithm (see Algorithm 2).

Let us compare Algorithms 1 and 2 with some existing methods and pinpoint the main differences. First, PGC is a proximal version of the conventional C-ADMM , where we have used the second order approximation $\langle G(x^{r}),x-x^{r}\rangle+\frac{1}{2}\|x-x^{r}\|^{2}_{\Omega}$ of the smooth function $g(x)$ in (21a) in the $x$-update, rather than exactly minimizing the augmented Lagrangian (as has been done in ). Moreover, a matrix penalty $\Gamma$ is used instead of a scalar one; the latter has been popular in the existing ADMM-based methods. Later we will show that by using a matrix penalty, the parameters $\rho_{ij}$’s can be chosen by only using local information, making the algorithm better suited to distributed implementation. Second, the DySPGC is a stochastic version of the algorithms proposed in , where we have used an iteration-dependent stochastic second order approximation

in the $x$-update step, rather than exactly minimizing the augmented Lagrangian

Detailed comparison with existing algorithms will be provided in Section V.

Algorithm 2. PGC Over Static Graphs At iteration , select $\lambda^{0},z^{0},x^{0}$ such that $B^{T}\lambda^{0}=0,\;z^{0}=\frac{1}{2}M^{T}_{+}x^{0}.$ At each iteration $r+1$, update the variable blocks by: $\displaystyle x^{r+1}$ \displaystyle=\arg\min_{x}\;\!\!\bigg{\{}\!\!\!\left\langle G(x^{r}),x-x^{r}\right\rangle\!+\!h(x)\!+\!\langle\lambda^{r},Ax+Bz^{r}\rangle \displaystyle~{}~{}~{}~{}~{}\quad+\frac{1}{2}\left\|Ax+Bz^{r}\right\|_{\Gamma}^{2}+\frac{1}{2}\|x-x^{r}\|^{2}_{\Omega}\bigg{\}} (21a) $\displaystyle z^{r+1}$ $\displaystyle=\arg\min_{z}\;\frac{1}{2}\left\|Ax^{r+1}+Bz+\Gamma^{-1}{\lambda^{r}}\right\|_{\Gamma}^{2}$ (21b) $\displaystyle\lambda^{r+1}$ $\displaystyle=\lambda^{r}+\Gamma\left(Ax^{r+1}+Bz^{r+1}\right)$ (21c)

III-B Distributed Implementation

Both algorithms proposed in the previous section can be implemented in a distributed manner, in which the information needed for updating each variable can be obtained from its immediate neighbors. To see this, note that in the original formulation (8) each node $i$ is only coupled with its neighboring links $\{e_{ij},e_{ji}\}_{j\in\mathcal{N}_{i}}$, and each link pair e_{ij},e_{ji}\in{\mbox{\mathcal{A}}} is only related to its two neighboring nodes \{i,j\}\in{\mbox{\mathcal{V}}}. Below we illustrate the distributed implementation of the PGC algorithm, as it takes a simple form.

To write the algorithm compactly, define the stepsize parameter $\beta_{i}$ as [with $\widehat{\rho}_{ij}:=1/2(\rho_{ij}+\rho_{ji})$]

Clearly $W$ is a row stochastic matrix satisfying the conditions in (7). However, generally $W$ constructed in this way is neither symmetric nor doubly stochastic, except when all $\beta_{i}$’s are identical. We note that each entry of $W[i,j]$ is directly related to how agent $i$ will combine agent $j$’s information (this point will be made clear shortly).

Surprisingly, Algorithm 2 admits a compact single-variable characterization, as we show in the following result.

The iteration (21a) – (21c) of Algorithm 2 (PGC) has the following compact characterization:

The proof for the above claim is relegated to Appendix -A.

Let us take a closer look at iterations (III.1) and (III.1). First, note that ${1}/{\beta_{i}}$ (or equivalently $\Upsilon^{-1}$) can be viewed as the stepsize for updating along the gradient direction. Second, for each node $i$, it is clear the penalty parameter $\hat{\rho}_{ij}=(\rho_{ij}+\rho_{ji})/2$ (or the $(i,j)$’s entry of the weight matrix $W$) is the weight that specifies how the user $j$’s information (i.e., $x^{r}_{j}$ and $x^{r-1}_{j}$) is combined with the user $i$’s information at each iteration. The larger the value of $\hat{\rho}_{ij}$ (or $W[i,j]$), the more emphasis that agent $i$ will put on agent $j$’s information.

Then we comment on how (III.1) and (III.1) can be carried out in practice.

If $h\equiv 0$, then $\zeta^{r}=0,\;\forall~{}r$. To perform (III.1) each agent $i$ needs its past iterate ($x^{r}_{i}$, $x^{r-1}_{i}$) the stepsize parameter $1/\beta_{i}$, the gradients ($\nabla g_{i}(x_{i}^{r}),\;\nabla g_{i}(x_{i}^{r-1})$), as well as the weighted sum of $x^{r}$ over its neighbors at the current and past iterations, i.e., $\sum_{j\in\mathcal{N}_{i}}\widehat{\rho}_{ij}x^{r}_{j}$ and $\sum_{j\in\mathcal{N}_{i}}\widehat{\rho}_{ij}x^{r-1}_{j}$, respectively. Also it is clear that the algorithm can be implemented in a fully distributed manner, since at iteration $r+1$, a given agent $i$ only communicates with its neighbors $\mathcal{N}_{i}$.

When $h\neq 0$, iterations (III.1) and (III.1) can be implemented in the following manner. Assume that $x^{0}=x^{-1}=0$ and $\zeta^{0}=0$ for initialization. Then according to (III.1) we have $x^{1}_{i}+\frac{1}{\beta_{i}}\zeta^{1}_{i}=0$, so $x^{1}_{i}$ and $\zeta^{1}_{i}$ can be obtained by solving the following problem

​​Then one can compute $c^{2}_{i}$ according to (III.1). To obtain $(x^{r+1}_{i},\zeta^{r+1}_{i})$, $r\geq 1$, suppose $\zeta^{r}_{i}$ and $c^{r+1}_{i}$ are available, then according to (III.1), we have

Finding $x^{r+1}_{i}$ is equivalent to solving the following

Once $x^{r+1}_{i}$ is obtained, we can compute $\zeta^{r+1}_{i}$ by

Clearly, as long as problem (30) can be solved easily, iteration (III.1) can be implemented efficiently in a distributed manner.

IV Convergence Analysis

We begin analyzing the (rate of) convergence of the proposed methods. Let us define a diagonal matrix of Lipschitz constants by

Let $w:=[x;z;\lambda]$ denote the vector of primal-dual iterates generated by PGC/DySPGC, and let $w^{*}:=[x^{*};z^{*};\lambda^{*}]$ denote a vector of optimal primal-dual solutions for problem (P). Our main convergence results are summarized in Table I. All the proofs of this section are relegated to the Appendix.

For the PGC algorithm which use static graph and exact gradients, we have the following convergence result.

Suppose that Assumption 1 holds, {\mbox{\mathcal{G}}}^{r}={\mbox{\mathcal{G}}} for all $r$, and $\mathcal{G}$ is connected. Then the following hold:

(a) Algorithm 2 converges to a primal-dual optimal solution of (P) if the following condition is satisfied

(b) Assume that $\mbox{dom}(h)$ is a bounded set, i.e., there exists a finite $C>0$ such that

Suppose that $w^{r}:=[x^{r};z^{r};\lambda^{r}]$ is generated by Algorithm 2 and the stepsize matrix satisfies

and $d_{\lambda}(\rho):=\sup_{\lambda\in\mathcal{B}_{\rho}}\|\lambda-\lambda^{0}\|_{\Gamma^{-1}}^{2}$ where $\mathcal{B}_{\rho}:=\{\lambda\mid\|\lambda\|\leq\rho\}$, for any $\rho>0$. Then for all $r>0$, we have

Let us briefly comment on the assumptions made in each part of the above statement. In part (a), the condition (33) imposes requirements on the parameters of the algorithm, such as the proximal matrix $\Omega$ and the matrix $\Xi$ (which contains all penalty parameters $\{\rho_{ij}\}$). Note that a sufficient condition for (33) is that $2\Omega\succ\widetilde{P}$, which is equivalent to $\omega_{i}>P_{i}/2$ for all i\in{\mbox{\mathcal{V}}}.

Also in part (b), $d_{x}$ represents the diameter of the feasible set $\mbox{dom}(h)$; $d_{z}$ can be viewed as the maximum size of any two $z$’s generated by the algorithm (see the first inequality in Appendix A-C); $d_{\lambda}(\rho)$ can be viewed as the distance between the initial solution $\lambda^{0}$ to the ball $\mathcal{B}_{\rho}$. Additionally, the boundedness of the set $\mbox{dom}(h)$ can be achieved for example when $h$ is the indicator function of certain bounded convex set.

The key novelty, as well as the main challenge, in the analysis of the proposed approach is a careful bounding of the proximal parameter $\omega_{i}$, which results in faster practical numerical convergence performance (to be shown in Section VII). Indeed, compared with the existing convergence results on proximal-based ADMM such as and , our bound for $\omega_{i}$ is reduced by at least a half. More importantly, no global information is needed at each agent to verify such condition, in contrast to . It is also interesting to note that the condition (33), which only guarantees convergence, is indeed weaker than the condition (34), which guarantees the global sublinear convergence rate.

Next we analyze the algorithm for static graph and stochastic gradient (i.e., Algorithm 1 applied to a static graph).

Suppose that Assumption 1 holds, and the graph is static and connected (with {\mbox{\mathcal{G}}}^{r}={\mbox{\mathcal{G}}} for all $r$). Suppose that $w^{r}$ is generated by Algorithm 1, and all the assumptions made Theorem IV.1(b) hold true. If additionally the penalty parameter sequence $\{\eta^{r}\}$ satisfies

In the previous two results, we have used $P(\bar{x}^{r},\bar{z}^{r})$ to measure the quality of the solution. This is a reasonable measure: according to [41, Lemma 2.4], when $\rho$ is large enough (in the sense that $\rho>\|\lambda^{*}\|$), $P(\bar{x}^{r},\bar{z}^{r})\leq\epsilon$ implies that

That is, both the constraint violation and the objective gap are in the same order as $\epsilon$. $\blacksquare$

We remark that the stochastic ADMM method for solving general linearly constrained problem has been discussed in several recent papers . However its application and the rate analysis in the context of distributed consensus based optimization appears to be new. In particular, compared with the SGADM proposed in , our scheme only linearizes the objective function $f_{i}$, but not the entire augmented Lagrangian. Further, the order of the updates of the two primal variables has been reversed. These key differences make the analysis in not directly applicable. $\blacksquare$

IV-B Analysis for Random Graphs

In this section we analyze the convergence properties of Algorithm 1 (DySPGC) for random graphs defined in Definition II.1. The convergence claims are similar to those given in the previous section, but in the sense of convergence in expectation or with probability 1 (w.p.1).

We first analyze the simple case with exact gradient. To proceed, define a new function $J(x,z,\lambda)$ as

where $\Psi$ and $\Phi$ are given in (14). These quantities can be viewed as matrix scaled versions of their counterparts {$d_{x},d_{\lambda}(\rho)$} in the statement of Theorem IV.1.

The derivation of the following result is mostly based on that of Theorem IV.1; the details can be found in our technical report .

Suppose that Assumption 1 holds, and $\widetilde{G}(x^{r},\xi^{r+1})=G(x^{r}),\;\forall~{}r$. Suppose that the graph \{{\mbox{\mathcal{G}}}^{r}\} is generated according to Definition II.1. Then the following two statements hold true.

(a) If the following condition is satisfied

then $w^{r}$ generated by Algorithm 1 converges w.p.1. to a primal-dual solution of problem (P).

(b) Define $\bar{w}^{r}$ similarly as in the statement of Theorem IV.1. Suppose the following holds true

then Algorithm 1 generates a sequence $\bar{w}^{r}$ that satisfies

where $d_{J}:=\sup_{\lambda\in\mathcal{B}_{\rho}}J(x^{0},z^{0},\lambda)$.

Let us briefly compare the assumptions made in each of the statement. The condition $2\Omega\succ\widetilde{P}$ is equivalent to the condition that $\omega_{i}>P_{i}/2,\;\forall~{}i$, which implies that each local agent’s proximal parameter should be chosen larger than $P_{i}/2$. Again, this condition is more relaxed compared with the one given in part (b), which requires a set of larger local proximal parameters $\{\omega_{i}\}$.

It is interesting to note that the stepsize rules (37) and (38) are both implied by their respective counterparts (33) and (34), but the new rules are no longer related to the network structure. Finally we analyze the case where the gradients are stochastic [i.e., Algorithm 1 (DySPGC) in its most general form].

Define $\bar{w}^{r}$ similarly as in the statement of Theorem IV.1. Suppose that

Then Algorithm 1 generates a sequence $\bar{w}^{r}$ that satisfies

where $d_{J}$ is defined in Theorem IV.3(b).

The detailed proof can be found in our technical report . We note that comparing with the existing analysis for random graphs in , our proof further takes into account inexact gradient information, and it does not assume any strong convexity on the objective functions.

V Comparison with Existing Algorithms

Our proposed DySPGC as well as its special case PGC is closely related to a few existing algorithms. In this section we provide a detailed account of such relations; see Table II for a summary.

Recently, an IC-ADMM algorithm is proposed in , which solves the following problem in a distributed manner

V-B Connection with the DLM algorithm

The Decentralized Linearized Alternating Direction Method of Multipliers (DLM) proposed in is closely related to IC-ADMM. The DLM solves (3) with $h_{i}\equiv 0$. Its basic iteration is again Algorithm 2 (PGC) with parameters $\rho_{ij}=\rho>0$ and $\omega_{i}=\omega\geq 0$ for all $i,j$. The convergence condition in [31, Theorem 1] is given by (described using our notation)

This condition is an immediate consequence of the condition (33) (with uniform $\rho_{ij}$’s and uniform $\omega_{i}$’s).

V-C Connection with EXTRA

We show that Algorithm 2 (PGC) can be viewed as a generalization of the EXTRA . Consider applying Algorithm 2 (PGC) to problem (P) with a smooth objective (i.e., $h_{i}\equiv 0$ for all $i$). According to Proposition III.1, one can write the iterates of Algorithm 2 as

Eq. (40) is precisely the EXTRA update developed in , except for the two relatively minor points:

In (40) a slightly more general matrix stepsize $\Upsilon^{-1}$ is used instead of the scalar stepsize used in EXTRA.

The EXTRA allows a slightly wider choice of $\widetilde{W}$, i.e., ${1}/{2}(I_{MN}+W\otimes I_{M})\succeq\widetilde{W}\succeq W$, $\mbox{null}\{W-\widetilde{W}\}=\mbox{span}\{\mathbf{1}\}$ and $\mbox{null}\{I_{MN}-\widetilde{W}\}\supseteq\mbox{span}\{\mathbf{1}\}$, where $\mathbf{1}$ is an all one vector of appropriate size. However except for the common choice (41), these conditions are difficult (if not impossible) to verify in a fully distributed manner.

When a single scalar stepsize is used (as was done in EXTRA), say $\beta=\beta_{i}=\beta_{j}>0$ for all $i,j$, then we can perform either one of the following procedures to identify the parameters of Algorithm 2 (PGC) (depending on whether the weight matrix $W$ is known a priori):

From Algorithm Parameters to Weight Matrix. Suppose the agents can select $\{\omega_{i}\}$ and $\{\rho_{ij}\}$. Then for any set of fixed $\{\rho_{ij}\}$’s, pick $\beta$ and $\omega_{i}$’s such that

To compare the convergence result in Theorem IV.1 and that of [14, Theorem 3.3], note that when the scalar stepsize is used, we have $\Upsilon=\beta I_{MN}$. Therefore a sufficient condition to guarantee the condition in Theorem IV.1 is that

This is precisely the condition set forth in [14, Theorem 3.3].

From the above expression it is clear that $\beta$ depends on all the local functions, therefore it has to be decided in a centralized manner. In contrast, the stepsize parameters in PGC can be chosen as: $\omega_{i}\geq P_{i}/2$ (cf. the remarks made after Theorem IV.1). The latter choice is simple, distributed implementable, and more importantly it results in improved convergence speed in practice, especially when the curvatures of $g_{i}$’s vary significantly, i.e., $\max_{i}{P_{i}}\gg\min_{i}{P_{i}}$. This will be demonstrated in Section VII.

We comment that in a couple of recent works and , the authors have established that EXTRA is also related (and in fact in most cases equivalent) to certain saddle point method, and certain proximal augmented Lagrangian method. Combining the observation made in this work, we can conclude that all these methods (i.e., the saddle point method , the proximal augemented Lagrangian method , the EXTRA and the PGC) are all closely connected We thank the anonymous reviewer for bringing these new developments to our attention..

V-D Connection with PG-EXTRA

One can also show that the proposed Algorithm 2 (PGC) generalizes the PG-EXTRA . According to the argument leading to (30), one can explicitly express (III.1) by

By the definition of $c_{i}$ in (III.1), we have

where $\widehat{W}_{i}$ and $\widetilde{W}_{i}$ denote the $i$th row of $\widehat{W}$ and $\widetilde{W}$ (as have been defined in (41)), respectively. Again by (III.1), and assume that $x^{0}=0$ and $\nabla g_{i}(x^{-1}_{i})=0$, we can check that

Combining the above three equalities we have

This is the PG-EXTRA proposed in [15, Algorithm 1].

V-E Connection with the DSG Method

Below we show that Algorithm 2 (PGC) is closely related to the DSG iteration (2). Assume for simplicity that $h_{i}\equiv 0$ for all $i$. Suppose that the $z$ and $\lambda$ steps of the PGC remain the same while the $x$-step (21a) is replaced by the following

That is, in the $x$-step we let $\lambda^{r}=0$. The claim is that by such modification one recovers the DSG iteration (2). To argue this, we write down the optimality condition of the modified iteration as

Following the derivation of Proposition III.1 until (52) in Appendix -A, we have

Note that compared with (52), the first equality above has an additional term $-\alpha^{r}$. Plugging the second equality into the first one, we obtain

By the definition of the matrices $M_{+}$ and $M_{-}$ in (19), one can verify the following identities

Utilizing (43), and by the definition of $\beta_{i}$ (22) and the definition of the weight matrix $W$ in (26), we can write the above iteration compactly as

After picking a uniform scalar stepsize $\beta_{i}=\beta_{j}=\beta>0$ (cf. Section V-C for how this can be done), we immediately get the DSG iteration (2) [with a weight matrix given by $\widetilde{W}=\frac{1}{2}(I_{MN}+W\otimes I_{M})$].

Obviously, our convergence analysis does not work for this variant, as the $x$-update is no longer related to the dual variable $\lambda$. Indeed, to prove convergence of the DSG, an iteration-dependent and increasing $\beta$ is needed, and such convergence is usually slower than $\mathcal{O}(1/r)$; see and the references therein. Nevertheless, the above observation reveals a fundamental connection between the ADMM-based method and the classical DSG method.

VI Extension to Accelerated DySPGC

The relationship identified between the DySPGC and the EXTRA, PG-EXTRA, IC-ADMM etc. provides a systematic way to analyze and generalize various existing algorithms. In this section, we provide one such generalization which accelerates the DySPGC (hence the EXTRA, PG-EXTRA, IC-ADMM, etc). The algorithm is inspired by .

For simplicity, we will restrict ourselves to the static graphs in this section. Let $\{\eta^{r},\theta^{r},\nu^{r}\geq 0\}$ denote a sequence of iteration-dependent parameters, whose values will be given shortly; Let $\{x^{r,{\rm md}},x^{r,{\rm ag}},z^{r,{\rm ag}},\lambda^{r,{\rm ag}}\}$ denote a sequence of auxiliary variables. The proposed accelerated algorithm is given in the table below.

Algorithm 3. Accelerated DySPGC Over Static Graphs At iteration , let $B^{T}\lambda^{0}=0$, $z^{0}=z^{0,{\rm ag}}=\frac{1}{2}M^{T}_{+}x^{0}$. At each iteration $r+1$, update the variable blocks by: $\displaystyle x^{r+1,{\rm md}}$ $\displaystyle=(1-\nu^{r})x^{r,{\rm ag}}+\nu^{r}x^{r}$ (44a) $\displaystyle x^{r+1}$ \displaystyle=\arg\min_{x}\;\bigg{\{}\left\langle\widetilde{G}(x^{r+1,{\rm md}},\xi^{r+1}),x-x^{r}\right\rangle+h(x) (44b) \displaystyle+\frac{1}{2}\left\|Ax+Bz^{r}+\Gamma^{-1}{\lambda^{r}}\right\|_{\Gamma}^{2}+\frac{1}{2}\|x-x^{r}\|^{2}_{\theta^{r}\Omega+{\eta^{r+1}}I_{MN}}\bigg{\}} $\displaystyle x^{r+1,{\rm ag}}$ $\displaystyle=(1-\nu^{r})x^{r,{\rm ag}}+\nu^{r}x^{r+1}$ (44c) $\displaystyle z^{r+1}$ $\displaystyle=\arg\min_{z}\;\frac{1}{2}\left\|Ax^{r+1}+Bz+\Gamma^{-1}{\lambda^{r}}\right\|_{\Gamma}^{2}$ (44d) $\displaystyle z^{r+1,{\rm ag}}$ $\displaystyle=(1-\nu^{r})z^{r,{\rm ag}}+\nu^{r}z^{r+1}$ (44e) $\displaystyle\lambda^{r+1}$ $\displaystyle=\lambda^{r}+\Gamma\left(Ax^{r+1}+Bz^{r+1}\right)$ (44f) $\displaystyle\lambda^{r+1,{\rm ag}}$ $\displaystyle=(1-\nu^{r})\lambda^{r,{\rm ag}}+\nu^{r}\lambda^{r+1}$ (44g)

First note that $x^{r,{\rm ag}},z^{r,{\rm ag}},\lambda^{r,{\rm ag}}$ are convex combinations of all previous iterates $\{x^{t}\}_{t=1}^{r}$, $\{z^{t}\}_{t=1}^{r}$, $\{\lambda^{t}\}_{t=1}^{r}$, respectively. Second, $x^{r+1,{\rm md}}$ is an intermediate point on which the stochastic gradient is evaluated. Therefore in total there are three sequences related to the $x$ update, resembling the Nesterov’s acceleration scheme .

The convergence rate of Algorithm 3 can be analyzed similarly as in , we include the proof in the Appendix for completeness. Compared with the bound given in Theorem IV.2, the accelerated version is able to significantly reduce the scaling with respect to $\max_{i}w_{i}$, which in turn depends on the network structure as well as the Lipschitz constants of the local gradients through (34).

Suppose that the assumptions made in Theorem IV.2 are true. Further let

Assume that the stepsize matrix satisfies

Then the iterates generated by Algorithm 3 satisfy

VII Numerical Results

In this section, we present some simulation results of the proposed algorithms by solving the following LASSO problem

We compare Algorithm 2 (PGC) with PG-EXTRA in , EXTRA in , DLM in , distributed gradient descent (DGD) algorithm in and the distributed Nesterov gradient descent (DNG) algorithm in . We also compare the static version of Algorithm 1 (i.e., SPGC) with a distributed stochastic gradient descent method (D-SGD) . The stepsize for the EXTRA/PG-EXTRA is chosen according to the sufficient condition suggested in [14, Theorem 3.3], and the weight matrix $W$ is the Metropolis constant edge weight matrix. For Algorithm 1 (resp. Algorithm 2), $\omega_{i}={P_{i}}/{2}$ (resp. $\omega_{i}={P_{i}}$) and $\rho_{ij}=10^{3}$ for all $i,j$. For DLM, the parameter $c$ in [31, Eqn. (21)] (which is equivalent to $\rho_{ij}$ here) is set to $10^{3}$, and $\rho$ (which corresponds to $\omega_{i}$ here) is set such that $\xi$ in [31, Eqn. (21)] equals zero. For the DGD, DNG and D-SGD algorithms, the Metropolis weight matrix is used and the stepsize is set as $0.01/(r+5000)$ The parameters are not searched in an exhaustive manner. Instead, for example, the parameter $\rho_{ij}$ in Algorithms 1 and 2 are tested for values 0.1, 1, 10, 50, 100, 500, 1000, 5000 and so on, and the value 1000 is chosen as the algorithm yields the fastest convergence performance among the others.. To measure the progress of different algorithms, we define the following two quantities

where $f^{*}$ is the optimal objective value of problem (47) and is obtained by the FISTA method , and $\hat{x}^{r}=\frac{1}{N}\sum_{i=1}^{N}x_{i}^{r}$.

The convergence curves of the proposed PGC (Algorithm 2) and the PG-EXTRA are shown in Figure 1, by assuming static networks and exact gradient information. Two problem settings are considered: Case 1). $M=1000,K=200,\nu=0.1$ and Case 2). $M=1000$, $K=50$, $\nu=50$. Given $N=16$, problem (47) is strongly convex for Case 1 and (non-strongly) convex for Case 2. One can see from Figure 1 the proposed algorithm outperforms the PG-EXTRA, in terms of both accuracy and consensus error. This is expected since compared with the PG-EXTRA, the PGC is able to use larger and more flexible stepsizes, as discussed in Section V-C.

As EXTRA, DLM, DGD and DNG are developed for smooth problems, we set $\nu=0$, $M=1000$ and $K=200$ to problem (47) and display the comparison results in Figure 2. Analogously, one can see from this figure that the proposed PGC performs the best and outperforms the DLM and EXTRA. Besides, the DLM, EXTRA and the proposed PGC all converge much faster than the DNG and DGD, which is consistent with the comparison results reported in .

In Figure 3, we present convergence curves of the proposed SPGC (Algorithm 1 with static graph), PG-EXTRA and D-SGD when the stochastic gradient information and the setting of Case 1 are used. The gradient noise power $\sigma^{2}$ is set to $0.1$ and $10$, respectively. Following Theorem IV.2, the stepsize $\eta^{r}$ of the SPGC is set to $\eta_{0}\sqrt{r}$, where $\eta_{0}\geq 0$. Note that when $\eta_{0}=0$, the SPGC reduces to the PGC in Algorithm 2. One can first observe from Figure 3(a) that, due to the stochastic gradients, both the SPGC ($\sigma^{2}=0.1,\eta_{0}=2500$) and the PG-EXTRA ($\sigma^{2}=0.1$) suffer higher error floors than their counterparts with exact gradients in Figure 1. Moreover, it can also be seen from Figure 3(b) that the PG-EXTRA ($\sigma^{2}=0.1$) achieves lower consensus error than the SPGC ($\sigma^{2}=0.1,\eta_{0}=2500$). However, as seen from Figure 3(a), not only the SPGC ($\sigma^{2}=0.1,\eta_{0}=2500$) converges faster than the PG-EXTRA ($\sigma^{2}=0.1$), but also the achieved solution accuracy keeps decreasing with the iteration number. This is contrast to the PG-EXTRA ($\sigma^{2}=0.1$) whose accuracy is limited by an error floor. We can observe similar convergence results for the case with $\sigma^{2}=10$. Finally, it can be seen that the D-SGD converges much slower than the other methods. Note that the convergence curves of the D-SGD for $\sigma^{2}=0.1$ and $\sigma^{2}=10$ overlap, implying that the noisy gradients have less impact on the D-SGD.

In the last example, we examine the convergence behavior of the proposed DySPGC (Algorithm 1) over a time-varying network. Following Definition II.1, we assume that each link $(i,j)\in\mathcal{E}$ has a probability $p_{ij}\in(0,1]$ being active (If $p_{ij}=1$ $\forall i,j$, then the DySPGC reduces to SPGC in static networks). The setting Case 1 is considered with gradient noise power $\sigma^{2}=0.1$, and the stepsize $\eta^{r}=2500\sqrt{r}$ of DySPGC is used. Figure 4 displays the convergence curves of the DySPGC for various values of $p_{ij}$. As seen, the DySPGC exhibits considerable robustness against the time-varying networks.

VIII Conclusion

In this paper we have proposed a dynamic stochastic proximal-gradient consensus (DySPGC) algorithm for solving a convex, possibly stochastic optimization problem over a randomly time-varying multi-agent network. We have analyzed the global convergence rate for DySPGC under various different scenarios, such as when the network is static/dynamic, or when the gradient is stochastic/deterministic. Our numerical results show that the proposed algorithms compare favorably with a few EXTRA based algorithms under various scenarios. Interestingly, Our algorithmic framework provides a unifying perspective for a few popular algorithms for distributed convex optimization. Such new perspective allows significant generalization of these methods based on existing theories such as the primal-dual methods. As an example, leveraging upon the recent work , we can develop an accelerated version of DySPGC, which is capable of further reducing certain constants in the convergence rate; see our technical report for details.

There are a few interesting directions that we would like to pursue in the future. For example, can we generalize the algorithms and their analysis to problems with nonconvex objective functions? Can we deal with a wider types of network dynamics such as the deterministic B-strongly connected networks? Is there a connection between distributed algorithms that we have studied in this work, with the optimization algorithms that minimizes a convex objective function consisting of a finite sum of components (which often arise in machine learning related applications), such as the SAG algorithm and the SVRG algorithm ? Some recent advancement in connecting distributed optimization methods with SAG/SVRG can be found in and .

The proof of this proposition is in fact straightforward. One only needs to start from the optimality condition of each step of Algorithm 2, and recognize that each $\lambda^{r}$ has some symmetric structure, and that the iterates $\{z^{r}\}$ can be expressed using $\{x^{r}\}$.

At iteration this is true due to the initialization $B^{T}\lambda^{0}=0$. At iteration $r\geq 0$, by (48b) and (48c) we have $B^{T}\lambda^{r+1}=0.$ This immediately implies \delta^{r+1}_{ij}=-\gamma^{r+1}_{ij},\;\forall~{}e_{ij}\in{\mbox{\mathcal{A}}}. Using this identity, we can define a new variable $\alpha$ as

Applying (49) to (48b) and the definition (13), we have

This fact combined with the update rule for $\lambda$ implies that

Utilizing the initial conditions $B^{T}\lambda^{0}=0$, $z^{0}=\frac{1}{2}M^{T}_{+}x^{0}$, and the fact that $A^{T}\lambda^{r+1}=\alpha^{r+1}$, the $x$-step optimality condition (48a) can be written as

Plugging in (43j) and utilizing (51), (52) becomes

Moreover, (51) can be expressed as, \forall~{}i\in{\mbox{\mathcal{V}}},

Next we remove the sequence $\{\alpha^{r}\}$ from the $x$ iterations. This is the key step towards obtaining a single-variable characterization. To this end, we subtract (53) by the same update for iteration $r$. By the definition of $\beta_{i}$ in (22) we can derive the following update rule for node $i$ at iteration $r+1$:

Note that by the definition of $W$ in (26), we have

​​which is simply a weighted average of $x^{r}$ over all the neighbors of node $i$ (including itself).

Writing in vector form and utilizing the definition of $W$ in (26), we have

Appendix A Preliminary Results

In this section we summarize a few preliminary results and identities that will be used later for proof of convergence of both Algorithm 1 and 2. We also provide a brief outline of the convergence proof.

First we discuss the optimality condition for problem (P). Suppose that Assumption 1 holds. Let $y^{*}\in X^{*}$ denote an optimal solution of (3). Let $z^{*}_{ij}=y^{*}$, (i,j)\in{\mbox{\mathcal{A}}} and $x^{*}_{i}=y^{*}$, for all $i$. Due to equivalence of problems (3) and (8), $(z^{*},x^{*})$ is an optimal solution of (P). From the assumed Slater condition we know that problem (P) has a saddle point $(x^{*},z^{*},\lambda^{*})$ satisfying the following condition $\forall~{}z,\lambda,$ and $\forall~{}x\in\mbox{dom}(h)$

where $L_{0}(\cdot)$ is given by (17) with $\Gamma\equiv 0$.

The second inequality in (55) is equivalent to the following

The second inequality in (55) also implies that $(x^{*},z^{*})$ is the optimizer for the following problem

The first-order optimality condition of the above problem is given by the following (for some $\zeta^{*}\in\partial h(x^{*})$)

For notational simplicity, let us define the left hand side by $U(w,w^{*})$.

It is easy to observe that for all $x\in\mbox{dom}(h)$ and all $z,\lambda$,

where in the second to the last equality we have used the fact that $Ax^{*}+Bz^{*}=0$. Using the above identity, (56)–(57) are equivalent to the following two inequalities, respectively

Let us characterize the optimality condition for the iterates. Define a block diagonal matrix ${H}^{r+1}(\eta)$:

Let ${H}(\eta)$ denote its time-invariant counterpart, that is, replacing $B^{r+1}$ in the above definition by $B$. Also define

Using the fact that $\lambda^{r+1}=\lambda^{r}+\Gamma(Ax^{r+1}+Bz^{r+1})$, the optimality conditions for the subproblems of Algorithm 1 are given by [for all $x\in\mbox{dom}(h)$ and all $z,\lambda$]

Adding these conditions we obtain $\forall~{}x\in\mbox{dom}(h),\;\forall~{}z,\lambda$,

Note that $(B^{r+1})^{T}\lambda^{r}=(B^{r})^{T}\lambda^{r}=0$ because $\lambda^{r}=[\delta^{r};-\delta^{r}]$, and each $B^{r+1}$ and $B^{r}$ stacks two identical matrices. Using this fact and the optimality condition of the $z$-step (60b), the following is true for any optimal solution $(z^{*},x^{*})$

where we have defined the gradient error as

Next let us briefly provide an outline of the proof of convergence. To show convergence (i.e., Theorem IV.1), we need to construct a potential function that decreases at each iteration. In our proof, the following quantity is used as the potential function

To show that such a measure decreases at each iteration, we need to utilize the optimality condition (A-A) that we have just derived from the execution of the algorithm, as well as the global optimality conditions (56) and (57).

A-B Proof of Theorem IV.1

We only prove the first part of the theorem. The second part is the consequence of Theorem IV.2. As $\eta^{r}=0$ for all $r$, and {\mbox{\mathcal{G}}}^{r}={\mbox{\mathcal{G}}} for all $r$, we denote $H:=H(0)$. Applying the static version of (A-A) and let $w^{*}:=(x^{*},z^{*},\lambda^{*})$, we have

where $\zeta^{*}\in\partial h(x^{*})$. Similarly, we have

​​where in $\rm(i)$ we have used the Young’s inequality: $\langle a,b\rangle\leq{\|a\|^{2}}/{(2\epsilon)}+{\epsilon\|b\|^{2}}/{2}$ for any $\epsilon>0$, and a key property due to Nesterov [47, Theorem 2.1.5]. Namely, if $g_{i}(x_{i})$ is convex with Lipschitzian gradient (constant $P_{i}$), then $\forall~{}x_{i},y_{i}\in X$,

Combining the optimality condition (59b) and the above two inequalities, we obtain

For time-invariant graph, (50) is true, which implies

Plugging this relation into (67) we obtain

Therefore, as long as $\Omega+\frac{1}{4}M_{+}B^{T}\Gamma BM^{T}_{+}-\frac{1}{2}\widetilde{P}\succ 0$ or equivalently $2\Omega+M_{+}(\Xi\otimes I_{M})M^{T}_{+}-\widetilde{P}\succ 0,$ we will have $x^{r+1}\to x^{r}$, $\lambda^{r+1}\to\lambda^{r}$. By a standard argument (cf. the derivation in [30, (A2.22)-(A2.25)]), we can argue that every limit point of the sequence $x^{r}$ and $\lambda^{r}$ is a primal dual optimal solution of problem (P). Finally, by noticing the identity $2\Omega+M_{+}(\Xi\otimes I_{M})M^{T}_{+}=\Upsilon(W\otimes I_{M}+I_{MN})$ by using the definitions of $W,\Upsilon$ and $M_{+}$, the theorem is proved.

A-C Proof of Theorem IV.2

Note that the assumption of boundedness of $x$ implies the boundedness of iterates $\{z^{r}\}$. This is because from the identity (50) we have $z^{r+1}_{ij}=\frac{1}{2}(x^{r+1}_{i}+x^{r+1}_{j})$ for all $r$. Therefore

​​By the convexity of $h$ and $g$ and the Lipschitz continuity of $\nabla g$, we have

​​where in $\rm(i)$ we have again used the Young’s inequality; in the last inequality we have used the assumption (34) (cf. the derivation in (69)). Evaluating the LHS based on the average of the iterates {\mbox{\bar{w}}}^{r+1}, and using convexity, we have

​​Further, we have the following series of inequalities

​​Taking the supreme of both sides of (71), we obtain

​​Taking the expectation on both sides of the above inequality and utilize the assumption made in (15) about the stochastic gradient, and the fact that

A-D Proof of Theorem IV.3

Our proof is motivated by . Suppose at iteration $r$ we have iterate $w^{r}=(x^{r},z^{r},\lambda^{r})$ and we are about to execute Algorithm 1. Consider the virtual sequence $(\hat{x}^{r+1},\hat{z}^{r+1},\hat{\lambda}^{r+1})$ generated by Algorithm 1 (based on $w^{r}$) with all nodes and edges being active (i.e., with {\mbox{\mathcal{A}}}^{r+1}={\mbox{\mathcal{A}}} and {\mbox{\mathcal{V}}}^{r+1}={\mbox{\mathcal{V}}}). Then from (67) in the proof of Theorem IV.1, we must have

​​where $\Psi$ and $\Phi$ are given in (14). Also define \mathcal{F}^{r}=\{x^{t},z^{t},\lambda^{t},{\mbox{\mathcal{G}}}^{t}_{d},t=1,\cdots,r\} as the filtration up to iteration $r$. The following is easy to verify

Summing up (75) – (77) and utilizing (74), we obtain

where in the last inequality we have removed the term $(\hat{z}^{r+1}-z^{r})^{T}B^{T}\Gamma B(\hat{z}^{r+1}-z^{r})\succeq 0$. Using the assumption (37), we conclude that the sequence $D_{x}(x^{r},x^{*})+D_{z}(z^{r},z^{*})+D_{\lambda}(\lambda^{r},\lambda^{*})$ is a nonnegative almost supermartingale, which is convergent by the nonnegative almost supermartigale convergence theorem [53, Theorem 1]:

Then again by a standard argument (cf. ) we conclude that $(x^{r},z^{r},\lambda^{r})$ as well as $(\hat{x}^{r},\hat{z}^{r},\hat{\lambda}^{r})$ converge with probability one to a primal-dual solution of problem (P).

A-E Proof of Theorem IV.4

​​Using (80) – (81) and the definition of $J(\cdot)$ in (36), the conditional expectation of $J(\cdot)$ can be expressed as below

then its conditional expectation is given by

Plugging (78) and (75) – (77) into (82), we obtain a bound on the conditional expectation of $J(\cdot)$, given below

Let us define $\bar{x}^{r+1}=\frac{1}{r+1}\sum_{t=0}^{r}x^{t}$ and $\bar{z}^{r+1}$ similarly. Taking expectation wrt $\mathcal{F}^{r}$ and summing over $t$, (83) becomes

The rest of the proof follows the last part of the proof of Theorem IV.2.

A-F Proof of Theorem VI.1

We first provide a lemma that bounds the quantity $Q(\cdot,\cdot)$ defined in (56).

for some $\zeta^{*}\in\partial h(x^{*})$. Proof. From the definition of $x^{r+1,{\rm ag}}$, $x^{r+1,{\rm md}}$ we have

First it is easy to show that for any feasible $x$, we have (cf. [42, Eq. (2.16)])

Using this result, we have the following series of inequalities

where the inequality uses (88), the convexity of $h(\cdot)$ and the update rule of $x^{r+1,{\rm ag}}$.

We then proceed to prove Theorem VI.1. Let us define $\varpi^{r}=\frac{2}{(r+1)r}$. From (45) one can check that the following two identities hold

Similarly as in (A-A), we can derive (for some $\zeta^{r+1}\in\partial h(x^{r+1})$)

From the assumption $\Omega$ should satisfy (46). Using such assumed bound, the definition of $\nu^{r}$ and $\theta^{r}$, and the fact that $\nu^{r}<1$ and $M_{+}(I_{M}\otimes\Xi)M^{T}_{+}\succeq 0$, we have

Applying the same derivation as in (70), and divide both sides of (89) by $\varpi^{r}$, we obtain

Let us then analyze the successive sum of the RHS of the above inequality. Note that the sequences $\{\frac{\nu^{r}}{2\varpi^{r}},\frac{\nu^{r}\eta^{r+1}}{2\varpi^{r}}\}$ are both increasing sequences, and the sequence $\frac{\nu^{r}\theta^{r}}{2\varpi^{r}}$ is non-increasing. Thus from [42, Lemam 2.4] we have

Combining these results, and noticing $(1-\nu^{1})/\varpi^{1}=0$, we obtain

Notice that from the derivation in (72) we have

Taking expectation on both sides of the above inequality and utilizing

References