Decentralized Primal-Dual Proximal Operator Algorithm for Constrained Nonsmooth Composite Optimization Problems over Networks

Abstract

In this paper, we focus on the nonsmooth composite optimization problems over networks, which consist of a smooth term and a nonsmooth term. Both equality constraints and box constraints for the decision variables are also considered. Based on the multi-agent networks, the objective problems are split into a series of agents on which the problems can be solved in a decentralized manner. By establishing the Lagrange function of the problems, the first-order optimal condition is obtained in the primal-dual domain. Then, we propose a decentralized algorithm with the proximal operators. The proposed algorithm has uncoordinated stepsizes with respect to agents or edges, where no global parameters are involved. By constructing the compact form of the algorithm with operators, we complete the convergence analysis with the fixed-point theory. With the constrained quadratic programming problem, simulations verify the effectiveness of the proposed algorithm.

1. Introduction

Recently, the distributed data processing methods based on multi-agent networks have received much attention. The traditional methods put all the data into one machine and perform the computation centrally. However, as the size of data continues to grow, this kind of centralized strategy is limited by the computing power of the hardware. In contrast to this, the distributed methods distribute computing tasks to agents over decentralized networks [1,2]. Each agent keeps an arithmetic unit and a memory unit. The agents interact with each other through communication links, and this communication occurs only among the neighboring agents. Under these conditions, the distributed methods can effectively solve the optimization problems common to sensor networks [3], economic dispatch [4,5,6], machine learning [7,8] and dynamic control [9].

The existing decentralized algorithms have included some successful results [10,11,12,13,14,15]. Previous works considered the problem models composed of a single function. With a fixed stepsize, Shi et al. designed EXTRA [10], which can exactly converge to the optimal solution. Lei et al. studied problems with bound constraints and proposed the primal-dual algorithm [11]. In addition, recent works [16,17,18] investigated a general distributed optimization with an objective function by designing decentralized subgradient-based algorithms, but diminishing or non-summable step-sizes are utilized, which may cause slow convergence rates [19].

In order to make full use of these special properties, some scholars have studied the nonsmooth composite optimization problems, which possess smooth and nonsmooth structures. By extending EXTRA to the nonsmooth combinational optimization, Shi et al. proposed PG-EXTRA [20]. Li et al. introduced the network-independent stepsize to PG-EXTRA and then developed NIDS [21]. In addition, Aybat et al. proposed DPDA-D [22] for time-varying networks. Considering the situation that the nonsmooth term cannot be split, Xu et al. proposed the [23]. PG-ADMM [24] was designed based on the distributed alternating direction multiplier method. Particularly, the nonsmooth combinational optimization problems also include a class of problems consisting of three functions and a linear operator. This structure is mainly discussed in the centralized optimization [25,26,27,28], and recently some distributed works also appear in [29,30]. In this paper, inspired by the constrained optimization problem [31], we study the constrained nonsmooth composite optimization problems over networks.

The contributions of this paper can be summarized as follows:

• This paper focuses on an optimization problem with partially smooth and nonsmooth objective functions, where the decision variable satisfies local equality and feasible constraints, unlike these works [10,16,18,19,20,21] without considering any constraints. Then, to solve this problem, we propose a novel decentralized algorithm by combining primal-dual frame with the proximal operators, which avoids the estimation of subgradients for nonsmooth terms.
• Different from existing node-based methods [16,17,18,19,20,21], the proposed algorithm adopts an edge-based communication pattern that explicitly highlights the process of information exchange among neighboring agents and further gets rid of the dependence on Laplacians [13]. Such a consideration also makes it possible to use uncoordinated stepsizes instead of commonly global or dynamic ones [10,12,16,18,19,21].
• By employing the first-order optimal conditions and fixed-point theory of operators, the convergence is proved, and its sublinear rate $O\left(1/k\right)$ (k is the number of iteration); i.e., at most, $O\left(1/ϵ\right)$ iterations in order to reach an accuracy of $ϵ$ is established.

Organization: The rest of this paper is organized as follows. In Section 2, the necessary notations and basic knowledge are first provided, and then we describe the optimization problem over the networks and necessary assumptions. Section 3 supplies the development of the proposed decentralized algorithm. In Section 4, the convergence analysis for the proposed algorithm is provided. In Section 5, we use the simulation experiments to verify the theoretical analysis. Finally, conclusions are given in Section 6.

2. Preliminaries

In this section, we introduce the notations involved in this paper. Meanwhile, the objective problem and its explanation are also supplied.

2.1. Graph Theory and Notations

The knowledge of graph theory is used to construct the mathematical model of the communication network. Let $\mathcal{G}=\left(\mathcal{V},\mathcal{E}\right)$ describe the network as a graph, where $\mathcal{V}$ is the set of vertices and $\mathcal{E}\subset \mathcal{V}\times \mathcal{V}$ is the set of edges. For an agent $i\in \mathcal{V}$, ${\mathcal{N}}_i$ denotes the set of its neighbors. Let the unordered pair $(i,j)\in \mathcal{E}$ represent the edge between agent i and agent j. However, $(i,j)$ or $(j,i)$ is still order, i.e., the variables with respect to them are different.

Next, we explain the notations that appear in this paper. Let $\mathbb{R}$ represent the set of real numbers. Therefore, ${\mathbb{R}}^n$ denotes the n-dimensional vector space, and ${\mathbb{R}}^{n\times m}$ denotes the set of all n-row and m-column real matrices. We define $I_n$ as the n-dimensional identity operator, ${\mathbf{0}}_n$ as the n-dimensional null vector, and ${\mathbf{0}}_{n\times n}$ as the null matrix. If their dimensions are clear from the context, we omit their subscript. Then, $\mathrm{blkdiag}\{P,Q\}$ is the block diagonal matrix grouped by matrices P and Q. For a matrix P, let M be its transpose. We denote ${∥x∥}_P=\sqrt{x^{\top }Px}$ as the induced norm with matrix P. The subdifferential of function f is $\partial f$, where $\partial f\left(x\right)=\left\{p|\forall q\in {\mathbb{R}}^n,f\left(x\right)+p^{\top }\left(q-x\right)\le f\left(q\right)\right\}$. The conjugate function $f^{*}$ is defined by $f^{*}\left(p\right)={sup}_{q\in {\mathbb{R}}^n}\left\{〈p,q〉-f\left(q\right)\right\}$. For a positive constant $\mu$, the resolvent of the proximal operator is $\mathrm{pro}{\mathrm{x}}_{\mu f}\left(y\right)={\left(I+\mu \partial f\right)}^{-1}\left(y\right)=arg{min}_{x}f\left(x\right)+\frac{1}{2\mu }{∥x-y∥}^2$, while the resolvent with respect to the matrix M is ${\mathrm{prox}}_f^M\left(y\right)={\left(I+M^{-1}\partial f\right)}^{-1}\left(y\right)=arg{min}_{x}f\left(x\right)+\frac{1}{2}{∥x-y∥}_M^2$. Moreover, let $\mathcal{S}$ represent the optimal solution set of a solvable optimization problem over networks.

2.2. Decentralized Optimization Problem

The constrained composite optimization problem over networks studied in this paper is based on the network $\mathcal{G}=\{\mathcal{V},\mathcal{E}\}$ with m agents. Specifically, the formulation of the problem is established as follows:

$$\begin{array}{cc}\hfill & \underset{\tilde{x}\in {\mathbb{R}}^n}{min}\sum _{i=1}^m{f}_i\left(\tilde{x}\right)+g_i\left(\tilde{x}\right),\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill & \mathrm{s}.\mathrm{t}.A_i\tilde{x}=b_i,i=1,\cdots ,m,\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill & \tilde{x}\in {\bigcap }_{i=1}^m{\mathsf{\Omega }}_i.\hfill \end{array}$$

(1c)

In problem (1), $\tilde{x}\in {\mathbb{R}}^n$ is the decision variable; $f_i:{\mathbb{R}}^n\to \mathbb{R}\cup \left\{+\infty \right\}$ and $g_i:{\mathbb{R}}^n\to \mathbb{R}\cup \left\{+\infty \right\}$ are two private cost functions to agent i, where the former has the Lipschitz continuous gradient, but the latter may be nonsmooth; $b_i\in {\mathbb{R}}^r$ is a vector and $A_i:{\mathbb{R}}^n\to {\mathbb{R}}^r$ is a linear operator. Convex set ${\mathsf{\Omega }}_i$ gives the box constraints to the decision variable of agent i.

To clarify the properties of problem (1), the following necessary assumption is given.

Assumption 1.

For any agent $i\in \mathcal{V}$:

(i) 

The cost function $f_i$ is Lipschitz continuous and convex; i.e., if we consider the positive Lipschitz constant ${\beta }_i$, then it holds the inequality for the gradient $\nabla f_i$:

$$\begin{array}{c}\hfill {∥\nabla f_i\left(\tilde{x}\right)-\nabla f_i\left(\tilde{y}\right)∥}^2\le {\beta }_i{\left(\tilde{x}-\tilde{y}\right)}^{\top }\left(\nabla f_i\left(\tilde{x}\right)-\nabla f\left(\tilde{y}\right)\right).\end{array}$$

(2)

(ii) 

The local cost function $g_i$ is a nonsmooth and convex function.

(iii) 

The optimal solution ${\tilde{x}}^{*}$ to objective problem (1) exists, which satisfies both the equality constraints and the box constraints.

(iv) 

The graph $\mathcal{G}$ is undirected and connected.

Note that the cost functions $f_i$ and $g_i$ are separable. Hence, we introduce the consensus constraint to transform problem (1) into the structure that can be computed in a decentralized manner:

$$\begin{array}{cc}\hfill & \underset{x_1,\cdots ,x_m}{min}\sum _{i=1}^m{f}_i\left(x_i\right)+g_i\left(x_i\right),\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill & \mathrm{s}.\mathrm{t}.A_i{x}_i=b_i,i=1,\cdots ,m,\hfill \end{array}$$

(3b)

$$\begin{array}{cc}\hfill & x_i=x_j,i=1,\cdots ,m,j\in {\mathcal{N}}_i,\hfill \end{array}$$

(3c)

$$\begin{array}{cc}\hfill & x_i\in {\mathsf{\Omega }}_i.\hfill \end{array}$$

(3d)

Define the set ${\mathcal{A}}_i=\left\{z\in {\mathbb{R}}^n\left|A_{i}z=b_i\right.\right\}$ and consider the indicator function

$$\begin{array}{c}\hfill {\delta }_C\left(e\right)=\left\{\begin{array}{c}0,e\in C,\hfill \\ +\infty ,e\notin C,\hfill \end{array}\right.\end{array}$$

such that Problem (3) can be processed by the penalty function method. For $i\in \mathcal{V}$ and $j\in {\mathcal{N}}_i$, let $C_{ij}=I$ if $i<j$ and $C_{ij}=-I$ otherwise. Thus, Problem (3) is equivalent to the following problem:

$$\begin{array}{cc}\hfill & \underset{x_1,\cdots ,x_m}{min}\sum _{i=1}^m{f}_i\left(x_i\right)+g_i\left(x_i\right)+{\delta }_{{\mathcal{A}}_i}\left(x_i\right)+{\delta }_{{\mathsf{\Omega }}_i}\left(x_i\right),\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill & \mathrm{s}.\mathrm{t}.C_{ij}{x}_i+C_{ji}{x}_j=0,i=1,\cdots ,m,j\in {\mathcal{N}}_i.\hfill \end{array}$$

(4b)

Then, let $x=col\left(x_1,\dots ,x_m\right)$ be the global variable. For $i\in \mathcal{V}$ and $j\in {\mathcal{N}}_i$, we introduce a linear operator $N_{\left(i,j\right)}:x\mapsto {({\left(C_{ij}{x}_i\right)}^{\top },{\left(C_{ji}{x}_j\right)}^{\top })}^{\top }$, which generates the edge-based variable from x. With the set ${\mathcal{C}}_{\left(i,j\right)}=\left\{{\left(z_1^{\top },z_2^{\top }\right)}^{\top }\left|z_1+z_2=0\right.\right\}$, the constraint in the problem (4) can be transformed into another penalty function. Therefore, the problem (1) is finally equivalent to the following problem:

$$\begin{array}{cc}\hfill \underset{x_1,\cdots ,x_m}{min}& \sum _{i=1}^m{f}_i\left(x_i\right)+g_i\left(x_i\right)+{\delta }_{{\mathcal{A}}_i}\left(x_i\right)+{\delta }_{{\mathsf{\Omega }}_i}\left(x_i\right)+\sum _{i=1}^m\sum _{\left(i,j\right)\in \mathcal{E}}{\delta }_{{\mathcal{C}}_{\left(i,j\right)}}\left(N_{\left(i,j\right)}x\right).\hfill \end{array}$$

(5)

Based on the problem (5), we design a novel decentralized algorithm to solve the constrained composite optimization problem over networks in the next section.

3. Algorithm Development

The introduction with respect to the design process of the proposed algorithm is provided in this section.

Notice that Problem (5) is an unconstrained problem. According to [32] (Proposition 19.20), we obtain the following Lagrangian function:

$$\begin{array}{c}\hfill \begin{array}{c}\mathcal{L}={\displaystyle \sum _{i=1}^m}\left(f_i\left(x_i\right)+g_i\left(x_i\right)+v_i^{\mathrm{T}}{x}_i-{\delta }_{{\mathcal{A}}_i}^{*}\left(v_i\right)+u_i^{\mathrm{T}}{x}_i-{\delta }_{{\mathsf{\Omega }}_i}^{*}\left(u_i\right)\right)\hfill \\ +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{\left(i,j\right)\in \mathcal{E}}}\left(w_{\left(i,j\right)}^{\mathrm{T}}{N}_{\left(i,j\right)}x-{\delta }_{{\mathcal{C}}_{\left(i,j\right)}}^{*}\left(w_{\left(i,j\right)}\right)\right),\hfill \end{array}\end{array}$$

(6)

where $v_i\in {\mathbb{R}}^n$, $u_i\in {\mathbb{R}}^n$ and $w_{(i,j)}\in {\mathbb{R}}^{2n}$ are dual variables, and ${\delta }_{{\mathcal{A}}_i}^{*}$, ${\delta }_{{\mathsf{\Omega }}_i}^{*}$ and ${\delta }_{{\mathcal{C}}_{\left(i,j\right)}}^{*}$ are the conjugate functions of ${\delta }_{{\mathcal{A}}_i}$, ${\delta }_{{\mathsf{\Omega }}_i}$, ${\delta }_{{\mathcal{C}}_{\left(i,j\right)}}$, respectively. Notice that $w_{\left(i,j\right)}={(w_{ij}^{\top },w_{ji}^{\top })}^{\top }\in {\mathbb{R}}^{2n}$ is an edge-based variable, where $w_{ij}\in {\mathbb{R}}^n$ is the local variable of agent i and $w_{ji}\in {\mathbb{R}}^n$ is for agent j. Then, the last term of the Lagrangian function (6) satisfies:

$$\begin{array}{c}\hfill \sum _{i=1}^m\sum _{\left(i,j\right)\in \mathcal{E}}\left(w_{\left(i,j\right)}^{\mathrm{T}}{N}_{\left(i,j\right)}x-{\delta }_{{\mathcal{C}}_{\left(i,j\right)}}^{*}\left(w_{\left(i,j\right)}\right)\right)=\sum _{i=1}^m\sum _{j\in {\mathcal{N}}_i}\left(w_{\left(i,j\right)}^{\mathrm{T}}{N}_{\left(i,j\right)}x-{\delta }_{{\mathcal{C}}_{\left(i,j\right)}}^{*}\left(w_{\left(i,j\right)}\right)\right),\end{array}$$

Thus, the Lagrangian function (6) can also be written as

$$\begin{array}{c}\hfill \begin{array}{c}\mathcal{L}={\displaystyle \sum _{i=1}^m}\left(f_i\left(x_i\right)+g_i\left(x_i\right)+v_i^{\mathrm{T}}{x}_i-{\delta }_{{\mathcal{A}}_i}^{*}\left(v_i\right)+u_i^{\mathrm{T}}{x}_i-{\delta }_{{\mathsf{\Omega }}_i}^{*}\left(u_i\right)\right)\hfill \\ +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j\in {\mathcal{N}}_i}}\left(w_{\left(i,j\right)}^{\mathrm{T}}{N}_{\left(i,j\right)}x-{\delta }_{{\mathcal{C}}_{\left(i,j\right)}}^{*}\left(w_{\left(i,j\right)}\right)\right),\hfill \end{array}\end{array}$$

(7)

Taking the partial derivatives of the Lagrangian function (7) and combining the operator splitting method [29], we propose a new update flow as follows:   

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\overline{w}}_{\left(i,j\right)}^k={\mathrm{prox}}_{{\omega }_{\left(i,j\right)}{\delta }_{{\mathcal{C}}_{\left(i,j\right)}}^{*}}\left(w_{\left(i,j\right)}^k+{\omega }_{\left(i,j\right)}\left(N_{\left(i,j\right)}{x}^k\right)\right),\hfill \\ {\overline{u}}_i^k={\mathrm{prox}}_{{\mu }_i{\delta }_{{\mathsf{\Omega }}_i}^{*}}\left(u_i^k+{\mu }_i{x}_i^k\right),\hfill \\ {\overline{v}}_i^k={\mathrm{prox}}_{{\sigma }_i{\delta }_{{\mathcal{A}}_i}^{*}}\left(v_i^k+{\sigma }_i{x}_i^k\right),\hfill \\ x_i^{k+1}={\mathrm{prox}}_{{\gamma }_i{g}_i}\left(\begin{array}{c}{x}_i^k-{\gamma }_i\nabla f_i\left(x_i^k\right)-{\gamma }_i{\overline{v}}_i^k-{\gamma }_i{\overline{u}}_i^k-{\gamma }_i{\sum }_{j\in {\mathcal{N}}_i}{C}_{ij}^{\top }{\overline{w}}_{ij}^k\hfill \end{array}\right),\hfill \\ w_{\left(i,j\right)}^{k+1}={\overline{w}}_{\left(i,j\right)}^k+{\omega }_{\left(i,j\right)}\left(N_{\left(i,j\right)}{x}^{k+1}-N_{\left(i,j\right)}{x}^k\right),\hfill \\ u_i^{k+1}={\overline{u}}_i^k+{\mu }_i\left(x_i^{k+1}-x_i^k\right),\hfill \\ v_i^{k+1}={\overline{v}}_i^k+{\sigma }_i\left(x_i^{k+1}-x_i^k\right),\hfill \end{array}\right.\end{array}$$

(8)

where ${\overline{w}}_{\left(i,j\right)}={({\overline{w}}_{(i,j)}^{\top },{\overline{w}}_{(j,i)}^{\top })}^{\top }\in {\mathbb{R}}^{2n}$, ${\overline{u}}_i\in {\mathbb{R}}^n$, ${\overline{v}}_i\in {\mathbb{R}}^n$ are the auxiliary variables, and ${\gamma }_i$, ${\sigma }_i$, and ${\mu }_i$ are positive stepsizes. Notice that the stepsizes are uncoordinated, which can be selected independently related to different agents and enjoy their own acceptable ranges. Additionally, the edge-based parameters ${\omega }_{(i,j)}$ can be seen as inherent parameters of the communication network, revealing the quality of the communication.

The steps related to the edge-based variables in update flow (8) cannot be conducted directly, so we next replace them with the agent-based variables. We apply the Moreau decomposition to the first step in update flow (8) such that for the second term on the right side, we have

$$\begin{array}{c}\hfill {\mathrm{prox}}_{{\omega }_{\left(i,j\right)}^{-1}{\delta }_{{\mathcal{C}}_{\left(i,j\right)}}}\left({\omega }_{\left(i,j\right)}^{-1}{w}_{\left(i,j\right)}^k+N_{\left(i,j\right)}{x}^k\right)=\underset{y\in {\mathcal{C}}_{\left(i,j\right)}}{argmin}\left\{∥y-\left({\omega }_{\left(i,j\right)}^{-1}{w}_{\left(i,j\right)}^k+N_{\left(i,j\right)}{x}^k\right)∥\right\}.\end{array}$$

(9)

Define (9) as the projection ${\mathcal{P}}_{{\mathcal{C}}_{\left(i,j\right)}}\left({\omega }_{\left(i,j\right)}^{-1}{w}_{\left(i,j\right)}^k+N_{\left(i,j\right)}{x}^k\right)$. Then, according to the definition of the set ${\mathcal{C}}_{\left(i,j\right)}$, the projection has the following explicit expression:

$${\mathcal{P}}_{{\mathcal{C}}_{\left(i,j\right)}}\left[\begin{array}{c}{a}_1\\ a_2\end{array}\right]=\frac{1}{2}\left[\begin{array}{c}{a}_1-a_2\\ a_2-a_1\end{array}\right].$$

Thus, for $i\in \mathcal{V}$, $j\in {\mathcal{N}}_i$, the update step for ${\overline{w}}_{\left(i,j\right)}$ can be decomposed into

$$\begin{array}{c}\hfill {\overline{w}}_{ij}^k=\frac{1}{2}\left(w_{ij}^k+w_{ji}^k\right)+\frac{{\omega }_{\left(i,j\right)}}{2}\left(C_{ij}{x}_i^k+C_{ji}{x}_j^k\right).\end{array}$$

(10)

Moreover, the update step for $w_{\left(i,j\right)}$ can be replaced by

$$\begin{array}{c}\hfill w_{ij}^{k+1}={\overline{w}}_{ij}^k+{\omega }_{\left(i,j\right)}{C}_{ij}\left(x_i^{k+1}-x_i^k\right).\end{array}$$

(11)

Combining the update flow (8), (10) and (11), we finally propose the decentralized algorithm for Problem (1) in Algorithm 1.

Here, we directly give the stepsize condition of Algorithm 1 in the following assumption. The specific theoretical origin of this condition can be found in the convergence analysis section.

Assumption 2.

(Stepsize conditions)

For any agent $i\in \mathcal{V}$ and $j\in {\mathcal{N}}_i$, the stepsizes ${\gamma }_i$, ${\mu }_i$, ${\sigma }_i$ and ${\omega }_{\left(i,j\right)}$ are positive. Let the following condition hold:

$${\gamma }_i<\frac{1}{\frac{{\beta }_i}{2}+{\mu }_i+{\sigma }_i+{\displaystyle \sum _{j\in {\mathcal{N}}_i}}{\omega }_{\left(i,j\right)}},$$

where ${\beta }_i$ is the Lipschitz constant for the gradient $\nabla f_i$.

Algorithm 1 The Decentralized Algorithm

$\mathbf{Initialization}:$ For each agent $i\in \mathcal{V}$ and all $j\in {\mathcal{N}}_i$, let ${\overline{w}}_{ij}^0\in {\mathbb{R}}^n$, ${\overline{u}}_i^0\in {\mathbb{R}}^n$, ${\overline{v}}_i^0\in {\mathbb{R}}^n$, $x_i^0\in {\mathbb{R}}^n$, $w_{ij}^0\in {\mathbb{R}}^n$, $u_i^0\in {\mathbb{R}}^n$ and $v_i^0\in {\mathbb{R}}^n$. $\mathbf{For}$  $k=0,1,2,\cdots$  $\mathbf{do}$ Each agent i repeats, for all $j\in {\mathcal{N}}_i$, $$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\overline{w}}_{ij}^k& =\frac{1}{2}\left(w_{ij}^k+w_{ji}^k\right)+\frac{{\omega }_{\left(i,j\right)}}{2}\left(C_{ij}{x}_i^k+C_{ji}{x}_j^k\right),\hfill \\ \hfill {\overline{u}}_i^k& =\mathrm{pro}{\mathrm{x}}_{{\mu }_i{\delta }_{{\mathsf{\Omega }}_i}^{*}}\left(u_i^k+{\mu }_i{x}_i^k\right),\hfill \\ \hfill {\overline{v}}_i^k& =\mathrm{pro}{\mathrm{x}}_{{\sigma }_i{\delta }_{{\mathcal{A}}_i}^{*}}\left(v_i^k+{\sigma }_i{x}_i^k\right),\hfill \\ \hfill x_i^{k+1}& \hfill =\mathrm{pro}{\mathrm{x}}_{{\tau }_i{g}_i}\left(\begin{array}{c}{x}_i^k-{\gamma }_i\nabla f_i\left(x_i^k\right)-{\gamma }_i{\overline{u}}_i^k-{\gamma }_i{\overline{v}}_i^k-{\gamma }_i{\sum }_{j\in {\mathcal{N}}_i}{C}_{ij}^{\top }{\overline{w}}_{ij}^k\end{array}\right),\\ \hfill w_{ij}^{k+1}& ={\overline{w}}_{ij}^k+{\omega }_{\left(i,j\right)}{C}_{ij}\left(x_i^{k+1}-x_i^k\right),\hfill \\ \hfill u_i^{k+1}& ={\overline{u}}_i^k+{\mu }_i\left(x_i^{k+1}-x_i^k\right),\hfill \\ \hfill v_i^{k+1}& ={\overline{v}}_i^k+{\sigma }_i\left(x_i^{k+1}-x_i^k\right).\hfill \end{array}$$ Agent i sends $w_{ij}^{k+1}$, ${\overline{w}}_{ij}^k$, $C_{ij}{x}_i^{k+1}$ to all of its neighbors. $\mathbf{End}$ $\mathbf{Output}:$ The sequence ${\left(x_i^k\right)}_{k=1}^{\infty }$ to estimate the optimal solution.

4. Convergence Analysis

In this section, we first establish the compact form with operators of the proposed algorithm. Then, the results of the theoretical analysis are provided.

For $i\in \mathcal{V}$, $j\in {\mathcal{N}}_i$, we make the following definitions. Let w and $\overline{w}$ represent the variables stacked by $w_{\left(i,j\right)}$ and ${\overline{w}}_{\left(i,j\right)}$, respectively. Define vectors $u=col\left(u_1,u_2,\dots ,u_m\right),v=col\left(v_1,v_2,\dots ,v_m\right),\overline{u}=col\left({\overline{u}}_1,{\overline{u}}_2,\dots ,{\overline{u}}_m\right)$ and $\overline{v}=col\left({\overline{v}}_1,{\overline{v}}_2,\dots ,{\overline{v}}_m\right)$. Then, we let $W=\mathrm{blkdiag}{\left\{{\omega }_{\left(i,j\right)}{I}_{2n}\right\}}_{(i,j)\in \mathcal{E}}$, $\Gamma =\mathrm{blkdiag}{\left\{{\gamma }_i{I}_n\right\}}_{i\in \mathcal{V}}$, $\Sigma =\mathrm{diag}{\left\{{\sigma }_i{I}_n\right\}}_{i\in \mathcal{V}}$ and $M=\mathrm{blkdiag}{\left\{{\mu }_i{I}_n\right\}}_{i\in \mathcal{V}}$ be the stepsize matrices. Then $\mathcal{C}={\prod }_{\left(i,j\right)\in \mathcal{E}}{\mathcal{C}}_{\left(i,j\right)}$, $\mathsf{\Omega }={\prod }_{i\in \mathcal{V}}{\mathsf{\Omega }}_i$ and $\mathcal{A}={\prod }_{i\in \mathcal{V}}{\mathcal{A}}_i$ hold such that there exist ${\delta }_{\mathcal{C}}={\sum }_{\left(i,j\right)\in \mathcal{E}}{\delta }_{{\mathcal{C}}_{\left(i,j\right)}}$, ${\delta }_{\mathsf{\Omega }}={\sum }_{i\in \mathcal{V}}{\delta }_{{\mathsf{\Omega }}_i}$ and ${\delta }_{\mathcal{A}}={\sum }_{i\in \mathcal{V}}{\delta }_{{\mathcal{A}}_i}$. The linear operator $N:x\mapsto {\left(N_{\left(i,j\right)}x\right)}_{\left(i,j\right)\in \mathcal{E}}$ is stacked by $N_{\left(i,j\right)}$. Considering the resolvent of the proximal operator, the update flow (8) leads to the following equalities:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\partial {\delta }_{\mathcal{C}}^{*}\left({\overline{w}}^k\right)+W^{-1}{\overline{w}}^k=W^{-1}{w}^k+N{x}^k,\hfill \\ \partial {\delta }_{\mathsf{\Omega }}^{*}\left({\overline{u}}^k\right)+M^{-1}{\overline{u}}^k=M^{-1}{u}^k+x^k,\hfill \\ \partial {\delta }_{\mathcal{A}}^{*}\left({\overline{v}}^k\right)+{\Sigma }^{-1}{\overline{v}}^k={\Sigma }^{-1}{v}^k+x^k,\hfill \\ \partial g\left({\overline{x}}^k\right)+{\Gamma }^{-1}{\overline{x}}^k+{\overline{u}}^k+{\overline{v}}^k+N^{\top }{\overline{w}}^k={\Gamma }^{-1}{x}^k-\nabla f\left(x^k\right),\hfill \\ w^{k+1}={\overline{w}}^k+WN{\overline{x}}^k-WN{x}^k,\hfill \\ u^{k+1}={\overline{u}}^k+M{\overline{x}}^k-M{x}^k,\hfill \\ v^{k+1}={\overline{v}}^k+\Sigma {\overline{x}}^k-\Sigma x^k,\hfill \end{array}\right.\end{array}$$

(12)

where ${\overline{x}}^k=x^{k+1}$ is the auxiliary variable.

Define two variables $U=col\left(w,u,v,x\right)$ and $\overline{U}=col\left(\overline{w},\overline{u},\overline{v},\overline{x}\right)$. Based on the equalities in (12), Algorithm 1 is equivalent to the following compact form described by the operators:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\overline{U}}^k={\left(T_H+T_A\right)}^{-1}\left(T_H-T_Q-T_D\right)U^k,\hfill \\ U^{k+1}=U^k+T_S^{-1}\left(T_H-T_Q\right)\left({\overline{U}}^k-U^k\right),\hfill \end{array}\right.\end{array}$$

(13)

where the operators are given as follows:

$$T_H:U\mapsto {\left(\begin{array}{c}{\left(W^{-1}w\right)}^{\top },{\left(M^{-1}u\right)}^{\top },{\left({\Sigma }^{-1}v\right)}^{\top },{\left(N^{\top }w+u+v+{\Gamma }^{-1}x\right)}^{\top }\hfill \end{array}\right)}^{\top },$$

$$T_A:U\mapsto {\left(\partial {\delta }_{\mathcal{C}}^{*}{\left(w\right)}^{\top },\partial {\delta }_{\mathsf{\Omega }}^{*}{\left(u\right)}^{\top },\partial {\delta }_{\mathcal{A}}^{*}{\left(v\right)}^{\top },\partial g{\left(x\right)}^{\top }\right)}^{\top },$$

$$T_Q:U\mapsto {\left({\left(-Nw\right)}^{\top },-u^{\top },-v^{\top },{\left(N^{\top }w+u+v\right)}^{\top }\right)}^{\top },$$

$$T_D:U\mapsto {\left({\mathbf{0}}_n,{\mathbf{0}}_n,{\mathbf{0}}_n,\nabla f{\left(x\right)}^{\top }\right)}^{\top },$$

$$T_S:U\mapsto {\left({\left(Ww\right)}^{\top },{\left(Mu\right)}^{\top },{\left(\Sigma v\right)}^{\top },{\left(\Gamma x\right)}^{\top }\right)}^{\top }.$$

Consider one iteration of the proposed algorithm as an operator T. Then we let $U^{*}=col\left(x^{*},w^{*},u^{*},v^{*}\right)$ be the fixed point of the operator T such that $U^{*}=T{U}^{*}$. Next, we conduct the convergence analysis.

Lemma 1.

(Optimal analysis) Let Assumption 1 be satisfied. The fixed point $U^{*}$ related to the operator T meets the first-order optimal conditions of the objective problem, and $x^{*}\in \mathcal{S}$ is an optimal solution.

Proof. 

Substituting the fixed point into (12), we have the following set of equalities:

$$\left\{\begin{array}{c}\partial {\delta }_{\mathcal{C}}^{*}\left(w^{*}\right)-N{x}^{*}=0,\hfill \\ \partial {\delta }_{\mathsf{\Omega }}^{*}\left(u^{*}\right)-x^{*}=0,\hfill \\ \partial {\delta }_{\mathcal{A}}^{*}\left(v^{*}\right)-x^{*}=0,\hfill \\ \partial g\left(x^{*}\right)+\nabla f\left(x^{*}\right)+u^{*}+v^{*}+N^{\top }{w}^{*}=0,\hfill \end{array}\right.$$

which is also the KKT condition of the Lagrangian function (6). Therefore, $x^{*}$ is an optimal solution to problem (1). □

The relationship between the fixed point and the optimal solution is ensured by Lemma 1. Split the operator $T_H$ as $T_H=T_P+T_K$, where we let

$$T_P=\left[\begin{array}{cccc}{W}^{-1}& 0& 0& \frac{1}{2}N\\ 0& M^{-1}& 0& \frac{1}{2}I\\ 0& 0& {\Sigma }^{-1}& \frac{1}{2}I\\ \frac{1}{2}{N}^{\top }& \frac{1}{2}I& \frac{1}{2}I& {\Gamma }^{-1}\end{array}\right],$$

$$T_K=\left[\begin{array}{cccc}0& 0& 0& -\frac{1}{2}N\\ 0& 0& 0& -\frac{1}{2}I\\ 0& 0& 0& -\frac{1}{2}I\\ \frac{1}{2}{N}^{\top }& \frac{1}{2}I& \frac{1}{2}I& 0\end{array}\right],$$

and further define another linear operator

$$T_{\tilde{P}}=\left[\begin{array}{cccc}{W}^{-1}& 0& 0& -\frac{1}{2}N\\ 0& M^{-1}& 0& -\frac{1}{2}I\\ 0& 0& {\Sigma }^{-1}& -\frac{1}{2}I\\ -\frac{1}{2}{N}^{\top }& -\frac{1}{2}I& -\frac{1}{2}I& {\Gamma }^{-1}\end{array}\right].$$

With these definitions above, the following lemma provides the property of the operator T for convergence analysis.

Lemma 2.

Under Assumption 1, there exists the following inequality for $U^{*}$:

$${\left({\overline{U}}^k-U^{*}\right)}^{\top }{T}_D\left(U^k-U^{*}\right)\ge -\frac{1}{4}{∥x^k-{\overline{x}}^k∥}_B^2,$$

where $B=\mathrm{blkdiag}\left\{{\beta }_i{I}_n\right\}$ for $i\in \mathcal{V}$ is the Lipschitz parameter matrix.

Proof. 

With the definition of operator $T_D$, we have the equality

$$\begin{array}{cc}\hfill & {\left({\overline{U}}^k-U^{*}\right)}^{\top }{T}_D\left(U^k-U^{*}\right)\hfill \\ \hfill & =-{\left(x^k-{\overline{x}}^k\right)}^{\top }\left(\nabla f\left(x^k\right)-\nabla f\left(x^{*}\right)\right)+{\left(x^k-x^{*}\right)}^{\top }\left(\nabla f\left(x^k\right)-\nabla f\left(x^{*}\right)\right).\hfill \end{array}$$

(14)

According to [32] (Theorem 18.16), for $i\in \mathcal{V}$, $\nabla f_i$ is cocoercive, i.e., it holds

$$\begin{array}{c}\hfill {∥\nabla f\left(x^k\right)-\nabla f\left(x^{*}\right)∥}_{B^{-1}}^2\le {\left(x^k-x^{*}\right)}^{\top }\left(\nabla f\left(x^k\right)-\nabla f\left(x^{*}\right)\right).\end{array}$$

(15)

Note that for any vector a and b in the same dimension and a diagonal positive definite matrix V, then there exists the inequality $x^{\top }y\le {∥x∥}_V^2+\frac{1}{4}{∥y∥}_{V^{-1}}^2$. Hence, we have

$$\begin{array}{c}\hfill {\left(x^k-{\overline{x}}^k\right)}^{\top }\left(\nabla f\left(x^k\right)-\nabla f\left(x^{*}\right)\right)\le {∥\nabla f\left(x^k\right)-\nabla f\left(x^{*}\right)∥}_{B^{-1}}^2+\frac{1}{4}{∥x_k-{\overline{x}}_k∥}_B^2.\end{array}$$

(16)

Combining (14)–(16), we can obtain the objective inequality and end the proof. □

Lemma 3.

Under Assumption 1, there exists the following inequality for $U^{*}$:

$$\begin{array}{c}\hfill {∥U^{k+1}-U^{*}∥}_{T_S}^2-{∥U^k-U^{*}∥}_{T_S}^2\le {∥U^{k+1}-U^k∥}_{T_S-2T_{\tilde{P}}}^2+\frac{1}{2}{∥x^k-{\overline{x}}^k∥}_B^2,\end{array}$$

(17)

where $T_{\tilde{P}}$ is defined before Lemma 2.

Proof. 

Considering the change of the optimal residual before and after one iteration, we have

$$\begin{array}{cc}\hfill & {∥U^{k+1}-U^{*}∥}_{T_S}^2-{∥U^k-U^{*}∥}_{T_S}^2\hfill \\ \hfill & =-\begin{array}{c}{∥U^{k+1}-U^k∥}_{T_S}^2+2{\left(U^{k+1}-U^{*}\right)}^{\top }{T}_S\left(U^{k+1}-U^k\right)\hfill \end{array}\hfill \\ \hfill & =\begin{array}{c}{∥U^{k+1}-U^k∥}_{T_S}^2+2{\left(U^k-U^{*}\right)}^{\top }{T}_S\left(U^{k+1}-U^k\right)\hfill \end{array}.\hfill \end{array}$$

(18)

From the second step of the update flow (13), there exists

$$T_S\left(U^{k+1}-U^k\right)=\left(T_H-T_Q\right)\left({\overline{U}}^k-U^k\right),$$

such that the equality (18) leads to

$$\begin{array}{cc}\hfill & {∥U^{k+1}-U^{*}∥}_{T_S}^2-{∥U^k-U^{*}∥}_{T_S}^2\hfill \\ \hfill & ={∥U^{k+1}-U^k∥}_{T_S}^2\hfill \\ \hfill & +2{\left({\overline{U}}^k-U^{*}\right)}^{\top }\left(T_H-T_Q\right)\left({\overline{U}}^k-U^k\right)\hfill \\ \hfill & -2{\left({\overline{U}}^k-U^k\right)}^{\top }\left(T_H-T_Q\right)\left({\overline{U}}^k-U^k\right)\hfill \\ \hfill & \le {∥U^{k+1}-U^k∥}_{T_S}^2\hfill \\ \hfill & +2{\left({\overline{U}}^k-U^{*}\right)}^{\top }\left(T_H-T_Q\right)\left({\overline{U}}^k-U^k\right)\hfill \\ \hfill & -2{\left({\overline{U}}^k-U^k\right)}^{\top }{T}_P\left({\overline{U}}^k-U^k\right).\hfill \end{array}$$

From the first step of the update flow (13), it holds that

$$\left(T_H-T_Q\right)\left({\overline{U}}^k-U^k\right)=-\left(T_D{U}^k+T_Q{\overline{U}}^k+T_A{\overline{U}}^k\right).$$

Thus, we further have

$$\begin{array}{cc}\hfill & {∥U^{k+1}-U^{*}∥}_{T_S}^2-{∥U^k-U^{*}∥}_{T_S}^2\hfill \\ \hfill & \le {∥U^{k+1}-U^k∥}_{T_S}^2-2{\left({\overline{U}}^k-U^k\right)}^{\top }{T}_P\left({\overline{U}}^k-U^k\right)\hfill \\ \hfill & -2{\left({\overline{U}}^k-U^{*}\right)}^{\top }\left(T_D{U}^k+T_Q{\overline{U}}^k+T_A{U}^{*}\right)\hfill \\ \hfill & -2{\left({\overline{U}}^k-U^{*}\right)}^{\top }\left(T_A{\overline{U}}^k-T_A{U}^{*}\right).\hfill \end{array}$$

(19)

Then, we discuss the right side of (19). Note that Lemma 1 proves the equivalence between the fixed point and the optimal solution. Substituting the property of fixed points into the update flow (13), we obtain $U^{*}={\overline{U}}^{*}$ and

$$T_A{U}^{*}=-T_Q{U}^{*}-T_D{U}^{*}.$$

Hence, the third term on the right side of (19) satisfies

$$\begin{array}{cc}\hfill & {\left({\overline{U}}^k-U^{*}\right)}^{\top }\left(T_D{U}^k+T_Q{\overline{U}}^k+T_A{U}^{*}\right)\hfill \\ \hfill & ={\left({\overline{U}}^k-U^{*}\right)}^{\top }{T}_D\left(U^k-U^{*}\right)+{\left({\overline{U}}^k-U^{*}\right)}^{\top }{T}_Q\left({\overline{U}}^k-U^{*}\right)\hfill \\ \hfill & \ge -\frac{1}{4}{∥x^k-{\overline{x}}^k∥}_B^2+{\left({\overline{U}}^k-U^{*}\right)}^{\top }{T}_Q\left({\overline{U}}^k-U^{*}\right),\hfill \end{array}$$

(20)

where the inequality is based on Lemma 2. Notice that the operator $T_A$ is monotone [32] (Theorem 21.2 and Proposition 20.23), i.e., it holds

$$\begin{array}{c}\hfill {\left({\overline{U}}^k-U^{*}\right)}^{\top }\left(T_A{\overline{U}}^k-T_A{U}^{*}\right)\ge 0.\end{array}$$

(21)

Since the linear operator $T_Q$ is a skew-symmetric matrix, it is monotone [29]. Combining (19)–(21), we obtain

$$\begin{array}{cc}\hfill & {∥U^{k+1}-U^{*}∥}_{T_S}^2-{∥U^k-U^{*}∥}_{T_S}^2\hfill \\ \hfill & \le {∥U^{k+1}-U^k∥}_{T_S}^2-2{\left({\overline{U}}^k-U^k\right)}^{\top }{T}_P\left({\overline{U}}^k-U^k\right)+\frac{1}{2}{∥x^k-{\overline{x}}^k∥}_B^2.\hfill \end{array}$$

(22)

From the second step of the update flow (13), it holds

$${\overline{U}}^k-U^k={\left(T_H-T_Q\right)}^{-1}{T}_S\left(U^{k+1}-U^k\right),$$

where $T_H$, $T_Q$ and $T_S$ are the linear operators. Considering that $T_P$ is also a linear operator, the second term on the right side of (22) has an equivalent form:

$$\begin{array}{c}\hfill {\left({\overline{U}}^k-U^k\right)}^{\top }{T}_P\left({\overline{U}}^k-U^k\right)={\left(U^{k+1}-U^k\right)}^{\top }{T}_{\tilde{P}}\left(U^{k+1}-U^k\right).\end{array}$$

(23)

Substituting (23) into (22), we complete the proof. □

Summarizing the above lemmas, the following theorem supplies the convergence results.

Theorem 1.

When Assumption 1 and 2 are satisfied, for the sequence ${\left(U^k\right)}_{k\ge 0}$ generated by the operator T, we have

$$\begin{array}{c}\hfill {∥U^{k+1}-U^{*}∥}_{T_s}^2-{∥U^k-U^{*}∥}_{T_s}^2\le -{∥U^{k+1}-U^k∥}_{2\left(T_{\tilde{P}}-\tilde{B}\right)-T_S}^2,\end{array}$$

(24)

where $\tilde{B}=\mathrm{blkdiag}\left\{{\mathbf{0}}_{n\times n},{\mathbf{0}}_{n\times n},{\mathbf{0}}_{n\times n},\frac{1}{4}B\right\}$. Then, the sequence ${\left(U^k\right)}_{k\ge 0}$ has sublinear rate $O\left(1/k\right)$, and the sequence ${\left(x^k\right)}_{k\ge 0}$ converges to an optimal solution $x^{*}\in \mathcal{S}$.

Proof. 

With the definition of $\tilde{B}$, we have the following equality:

$$\begin{array}{c}\hfill \frac{1}{4}{∥x^k-{\overline{x}}^k∥}_B^2={∥U^{k+1}-U^k∥}_{\tilde{B}}^2.\end{array}$$

(25)

Substituting (25) into (17), we obtain the inequality (24). Note that under Assumption 2, the matrix $2\left(T_{\tilde{P}}-\tilde{B}\right)-T_S$ is positive definite. Hence, the sequence ${\left(U^k\right)}_{k\ge 0}$ converges to the fixed point $U^{*}$. Meanwhile, utilizing [10] (Proposition 1) results in the $O\left(1/k\right)$ rate, and based on Lemma 1, the convergence of ${\left(x^k\right)}_{k\ge 0}$ holds. □

In Theorem 1, the positive definite property is needed for the induced matrices, which leads to the stepsize conditions in Assumption 2.

5. Numerical Simulation

The correctness of the theoretical analysis is verified through numerical simulation on a constrained optimization problem over networks in this section.

The constrained quadratic programming problem [33] is considered in the experiments, which has the formulation as follows:

$$\begin{array}{cc}\hfill & \underset{\tilde{x}}{min}\sum _{i=1}^m{∥\tilde{x}∥}_{E_i}^2+e_i^{\top }\tilde{x}+{\rho }_i{∥\tilde{x}∥}_1,\hfill \end{array}$$

(26a)

$$\begin{array}{cc}\hfill & s.t.A_i\tilde{x}=b_i,i\in \mathcal{V},\hfill \end{array}$$

(26b)

$$\begin{array}{cc}\hfill & {\tilde{x}}^{\min }\le \tilde{x}\le {\tilde{x}}^{\max },i\in \mathcal{V},\hfill \end{array}$$

(26c)

where matrix $E_i\in {\mathbb{R}}^{n\times n}$ is diagonal and positive definite, $e_i\in {\mathbb{R}}^n$ is a vector, and ${\rho }_i$ is the penalty factor. Both $x_i^{\min }$ and $x_i^{\max }$ are vectors with constants, which give the bounds of the decision variable $\tilde{x}$. In the light of (1), we can set $f_i\left(\tilde{x}\right)={∥\tilde{x}∥}_{E_i}^2+e_i^{\mathrm{T}}\tilde{x}$ and $g_i\left(\tilde{x}\right)={\rho }_1{∥\tilde{x}∥}_1$.

In this case, the dimension of the decision variable is set as $n=4$, and we let $r=1$. For $i\in \mathcal{V}$, the paramount data of Problem (26) are selected randomly. The elements of matrix $E_i$ are in $[1,2]$, and the elements of the linear operator $A_i$ are in $[1,15]$. Both vectors $e_i$ and $b_i$ take values in $[-5,5]$. The box constraints are considered as $[-2.5,2.5]$. Then, we set the uncoordinated setpsizes randomly as ${\gamma }_i\in [0.005,0.006]$, while ${\sigma }_i$, ${\mu }_i$ and ${\omega }_{(i,j)}$ are in $[5,6]$. The numerical experiments are performed over the generated network with eight agents, which is displayed in Figure 1. The simulations are carried by running the distributed algorithms on a laptop with Intel(R) Core i5-5500U CPU @ 2.40 GHz, 8.0 GB of RAM, and Matlab R2016a on Windows 10 operating system.

The simulation results are shown in Figure 2 and Figure 3. The transient behaviors of each component of $x_i^k$ is displayed in Figure 2, in which a node-based consensus algorithm [34] is introduced as a comparative profile. Note that the obtained optimal solution from the proposed algorithm is in line with that of the node-based consensus one, i.e., $x^{*}={\left[0.6900,0.6270,0.8046,0.4400\right]}^{\mathrm{T}}$, but the latter achieves a stable consensus after 15,000 iterations. Figure 3 shows that our proposed algorithm outperforms the node-based and subgradient algorithms [35] in terms of convergence performance by evaluating the relative errors ${\sum }_{i=1}^m\parallel x_i^k-{\tilde{x}}^{*}\parallel /m\parallel {\tilde{x}}^{*}\parallel$.

6. Conclusions

In this paper, a distributed algorithm based on proximal operators has been designed to deal with discussed a class of distributed composite optimization problems, in which the local function has a smooth and nonsmooth structure and the decision variable abides by both affine and feasible constraints. Distinguishing attributes of the proposed algorithm include the use of uncoordinated stepsizes and the edge-based communication that avoids the dependency on Laplacian weight matrices. Meanwhile, the algorithm has been verified in theory and simulation. However, there are still some aspects worthy of improvement in this paper. For example, it is worth adopting efficient accelerated protocols (such as the Nesterov-based method and heavy ball method) to improve the convergence rate and developing asynchronous distributed algorithms to deal with the issue of communication latency. In addition, more general optimization models and more efficient algorithms should be investigated in order to address potential applications, e.g., [36,37,38] with nonconvex objectives, coupled and nonlinear constraints.