Accelerated Primal-Dual Algorithms for Distributed Smooth Convex Optimization over Networks

Jinming Xu, Ye Tian, Ying Sun, Gesualdo Scutari

Introduction

We study distributed (smooth) convex optimization over multi-agent networks, modeled as a fixed, undirected graph. Agents aim to cooperatively solve

The focus of this paper is on optimal rate decentralized algorithms for Problem (1) that use only gradient information and gossip communications. By optimal we mean that these algorithms provably achieve lower complexity bounds for such a class of problems and oracle decentralized algorithms. Primal (Duchi et al., 2012; Yuan et al., 2016; Jakovetic et al., 2014; Nedic and Olshevsky, 2014; Di Lorenzo and Scutari, 2016; Nedich et al., 2017; Qu and Li, 2017b; Xu et al., 2015; Sun et al., 2019) and primal-dual distributed methods (Shi et al., 2015, 2014; Ling et al., 2015; Wei and Ozdaglar, 2012; Chang et al., 2015) applicable to Problem (1) have been extensively studied in the literature, enjoying different convergence rates. In general, these rates are not optimal for several reasons: i) the schemes do not employ any acceleration on the local optimization step and/or communications; or ii) they do not balance optimally the number of optimization and communications steps. Optimal rates of first-order distributed algorithms have been recently studied in Scaman et al. (2017, 2018); Sun and Hong (2018); Uribe et al. (2018); Lan et al. (2017); Shamir (2014); Arjevani and Shamir (2015) for different classes of optimization problems and network topologies; they however are not optimal or applicable to the formulation considered in this paper.

Related works. Optimal lower complexity bounds and matching distributed algorithms have been recently investigated in Scaman et al. (2017) for smooth strongly convex functions, in Scaman et al. (2018) for nonsmooth convex functions, and in Sun and Hong (2018) for smooth nonconvex functions. Fully connected networks have been considered in Shamir (2014); Arjevani and Shamir (2015). However, to our knowledge, no first-order gossip algorithm is known that achieves both computation and communication lower complexity bound for the minimization of smooth convex functions over graphs. Attempts of designing accelerated distributed algorithms for Problem (1) can be found in Li et al. (2018); Qu and Li (2017a); Uribe et al. (2018) and are briefly discussed next. The scheme in Qu and Li (2017a) combines the technique of gradient tracking (Di Lorenzo and Scutari, 2016; Xu et al., 2015; Nedich et al., 2017) with Nesterov acceleration of local computations and achieves an ϵ>0\epsilon>0 solution in O(1/ϵ5/7)O\left({1}/{\epsilon^{5/7}}\right) gradient and communication steps, under the assumption that the solution set of the optimization problem (1) is compact. Algorithm 7 in Uribe et al. (2018) is designed for general smooth convex objectives; it reaches an ϵ\epsilon solution in O(Lf/(η ϵ) log⁡1/ϵ)O\left(\sqrt{L_{f}/(\eta\,\epsilon)}\,\log 1/\epsilon\right) outer loops of communications and O(Lf/ϵlog⁡1/ϵ)O\left(\sqrt{L_{f}/\epsilon}\log 1/\epsilon\right) inner loops of computations (per communication), resulting in an overall gradient evaluations of O(Lf/(ϵη) log⁡21/ϵ)O\left(L_{f}/(\epsilon\sqrt{\eta})\,\log^{2}1/\epsilon\right), which do not match existing lower bounds. The subsequent work (Li et al., 2018) proposes an accelerated penalty-based method with increasing penalty values; the algorithm achieves the lower bound of O(Lf/ϵ)O\left(\sqrt{{L_{f}}/{\epsilon}}\right) gradient evaluations but at the cost of an increasing number of communications per gradient evaluation (iteration)–namely: O(Lf/(ηϵ)log⁡1/ϵ)O\left(\sqrt{{L_{f}}/\left(\eta\epsilon\right)}\log{1}/{\epsilon}\right), making it not optimal in terms of communication steps.

Summary of the contributions. We propose a novel family of primal-dual-based distributed algorithms for Problem (1) that use only gradient information and gossip communications. The algorithms can also employ acceleration on the computation and communications. We provide a unified analysis of their convergence rate, measured in terms of the Bregman distance associated to the saddle point reformation of (1). When acceleration on both computation and communications is properly designed, the proposed algorithms are shown to be optimal, in the sense that they match existing complexity lower bounds (Li et al., 2018), rewritten in terms of the Bregman distance metric. Furthermore, differently from Scaman et al. (2017); Uribe et al. (2018), our algorithms do not require any information on the Fenchel conjugate of the agents’ functions, which significantly enlarge the class of functions to which provably optimal rate algorithms can be applied to. Hence, we termed our algorithms OPTRA (optimal conjugate-free distributed primal-dual methods). Our preliminary numerical results show that OPTRA compares favorably with existing distributed accelerated methods (Li et al., 2018; Qu and Li, 2017a; Uribe et al., 2018) proposed for Problem (1), which supports our theoretical findings.

Technical novelties. While the genesis of OPTRA finds routs in the primal-dual algorithm (Chambolle and Pock, 2011) and employs Nesterov acceleration similarly to Chen et al. (2014) (which also builds on Chambolle and Pock (2011)), there are some substantial differences between the proposed distributed algorithms and the aforementioned schemes (Chambolle and Pock, 2011; Chen et al., 2014), which are briefly discussed next. The scheme in Chambolle and Pock (2011) is meant for abstract saddle-point problems and so Chen et al. (2014) does; the focus therein is not on distributed optimization. Hence communications over networks are not explicitly accounted. Furthermore, Chambolle and Pock (2011) does not employ any acceleration while Chen et al. (2014) accelerates the computation but lacks of the communication (networking) component (no gossip-based updates are present in (Chen et al., 2014, Alg. 2)). On the other hand, OPTRA adopts Nesterov and Chebyshev acceleration to balance computation and communication, so that lower complexity bounds on both are achieved (in terms of Bregman distance). This is a major novelty with respect to Chambolle and Pock (2011); Chen et al. (2014). Because of these differences, the convergence analysis of OPTRA can not be deducted by that of Chambolle and Pock (2011); Chen et al. (2014); a novel convergence proof is provided, which shows an explicit dependence of the rate on key network parameters.

Notations: We use null(⋅){\bf null}(\cdot) (resp. span(⋅){\bf span}(\cdot)) to denote the null space (resp. range space) of the matrix argument. The vector or matrix (with proper dimension) of all ones (resp. all zeros) is denoted by 1{\bf 1} (resp. 0{\bf 0}); eie_{i} denotes the ii-th canonical vector; and the identity matrix is denoted by I\mathbf{I}; the dimensions of these vector and matrices will be clear from the context. The inner product between two matrices x,y\mathbf{x},\mathbf{y} is defined as ⟨x,y⟩:=trace(x,y)\left\langle\mathbf{x},\mathbf{y}\right\rangle:=\text{trace}(\mathbf{x},\mathbf{y}) while the induced norm is ∥x∥:=∥x∥F\left\|\mathbf{x}\right\|:=\left\|\mathbf{x}\right\|_{F}; we will use the same notation for vectors, treated as special cases. Given a positive semidefinite matrix G\mathbf{G}, we define ⟨x,x′⟩G=⟨Gx,x′⟩ and ∥x∥G=⟨Gx,x⟩.\left\langle\mathbf{x},\mathbf{x}^{\prime}\right\rangle_{\mathbf{G}}=\left\langle\mathbf{G}\mathbf{x},\mathbf{x}^{\prime}\right\rangle~{}\text{and}~{}\left\|\mathbf{x}\right\|_{\mathbf{G}}=\sqrt{\left\langle\mathbf{G}\mathbf{x},\mathbf{x}\right\rangle}.

Problem formulation

We study Problem (1) under the following assumptions.

Network model Agents are embedded in a communication network, modeled as an undirected graph G=(E,V)\mathcal{G}=(\mathcal{E},\mathcal{V}), where V\mathcal{V} is the set of vertices–the agents–and E\mathcal{E} is the set of edges; {i,j}∈E\{i,j\}\in\mathcal{E} if there is a communication link between agent ii and agent jj. We assume that the graph has no self-loops, i.e., {i,i}∉E\{i,i\}\notin\mathcal{E}. We use Ni:={j∣{i,j}∈E}\mathcal{N}_{i}:=\{j|\{i,j\}\in\mathcal{E}\} to denote the set of neighbors of agent ii.

Since we are interested in optimization over networks with no centralized nodes, we will focus on distributed algorithms whereby agents communicate with their neighbors using a suitably designed gossip matrix. Standard assumptions on such matrices are the following.

L∈WG\mathbf{L}\in\mathcal{W}_{\mathcal{G}};

Positive semi-definiteness: L⪰0\mathbf{L}\succeq 0, with 0=λ1≤λ2≤λ3≤...≤λm0=\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}\leq...\leq\lambda_{m};

Connectivity: null(L)=span(1){\bf null}(\mathbf{L})={\bf span}({\bf 1});

where {λi}i=1m\{\lambda_{i}\}_{i=1}^{m} are the eigenvalues of L\mathbf{L}.

It is not difficult to check a gossip matrix satisfying Assumption 2 always exists if the associated graph is connected; see, e.g., Olfati-Saber et al. (2007). Several gossip matrices have been considered in the literature; we refer the reader to Xiao and Boyd (2004); Nedić et al. (2018) and references therein for specific examples.

2 Saddle-point reformulation

A standard approach for solving (1) consists in rewriting the optimization problem in the so-called consensus optimization form, that is

To solve Problem (2), we consider the following closely related saddle point formulation

and the saddle-point property Φ(x⋆,y)≤Φ(x⋆,y⋆)≤Φ(x,y⋆),\Phi(\mathbf{x}^{\star},\mathbf{y})\leq\Phi(\mathbf{x}^{\star},\mathbf{y}^{\star})\leq\Phi(\mathbf{x},\mathbf{y}^{\star}), for all (x,y)∈D(\mathbf{x},\mathbf{y})\in\mathcal{D}. Note that x⋆\mathbf{x}^{\star} solves Problem (2) and thus it is also a solution of the original formulation (1) (Bertsekas et al., 2003).

where G(x,x⋆)G(\mathbf{x},\mathbf{x}^{\star}) is the Bregman distance. The following properties of GG are instrumental for our develoments (the proof is provided in the supporting material).

Let x⋆\mathbf{x}^{\star} be any optimal solution of (2); the following hold for GG defined in (5):

xˉ\bar{\mathbf{x}} is an optimal solution of (2) if and only if xˉ∈C\bar{\mathbf{x}}\in\mathcal{C} and G(xˉ,x⋆)=0G(\bar{\mathbf{x}},\mathbf{x}^{\star})=0;

G(x,∙)G(\mathbf{x},\bullet) is constant over the solution set of (2).

Due to (b), for notational simplicity, in what follows, we will write G(x)G(\mathbf{x}) for G(x,x⋆)G(\mathbf{x},\mathbf{x}^{\star}).

In this paper we will use GG as metric to assess the (worst-case) convergence rate of the proposed algorithms as well as to state lower complexity bounds. Note that, since ff is not assumed to be strictly convex, G(x)=0G(\mathbf{x})=0 does not imply x=x⋆\mathbf{x}=\mathbf{x}^{\star}, but it is only a necessary condition for x\mathbf{x} to be optimal (cf. Proposition 1(a)). Still, GG is a valid merit function for both purposes above, as explained next. First, G(x)>ϵG(\mathbf{x})>\epsilon implies that x\mathbf{x} is ϵ\epsilon “far” away (in the GG-measure) from any optimal solution of (2); hence, a lower bound in terms of GG is an informative measure. Furthermore, when it comes to the convergence rate analysis of distributed algorithms, Proposition 1-(a) legitimates the use of (the decay rate of) GG along the agents’ iterates {xk}k=0∞\{\mathbf{x}^{k}\}_{k=0}^{\infty}, as the distance of xk\mathbf{x}^{k} from C\mathcal{C} is proved to be vanishing–see Sec. 4.

Lower Complexity Bounds

We recall here existing lower complexity bounds for decentralized first-order schemes belonging to the same oracle class of the distributed algorithms we are going to introduce. The difference from the literature is that we will write such bounds in terms of the Bregman distance GG. We begin introducing the distributed oracle model, followed by the lower complexity bound.

Distributed first-order oracle A\mathcal{A}: A distributed first order iterative method generates a sequence \big{\{}\mathbf{x}^{(t)}\big{\}}_{t\geq 0}, with x(t)≜[x1(t),…,xm(t)]\mathbf{x}^{(t)}\triangleq[x_{1}^{(t)},\ldots,x_{m}^{(t)}], such that

for all i∈Vi\in\mathcal{V}. We made the blanket assumption that each xi0=0x_{i}^{0}=0, without loss of generality.

The oracle (6) allows each agent to use all the historical values of its local gradients (local computations) as well as that of the decision variables received from its neighbors (local communications). Furthermore, (6) also captures algorithms employing multiple rounds of communications (resp. gradient computations) per gradient evaluation (resp. communication). In the supporting material (Appendix A), we show that the above oracle indeed accounts for most existing distributed algorithms, such as primal-dual methods (Shi et al., 2015) and gradient tracking methods (Di Lorenzo and Scutari, 2016; Nedich et al., 2017; Qu and Li, 2017b; Xu et al., 2015).

A similar black-box procedure has been introduced in Scaman et al. (2017) for strongly convex instances of (1). The difference here is that the oracle in (6) cannot return the gradient of the conjugate of the fif_{i}’s. The reason of considering such “less powerful” methods is that, in practice, it is hard to compute the gradient of conjugate functions. This means that the gossip (dual-based) methods in Scaman et al. (2017) do not belong to the oracle considered in this paper.

2 Lower complexity bounds

We state now lower complexity bounds in the GG-metric for the class of algorithms A\mathcal{A} applied to Problem (2) [and thus (1)] over a connected graph G\mathcal{G}. In Section 4 we will introduce a primal-dual distributed algorithm that indeed converges to an optimal solution of (2) driving GG to zero at a rate that matches the lower complexity bound (see the proofs in the supporting material).

for all t∈[0,d−12(1+⌈15η⌉τc)]t\in\left[0,\frac{d-1}{2}\left({1+\left\lceil\frac{1}{5\sqrt{\eta}}\right\rceil\tau_{c}}\right)\right], where R≜∥x0−x⋆∥R\triangleq\|\mathbf{x}^{0}-\mathbf{x}^{\star}\|. Furthermore,

In the setting of Theorem 2, the overall time needed by any first-order algorithm in A\mathcal{A} using the gossip matrix L\mathbf{L} to drive GG below ϵ>0\epsilon>0, with ff given in Theorem 2, is

Notice that, because of (8), the lower bound (9) can be equivalently stated as

It is not difficult to check that the lower bound in terms of the traditional function-error-based metric (FEM):

has the same expression as (7) [and thus (9) and (10)] up to some constants. This observation is also reported in Li et al. (2018) without proof, and stated formally below for completeness (see the proofs the supporting material).

The above lower bounds tell us that one cannot reach an ϵ\epsilon-solution of (2) (measured either in terms of the GG or FEM-metrics) in less than O(LfR2/ϵ+R∥∇f(x⋆)∥/ϵ)O\left(\sqrt{L_{f}R^{2}/{\epsilon}}+R\|\nabla f(\mathbf{x}^{\star})\|/{\epsilon}\right) computing time and O(τc/η⋅(LfR2/ϵ+R∥∇f(x⋆)∥/ϵ))O\left(\tau_{c}/\sqrt{\eta}\cdot\left(\sqrt{L_{f}R^{2}/\epsilon}+R\|\nabla f(\mathbf{x}^{\star})\|/\epsilon\right)\right) communication time for the worst-case problem as stated in Theorem 2 (see Eq. (28) in the supporting material for a concrete example). Since the time for a single gradient evaluation has been normalized to one, the former lower bound corresponds also to the overall number of gradient evaluations while the overall communication steps read Ω(1/η⋅(LfR2/ϵ+R∥∇f(x⋆)∥/ϵ))\Omega\left(1/\sqrt{\eta}\cdot\left(\sqrt{L_{f}R^{2}/{\epsilon}}+R\|\nabla f(\mathbf{x}^{\star})\|/{\epsilon}\right)\right). This sheds light also on the optimal balance between computation and communication: the optimal number of communication steps per gradient evaluations is ⌈1/η⌉\lceil 1/\sqrt{\eta}\rceil (in the worst case). In the next section, we introduce a distributed, gossip-based algorithm that achieves lower complexity bounds in the GG-metric.

Distributed primal-dual algorithms

A gamut of primal-dual algorithms has been proposed in the literature to solve Problem (2) in a centralized setting; see, e.g., Condat (2013); Chambolle and Pock (2011) and references therein for details. Building on Condat (2013); Chambolle and Pock (2011), here, we propose a general primal-dual algorithm to solve the saddle point problem (3) in a distributed manner. The algorithm reads: given xk\mathbf{x}^{k} and yk\mathbf{y}^{k} at iteration kk,

A=A⊤\mathbf{A}=\mathbf{A}^{\top}, 0⪯A⪯I{\bf 0}\preceq\mathbf{A}\preceq\mathbf{I}, and null(I−A)⊇span(1){\bf null}(\mathbf{I}-\mathbf{A})\supseteq{\bf span}(\mathbf{1});

B=B⊤\mathbf{B}=\mathbf{B}^{\top}, B⪰0\mathbf{B}\succeq{\bf 0}, and null(B)=span(1){\bf null}(\mathbf{B})={\bf span}(\mathbf{1}).

Several choices for A\mathbf{A} and B\mathbf{B} satisfying Assumption 3 are possible, resulting in a gamut of specific algorithms, obtained as instances of (12). Note that, when A\mathbf{A} and B\mathbf{B} satisfy also Assumption 2, all these algorithms are implementable over the graph G\mathcal{G}. Several examples of such distributed algorithms are discussed in details in Appendix A. Here, we only mention that the gradient tracking methods (Di Lorenzo and Scutari, 2016; Nedich et al., 2017; Qu and Li, 2017b; Xu et al., 2015) and primal-dual methods, such as EXTRA (Shi et al., 2015), are all special cases of (12); the former schemes are obtained setting A=W2\mathbf{A}=\mathbf{W}^{2} and B=(I−W)2\mathbf{B}=(\mathbf{\mathbf{I}-\mathbf{W}})^{2}, where W∈WG\mathbf{W}\in\mathcal{W}_{\mathcal{G}} is the weight matrix used by the agents to employ the consensus step; and EXTRA is obtained setting A=W\mathbf{A}=\mathbf{W} and B=I−W\mathbf{B}=\mathbf{I}-\mathbf{W}.

We begin studying convergence of the general primal-dual algorithm (12), under the following tuning of the free parameters:

where λm(B)\lambda_{m}(\mathbf{B}) is the largest eigenvalue of B\mathbf{B}.

Consider Problem (1) under Assumption 1. Given (x1,y1)(\mathbf{x}^{1},\mathbf{y}^{1}), let {(xk,yk)}k=1∞\{(\mathbf{x}^{k},\mathbf{y}^{k})\}_{k=1}^{\infty} be the sequence generated by Algorithm (12), under Assumption 3 and the setting in (13). Define xˉk:=1k−1∑t=2kxt\bar{\mathbf{x}}^{k}:=\frac{1}{k-1}\sum_{t=2}^{k}\mathbf{x}_{t} and R≜∥x1−x⋆∥R\triangleq\|\mathbf{x}^{1}-\mathbf{x}^{\star}\|. Then, the following hold: (i) {xk}k=0∞\{\mathbf{x}^{k}\}_{k=0}^{\infty} converges to an optimal solution x⋆\mathbf{x}^{\star} of (2) [thus x⋆=1x⋆\mathbf{x}^{\star}=\mathbf{1}x^{\star}, for some solution x⋆x^{\star} of (1)]; therefore lim⁡k→∞G(xk)=0\lim_{k\to\infty}G(\mathbf{x}^{k})=0; and (ii) the number of iterations needed for G(xˉk)G(\bar{\mathbf{x}}^{k}) to go below ϵ>0\epsilon>0 isWe use η(B)\eta(\mathbf{B}) to denote the eigengap of B\mathbf{B}.

To match the lower lower bound given in Theorem 2, our next step is accelerating the algorithm, both the computational part and the communication step; we leverage Nesterov acceleration (Nesterov, 2013) for the optimization step while employ Chebyshev polynomials (Wien, 2011) to accelerate communications. To provide some insight of our construction, we begin with the former acceleration; the latter is added in Section 4.3.

2 Accelerated primal-dual algorithms

We accelerate the primal-dual algorithm (12) as follows:

where uk,\mathbf{u}^{k}, x^k,\hat{\mathbf{x}}^{k}, and y^k\hat{\mathbf{y}}^{k} are auxiliary variables and αk,σk,τk,βk\alpha_{k},\sigma_{k},\tau_{k},\beta_{k} are parameters to be properly chosen. Roughly speaking, (15a), (15d) and (15e) are the standard primal-dual steps while (15b) and (15c) are the extra steps meant for the acceleration, with (15b) being the standard Nesterov momentum step and (15c) being a correction step. Note that setting αk≡0,σk≡1,τk≡τ,βk≡1\alpha_{k}\equiv 0,\sigma_{k}\equiv 1,\tau_{k}\equiv\tau,\beta_{k}\equiv 1, the algorithm reduces to the primal-dual method (12). We provide next an instance of (15) that is suitable for a distributed implementation.

Let TT be the overall number of iterations being carried out. The free parameters in (15) is chosen as follows:

The resulting scheme is summarized in Algorithm 1, and its convergence properties are stated in Theorem 6. We point out that Theorem 6, although stated for Algorithm 1, can be readily extended to the more general accelerated primal-dual scheme (15), with other choices of A\mathbf{A} and B\mathbf{B} just satisfying Assumption 3.

If one can set ν=O(ηR/∥∇f(x⋆)∥),\nu=O\left(\sqrt{\eta}R/\left\|\nabla f(\mathbf{x}^{\star})\right\|\right), the above bound can be improved to

Furthermore, the consensus error decays as

While the convergence time of Algorithm 1 benefits from the Nesterov acceleration of the computation step, it is not optimal in terms of communications (optimal dependence on η\eta). In fact, when the network is poorly connected, the second term on the RHS of (17) becomes dominant with respect to the first one, and (17) overall will be larger than (7). This is due to the fact that Algorithm 1 performs a one-consensus-one-gradient update while the lower bound shows an optimal ratio of ⌈1/η⌉\lceil 1/\sqrt{\eta}\rceil in the worst case (cf. Remark 2). This optimal ratio can be achieved accelerating also the communication step, as described in the next section.

3 Optimal primal-dual algorithms with Chebyshev acceleration

We employ the acceleration of the communication step in Algorithm 1 by replacing the gossip matrix L\mathbf{L} with PK(L)P_{K}(\mathbf{L}), where PK(⋅)P_{K}(\cdot) is a polynomial of at most KK degree that maximizes the eigengap of PK(L)P_{K}(\mathbf{L}). This leads to a widely used acceleration scheme known as Chebyshev acceleration and the choice PK(x)=1−TK(c1(1−x))/TK(c1)P_{K}(x)=1-T_{K}(c_{1}(1-x))/T_{K}(c_{1}), with c1=(1+η(L))/(1−η(L))c_{1}=(1+\eta(\mathbf{L}))/(1-\eta(\mathbf{L})) and TK(⋅)T_{K}(\cdot), are the Chebyshev polynomials (Wien, 2011). It is not difficult to check that such a PK(L)P_{K}(\mathbf{L}) is still a gossip matrix. Using in (15) the following setting:

leads to the distributed scheme described in Algorithm 2, whose convergence rate achieves the lower bound (9), as proved in Theorem 7 below. Although the idea of using Chebyshev polynomial has been already used in some (centralized and distributed) algorithms in the literature (Wien, 2011; Scaman et al., 2017), Algorithm 2 substantially differs from that of Scaman et al. (2017), which assumes strongly-convex cost functions and is not rate-optimal in the setting considered in this paper (cf. Sec. G in the supporting material for more details).

If one can set ν=O(R/∥∇f(x⋆)∥),\nu=O\left({R}/{\left\|\nabla f(\mathbf{x}^{\star})\right\|}\right), the above bound can be improved to

Furthermore, the consensus error ∥(I−11Tm)u(t)∥\left\|(\mathbf{I}-\frac{{\bf 1}{\bf 1}^{T}}{m})\mathbf{u}^{(t)}\right\| decays as

According to Theorem 7, given ϵ>0,\epsilon>0, the time needed by the algorithm to drive GG below ϵ>0\epsilon>0 is

matching the lower complexity bound given in (9).

Note that the optimality is stated in terms of the G-metric and does not imply that the algorithm is rate optimal also in the FEM-metric (11), which to date remains an open question. In our experiments (cf. Sec. 5) we observed i) the same behavior of the two errors as a function of the total number of computations and communications; and ii) that Algorithm 2 in fact outperforms existing distributed schemes.

Numerical Results

We report here some preliminary numerical resultsCode: https://github.com/YeTian-93/OPTRA. validating our theoretical findings. We compare the proposed rate-optimal algorithm–OPTRA–with existing accelerated ones designed for convex smooth problems, i.e., Acc-DNGD-NSC (Qu and Li, 2017a) and APM-C (Li et al., 2018). We also included non-accelerated schemes that perform quite well in practice, i.e., i) the gradient tracking method, NEXT/DIGing (Di Lorenzo and Scutari, 2016; Nedich et al., 2017); ii) the primal-dual method, EXTRA (Shi et al., 2015); and iii) the decentralized stochastic gradient method, DPSGD (Lian et al., 2017).

Our experiments are reported in Figure 1, where we plot the Bregman distance versus the overall number of communications and computations performed by each agent (left panel), the number of communications (middle panel), and the number of computations (right panel). The time for local communications and gradient computations using all the local data samples is normalized to one; for DPSGD, the computation time unit is scaled proportionally to the size of the local mini-batch. The plots in terms of the more traditional FEM-metric are reported in the supporting material, the behavior is consistent with the results in Figure 1.

The following comments are in order. The accelerated schemes and the stochastic algorithm–DPSGD–converge faster than the non-accelerated schemes–NEXT/DIGing, EXTRA (the curves of EXTRA and NEXT/DIGing coincide in all the panels). In our experiments, we observed that this gap is quite evident when problems are ill-conditioned. From the right panel, one can see that APM-C performs better than OPTRA and Acc-DNGD-NSC in terms of the number of gradient evaluations, which is expected since APM-C employs an increasing number of communication steps per gradient evaluation. On the other hand, APM-C suffers from high communication cost (which is evident from the middle panel), making it not competitive with respect to OPTRA in terms of communications. When both communication and computation costs are considered (left panel), OPTRA outperforms all the other simulated schemes, which support our theoretical findings.

Conclusion

We studied distributed gossip first-order methods for smooth convex optimization over networks. We provided a novel primal-dual distributed algorithm that employs Nesterov acceleration on the optimization step and acceleration of the communication step via Chebyshev polynomials, balancing thus computation and communication. We also proved that the algorithm achieves the lower complexity bound in the Bregman distance-metric. Preliminary numerical results showed that the proposed scheme outperforms existing distributed algorithms proposed for the same class of problems. An open question, currently under investigation, is whether the proposed distributed algorithms are rate optimal also in terms of the FEM metric. No such an algorithm is known so far in the literature.

Acknowledgments

This work has been supported by the following grants: NSF of USA under Grants CIF 1719205 and CMMI 1832688; in part by the Army Research Office under Grant W911NF1810238; and in part by NSF of China under Grants U1909207 and 61922058.

References

Appendix A Review of existing distributed algorithms and their connections

This section shows the generality of the first-order oracle A\mathcal{A} in (6) and the proposed distributed primal-dual algorithmic framework (12) by casting several existing distributed algorithms in the oracle form (6) and algorithmic form (12).

One of the first distributed algorithms for Problem (1) was proposed in the seminal work Nedic and Ozdaglar (2009) and called Distributed Gradient Algorithm (DGD). DGD employing constant step-size can be written in compact form as:

where W∈WG\mathbf{W}\in\mathcal{W}_{\mathcal{G}}. Defining x(tk)=xk\mathbf{x}^{(t_{k})}=\mathbf{x}^{k}, DGD can be rewritten in a piece-wise continuous form as

which is an instance of the oracle A\mathcal{A}.

The distributed gradient tracking algorithm, first proposed in Di Lorenzo and Scutari (2016); Xu et al. (2015) and further analyzed in Nedich et al. (2017); Qu and Li (2017b), reads

where yk\mathbf{y}_{k} is an auxiliary variable aiming at tracking the gradient of the sum-cost function. The above algorithm is proved to converge at linear rate to a solution of Problem (2), under proper conditions on the stepsize γ\gamma. To show its relationship to the oracle, we first rewrite (22) absorbing the tracking variable y\mathbf{y}, which yields

with x1=Wx0−γ∇f(x0)\mathbf{x}^{1}=\mathbf{W}\mathbf{x}^{0}-\gamma\nabla f(\mathbf{x}^{0}). It is clear that the gradient tracking algorithm belongs to the oracle S\mathcal{S}, as each iteration kk only involves the historical neighboring information and local gradients at k−1k-1 and k−2k-2.

Distributed primal-dual algorithms can be generally written in the following form Shi et al. (2015)

where yk\mathbf{y}_{k} is the dual variable. When y0=0\mathbf{y}^{0}={\bf 0}, the algorithm 23 can solve problem (2). Evaluating (23a) at k+1k+1 and substituting it into (23b) yields

with x1=Wx0−γW∇f(x0)\mathbf{x}^{1}=\mathbf{W}\mathbf{x}^{0}-\gamma\mathbf{W}\nabla f(\mathbf{x}^{0}). It is easy to check that (24) belongs to the oracle A\mathcal{A}.

There are some other distributed algorithms that do not belong to the categories above such as Chen and Sayed (2012). However, using similar arguments as above, one can show that they are instances of the oracle A\mathcal{A}.

A.2 Connections between gradient tracking and primal-dual methods

We reveal here an interesting connection between primal-dual methods and gradient tracking based methods. More specifically, setting in (12a) A=W2\mathbf{A}=\mathbf{W}^{2} and B=(I−W)2\mathbf{B}=(\mathbf{\mathbf{I}-\mathbf{W}})^{2}, one can easily recover gradient tracking methods from the primal-dual ones. To simplify the presentation, we consider a slightly different form of (12a), i.e.,

Then, from (25a), we have at iteration k+1k+1

Subtracting (25a) from the above equation we have

Let −γWyk=xk+1−Wxk-\gamma\mathbf{W}\mathbf{y}^{k}=\mathbf{x}^{k+1}-\mathbf{W}\mathbf{x}^{k} and suppose W\mathbf{W} is invertible. Then, we have

which is exactly the standard gradient tracking method in the ATC form (Di Lorenzo and Scutari, 2016; Xu et al., 2015).

Appendix B Proof of Proposition 1

Statement (a) is a direct result of (Bertsekas et al., 2003, Prop. 6.1.1). We prove next statement (b). Suppose that there are two optimal solutions x⋆\mathbf{x}^{\star} and x~⋆\widetilde{\mathbf{x}}^{\star} such that

where we have used the fact that ⟨∇f(z),z⟩=0\left\langle\nabla f(\mathbf{z}),\mathbf{z}\right\rangle=0 for any optimal solution z\mathbf{z}.

Appendix C Proof of Theorem 2

The proof is based on building a worst-case objective function in (27) and network graph for which the lower bound is achieved by the best available gossip, distributed algorithm in the oracle A\mathcal{A}. To do so we build on the cost function first introduced in Arjevani and Shamir (2015) for a fully connected network and later used for a peer-to-peer network in Scaman et al. (2017), both for smooth strongly convex problems. Since we use a different metric (the Bregman distance) to define the lower bound and consider smooth convex problems (not necessarily strongly-convex), the analysis in Scaman et al. (2017) cannot be readily applied to our setting and an ad-hoc proof of the theorem is needed.

The path of our proof is the following: i) We start with a simple network consisting of two agents such that the diameter of the network will not come into play–see Sec. C.1; and ii) then we extend our results to a general network composed by an arbitrary number of agents–see Sec. C.2.

We state the result on the simple two-agent network as the following.

Consider a two-agent network with cost functions given in (28). Let {xk}k=0∞\{\mathbf{x}^{k}\}_{k=0}^{\infty} be the sequence generated by any first-order algorithm A\mathcal{A}. Suppose 0≤k≤d−120\leq k\leq\frac{d-1}{2}. Then, we have

We prove the above result in three steps: i) we construct the hard function in Sec. C.1.1, which is the worst-case function for all methods belonging to the oracle A\mathcal{A}; ii) we introduce some intermediate result in Sec. C.1.2, which is related to our specific metric–the Bregman distance GG, and iii) building on step i-ii, we derive the lower bound in Sec. C.1.3.

Consider a network composed of two agents. The idea of the proof of the lower complexity bound relies on splitting the “hard” function used by Nesterov to prove the iteration complexity of first-order gradient methods for (centralized) smooth convex problems across the agents(Nesterov, 2013, Chapter 2). More specifically, consider the following cost functions for the two agents:

are two d×dd\times d matrices with their leading principal minors of order k∈[1,d]k\in[1,d] having non-zero block diagonals while the rest being zero.

The key idea of Nesterov proof for the lower complexity bound of centralized first-order gradient methods consists in designing the “hardest” function to be minimized by any method belonging to the oracle. This function was shown to be such that, at iteration kk, all these methods produce a new iterate whereby only the kkth component is updated. The choice of the two agents’ cost functions in (28) follows the same rationale: the structure of A1,[k]\mathbf{A}_{1,[k]} and A2,[k]\mathbf{A}_{2,[k]} is such that none of the two agents is able to make progresses towards optimality, i.e., updating the next component in their local optimization vector (with odd index for agent 11 and even index for agent 22) just performing local gradient updates and without communication with each other. This means that at certain stages a communication between the two agents is necessary for the algorithm to make progresses towards optimality. Building on the above idea, we begin establishing the lower complexity bound for the two-agent network problem in terms of gradient evaluations.

C.1.2 Intermediate results

Now substituting f(x)=f1,[k](x1)+f2,[k](x2)f(\mathbf{x})=f_{1,[k]}(x_{1})+f_{2,[k]}(x_{2}) in (27) and ignoring constants, we obtain

We denote the optimal function value of the above problem as f[k]⋆.f_{[k]}^{\star}. It is obvious that, when agents reach consensus, i.e., x1=x2x_{1}=x_{2}, the function f[k](x)f_{[k]}(\mathbf{x}) will reduce to the Nesterov’s “hard” function (Nesterov, 2013, Section 2.1.2), for which we have the optimal solution

and f[k]⋆=Lf8(−1+1k+1)f_{[k]}^{\star}=\frac{L_{f}}{8}(-1+\frac{1}{k+1}). Also, we have

Note that quantities (31) and (33) will be useful later to relate the complexities with ∥x0−x⋆∥\left\|\mathbf{x}^{0}-\mathbf{x}^{\star}\right\| and ∥∇f(x⋆)∥.\left\|\nabla f(\mathbf{x}^{\star})\right\|. According to (32), Problem (30) further becomes

Let {xk}k=0∞\{\mathbf{x}^{k}\}_{k=0}^{\infty} be the sequence generated by any distributed first-order algorithm A\mathcal{A} with x0=0\mathbf{x}^{0}={\bf 0}. Then, xik∈Lkx^{k}_{i}\in\mathcal{L}^{k} for all k≥0k\geq 0 and all i∈Vi\in\mathcal{V}.

The proof of the above lemma is straightforward, since local communication steps do not change the space spanned by the historical gradient vectors generated over the network.

which attains the optimum f[k,1]⋆=Lf8(−k2(k+1)2−1(k+1)2)f^{\star}_{[k,1]}=\frac{L_{f}}{8}(-\frac{k^{2}}{(k+1)^{2}}-\frac{1}{(k+1)^{2}}).

which gives f[k,3]⋆=Lf8(−k2(k+1)2−3(k+1)2)f^{\star}_{[k,3]}=\frac{L_{f}}{8}(-\frac{k^{2}}{(k+1)^{2}}-\frac{3}{(k+1)^{2}}).

The proof is completed by combining the two cases above. ∎

C.1.3 Proof of Theorem 8

We can now prove the theorem. Let us fix kk and apply the first-order gossip algorithm A\mathcal{A} to minimize f[2k+1]f_{[2k+1]}. Since x0=0\mathbf{x}^{0}={\bf 0}, invoking Lemma 11, we have

where the last inequality comes from the previously developed facts ∥x⋆∥2=Θ(k+1)\left\|\mathbf{x}^{\star}\right\|^{2}=\Theta\left(k+1\right), ∥∇f(x⋆)∥=Θ(Lfk+1)\left\|\nabla f(\mathbf{x}^{\star})\right\|=\Theta\left(\frac{L_{f}}{\sqrt{k+1}}\right) and thus \left\|\mathbf{x}^{0}-\mathbf{x}^{\star}\right\|^{2}=\Theta\big{(}\frac{k+1}{L_{f}}\left\|\mathbf{x}^{0}-\mathbf{x}^{\star}\right\|\left\|\nabla f(\mathbf{x}^{\star})\right\|\big{)}. This completes the proof for the two-agent network. □\square

The lower bound we develop in Theorem 8 for distributed scenarios has similar structure of that of the recent paper Ouyang and Xu (2018), where the lower bound is derived for general equality-constrained problems in centralized scenarios (i.e., Ax=b\mathbf{Ax=b}). Notice that the results and techniques therein can not apply to our distributed setting, as we require b=0\mathbf{b=0} and A∈WG\mathbf{A}\in\mathcal{W}_{\mathcal{G}} while the lower bound in Ouyang and Xu (2018) is determined by a choice of b\mathbf{b} and A\mathbf{A} that does not meet our requirement.

C.2 Proof of Theorem 2

Following the same path of Scaman et al. (2017), we now extend the above analysis to the general network setting (arbitrary number of agents) by employing a line graph and constructing certain number of pairwise two-agent networks as in (28) from the left and the right of the line graph, respectively, yielding two subgroups. Between these two subgroups, we place a number (proportional to the diameter of the network) of agents with zero cost functions to ensure the necessity of communications between the agents in the two subgroups. To prove the time complexity lower bound, we then leverage the effect of the network by establishing the connection between the diameter of the network and the eigengap of the gossip matrix.

Let ηn=1−cos⁡(πT)1+cos⁡(πT).\eta_{n}=\frac{1-\cos\left(\frac{\pi}{T}\right)}{1+\cos\left(\frac{\pi}{T}\right)}. For a given η∈(0,1],\eta\in(0,1], there exists n≥2n\geq 2 such that ηn≥η>ηn+1.\eta_{n}\geq\eta>\eta_{n+1}. We treat the cases n=2n=2 and n≥3n\geq 3 separately. Let us first consider the case n≥3n\geq 3. There exists a line graph of m=nm=n agents and associated Laplacian weight matrix with eigengap η.\eta. Now, let us define two subsets of agents as \mathcal{A}_{l}=\left\{i\big{|}1\leq i\leq\lceil\zeta m\rceil\right\} and \mathcal{A}_{r}=\left\{i\big{|}\lfloor(1-\zeta)m\rfloor+1\leq i\leq m\right\}, which lie on the left and the right of the line graph, respectively; the parameter ζ∈(0,12)\zeta\in(0,\frac{1}{2}) is to be determined. The distance between the two subsets is thus dc≜⌊(1−ζ)m⌋+1−⌈ζm⌉.d_{c}\triangleq\lfloor(1-\zeta)m\rfloor+1-\lceil\zeta m\rceil. The class of local functions is defined as follows

where A1,[k],A2,[k]\mathbf{A}_{1,[k]},\mathbf{A}_{2,[k]} are the two matrices defined in (29). Similarly to the two-agent network case (cf. Sec. C.1), we have

Similarly as the two-agent case, one can verify that \left\|\mathbf{x}^{0}-\mathbf{x}^{\star}\right\|^{2}=\Theta\big{(}\frac{k+1}{L_{f}}\left\|\mathbf{x}^{0}-\mathbf{x}^{\star}\right\|\left\|\nabla f(\mathbf{x}^{\star})\right\|\big{)}.

To have at least one non-zero element at the kkth component among the local copies of agents in both of the above two subsets, one must perform at least kk local computation steps and (k−1)dc(k-1)d_{c} communication steps. Thus, we have

where (a) is due to η>ηm+1>2(m+1)2\eta>\eta_{m+1}>\frac{2}{(m+1)^{2}} and (b) is due to η≤η3=13.\eta\leq\eta_{3}=\frac{1}{3}. Further, since dcd_{c} is an integer, we have dc≥⌈15η⌉.d_{c}\geq\left\lceil\frac{1}{5\sqrt{\eta}}\right\rceil. Combining (37) and (38) leads to

We focus now on the case n=2n=2. Consider a complete graph of 33 agents with associated Laplacian matrix having eigengap equal to η\eta. The agents’ cost functions are

Following similar steps as above, one can show that

which leads to the same expression of the lower bound as in (39). This concludes the proofs.

Appendix D Proof of Theorem 4

For the cost functions as mentioned above, one can also verify that (cf. Appendix C.1.3)

and thus Lf∥x⋆−x0∥2k+1=Θ(∥x⋆−x0∥∥∇f(x⋆)∥)\frac{L_{f}\left\|\mathbf{x}^{\star}-\mathbf{x}^{0}\right\|^{2}}{k+1}=\Theta\left(\left\|\mathbf{x}^{\star}-\mathbf{x}^{0}\right\|\left\|\nabla f(\mathbf{x}^{\star})\right\|\right). As a result, the RHS of (40) can be rewritten as:

which translate to the following lower bounds in terms of number of iterations, respectively:

The rest of proof follows by the same argument as in Section C.2 to relate kk to the absolute time tt as well as the eigengap η\eta of the network.

Appendix E Proofs for the Upper Complexity Bounds

This section is devoted to the proofs of the upper complexity bounds of the proposed algorithms. We begin in Sec. E.1 establishing two fundamental inequalities that are valid for all feasible primal-dual solutions of Problem (3); see Lemma 12 and Lemma 13. Then, applying these inequalities to a saddle point solution of Problem (3), we obtain the convergence rate of Algorithms 1 and 2 in terms of the Bregman distance, see Sec. E.2and Sec. E.3 respectively. Finally in Sec. E.4, we apply the analysis of Chebyshev polynomials to show that the eigengap of the communication matrix B\mathbf{B}, as a polynomial of the gossip matrix, can be upper bounded by a constant, leading to the upper complexity bound that matches the established lower bound.

Consider Algorithm (15). We define τ=1νTλm(B).\tau=\frac{1}{\nu T\lambda_{m}(\mathbf{B})}. Then we have

where h(⋅)=12γ∥⋅∥A−A22h(\cdot)=\frac{1}{2\gamma}\left\|\cdot\right\|^{2}_{\mathbf{A}-\mathbf{A}^{2}}, uk+12=xk−γ(∇f(xk)+y^k)\mathbf{u}^{k+\frac{1}{2}}=\mathbf{x}^{k}-\gamma(\nabla f(\mathbf{x}^{k})+\hat{\mathbf{y}}^{k}) and Lf=max⁡i{Lfi}L_{f}=\max_{i}\{L_{f_{i}}\}.

Since ff is LfL_{f}-smooth by Assumption 1, we have

and using f(Ax)≥f(xk)+⟨∇f(xk),Ax−xk⟩f(\mathbf{A}\mathbf{x})\geq f(\mathbf{x}^{k})+\left\langle\nabla f(\mathbf{x}^{k}),\mathbf{A}\mathbf{x}-\mathbf{x}^{k}\right\rangle, further gives

Also, subtracting Auk+1\mathbf{A}\mathbf{u}^{k+1} from both sides of (15a), multiplying (15d) by γA\gamma\mathbf{A}, and adding the obtained two equations while using (15e) lead to

where in (∗)(*) we used βk−1=τkτk−1,τk=τθk\beta_{k-1}=\frac{\tau_{k}}{\tau_{k-1}},\tau_{k}=\frac{\tau}{\theta_{k}}. Notice that, for the above derivation, we implicitly assume that k≥2.k\geq 2. However, with the definition of x^1:=x1\hat{\mathbf{x}}^{1}:=\mathbf{x}^{1} and the fact that y^1=τ1Bx1\hat{\mathbf{y}}^{1}=\tau_{1}\mathbf{B}\mathbf{x}^{1}, we still have (I−A)u2=−A(u2−x1+γ(∇f(x1)+y2))−γτθ1AB(x^1−x^2)(\mathbf{I}-\mathbf{A})\mathbf{u}^{2}=-\mathbf{A}\left(\mathbf{u}^{2}-\mathbf{x}^{1}+\gamma\left(\nabla f(\mathbf{x}^{1})+\mathbf{y}^{2}\right)\right)-\frac{\gamma\tau}{\theta_{1}}\mathbf{A}\mathbf{B}(\hat{\mathbf{x}}^{1}-\hat{\mathbf{x}}^{2}).

Multiplying uk+12−x\mathbf{u}^{k+\frac{1}{2}}-\mathbf{x} from both sides of the above equation and using the convexity of h(⋅)h(\cdot) and the fact that uk+1=Auk+12\mathbf{u}^{k+1}=\mathbf{A}\mathbf{u}^{k+\frac{1}{2}} we obtain

Since σk=1θk+1\sigma_{k}=\frac{1}{\theta_{k+1}} and αk=θk+1θk−θk+1\alpha_{k}=\frac{\theta_{k+1}}{\theta_{k}}-\theta_{k+1}, using (15b) and (15c) leads to

We implicitly assumed k≥2k\geq 2; still we have x^2−x^1=1θk(u2−x1)\hat{\mathbf{x}}^{2}-\hat{\mathbf{x}}^{1}=\frac{1}{\theta_{k}}\left(\mathbf{u}^{2}-\mathbf{x}^{1}\right), recalling that x^1=x1\hat{\mathbf{x}}^{1}=\mathbf{x}^{1}. Thus, (42) becomes

which, recalling that Φ(x,y)=f(x)+⟨y,x⟩\Phi(\mathbf{x},\mathbf{y})=f(\mathbf{x})+\left\langle\mathbf{y},\mathbf{x}\right\rangle, completes the proof. ∎

In the setting of Lemma 12, let 1θk−12−1−θkθk2=0\frac{1}{\theta_{k-1}^{2}}-\frac{1-\theta_{k}}{\theta_{k}^{2}}=0, that is, 1θk=1+1+4(1θk−1)22\frac{1}{\theta_{k}}=\frac{1+\sqrt{1+4(\frac{1}{\theta_{k-1}})^{2}}}{2}, with θ1=1\theta_{1}=1 and (1−γLf)I−γτθk2B⪰0(1-\gamma L_{f})\mathbf{I}-\frac{\gamma\tau}{\theta_{k}^{2}}\mathbf{B}\succeq{\bf 0}, for all 1≤k≤T−11\leq k\leq T-1. Suppose Assumptions 1 and 3 hold. Then, for any x∈C,y∈C⊥\mathbf{x}\in\mathcal{C},\mathbf{y}\in\mathcal{C}^{\perp}, we have

where NN is the overall number of iterations.

Applying Lemma 12 with x∈C\mathbf{x}\in\mathcal{C}, we have (note that Ax=x\mathbf{A}\mathbf{x}=\mathbf{x} by Assumption 3)

Likewise, with x=uk−12\mathbf{x}=\mathbf{u}^{k-\frac{1}{2}} we have

Let Vk=Φ(uk,y)−Φ(x,y)+h(uk+12)−h(x)V_{k}=\Phi(\mathbf{u}^{k},\mathbf{y})-\Phi(\mathbf{x},\mathbf{y})+h(\mathbf{u}^{k+\frac{1}{2}})-h(\mathbf{x}). Then, multiplying (46) by 1−θk1-\theta_{k} and (45) by θk\theta_{k}, and combing the obtained equations yield

where in the last equality we used 1⊤yk=0,∀k≥1{\bf 1}^{\top}\mathbf{y}^{k}={\bf 0},\forall k\geq 1 and the following result (recall BJ=JB=0\mathbf{B}\mathbf{J}=\mathbf{J}\mathbf{B}={\bf 0} and y∈C⊥\mathbf{y}\in\mathcal{C}^{\perp}):

Dividing θk2\theta_{k}^{2} from both sides of (47) leads to

Summing (48) over kk from 11 to T−1T-1 yields

Recalling that 1θk−12−1−θkθk2=0\frac{1}{\theta_{k-1}^{2}}-\frac{1-\theta_{k}}{\theta_{k}^{2}}=0 and θ1=1\theta_{1}=1, by induction it is easy to see that k+1>1θk≥k+12k+1>\frac{1}{\theta_{k}}\geq\frac{k+1}{2} and thus 1θk2−1θk−12=1θk>0\frac{1}{\theta^{2}_{k}}-\frac{1}{\theta^{2}_{k-1}}=\frac{1}{\theta_{k}}>0. Then, with x^1=x1:=u1,y1:=0\hat{\mathbf{x}}^{1}=\mathbf{x}^{1}:=\mathbf{u}^{1},\mathbf{y}^{1}:={\bf 0}, (49) can be simplified as

Since ρ((B+J)−1)=1λmin(B+J)=1λ2(B)\rho\left(\left(\mathbf{B}+\mathbf{J}\right)^{-1}\right)=\frac{1}{\lambda_{\text{min}}\left(\mathbf{B}+\mathbf{J}\right)}=\frac{1}{\lambda_{2}(\mathbf{B})}, B⪰0\mathbf{B}\succeq{\bf 0} and (1−γLf)I−γτθk2B⪰0(1-\gamma L_{f})\mathbf{I}-\frac{\gamma\tau}{\theta_{k}^{2}}\mathbf{B}\succeq{\bf 0}, we further have

which, together with the fact that Vk≥Φ(uk,y)−Φ(x,y)V_{k}\geq\Phi(\mathbf{u}^{k},\mathbf{y})-\Phi(\mathbf{x},\mathbf{y}), completes the proof. ∎

E.2 Proof of Theorem 5

Note that the primal-dual method (12) is a special case of the update (15) with the setting θk≡1,αk≡0,σk≡1,τk≡τ,βk≡1\theta_{k}\equiv 1,\alpha_{k}\equiv 0,\sigma_{k}\equiv 1,\tau_{k}\equiv\tau,\beta_{k}\equiv 1. Furthermore, γ\gamma and τ\tau defined in (13) satisfy (1−γLf)I−γτB≥0(1-\gamma L_{f})\mathbf{I}-\gamma\tau\mathbf{B}\geq 0 and xk≡uk\mathbf{x}^{k}\equiv\mathbf{u}^{k}. Invoking (49) with these parameter settings and x:=x⋆, y:=y⋆=−∇f(x⋆), x^1:=x1,y1:=0,\mathbf{x}:=\mathbf{x}^{\star},\,\mathbf{y}:=\mathbf{y}^{\star}=-\nabla f(\mathbf{x}^{\star}),\,\hat{\mathbf{x}}^{1}:=\mathbf{x}^{1},\mathbf{y}^{1}:={\bf 0}, we have

Let xˉT:=1T−1∑k=2Txk\bar{\mathbf{x}}^{T}:=\frac{1}{T-1}\sum_{k=2}^{T}\mathbf{x}_{k}. Using the convexity of Φ\Phi, we furhter have

E.3 Proof of Theorem 6

Since γ=ννLf+T,τ=1νTλm(B)\gamma=\frac{\nu}{\nu L_{f}+T},\tau=\frac{1}{\nu T\lambda_{m}(\mathbf{B})} and 1k+1<θk<2k+1\frac{1}{k+1}<\theta_{k}<\frac{2}{k+1}, we have

Then, invoking Lemma 13 with x=x⋆,y=y⋆=−∇f(x⋆)\mathbf{x}=\mathbf{x}^{\star},\mathbf{y}=\mathbf{y}^{\star}=-\nabla f(\mathbf{x}^{\star}) and knowing that Φ(uk,y⋆)−Φ(x⋆,y⋆)=G(uk)≥0\Phi(\mathbf{u}^{k},\mathbf{y}^{\star})-\Phi(\mathbf{x}^{\star},\mathbf{y}^{\star})=G(\mathbf{u}^{k})\geq 0 (cf., the relation (5) in the main text), we obtain

where Rx=∥u1−x⋆∥2,Ry=∥∇f(x⋆)∥2R_{x}=\left\|\mathbf{u}^{1}-\mathbf{x}^{\star}\right\|^{2},R_{y}=\left\|\nabla f(\mathbf{x}^{\star})\right\|^{2}. Setting ν=η(B)\nu=\sqrt{\eta(\mathbf{B})} we have

which, together with the time (1+tc1+t_{c}) needed at each iteration, gives the overall time complexity.

If we setNote that this requires accurate estimates on the ratio of Rx/RyR_{x}/R_{y}, which, indeed, plays a key role of trade-off parameter balancing gradient computation steps and communication steps. ν=η(B)RxRy\nu=\sqrt{\frac{\eta(\mathbf{B})R_{x}}{R_{y}}}, we have

which matches the lower bound also with respect to Rx,RyR_{x},R_{y}.

In the following, we show that both consensus error and the absolute value of the objective error will converge at the same rate as the Bregman distance. Invoking Lemma 13 with x=x⋆\mathbf{x}=\mathbf{x}^{\star}, γ=ννLf+T\gamma=\frac{\nu}{\nu L_{f}+T}, τ=1νTλm(B)\tau=\frac{1}{\nu T\lambda_{m}(\mathbf{B})} and ν=η(B)\nu=\sqrt{\eta(\mathbf{B})}, we have

where ϕ(⋅):=2LfRxT2+2η(Rx+(⋅)2)T\phi(\cdot):=\frac{2L_{f}{R}_{x}}{T^{2}}+\frac{\frac{2}{\sqrt{\eta}}({R}_{x}+\left(\cdot\right)^{2})}{T}.

E.4 Proof of Theorem 7

Following the similar lines in (Scaman et al., 2017, Theorem 4), we first consider the normalized Laplacian L\mathbf{L} has a spectrum in [1−c1−1,1+c1−1].[1-c_{1}^{-1},1+c_{1}^{-1}]. According to Scaman et al. (2017); Wien (2011), the Chebyshev polynomail PK(x)=1−TK(c1(1−x))TK(c1)P_{K}(x)=1-\frac{T_{K}(c_{1}(1-x))}{T_{K}(c_{1})} is the solution of the following problem

Define δ=2c0K1+c02K.\delta=2\frac{c_{0}^{K}}{1+c_{0}^{2K}}. Since Algorithm 2 amounts to an instance of Procedure (15) with A=I−c2⋅PK(L)\mathbf{A}=\mathbf{I}-c_{2}\cdot P_{K}(\mathbf{L}) and B=PK(L)\mathbf{B}=P_{K}(\mathbf{L}), its convergence proof follows the same lines as that of Theorem 6 with the following properties of PK(L)P_{K}(\mathbf{L}): i) PK(L)P_{K}(\mathbf{L}) is symmetric; ii) according to (52), 0⪯I−c2⋅PK(L)⪯I{\bf 0}\preceq\mathbf{I}-c_{2}\cdot P_{K}(\mathbf{L})\preceq\mathbf{I} and PK(L)⪰0P_{K}(\mathbf{L})\succeq{\bf 0}, and null(PK(L))=C;{\bf null}(P_{K}(\mathbf{L}))=\mathcal{C}; iii) The values given for γ\gamma and τ\tau in Algorithm 2 ensures that (1−γLf)I−γτθk2PK(L)⪰0(1-\gamma L_{f})\mathbf{I}-\frac{\gamma\tau}{\theta_{k}^{2}}P_{K}(\mathbf{L})\succeq{\bf 0}, analogus to (51). Therefore we have

where (*) requires a specified ν\nu. Finally, we have

Taking K=⌈1η⌉K=\left\lceil\frac{1}{\sqrt{\eta}}\right\rceil, we have

where (*) is due to the fact that ⌈1η⌉≥1\left\lceil\frac{1}{\sqrt{\eta}}\right\rceil\geq 1. Thus, we have 1+δ1−δ≤1+e−11−e−1≤2.5\sqrt{\frac{1+\delta}{1-\delta}}\leq\frac{1+e^{-1}}{1-e^{-1}}\leq 2.5, which, together with the time (1+Kτc)(1+K\tau_{c}) needed at each iteration, gives the time complexity as announced. □\square

Appendix F Additional numerical results

This section provides additional numerical results, complementing those reported in the paper (cf. Section 5). In Figure 2 we plot the FEM-metric (11) versus the overall number of communications and computations performed by each agent (left panel), the number of communications (middle panel), and the number of computations (right panel). The comparison of the different schemes suggests to the same conclusions as in Sec.5; the only exception is that in the FEM-metric, the stochastic algorithm–DPSGD–does not present a significant advantage with respect to the non-accelerated algorithms–NEXT, DIGing and EXTRA.

F.2 Decentralized logistic regression

We estimated LfL_{f} for the problem as Lf=24L_{f}=24 and tuned the free parameters of the simulated algorithms manually to achieve the best practical performance for each algorithm. This leads to the following choices: i) the step size of NEXT/DIGing is set to 0.010.01; ii) the step size of EXTRA is set to 0.005;0.005; iii) for Acc-DNGD-NSC, we used the fixed step-size rule, with η=0.01/Lf\eta=0.01/L_{f}; iv) for APM-C, we set (see notation therein) Tk=⌈c⋅(log⁡k/1−σ2(W))⌉T_{k}=\lceil c\cdot({\log k}/{\sqrt{1-\sigma_{2}(\mathbf{W}))}}\rceil, with c=0.2c=0.2 and β0=104;\beta_{0}=10^{4}; v) for DPSGD, we set its step size as 0.0010.001 and the portion of batch size to the full local data set as 20%20\% and for vi) for our algorithm, we set ν=1500\nu=1500 and K=2.K=2.

The experiment result is reported in Figure 3. The first row of panels shows the Bregman distance versus the total cost (left panel), the communication cost (middle panel), and the gradient computation cost (right panel). The second row plots the FEM-metric versus the same quantities as in the first row. Both the communication time unit and the computation time unit for a full epoch of local data is set as 1. For DPSGD, the computation time unit is scaled in proportion to the local batch size. The only existing algorithm that has a comparable performance with the proposed OPTRA is APM-C. As discussed in the task of decentralized linear regression, APM-C performs better than OPTRA in terms of the number of gradient computations, while suffers from high communication cost. In terms of the overall number of communications and computations, OPTRA outperforms all the other simulated schemes under the above setting.

F.3 Different ratio of communication time versus computation time

In all the previous experiments, we set both the communication time unit and the computation time unit for a full epoch of local data as 1. To incorporate scenarios where a full epoch computation of local gradient is much more expensive than one communication process, we re-conducted the previous experiments in the setting where the communication time unit is 1 while the computation time unit for a full epoch of local data is 5. Note that all the process of data generation and parameters tunings are the same as in the Sec. F.2. The results are reported in Figure 4 and Figure 5 respectively for decentralized linear regression problem and the decentralized logistic regression problem. It can be seen that OPTRA outperforms all the other simulated schemes in terms of the overall number of communications and computations, especially when the communication cost is not negligible.

Appendix G Additional Comments

Scaman et al. (2017) presented optimal algorithms for decentralized optimization of strongly convex smooth functions. The functions considered in our paper are just convex (smooth). However, adding a small regularization (of the order of ϵ∥x∥2/R2\epsilon\left\|\mathbf{x}\right\|^{2}/R^{2}), the problem becomes strongly convex and the results of (Scaman et al., 2017) apply. In so doing, one can show that the method in (Scaman et al., 2017) achieves an ϵ\epsilon-solution in O((1+1η τc)LfR2ϵlog⁡(1ϵ))O\left(\left(1+\frac{1}{\sqrt{\eta}}\,\tau_{c}\right)\sqrt{\frac{L_{f}R^{2}}{\epsilon}}\log\left(\frac{1}{\epsilon}\right)\right). However, this requires an accurate estimate of RR and the resulting rate differs from the lower bound for smooth convex functions by an extra log-factor “log⁡1/ϵ\log{1/\epsilon}”, meaning that this is not optimal as our scheme. Also, more importantly, i) Scaman et al. (2017) requires the computation in closed form of the gradient of the Fenchel conjugate while our scheme does not have this limitation; ii) the Lipschitz constant of the gradient of the Fenchel conjugate in the setting above will scale as 1/ϵ1/\epsilon and thus the condition number of the dual problem is O(1/ϵ)O(1/\epsilon), which becomes arbitrarily large as ϵ\epsilon decreases. As a consequence, (Scaman et al., 2017) significantly slow-downs in practice. This motivates our design of distributed algorithms specifically for convex (but not strongly convex) functions.