A General Framework of Exact Primal-Dual First Order Algorithms for Distributed Optimization

TL;DR

This paper introduces a flexible distributed primal-dual optimization framework that achieves exact convergence with fixed stepsizes, balancing speed and computational effort, and demonstrates linear convergence for strongly convex problems.

Contribution

It proposes a general class of distributed primal-dual algorithms with multiple primal steps per iteration, enabling controlled trade-offs and exact convergence with fixed stepsizes.

Findings

Achieves linear convergence for strongly convex, Lipschitz gradient functions.

Outperforms existing methods in simulation tests.

Allows multiple primal steps to balance performance and complexity.

Abstract

We study the problem of minimizing a sum of local objective convex functions over a network of processors/agents. This problem naturally calls for distributed optimization algorithms, in which the agents cooperatively solve the problem through local computations and communications with neighbors. While many of the existing distributed algorithms with constant stepsize can only converge to a neighborhood of optimal solution, some recent methods based on augmented Lagrangian and method of multipliers can achieve exact convergence with a fixed stepsize. However, these methods either suffer from slow convergence speed or require minimization at each iteration. In this work, we develop a class of distributed first order primal-dual methods, which allows for multiple primal steps per iteration. This general framework makes it possible to control the trade off between the performance and the execution complexity in primal-dual algorithms. We show that for strongly convex and Lipschitz gradient objective functions, this class of algorithms converges linearly to the optimal solution under appropriate constant stepsize choices. Simulation results confirm the superior performance of our algorithm compared to existing methods.

Equations

$${\tilde{x}}^{\ast }:={\mathrm{argmin}}_{\tilde{x}\in \mathbb{R}}\sum _{i=1}^n{f}_i(\tilde{x})$$

$$\underset{x_i\in \mathbb{R}}{\min }\sum _{i=1}^n{f}_i(x_i)\text{mbox}s.t{x}_i=x_j\forall (i,j)\in \mathcal{E}.$$

$$\underset{x}{\min }f(x)=\sum _{i=1}^n{f}_i(x_i)\text{mbox}s.tAx=0,$$

$$L(x,\lambda )=f(x)+{\lambda }^{\prime }Ax+\frac{1}{2}{x}^{\prime }Bx,$$

$$x^{k+1}=x^k-\alpha {\nabla }_{x}L(x^k,{\lambda }^k)$$

$$=x^k-\alpha \nabla f(x^k)-\alpha A^{\prime }{\lambda }^k-\alpha B{x}^k,$$

$${\lambda }^{k+1}={\lambda }^k+\beta {\nabla }_{\lambda }L(x^{k+1},{\lambda }^k)={\lambda }^k+\beta A{x}^{k+1},$$

$$x_i^{k+1,0}=x_i^k$$

$$x_i^{k+1,t}=x_i^{k+1,t-1}-\alpha \nabla f_i(x_i^k)-\alpha \sum _{l=1}^e{A}_{il}^{\prime }{\lambda }_l^k$$

$$-\alpha \sum _{j=1}^n{B}_{ij}{x}_j^{k+1,t-1}$$

$$x_i^{k+1}=x_i^{k+1,T}$$

$${\lambda }_{l_i}^{k+1}={\lambda }_{l_i}^k+\beta \sum _{j=1}^n{A}_{l_{i}j}{x}_j^{k+1}\forall {\lambda }_{l_i}\in {\Lambda }_i.$$

$$x^{k+1}=(I-\alpha B{)}^T{x}^k-\alpha C\nabla f(x^k)-\alpha C{A}^{\prime }{\lambda }^k,$$

$${\lambda }^{k+1}={\lambda }^k+\beta A{x}^{k+1},$$

$$M=C^{-1}(I-\alpha B{)}^T\text{mbox}andN=\frac{1}{\alpha }(C^{-1}-M).$$

$$C=\sum_{t=0}^{T-1}(I-\alpha B)^{t}=V\big{(}\sum_{t=0}^{T-1}Z^{t}\big{)}V^{\prime}=VQV^{\prime},$$

$$\frac{(1-\alpha \rho (B){)}^T}{{\sum }_{t=0}^{T-1}(1-\alpha \rho (B){)}^t}I⪯M⪯\frac{1}{T}I.$$

$$C^{-1}-M=V{Q}^{-1}{V}^{\prime }-VP{V}^{\prime }=V(Q^{-1}-P)V^{\prime },$$

$$\alpha (\nabla f(x^k)-\nabla f(x^{\ast }))=M(x^k-x^{k+1})+$$

$$\alpha (\beta A^{\prime }A-N)(x^{k+1}-x^{\ast })-\alpha A^{\prime }({\lambda }^{k+1}-{\lambda }^{\ast }),$$

$$\alpha C\nabla f(x^k)=(I-\alpha B{)}^T{x}^k-x^{k+1}-\alpha C{A}^{\prime }{\lambda }^k.$$

$$\alpha C\nabla f(x^k)=(I-\alpha B{)}^T{x}^k-x^{k+1}$$

$$-\alpha C{A}^{\prime }({\lambda }^{k+1}-\beta A{x}^{k+1})=(I-\alpha B{)}^T(x^k-x^{k+1})+$$

$$\displaystyle\big{(}\alpha\beta CA^{\prime}A-I+(I-\alpha B)^{T}\big{)}x^{k+1}-\alpha CA^{\prime}\lambda^{k+1},$$

$$\alpha C\nabla f(x^{*})=\big{(}\alpha\beta CA^{\prime}A-I+(I-\alpha B)^{T}\big{)}x^{*}-\alpha CA^{\prime}\lambda^{{}^{*}}.$$

$$\alpha (\nabla f(x^k)-\nabla f(x^{\ast }))=C^{-1}(I-\alpha B{)}^T(x^k-x^{k+1})$$

$$\displaystyle+\alpha\Big{(}\beta A^{\prime}A-\frac{1}{\alpha}\big{(}C^{-1}-C^{1}(I-\alpha B)^{T}\big{)}\Big{)}(x^{k+1}-x^{*})$$

$$-\alpha A^{\prime }({\lambda }^{k+1}-{\lambda }^{\ast }).$$

$$(x^{k+1}-x^{\ast }{)}^{\prime }(2\alpha mI-\alpha \eta I+2\alpha (N-\beta A^{\prime }A))(x^{k+1}-x^{\ast })+$$

$$\frac{\alpha }{\beta }{\left|\left|{\lambda }^{k+1}-{\lambda }^k\right|\right|}^2+(x^{k+1}-x^k{)}^{\prime }(M-\frac{\alpha L^2}{\eta }I)(x^{k+1}-x^k)$$

$$\le (x^k-x^{\ast }{)}^{\prime }M(x^k-x^{\ast })-(x^{k+1}-x^{\ast }{)}^{\prime }M(x^{k+1}-x^{\ast })$$

$$+\frac{\alpha }{\beta }\left({\left|\left|{\lambda }^k-{\lambda }^{\ast }\right|\right|}^2-{\left|\left|{\lambda }^{k+1}-{\lambda }^{\ast }\right|\right|}^2\right).$$

$$2\alpha m{\left|\left|x^{k+1}-x^{\ast }\right|\right|}^2\le$$

$$2\alpha (x^{k+1}-x^{\ast }{)}^{\prime }(\nabla f(x^{k+1})-\nabla f(x^{\ast }))$$

$$=2\alpha (x^{k+1}-x^{\ast }{)}^{\prime }(\nabla f(x^{k+1})-\nabla f(x^k))$$

$$+2\alpha (x^{k+1}-x^{\ast }{)}^{\prime }(\nabla f(x^k)-\nabla f(x^{\ast })),$$

$$2\alpha m{\left|\left|x^{k+1}-x^{\ast }\right|\right|}^2\le 2\alpha (x^{k+1}-x^{\ast }{)}^{\prime }\times$$

Full text

I Introduction

This paper focuses on solving the following optimization problem

[TABLE]

over a network of $n$ agents (processors)111For representation simplicity, we consider the case with $\tilde{x}\in\mathbb{R}$. Our analysis can be easily generalized to the multidimensional case.. The agents are connected with an undirected static graph $\mathcal{G(V,E)}$, with $\mathcal{V}$ and $\mathcal{E}$ being the set of vertices and edges respectively. Set $\mathcal{N}_{i}$ denotes the neighbors of agent $i$. Each agent $i$ in the network has access to a real-valued convex local objective function $f_{i}$, and can only communicate with its neighbors defined by the graph.

Problems of form (1) that require distributed optimization algorithms are applicable to many important areas such as sensor networks, smart grids, robotics, and machine learning [26, 29, 13, 33]. The distributed equivalent of problem (1) can be formulated by defining local copies of the decision variable at each agent and rewriting the problem as the following consensus problem

[TABLE]

In this distributed setting, each agent balances the goal of minimizing its local objective function and ensuring that its decision variable is equal to all of its neighbors.

One main category of distributed algorithms for solving problem (2) is first order primal methods that can parallelize computations across multiple processors. These methods include different variations of distributed (sub)gradient descent (DGD) [22, 27, 25, 28, 21, 7, 12]. In DGD-based methods, each agent updates its local estimate of the solution through a combination of a local gradient decent step and a consensus step (weighted average with neighbors variables). In order to converge to the exact solution, these algorithms need to use a diminishing stepsize, which results in a slow rate of convergence. A faster convergence to a neighborhood of the exact solution can be achieved by using a fixed stepsize. Recently, a class of algorithms based on DGD has been developed in [1], which uses increasing number of consensus steps (linear increase with iteration number) to achieve exact linear convergence with fixed stepsize.

Another line of distributed optimization is inspired by Lagrange multiplier methods for constrained optimization. Specifically, Method of Multipliers (MM), based on augmented Lagrangian, has nice convergence guarantees in centralized setting [9, 3]. Each iteration of MM includes minimizing the augmented Lagrangian in the primal space, followed by a dual variable update. This method is not implementable in a distributed setting due to the nonseparable augmentation term. Decentralized versions of augmented Lagrangian methods and alternating direction method of multiplier (ADMM) are studied in [2, 5, 35, 20, 36, 38, 32, 14, 18, 10, 11, 37]. Although these algorithms involve more computational complexity compared to primal methods, they guarantee convergence to the exact solution. Specifically algorithms in [32, 14, 18, 24] are shown to converge linearly under the assumptions of strong convexity and Lipschitz gradient. Some other algorithms like EXTRA [31], DIGing [23], and DSA [17] can achieve exact convergence with a fixed stepsize without explicitly using dual variables. These algorithms, however, can be viewed as augmented Lagrangian primal-dual methods with a single gradient step in the primal space. These methods enjoy less computational complexity compared to ADMM-based methods and are proved to converge with a linear rate, under standard assumptions.

Besides these first order methods, second order distributed algorithms are studied in [16, 19, 34, 8, 15]. These methods use Newton-type updates to achieve faster convergence and thus involve more computational complexity than first order methods.

While existing distributed first order primal-dual algorithms can achieve exact convergence, they either suffer from slower convergence compared to MM or they require minimization at each iteration. In this paper we develop a class of distributed first order primal-dual methods, which does not require exact minimization at each iteration but allows for multiple primal steps per iteration. This general framework makes it possible to control the trade off between the performance and the execution complexity in primal-dual algorithms.

We develop our algorithm based on a general form of augmented Lagrangian, which is flexible in the augmentation term, and prove global linear rate of convergence with a fixed stepsize, under the assumptions of strong convexity and Lipschitz gradient. Our framework includes EXTRA and DIGing as special cases for specific choices of augmented Lagrangian function and one primal step per iteration.

The rest of this paper is organized as follows: Section II describes the development of our general framework. Section III contains the convergence analysis. Section IV presents the numerical experiments and Section V contains the concluding remarks.

II Algorithm Development

To develop our algorithm, we express problem (2) in the following compact form

[TABLE]

where $x=[x_{1},x_{2},...,x_{n}]^{\prime}\in\mathbb{R}^{n}$, 222For a vector $x$ we denote its transpose by $x^{\prime}$ and for a matrix $A$ we denote its transpose by $A^{\prime}$. and $Ax=0$ represents all equality constraints. Matrix $A\in\mathbb{R}^{e\times n}$, $e=|\mathcal{E}|$, is the edge-node incidence matrix of the network graph and its null space is spanned by a vector of all ones. Row $l$ of matrix $A$ corresponds to edge $l$, connecting vertices $i$ and $j$, and has $+1$ in column $i$ and $-1$ in column $j$ (or vice versa) and [math] in all other columns [4].

Remark 1

Other choices for matrix $A$ include weighted incidence matrix [37], graph Laplacian [34] matrix, and weighted Laplacian matrix [19, 1].

We denote by $x^{*}=[\tilde{x}^{*},\tilde{x}^{*},...,\tilde{x}^{*}]^{\prime}$ the minimizer of problem (3) and we focus on developing a distributed algorithm which converges to $x^{*}$. In order to achieve exact convergence, we develop our algorithm based on the Lagrange multiplier methods. We form the following augmented Lagrangian

[TABLE]

where $\lambda\in\mathbb{R}^{e}$ is the vector of Lagrange multipliers. Every dual variable $\lambda_{l}$ is associated with an edge $l=(i,j)$ and thus coupled between two agents. We assume that one of the two agents $i$ or $j$ choose to update the dual variable $\lambda_{l}$ through negotiation. The set of dual variables that agent $i$ updates is denoted by $\Lambda_{i}$ . We adopt the following assumptions on our problem.

Assumption 1

The local objective functions $f_{i}(x)$ are $m-$ strongly convex, twice differentiable, and $L-$ Lipschitz gradient.

Assumption 2

Matrix $B\in\mathbb{R}^{n\times n}$ is a symmetric positive semidefinite matrix, has the same null space as matrix $A$, i.e., $Bx=0$ only if $Ax=0$, and is compatible with network topology, i.e., $B_{ij}\neq 0$ only if $(i,j)\in\mathcal{E}$.

We assume these conditions hold for the rest of the paper. The first assumption requires the eigenvalues of the Hessian matrix of local objective functions to be bounded, i.e., $mI\preceq\nabla^{2}f_{i}(x)\preceq LI$ for all $i$. This is a standard assumption in proving global linear rate of convergence. The second assumption requires matrix $B$ to represent the network topology, which is needed for distributed implementation. The other assumptions on matrix $B$ are required for convergence guarantees. Some examples of matrix $B$ are $B=A^{\prime}A=L$, with $L$ being the graph Laplacian matrix, and weighted Laplacian matrix.

One important strand of Lagrange multiplier methods is based on the Method of Multipliers (MM), which requires solving a minimization problem in the primal space, followed by a dual update at each iteration. Despite the nice convergence properties of MM (global linear convergence for strongly convex objective functions), its primal update is computationally costly (because of the exact minimization) and cannot be implemented in a distributed manner (because $x^{\prime}Bx$ is not separable). An alternative is to replace the exact minimization in primal space with a single gradient step, which results in the following primal-dual update

[TABLE]

[TABLE]

where $\alpha$ and $\beta$ are constant stepsize parameters. The above iterate is implementable in a distributed manner because $\nabla f(x)$ is separable and thus locally available at agents and the terms $\lambda^{\prime}A$, $Bx$ and $Ax$ can be computed at each agent through communication with neighbors (because matrices $A$ and $B$ represent the graph topology). Although different variations of the above iteration have linear convergence guarantees [31, 23], the linear convergence of MM is shown to be faster [19]. One interesting question is whether by having multiple updates in the primal space, the performance of the algorithm gets closer to MM. By looking carefully through Eq. (5), we notice that each primal update at each agent requires one gradient evaluation and one round of communication with the neighbors.

A follow up question is whether in an algorithm with multiple primal updates per iteration, it is necessary or efficient to recompute the gradient or recommunicate with neighbors at each update.

Motivated by these questions, we present an improved class of distributed primal-dual algorithms for solving problem (3), which allows for multiple primal updates at each iteration. To avoid computational complexity, in our algorithm, the gradient is evaluated once per iteration and is used for all primal updates in that iteration. In our proposed framework in Algorithm 1, at each iteration $k$, agent $i$ computes its local gradient $\nabla f_{i}(x_{i}^{k})$, and performs a predetermined number ($T$) of primal updates by repeatedly communicating with neighbors without recomputing its gradient. Each agent $i$ then updates its corresponding dual variables $\lambda_{l_{i}}^{k}$ by using local $x_{i}^{k+1,T}$ and $x_{j}^{k+1,T}$ from its neighbors.

Under Assumption 1, there exists a unique optimal solution $\tilde{x}^{*}$ for problem (1) and thus a unique $x^{*}$ exists, at which the function value is bounded. Moreover, since $Null(A)\neq\emptyset$, the Slater’s condition is satisfied. Consequently, strong duality holds and a dual optimal solution $\lambda^{*}$ exists. We note that the projection of $\lambda^{k}$ in the null space of matrix $A^{\prime}$ would not affect the performance of algorithm, and if the algorithm starts at $\lambda=0$, then all the iterates $\lambda^{k}$ are in the column space of $A$ and hence orthogonal to null space of $A^{\prime}$. Hence, the optimal dual solution is not uniquely defined, since for any optimal dual solution $\lambda^{*}$, the dual solution $\lambda^{*}+u$, where $u$ is in the null space of $A^{\prime}$ is also optimal. Without loss of generality, we assume that in Algorithm 1 $\lambda^{0}=0$, and when we refer to an optimal dual solution $\lambda^{*}$, we assume its projection onto the null space of $A^{\prime}$ is [math]. We note that $(x^{*},\lambda^{*})$ is a fixed point of Algorithm 1.

III Convergence Analysis

In order to analyze the convergence properties of our algorithm, we first express the primal-dual update in Algorithm 1 in the following compact form

[TABLE]

[TABLE]

where $I$ denotes the identity matrix and $C=\sum_{t=0}^{T-1}(I-\alpha B)^{t}.$ We next proceed to prove the linear convergence rate for our proposed framework. In Lemmas III.1 and III.2 we establish some key relations which we use to derive two fundamental inequalities in Lemmas III.3 and III.4. Finally we use these key inequalities to prove the global linear rate of convergence in Theorem III.5. In the following analysis, we denote by $\rho(\cdot)$ the largest eigenvalue of a symmetric positive semidefinite matrix and we define matrices $M$ and $N$ as follows

[TABLE]

In the next lemma we show that matrix $C$ is invertible and thus matrices $M$ and $N$ are well-defined.

Lemma III.1

Consider the symmetric positive semidefinite matrix $B$ and matrices $C$, $M$, and $N$. If we choose $\alpha$ such that $I-\alpha B$ is positive definite, i.e., $\alpha<\frac{1}{\rho(B)}$, then matrix $C$ is invertible and symmetric, matrix $N$ is symmetric positive semidefinite, and matrix $M$ is symmetric positive definite with $\frac{(1-\alpha\rho(B))^{T}}{\sum_{t=0}^{T-1}(1-\alpha\rho(B))^{t}}I\preceq M\preceq\frac{1}{T}I.$

Proof:

Since $I-\alpha B$ is symmetric, it can be written as $I-\alpha B=VZV^{\prime}$, where $V\in\mathbb{R}^{n\times n}$ is an orthonormal matrix, i.e., $VV^{\prime}=I$, whose $i^{th}$ column $v_{i}$ is the eigenvector of $(I-\alpha B)$ and $v_{i}^{\prime}v_{t}=0$ for $i\neq t$ and $Z$ is the diagonal matrix whose diagonal elements, $Z_{ii}=\mu_{i}>0$, are the corresponding eigenvalues. We also note that since $V$ is an orthonormal matrix, $(I-\alpha B)^{t}=VZ^{t}V^{\prime}$. We have

[TABLE]

Hence matrix $C$ is symmetric. We note that matrix $Q$ is a diagonal matrix with $Q_{ii}=1+\sum_{t=1}^{T-1}\mu_{i}^{t}$. Since $\mu_{i}>0$ for all $i$, $Q_{ii}\neq 0$ and thus $Q$ is invertible and we have $C^{-1}=VQ^{-1}V^{\prime}.$ We also have $M=C^{-1}(I-\alpha B)^{T}=VQ^{-1}V^{\prime}VZ^{T}V^{\prime}=VQ^{-1}Z^{T}V^{\prime}=VPV^{\prime}$ where $P$ is a diagonal matrix with $P_{ii}=\frac{\mu_{i}^{T}}{1+\sum_{t=1}^{T-1}\mu_{i}^{t}}$, consequently, matrix $M$ is symmetric. We next find the smallest and largest eigenvalues of matrix $M$. We note that since $P_{ii}$ is increasing in $\mu_{i}$, the smallest and largest eigenvalues of $M$ can be computed using the smallest and largest eigenvalues of $I-\alpha B$. We have $0\preceq B\preceq\rho(B)I$, where $\rho(B)$ is the largest eigenvalue of matrix $B$. Therefore, the largest and smallest eigenvalues of $I-\alpha B$ are $1$ and $1-\alpha\rho(B)$ respectively. Hence,

[TABLE]

We next use the eigenvalue decomposition of matrices $C^{-1}$ and $M$ to obtain

[TABLE]

where $Q^{-1}-P$ is a diagonal matrix, and its $i^{th}$ diagonal element is equal to $\frac{1-\mu_{i}^{T}}{1+\sum_{t=1}^{T-1}\mu_{i}^{t}}$. Since $0<\mu_{i}\leq 1$ for all $i$, $\frac{1-\mu_{i}^{T}}{1+\sum_{t=1}^{T-1}\mu_{i}^{t}}\geq 0$ and hence $N$ is symmetric positive semidefinite.∎

Lemma III.2

Consider the primal-dual iterates as in Algorithm 1 and recall the definitions of matrices $M$ and $N$ from Eq. (8), if $\alpha<\frac{1}{\rho(B)}$, then

[TABLE]

Proof:

At each iteration, from Eq. (6) we have

[TABLE]

Moreover, from Eq. (7) we have $\lambda^{k}=\lambda^{k+1}-\beta Ax^{k+1}.$ We can substitute this expression for $\lambda^{k}$ into the previous equation and have

[TABLE]

where we added and subtracted a term of $(I-\alpha B)^{T}x^{k+1}$. Since an optimal solution pair $(x^{*},\lambda^{*})$ is a fixed point of the algorithm update, we also have

[TABLE]

We then subtract the above inequality from Eq. (9) and multiply both sides by $C^{-1}$ [c.f. Lemma III.1], to obtain

[TABLE]

∎

Lemma III.3

Consider the primal-dual iterates as in Algorithm 1 and recall the definition of matrices $M$ and $N$ from Eq. (8). If $\alpha<\frac{1}{\rho(B)}$, we have for any $\eta>0$,

[TABLE]

Proof:

From strong convexity of function $f(x)$, we have

[TABLE]

where we add and subtract a term $(x^{k+1}-x^{*})^{\prime}\nabla f(x^{k})$. We can substitute the equivalent expression of $\alpha(\nabla f(x^{k})-\nabla f(x^{*}))$ from Lemma III.2 and have

[TABLE]

From Young’s inequality we have $2\alpha(x^{k+1}-x^{*})^{\prime}(\nabla f(x^{k+1})-\nabla f(x^{k}))\leq\alpha\eta\left|\left|x^{k+1}-x^{*}\right|\right|^{2}+\frac{\alpha}{\eta}\left|\left|\nabla f(x^{k+1})-\nabla f(x^{k})\right|\right|^{2}$ for all $\eta>0$. This can be further reduced to $2\alpha(x^{k+1}-x^{*})^{\prime}(\nabla f(x^{k+1})-\nabla f(x^{k}))\leq\alpha\eta\left|\left|x^{k+1}-x^{*}\right|\right|^{2}+\frac{\alpha L^{2}}{\eta}\left|\left|x^{k+1}-x^{k}\right|\right|^{2}$, by Lipschitz gradient property of function $f$. By dual update Eq. (7) and feasibility of $x^{*}$, we have

[TABLE]

These two equations combined yields

[TABLE]

Hence we can rewrite the right hand side of Eq. (11) as

[TABLE]

We now focus on the last two terms. First, since matrix $M$ is symmetric, we have

[TABLE]

Similarly, we have $-2\frac{\alpha}{\beta}(\lambda^{k+1}-\lambda^{k})^{\prime}(\lambda^{k+1}-\lambda^{*})=$

[TABLE]

Now we combine the terms in the preceding three relations and have

[TABLE]

By rearranging the terms in the above inequality, we complete the proof. ∎

Lemma III.4

Consider the primal-dual iterates as in Algorithm 1 and recall the definition of symmetric matrices $M$ and $N$ from Eq.(8) then if $\alpha<\frac{1}{\rho(B)}$, for $d,g,e>1$ we have

[TABLE]

with $s(AA^{\prime})$ being the smallest nonzero eigenvalue of the positive semidefinite matrix $AA^{\prime}$.

Proof:

We recall the following relation from Lemma III.2

[TABLE]

Thus we have

[TABLE]

We can add and subtract a term of $\nabla f(x^{k+1})$ and equivalently express the above relation as

[TABLE]

By using the fact that $(a+b)^{\prime}(a+b)\leq\frac{\xi}{\xi-1}a^{\prime}a+\xi b^{\prime}b$, for any vectors $a$, $b$, and scalar $\xi>1$, we have for any scalars $d,e,g>1$,

[TABLE]

Since $\lambda^{0}=0$ and $\lambda^{k+1}=\lambda^{k}+\beta Ax^{k+1}$, we have that $\lambda^{k}$ is in the column space of $A$ and hence orthogonal to the null space of $A^{\prime}$, therefore, we have $\left|\left|\alpha A^{\prime}(\lambda^{k+1}-\lambda^{*})\right|\right|^{2}\geq\alpha^{2}s(AA^{\prime})\left|\left|\lambda^{k+1}-\lambda^{*}\right|\right|^{2}.$ Using the above two relations and Lipschitz gradient property of function $f$ , we have

[TABLE]

By using the facts that $M^{2}\preceq\rho(M^{2})I$ and (\beta A^{\prime}A-N)^{2}\preceq\rho\big{(}(\beta A^{\prime}A-N)^{2}\big{)}I, we complete the proof. ∎

Theorem III.5

Consider the primal-dual iterates as in Algorithm 1, recall the definition of matrix $M$ from Eq. (8), and define $z^{k}=\begin{bmatrix}x^{k}\\ \lambda^{k}\end{bmatrix},\quad G=\begin{bmatrix}M&\textbf{0}\\ \textbf{0}&\frac{\alpha}{\beta}I\end{bmatrix}.$ If Assumption 1 holds true, and the primal and dual stepsizes satisfy

[TABLE]

[TABLE]

with $0<\eta<2m-\beta\rho(A^{\prime}A)$, then there exists a $\delta>0$ such that

[TABLE]

that is $\left|\left|z^{k}-z^{*}\right|\right|_{G}$ converges Q-linearly to [math] and consequently $\left|\left|x^{k}-x^{*}\right|\right|_{M}$ converges R-linearly to [math].

Proof:

To show linear convergence, we will show that

[TABLE]

for some $\delta>0$. By using the result of Lemma III.3 [c.f. Eq. (10)], it suffices to show that there exists a $\delta>0$ such that

[TABLE]

We now use Eq. (7) together with the fact that $Ax^{*}=0$ to obtain $\left|\left|\lambda^{k+1}-\lambda^{k}\right|\right|^{2}=\beta^{2}(x^{k+1}-x^{*})^{\prime}(A^{\prime}A)(x^{k+1}-x^{*}),$ which we can substitute into the previous inequality and have that we need to show

[TABLE]

We collect the terms and we will focus on showing

[TABLE]

We compare Eq. (LABEL:ineq:linearRateToShow) with Eq. (LABEL:ineq:toCompare) [c.f. Lemma III.4], and we need to have for some $\delta>0$

[TABLE]

This is satisfied if

[TABLE]

Since the previous two inequalities holds for all $e,d,g>1$, we can find the parameters $e$ and $g$ to make the right hand sides the smallest, which would give us the most freedom to choose algorithm parameters. The term $\frac{e}{e-1}\rho(M)^{2}+e\alpha^{2}L^{2}$ is convex in $e$ and to minimize it we set derivative to [math] and have $e=1+\frac{\rho(M)}{\alpha L}.$ Similarly, we choose $g$ to be g=1+\frac{\sqrt{\rho\big{(}(\beta A^{\prime}A-N)^{2}\big{)}}}{L}. With these parameter choices, we have $\left(\frac{e}{e-1}\rho(M)^{2}+e\alpha^{2}L^{2}\right)=(\rho(M)+\alpha L)^{2},$ and \left(\frac{g}{g-1}\alpha^{2}\rho\big{(}(\beta A^{\prime}A-N)^{2}\big{)}+g\alpha^{2}L^{2}\right)=\alpha^{2}\Big{(}\sqrt{\rho\big{(}(\beta A^{\prime}A-N)^{2}\big{)}}+L\Big{)}^{2}. The desired relations can now be expressed as \frac{\beta}{\delta\alpha}(M-\frac{\alpha L^{2}}{\eta}I)\succcurlyeq\frac{d}{\alpha^{2}s(AA^{\prime})}\big{(}\rho(M)+\alpha L\big{)}^{2}I and \frac{\beta}{\delta\alpha}(2\alpha m-\alpha\eta I+2\alpha N-\alpha\beta A^{\prime}A-\delta M)\succcurlyeq\frac{d}{(d-1)s(AA^{\prime})}\Big{(}\sqrt{\rho\big{(}(\beta A^{\prime}A-N)^{2}\big{)}}+L\Big{)}^{2}I. By using the fact that $N$ and $A^{\prime}A$ are positive semidefinite matrices, and the result of Lemma III.1 to bound eigenvalues of matrix $M$, the desired relations can be satisfied if

[TABLE]

For the first inequality, we can multiply both sides by $\delta\alpha$ and rearrange the terms to have

[TABLE]

We can similarly solve for the second inequality,

[TABLE]

We next show that the upper bounds on $\delta$ given in Eq. (16) and Eq. (17) are both positive. Since the dual step size $\beta$ satisfies Eq. (13) and $0<\eta<2m-\beta\rho(A^{\prime}A)$, the right hand side of Eq. (17) is positive.

Moreover, in order for the right hand side of Eq. (16) to be positive we need

[TABLE]

Since $1-\alpha\rho(B)\neq 1$ [c.f. Eq. (14)], we have

[TABLE]

Therefore inequality (18) can be written as

[TABLE]

which holds true for $\alpha$ satisfying Eq. (14).

Hence, the parameter set is nonempty and thus we can find

[TABLE]

which establishes linear rate of convergence. ∎

Remark 2

If we choose $B=\beta A^{\prime}A$, we have $N-\beta A^{\prime}A=0$, and from the analysis of the above theorem we can see that the upper bound on $\beta$ can be removed. In this case, by choosing $T=1$, the stepsize upper bound in Eq. (14) and the linear rate parameter recover those of EXTRA[31].

Remark 3

To find an optimal value for the number of primal updates per iteration, $T$, leading to the best convergence rate, we can optimize over various parameters in the analysis. While a general result is quite messy, we can show that the upper bound on $T\alpha$ [c.f. Eq. (14)] is increasing in $T$ and approaches $-\ln{\frac{L^{2}}{L^{2}+\eta\rho(B)}}$ for large values of $T$. This suggests that the benefits in improving convergence speed from increasing $T$ diminishes for large $T$.

Remark 4

In our algorithm, the stepsize parameters are common among all agents and computing them requires global variables across the network. This global variables can be obtained by applying a consensus algorithm prior to the main algorithm [31, 37, 16].

IV Simulation Results

In this section we present some numerical experiments, where we compare the performance of our proposed algorithms with other first order exact methods. We also study the performance of our algorithms on networks with different sizes and topologies. In all experiments we set $B=\beta A^{\prime}A$ for our algorithm.

We first consider solving a classification problem using regularized logistic regression. We consider a problem with $K$ training samples that are uniformly distributed over $n=10$ agents in a network with $4-$regular graph, in which agents first form a ring and then each agent gets connected to two other neighbors (one from each side). Each agent $i$ has access to a batch of data with size $k_{i}=\lfloor\frac{K}{n}\rfloor$. This problem can be formulated as follows

[TABLE]

where $u_{ij}$ and $v_{ij}$, $j\in\{1,2,...,k_{i}\}$ are the feature vector and the label for the data point $j$ associated with agent $i$ and the regularizer \frac{\nu}{2}\big{\|}x\big{\|}^{2} is added to avoid overfitting. We can write this objective function in the form of $f(x)=\sum_{i=1}^{n}f_{i}(x)$, where $f_{i}(x)$ is defined as

[TABLE]

In order to study the performance of our proposed framework on different networks, we consider $5-30$ agents which are connected with random $4-$regular graphs (agents first form a ring and then each agent gets connected to two other random agents). The objective function at each agent $i$ is of the form $f_{i}(x)=c_{i}(x_{i}-b_{i})^{2}$ with $c_{i}$ and $b_{i}$ being integers that are randomly chosen from $[1,10^{4}]$ and $[1,100]$. We run the simulation for $1000$ random seeds and we plot the average number of steps and communications until the relative error is less than $\epsilon=0.01$, i.e., $\frac{\left|\left|x^{k}-x^{*}\right|\right|}{\left|\left|x^{0}-x^{*}\right|\right|}<0.01$ in Figure 2. The method of multipliers (centralized implementation) is also included as a benchmark. The primal stepsize parameter $\alpha$ at each seed is chosen based on the theoretical bound given in Eq. (14) and the dual stepsize is $\beta=T$. We observe that regardless the size of network, increasing the number of primal updates per iteration improves the performance of the algorithm. We also observe that as the network size grows, the number of steps to optimality of our proposed method grows sublinearly and the communication needed increases almost linearly. In our simulations, we use the mushrooms dataset [6], with $8124$ data points, distributed uniformly over $10$ agents. Each data point has a feature vector of size $112$ and a label which is either $1$ or $-1$. In Figure 1 we compare the performance of our primal-dual algorithms [c.f. Algorithm 1] with $T=1,...,4$ (represented by PD - 1 to PD - 4), with EXTRA [31], DIGing [23], and NEAR-DGD+ [1], in terms of relative error, $\frac{\left|\left|x^{k}-x^{*}\right|\right|}{\left|\left|x^{0}-x^{*}\right|\right|}$, with respect to number of iterations and number of communications. The benchmark $x^{*}$ is computed using minFunc software [30] and the stepsize parameters are manually tuned for each algorithm. We can clearly see that increasing the number of primal updates improves the performance of the algorithm, while incurring a higher communication cost. In our experiments, we observed that by increasing $T$ to larger numbers, the performance of the algorithms does not improve much, which can be explained through Remark 3 and also the effect of the outdated gradients. EXTRA and DIGing algorithms are special cases of our algorithm for specific choices of matrices $A$ and $B$ and one primal update per iteration, and thus have close performance to PD - 1. In the NEAR-DGD*+* the number of communications increases linearly with the iteration number, which results in slow convergence with respect to the number of communications. We obtained similar results for other standard machine learning datasets, including diabetes-scale, heart-scale, a1a, australian-scale, and german [6].

V Final Remarks

In this paper we propose a general framework of first order decentralized primal-dual optimization algorithms. Our general class of algorithms allows for multiple primal updates per iteration, which results in design flexibility to control the trade off between execution complexity and performance of the algorithm. We show that the proposed class of algorithms converges to the exact solution with a global linear rate. The numerical experiments show the convergence speed improvement of primal-dual algorithms with multiple primal updates per iteration compared to other known first order methods like EXTRA, DIGing, and NEAR-DGD+. Possible future work includes analysis of the convergence properties for non-convex objective functions and extending the framework to second order primal-dual algorithms.