Achieving Geometric Convergence for Distributed Optimization over Time-Varying Graphs

Angelia Nedich, Alex Olshevsky, Wei Shi

Introduction

This paper focuses on the following distributed convex optimization problem:

We assume that the functions fif_{i} in problem (1) are convex and continuously differentiable. For such a problem, we propose a class of distributed algorithms that solve the problem over time-varying connectivity graphs for two different cases, namely, the case when the graphs are undirected and the case when they are directed. The algorithms employ consensus ideas for estimating the gradient of the global objective function in (1). When at least one of the objective functions is strongly convex, we show that the algorithms achieve R-linear convergence rates Suppose that a sequence {x(k)}\{x(k)\} converges to xx^{*} in some norm \|\cdot\|. We say that the convergence is: (i) Q-linear if there exists λ(0,1)\lambda\in(0,1) such that x(k+1)xx(k)xλ\frac{\|x(k+1)-x^{*}\|}{\|x(k)-x^{*}\|}\leq\lambda for all kk; (ii) R-linear if there exists λ(0,1)\lambda\in(0,1) and some positive constant CC such that x(k)xCλk\|x(k)-x^{*}\|\leq C\lambda^{k} for all kk. Both of these rates are geometric, and they are often referred to as global rates to be distinguished from the case when the given relations are valid for some sufficiently large indices kk. The difference between these two types of geometric rate is in that Q-linear rate implies monotonic decrease of x(k)x\|x(k)-x^{*}\|, while R-linear rate does not..

The research on distributed optimization dates back to 1980s . Since then, due to the emergence of large-scale networks, the development of distributed algorithms for solving problem in (1) has received significant attention recently. Besides the decentralized and distributed approaches we are going to discuss below, many efforts have been made to solve (1) in a master-slave structured isotropic network. The distributed algorithms designed over such special structure are usually fast in practice and mostly used in machine learning to handle big-data in a cluster of computers . Such scheme is “centralized” due to the use of a “master”. In this paper, we focus on solving (1) in a decentralized fashion motivated by the applications mentioned above.

Some earlier methods include distributed incremental (sub)gradient methods and incremental proximal methods , while a more recent work includes incremental aggregated gradient methods and its proximal gradient variants . All of the incremental methods require a special ring networks due to the nature of these methods. To handle a more general (possibly time-varying) networks, distributed subgradient algorithm was proposed in , while its stochastic variant was studied in and its asynchronous variant in with provable convergence rates. These algorithms are intuitive and simple but usually slow due to the fact that even if the objective functions are differentiable and strongly convex, these methods still need to use diminishing step-size to converge to a consensual solution. Other works on distributed algorithms that also require the use of diminishing step-sizes include . With a fixed step-size, these distributed methods can be fast, but they only converge to a neighborhood of the solution set. This phenomenon creates an exactness-speed dilemma. A different class of distributed approaches that bypasses this dilemma is based on introducing Lagrangian dual variables and working with the Lagrangian function. The resulting algorithms include distributed dual decomposition and decentralized alternating direction method of multipliers (ADMM) . Specifically, the decentralized ADMM can employ a fixed step-size and it has nice provable rates . Under the strong convexity assumption, the decentralized ADMM has been shown to have linear convergence for time-invariant undirected graphs . Building on (augmented) Lagrangian, a few improvements have been made via proximal-gradient , stochastic gradient , and second-order approximation . In particular, ADMM over a random undirected network has been shown to have O(1/k)O(1/k) rate for convex functions . However, de-synchronization and extensions of these methods to time-varying undirected graphs are more involved , while their extensions to directed graphs are non-existent in the current literature.

Some distributed methods exist that do not (explicitly) use dual variables but can still converge to an exact consensual solution while using fixed step-sizes. In particular, work in employs multi-consensus inner loop and Nesterov’s acceleration method, which gives a proximal-gradient algorithm with a rate at least O(1/k)O(1/k). By utilizing multi-consensus inner loop, adapt-then-combine (ATC) strategy, and Nesterov’s acceleration, the algorithm proposed in is shown to have O(ln(k)/k2)O\left(\ln(k)/k^{2}\right) rate under the assumption of bounded and Lipschitz gradients. For least squares, the general diffusion strategy (a generalization of ATC) can converge to the global minimizer‘. Although it is unknown to the literature, the above algorithms that do not use dual variable but use fixed step-size are not likely to reach linear convergence even under the strong convexity assumption. References use a difference structure to cancel the steady state error in decentralized gradient descent , thereby developing the algorithm EXTRA and its proximal-gradient variant. EXTRA converges at an o(1/k)o(1/k) rate when the objective function in (1) is convex, and it has a Q-linear rate when the objective function is strongly convex.

Aside from the diminishing step-size issue, another topic of interest is distributed optimization over time-varying directed graphs. The distributed algorithms over time-varying graphs require the use of doubly stochastic weight matrices, which are not easily constructed in a distributed fashion when the graphs are directed. To overcome this issue, reference is the first to propose a different distributed approach, namely a subgradient-push algorithm that combines the distributed subgradient method with the push-sum protocol . While the subgradient-push eliminates the requirement of graph balancing , it suffers from a slow sublinearWhen an algorithm has convergence rate of O(θ(k))O(\theta(k)), we say that the rate is sublinear if limk+λkθ(k)=0\lim_{k\rightarrow+\infty}\frac{\lambda^{k}}{\theta(k)}=0 for any constant λ(0,1)\lambda\in(0,1). A typical sublinear rates include O(1/kp)O(1/k^{p}) with p>0p>0. convergence rate even for strongly convex smooth functions due to its employment of diminishing step-size . On the other hand, noticing that EXTRA has satisfactory convergence rates for undirected graphs, references combine EXTRA with the push-sum protocol to produce DEXTRA (ExtraPush) algorithm in hope of making it work over directed graph. It turns out that for a time-invariant strongly connected directed graph, DEXTRA converges at an R-linear rate under strong convexity assumption but the step-size has to be carefully chosen in some interval. However, the feasible set of step-sizes for DEXTRA can even be empty in some situations . Finally, the paper proposed an algorithm with diminishing step-size for nonconvex optimization over directed graphs based on the push-sum method and showed convergence to a stationary point.

The algorithm we will study in this paper crucially relies on tracking differences of gradients and is a minor variation of algorithms that have appeared in and . To be specific, references utilize an Adapt-then-Combine variation of the dynamic average consensus approach , and thereby develop an Aug-DGM algorithm which is capable of employing uncoordinated step-sizes for multi-agent optimization. Independently, a scheme based on difference of gradients was proposed for the more general class of non-convex functions in where a large class of distributed algorithms is developed. Finally, we note that reference , appearing simultaneously with this work, also proposed a method for distributed optimization based on gradient differences. However, none of the papers mentioned in this paragraph provide a theoretical analysis for strongly convex optimization problems over time-varying graphs, which is the goal of this paper.

2 Summary of Contributions

Prior to this work, an open question was how to construct a linearly convergent method for distributed optimization over time-varying (undirected or directed) graphs.

The present paper resolves the issue. Specifically, we construct distributed methods which are linearly convergent over graphs which are time-varying and directed. Furthermore, we show that when the graphs are time-varying and undirected, a particular (distributed) choice of weights results in a polynomial iteration complexity (meaning that the number of iterations until the protocol reaches any fixed accuracy is polynomial in the total number of nodes).

Our protocols work for all step-sizes small enough. To set the step-size to achieve the convergence rate guaranteed by our theorems, nodes need to know (i) an upper bound on the total number of agents and (ii) an upper bound on the number of time steps needed to achieve long-term connectivity. This compares favorably to some of the existing literature, e.g., , which require detailed spectral information about the network for step-size selection.

Moreover the technical tools we use are of independent interest. Although linearly convergent distributed optimization methods over fixed graphs were first developed in , extending the proof of to time-varying graphs does not appear to be possible. The current paper develops a new approach to the problem based on the small-gain theorem, a standard tool for proving stability of interconnected dynamical systems in control theory. In fact, to our knowledge, our work is the first to use a small-gain based analysis to show convergence (and to bound convergence time) of an optimization protocol.

3 Paper Organization

The rest of this paper is organized as follows. To facilitate the description of the technical ideas, the algorithms, and the analysis, we first introduce the notation in Subsection 1.4. In Section 2 we consider the case of undirected time-varying graphs, and we introduce a distributed consensus-based algorithm in Section 2.1. The algorithm uses “distributed inexact gradients” and, also, employs a “gradient tracking” technique, thus we term the algorithm as DIGing to account for its main design features. In Section 3 we establish that the DIGing algorithm converges at an R-linear rate under standard assumptions including uniform joint strong connectivity of the graphs, the strong convexity of the objective function, and the Lipschitz continuity of the gradients. Moreover, we show that the convergence rate of DIGing scales polynomially in the total number of agents in the network. After this, we consider the case of time-varying directed graphs, and we propose a push-sum consensus-based variant of DIGing in Section 4. We establish its R-linear rate in Section 5. Finally, some numerical simulations are given in Section 6, and the paper concludes with some final remarks in Section 7.

4 Notation

respectively. Each row ii of x\mathbf{x} and f(x)\nabla\mathbf{f}(\mathbf{x}) is associated with agent ii. We say that x\mathbf{x} is consensual if all of its rows are identical, i.e., x1=x2==xnx_{1}=x_{2}=\cdots=x_{n}. The analysis and results of this paper hold for all p1p\geq 1. The reader can assume p=1p=1 for convenience (so x\mathbf{x} and f\nabla\mathbf{f} become vectors) without loss of generality. The notation is not standard but it enables us to present our algorithm and analysis in a compact form.

Distributed Optimization over Undirected Graphs

In this section, we consider the case when the agents want to jointly solve the problem (1) while interacting over time-varying undirected graphs. We describe our algorithm, provide an interpretation of its steps and discuss its connection to some of the existing approaches.

We introduce the algorithm and provide some insights into its iterations. In what follows we will make use of the following proposition, which provides the optimality conditions for problem (1).

To illustrate the idea of the DIGing algorithm, let us focus on the case of a static graph for the moment. Consider the distributed gradient descent (DGD), given as follows:

where W\mathbf{W} is a doubly stochastic mixing matrix and α>0\alpha>0 is a fixed step-size. The mixing part “Wx(k)\mathbf{W}\mathbf{x}(k)” is necessary for reaching consensus while DGD exhibits undesirable behavior due to its use of the gradient direction, “αf(x(k))-\alpha\nabla\mathbf{f}(\mathbf{x}(k))”. To see this, let us break the update into steps per agent: for every agent ii, we have xi(k+1)=Wiixi(k)+jNiWijxj(k)αfi(xi(k))x_{i}(k+1)=W_{ii}x_{i}(k)+\sum_{j\in\mathcal{N}_{i}}W_{ij}x_{j}(k)-\alpha\nabla f_{i}(x_{i}(k)), where Ni\mathcal{N}_{i} is the set of the neighbors of agent ii in the given graph. Thus, each agent is updating using only the gradient of its local objective function fif_{i}. Suppose now that the values xi(k)x_{i}(k) have reached consensus and that xi(k)=xx_{i}(k)=x^{*} for all ii and some solution xx^{*} of the problem (1). Then, the mixing part gives Wiixi(k)+jNiWijxj(k)=xW_{ii}x_{i}(k)+\sum_{j\in\mathcal{N}_{i}}W_{ij}x_{j}(k)=x^{*} for all ii. However, the gradient-based term gives αfi(xi(k))-\alpha\nabla f_{i}(x_{i}(k)) for all ii, which need not be zero in general, thus resulting in xi(k+1)x_{i}(k+1) that will move away from the solution xx^{*} (recall that a solution to the problem (1) is at a point xx where j=1nfj(x)=0 i\sum_{j=1}^{n}\nabla f_{j}(x)=0\ \forall i and not necessarily a point where fi(x)=0 i\nabla f_{i}(x)=0\ \forall i).

Conceptually, one (non-distributed) scheme that bypasses this limitation is the update

which can be implemented if every agent has access to the average of all the gradients fj(xj(k))\nabla f_{j}(x_{j}(k)), j=1,,nj=1,\ldots,n (evaluated at each agent’s local copy). One can verify that if (2) converges, its limit point x()\mathbf{x}(\infty) satisfies the optimality conditions as given in Proposition 2.1. However, the update in (2) is not distributed among the agents as it requires a central entity to provide the average of the gradients.

Nevertheless, one may approximate the update in (2) through a surrogate direction that tracks the gradient average. To track the average of the gradients, namely, 1n11f(x(k))\frac{1}{n}\mathbf{1}\mathbf{1}^{\top}\nabla\mathbf{f}(\mathbf{x}(k)), we introduce a variable y(k)\mathbf{y}(k) that is updated as follows:

which is exactly what we need in view of (2). Replacing 1n11f(x(k))\frac{1}{n}\mathbf{1}\mathbf{1}^{\top}\nabla\mathbf{f}(\mathbf{x}(k)) in (2) by its dynamic approximation y(k)\mathbf{y}(k) is exactly what we use to construct the DIGing algorithm. Furthermore, to accommodate time-varying graphs, the static weight matrix W\mathbf{W} is replaced by a time varying matrix W(k)\mathbf{W}(k), thus resulting in the DIGing algorithm, as given below.

where Ni(k)\mathcal{N}_{i}(k) is the set of all neighbors of agent ii at time kk. At every iteration kk, each agent ii sends its current solution estimate xi(k)x_{i}(k) and average gradient estimate yi(k)y_{i}(k) to all of its neighbors Ni(k)\mathcal{N}_{i}(k) while receiving all of its neighbors’ solution estimates xj(k)x_{j}(k) and average gradient estimates yj(k)y_{j}(k), jNi(k)\forall j\in\mathcal{N}_{i}(k). Then, each agent ii updates its vector xi(k+1)x_{i}(k+1) by mixing its own xi(k)x_{i}(k) and the neighbors’ copies xj(k),x_{j}(k), for jNi(k)j\in\mathcal{N}_{i}(k), with specific weights, and adjusting along the direction of yi(k)-y_{i}(k). Also, each agent ii updates its direction yi(k+1)y_{i}(k+1) by mixing its own yi(k)y_{i}(k) and the neighbors’ directions yj(k)y_{j}(k), for jNi(k)j\in\mathcal{N}_{i}(k) with specific weights, and by taking into account only the new information contained in the most recent gradient evaluation, as captured in the gradient-difference term fi(xi(k+1))fi(xi(k))\nabla f_{i}(x_{i}(k+1))-\nabla f_{i}(x_{i}(k)).

2 Relation of DIGing to some of the existing approaches

This section explains how the introduced DIGing algorithm is related to some of the other distributed algorithms.

When W(k)=W\mathbf{W}(k)=\mathbf{W} (time-invariant case) and W\mathbf{W} is symmetric, the DIGing algorithm shares some similarity with EXTRA. If we eliminate the variables y(k)\mathbf{y}(k) in the recursion of DIGing, we will obtain

In this form, (2WI)(2\mathbf{W}-I) and W2\mathbf{W}^{2} will be the two mixing matrices in EXTRA. As long as we have

the convergence properties of DIGing will follow from the results in . It can be seen that, when W0\mathbf{W}\succ 0, the convergence of DIGing follows from the convergence of EXTRA immediately. In this paper, we conduct the convergence analysis with more general choices of time-varying W(k)\mathbf{W}(k).

2.2 Connection with primal-dual approaches

DIGing has a primal-dual interpretation when the mixing matrices are static and symmetric. Indeed, suppose that W(k)=W\mathbf{W}(k)=\mathbf{W} where W\mathbf{W} is a symmetric doubly stochastic matrix, and consider the augmented Lagrangian function

If we apply the basic gradient method with a step-size α\alpha in Gauss-Seidel-like order for computing the saddle point of the augmented Lagrangian function (4), we will have

By eliminating the dual variable r\mathbf{r}, we will obtain the same updates as in the DIGing algorithm (where the variable y\mathbf{y} is eliminated). The same will happen if, alternatively, we consider the augmented Lagrangian function

and apply the basic gradient method with a step-size α\alpha in Jacobi-like order for seeking the saddle point of the augmented Lagrangian function (5). Similar connections between the EXTRA algorithm and a general primal-dual have been made in a few recent references, , , and .

The above discussion assumes a time-invariant matrix W\mathbf{W} which is symmetric. Even in the time-invariant graph case, but for asymmetric W\mathbf{W}, it appears to be difficult to adapt the classical primal-dual analysis for the recursion of DIGing. Our arguments in this paper are from a pure primal perspective and do not assume any symmetry property of W(k)\mathbf{W}(k). This suggests the possibility of its extension to the case of directed graphs, which will be addressed in Section 4.

2.3 Connection with Aug-DGM [74]

A very recent reference that comes to our attention is that proposes a distributed consensus optimization algorithm, Aug-DGM, which is applicable to general time-invariant graphs and is quite similar to DIGing. The proposed algorithm is based on the combination of Adapt-then-Combine (ATC) strategy of and the dynamic average consensus for gradient tracking of. It differs from DIGing only in the dynamic average gradient-consensus update, which uses an ATC variant. The updates of Aug-DGM are given by

where W\mathbf{W} is a doubly stochastic matrix and D\mathbf{D} is a diagonal step-size matrix. When D\mathbf{D} is chosen as αI\alpha I, it turns into an ATC variant of DIGing. With a general (positive) diagonal matrix D\mathbf{D}, Aug-DGM allows different agents to use different step-sizes and it still drives all the agent to reach a consensus on a global minimizer. The convergence of Aug-DGM is provided under general convexity and Lipschitz gradient assumptions.

Convergence Analysis for DIGing over Undirected Graphs

In the sequel, we use the following notation:

with the convention that Wb(k)=I\mathbf{W}_{b}(k)=I for any needed k<0k<0 and W0(k)=I\mathbf{W}_{0}(k)=I for any kk.

Next, we give the basic assumption that we impose on the weight matrices.

(Decentralized property) If iji\neq j and the edge (j,i)E(k)(j,i)\notin\mathcal{E}(k), then Wij(k)=0W_{ij}(k)=0;

(Double stochasticity) W(k)1=1\mathbf{W}(k)\mathbf{1}=\mathbf{1}, 1W(k)=1\mathbf{1}^{\top}\mathbf{W}(k)=\mathbf{1}^{\top};

(Joint spectrum property) There exists a positive integer BB such that

In Assumption 1, item (i) is due to the physical restriction of the network. Properties (ii) and (iii) are commonly used in the analysis of the rate of consensus algorithms. Several different mixing rules exist that yield the matrix sequences which have property (iii) (see subsection 2.4 of reference ).

In particular, the following two assumptions taken together imply Assumption 1 .

(Double stochasticity) W(k)1=1\mathbf{W}(k)\mathbf{1}=\mathbf{1}, 1W(k)=1\mathbf{1}^{\top}\mathbf{W}(k)=\mathbf{1}^{\top};

(Positive diagonal) For all ii, Wii(k)>0\mathbf{W}_{ii}(k)>0;

(Edge utilization) If (j,i)E(k)(j,i)\in\mathcal{E}(k), then Wij(k)>0W_{ij}(k)>0; otherwise Wij(k)=0W_{ij}(k)=0;

(Non-vanishing weights) There exists some τ>0\tau>0 such that if Wij(k)>0W_{ij}(k)>0, then Wij(k)τW_{ij}(k)\geq\tau;

Assumption 3 is strong but typical for multi-agent coordination and optimization. For undirected graph it can be fulfilled, for example, by using Metropolis weights:

where di(k)=Ni(k)d_{i}(k)=|\mathcal{N}_{i}(k)| is the degree of agent ii at time kk. In this case, Assumption 3 will be satisfied with the choice of τ=1/n\tau=1/n.

in the context wherever Assumption 2 is used.

The following lemma provides an important relation for later use.

Under Assumption 1, for any k=B1,B,k=B-1,B,\ldots, and any matrix b\mathbf{b} with appropriate dimensions, if a=WB(k)b\mathbf{a}=\mathbf{W}_{B}(k)\mathbf{b}, then we have a\Lδ(k)b\L\|\mathbf{a}\|_{\emph{{{\L}}}}\leq\delta(k)\|\mathbf{b}\|_{\emph{{{\L}}}}, where δ(k)\delta(k) is as given in Assumption 1(iii).

We do not claim any originality of this lemma. This lemma is fairly standard in consensus theory and it is a direct consequence of Assumption 1 due to the fact that WB(k)\mathbf{W}_{B}(k) is doubly stochastic:

To make our arguments more concise, we will use δ=supkB1{δ(k)}\delta=\sup_{k\geq B-1}\{\delta(k)\} in our analysis of the algorithm. An explicit expression of δ\delta in terms of nn can be found in if the more specific Assumption 3 is made.

We also need the following two assumptions on the objective functions, which are standard for deriving linear (geometric) rate of gradient algorithms for minimizing strongly convex smooth functions.

When Assumption 4 holds, we will also say that each fi\nabla f_{i} is LiL_{i}-Lipschitz (continuous). In the forthcoming analysis, we will use Lmaxi{Li}L\triangleq\max_{i}\{L_{i}\}, which is the Lipschitz constant of f(x)\nabla\mathbf{f}(\mathbf{x}), and Lˉ(1/n)i=1nLi\bar{L}\triangleq(1/n)\sum_{i=1}^{n}L_{i} which is the Lipschitz constant of f(x)\nabla f(x).

where μi[0,+)\mu_{i}\in[0,+\infty) and at least one μi\mu_{i} is nonzero.

When μi>0\mu_{i}>0, we will say that fif_{i} is μi\mu_{i}-strongly convex. In the analysis we will use μ^maxi{μi}\hat{\mu}\triangleq\max_{i}\{\mu_{i}\} and μˉ(1/n)i=1nμi\bar{\mu}\triangleq(1/n)\sum_{i=1}^{n}\mu_{i}. Assumption 5 implies the μˉ\bar{\mu}-strong convexity of f(x)f(x). Under this assumption, the optimal solution to problem (1) is guaranteed to exist and to be unique since μˉ>0\bar{\mu}>0. We note that all the convergence results in our analysis are achieved under Assumption 5. We will also use κˉL/μˉ\bar{\kappa}\triangleq L/\bar{\mu}.

To establish the R-linear rate of the algorithm, one of our technical innovations will be to resort to a somewhat unusual version of small gain theorem under a well-chosen metric, whose original version has received an extensive research and been widely applied in control theory . We will give an intuition of the whole analytical approach shortly, after stating the small gain theorem at first.

Suppose that s1,,sm\mathbf{s}^{1},\ldots,\mathbf{s}^{m} are sequences such that for all positive integers KK and for each i=1,,mi=1,\ldots,m, we have an arrow sis(imodm)+1\mathbf{s}^{i}\rightarrow\mathbf{s}^{(i\mod m)+1}, that is

where the constants (gains) γ1,,γm\gamma_{1},\ldots,\gamma_{m} are nonnegative and satisfy γ1γ2γm<1\gamma_{1}\gamma_{2}\cdots\gamma_{m}<1. Then

By iterating inequality (7) for ii from mm down to 11, we obtain

Since (8) holds for all KK and its right-hand side does not depend on KK, taking KK\rightarrow\infty implies the desired relation.

2 Sketch of the Main Idea

Before summarizing our main proof idea, let us define some quantities which we will use frequently in our analysis. We define x1(x)\mathbf{x}^{*}\triangleq\mathbf{1}(x^{*})^{\top} where xx^{*} is the optimal solution of problem (1). Also, define

which is the optimality residual of the iterates x(k)\mathbf{x}(k) (at the kk-th iteration). Moreover, let us adopt the notation

and with the convention that z(0)0\mathbf{z}(0)\triangleq{\bf 0}.

where, recall, q\mathbf{q} is the difference between local copies and the global optimizer, z\mathbf{z} is the successive difference of gradients, yˇ\check{\mathbf{y}} is the consensus violation of the estimation of gradient average across agents, and xˇ\check{\mathbf{x}} is the consensus violation of local copies (see Subsection 1.4 for the definition of operator “ ˇ\check{\ }”).

Intuitively: as long q\mathbf{q} is small, the successive difference of the gradients z\mathbf{z} is small since the gradients are close to zero in the neighborhood of the optimal point; as long as the successive difference of the gradients z\mathbf{z} is small, the structure of DIGing implies that y\mathbf{y} is close to consensual; as long as y\mathbf{y} is close to consensual, then by the structure of DIGing so is x\mathbf{x}; and, finally, as long as x\mathbf{x} is close to consensual, DIGing is very similar to gradient descent and drives the distance to the optimal point q\mathbf{q} to zero and thus completes the cycle.

Note that to apply the small gain theorem, we would need to have gains (γi\gamma_{i}) that multiply to less than one. This is achieved by choosing an appropriate step-size α\alpha. Indeed, by looking at the algorithm, we can see that the step-size appears only in one place – the third arrow (i.e., the arrow yˇxˇ\check{\mathbf{y}}\rightarrow\check{\mathbf{x}}), and the dependence of the corresponding gain in that arrow is linear in α\alpha. Thus we should be able to apply the small gain theorem after choosing small enough α\alpha.

3 The Establishment of Each Arrow

We now discuss the establishment of each arrow/relation in the sketch above [cf. (9)].

The first arrow demonstrated in Lemma 3.6 is a simple consequence of Assumption 4 (namely, it is a consequence of the fact that the gradient of f{\bf f} is LL-Lipschitz).

Under Assumption 4, we have that for all K=0,1,K=0,1,\ldots and any λ(0,1)\lambda\in(0,1),

By Assumption 4, f(x)\nabla\mathbf{f}(\mathbf{x}) is LL-Lipschitz and we have

By the definition of z\mathbf{z} and q\mathbf{q}, it follows from (10) that

Taking maxk=0,1,,K1{}\max_{k=0,1,\ldots,K-1}\{\cdot\} on both sides of (11) gives

Next we provide the lemmata for the second and third arrows in the cycle (9). They are proved by an almost identical analysis based on Lemma 3.1: indeed, a glance at the structure of DIGing implies that some (semi)-norm of x\mathbf{x} can be bounded in terms of some (semi)-norm of y\mathbf{y}, while some (semi)-norm of y\mathbf{y} can be bounded in terms of some (semi)-norm of z\mathbf{z}. This is a fairly straightforward application of Lemma 3.1, which shows how multiplication by W(k)\mathbf{W}(k) shrinks the distance toward the consensus subspace.

Let Assumption 1 hold, and let δ=supkB1{δ(k)}\delta=\sup_{k\geq B-1}\{\delta(k)\}, where δ(k)\delta(k) is as given in Assumption 1(iii). Also, let λ\lambda be such that δ<λB<1\delta<\lambda^{B}<1. Then, we have for all K=0,1,K=0,1,\ldots,

The equivalent relation in DIGing involving y\mathbf{y} and z\mathbf{z} is

From (13), using Lemma 3.1, for all kB1k\geq B-1, it follows that

for k=1,,B2k=-1,\ldots,B-2. Taking the maximum over k=1,0,,B2k=-1,0,\ldots,B-2 on both sides of (16) and the maximum over k=B1,,Kk=B-1,\ldots,K in (15), and then by combining the obtained relations, we obtain

Let Assumption 1 hold, and let δ=supkB1{δ(k)}\delta=\sup_{k\geq B-1}\{\delta(k)\}, where δ(k)\delta(k) is as given in Assumption 1(iii). Furthermore, let λ\lambda be such that δ<λB<1\delta<\lambda^{B}<1. Then, we have for all K=0,1,K=0,1,\ldots,

The relation in DIGing involving x\mathbf{x} and y\mathbf{y} is given by

Noticing the similarity between (17) and (13), we omit the proof of Lemma 3.10 since it is almost identical to that of Lemma 3.8.

With all the above lemmata in place, the last arrow of our proof sketch remains to be addressed. For this, we need an interlude on gradient descent with errors in the gradient. Since this part is relatively independent from the preceding development, we provide it in the next subsection.

4 The Inexact Gradient Descent on a Sum of Strongly Convex Functions

In this subsection, we consider the basic (centralized) first-order method for problem (1) under inexact first-order oracle. To distinguish from the notation used for our distributed optimization problem/algorithm/analysis, let us make some definitions that are only used in this subsection. Problem (1) is restated as follows with different notation,

where all gig_{i}’s satisfy Assumptions 4 and 5 with fif_{i} being replaced by gig_{i}. Let us consider the inexact gradient descent (IGD) on the function gg:

where θ\theta is the step-size. Note that since this subsection has nothing to do with time-varying setup, to avoid heavy notation, we use the upper right corner pkp^{k} instead of p(k)p(k) to denote the value of pp at iteration kk. In particular, we use (p)a(p)^{a} instead of pap^{a} to denote the aa-th power of pp when it may cause confusion. Let pp^{*} be the global minimum of gg, and define

The main lemma of this subsection is stated next; it is basically obtained by following the ideas in .

where β>0\beta>0 and η>0\eta>0. Then under Assumptions 4 and 5 with fif_{i}’s replaced by gig_{i}’s, the tuple sequence {rk,pk;s1k,s2k,,snk}\{r^{k},p^{k};s_{1}^{k},s_{2}^{k},\ldots,s_{n}^{k}\} generated by the inexact gradient method (18) obeys

By assumptions, for each i{1,2,,n}i\in\{1,2,\ldots,n\} and k=0,1,k=0,1,\ldots, we have

Averaging (21) over ii through 11 to nn gives

On the other hand, we also have that for any vector Δ\Delta,

where η>0\eta>0 is some tunable parameter, and therefore

Averaging (23) over ii through 11 to nn gives

Next, in (25), we substitute (22) for the second term, and we substitute (24) with Δ=pk+1pk\Delta=p^{k+1}-p^{k} for the third term. Thus, we obtain that

Let us look into the last two terms of (27). Noticing that μˉ=(1/n)i=1nμi\bar{\mu}=(1/n)\sum_{i=1}^{n}\mu_{i} is a strong convexity constant of g(p)g(p), there are two possibilities that could happen at time kk. Possibility A is that

while possibility B is the opposite, namely that

Considering both possibilities A and B, it follows that

Recursively using the inequality (28) we can see that

Taking square root on both sides of (29) gives us

Choose λ\lambda that satisfies (19) so to have c(λ)2(1θμˉββ+1)1c\triangleq(\lambda)^{-2}\left(1-\frac{\theta\bar{\mu}\beta}{\beta+1}\right)\leq 1, then from (30) we get

and combining it with (31), it follows that

Taking maxk=0,1,,K1{}\max_{k=0,1,\ldots,K-1}\{\cdot\} on both sides of (32) gives

5 The Last Arrow

Now we prove the last arrow of our proof sketch [cf. (9)] in the following lemma. Its establishment will use the error bound on the IGD of Lemma 3.11, as a key ingredient.

Let Assumptions 1, 4, and 5 hold. In addition, assume that the stepsize α>0\alpha>0 and the parameter λ\lambda are such that

where β>0\beta>0 and η>0\eta>0 are some tunable parameters. Then, we have

First, let us consider the evolution of x(k)\overline{x}(k). Noticing that

Applying Lemma 3.11 to the recursion relation of x\overline{x}, namely (34), we obtain

6 Linear Convergence of DIGing

We now state our main results on the convergence rates of DIGing. Our first theorem gives an explicit convergence rate for DIGing in terms of the network parameters (BB, nn, and δ\delta), objective parameters (μˉ\bar{\mu} and κˉ=Lμˉ\bar{\kappa}=\frac{L}{\bar{\mu}}), and the algorithmic step-size (α\alpha).

Suppose that Assumptions 1, 4, and 5 hold. Let

Then, for any step-size α(0,1.5(1δ)2μˉJ1]\alpha\in\left(0,\frac{1.5(1-\delta)^{2}}{\bar{\mu}J_{1}}\right], the sequence {x(k)}\{\mathbf{x}(k)\} generated by DIGing algorithm converges to the matrix x=1(x)\mathbf{x}^{*}=\mathbf{1}(x^{*})^{\top}, where xx^{*} is the unique optimal solution of problem (1) at a global R-linear (geometric) rate O(λk)O(\lambda^{k}), where the parameter λ\lambda is given by

Let us collect all the relations/arrows at hand [cf. Lemmata 3.6, 3.8, 3.10, and 3.13]:

along with other restrictions on parameters that appear in Lemmata 3.8, 3.10, and 3.13:

We next use relations (37)–(41) with a specific values for the parameters β,η\beta,\eta and the stepsize α\alpha, which yields the desired result. Specifically, let β=2L/μ^\beta=2L/\hat{\mu} and η=1\eta=1 in relation (38). By further using 0.5λ<10.5\leq\lambda<1 and (1λB)/(1λ)B(1-\lambda^{B})/(1-\lambda)\leq B, (37) and (41) together yield

On the other hand, since 1.51+1/β1.5\geq 1+1/\beta, relation (40) implies that

Using (42) and (43), it remains to show that there exists λ(δB,1)\lambda\in(\sqrt[B]{\delta},1) [cf. (39)] such that

where J1=3κˉB2(1+4nκˉ)J_{1}=3\bar{\kappa}B^{2}\left(1+4\sqrt{n}\sqrt{\bar{\kappa}}\right). We consider a smaller interval by enlarging the left-bound of the interval in (44), i.e., we will prove that

When λ\lambda varies from δB\sqrt[B]{\delta} to 11, the left-bound of the interval in (45) is monotonically decreasing from 1.5(1δ2)μˉ\frac{1.5(1-\delta^{2})}{\bar{\mu}} to , while its right-bound is monotonically increasing from to 1.5(1δ)2μˉJ1\frac{1.5(1-\delta)^{2}}{\bar{\mu}J_{1}}. In particular, the relation (45) is valid when λ\lambda (as small as the current choice of all parameters can give) is given by

we can set λ=1αμˉ1.52B\lambda=\sqrt[2B]{1-\frac{\alpha\bar{\mu}}{1.5}}, while for

we can use λ=αμˉJ11.5+δB\lambda=\sqrt[B]{\sqrt{\frac{\alpha\bar{\mu}J_{1}}{1.5}}+\delta}. The rest of the statements follow from Theorem 3.2 and Lemma 3.4.

Other possible choices of β\beta, η\eta, α\alpha, and λ\lambda exist and may give tighter bounds but here we only aim to give an explicit estimation on the rate.

To see how the geometric rate scales with the number of agents, we further have the following corollary.

Under Assumptions 2, 3, 4, and 5, if the agents choose the step-size to be

where τ\tau is the smallest nonzero positive element of the nonnegative matrices W(k)\mathbf{W}(k) for all kk [cf. Assumption 3]. Then, the sequence {x(k)}\{{\bf x}(k)\} generated by DIGing converges to the unique optimal solution x{\bf x}^{*} at a global R-linear (geometric) rate of O((λ(τ))k)O\left((\lambda(\tau))^{k}\right) where

Define φ=1λB\varphi=1-\lambda^{B}, then requiring that the interval in (45) be nonempty is equivalent to showing that the following inequality has a solution:

Therefore, we can show that an achievable φ\varphi is

By Assumption 3, from Lemma 9 of reference , we have that δ1τ2n2\delta\leq 1-\frac{\tau}{2n^{2}}. Substituting J1=3B2κˉ(1+4nκˉ)J_{1}=3B^{2}\bar{\kappa}\left(1+4\sqrt{n}\sqrt{\bar{\kappa}}\right) into (47) gives us

Thus the final rate is λ=1τ2128B2n4.5κˉ1.5B\lambda=\sqrt[B]{1-\frac{\tau^{2}}{128B^{2}n^{4.5}\bar{\kappa}^{1.5}}}, and a step-size to reach this rate is α=1.5(1λ2B)μˉ=3τ2128B2n4.5Lκˉ1.5μˉ(τ2128B2n4.5κˉ1.5)2\alpha=\frac{1.5(1-\lambda^{2B})}{\bar{\mu}}=\frac{3\tau^{2}}{128B^{2}n^{4.5}L\sqrt{\bar{\kappa}}}-\frac{1.5}{\bar{\mu}}\left(\frac{\tau^{2}}{128B^{2}n^{4.5}\bar{\kappa}^{1.5}}\right)^{2}.

Corollary 3.17 explicitly shows how the linear convergence rate of the DIGing algorithm depends on the condition number κˉ\bar{\kappa}, time-varying graph connectivity constant BB, and the network size nn. To reach ε\varepsilon-accuracy, the iteration complexity under conditions of Corollary 3.17 is O(τ2B3n4.5κˉ1.5ln1ε)O\left(\tau^{-2}B^{3}n^{4.5}\bar{\kappa}^{1.5}\ln{\frac{1}{\varepsilon}}\right) which is polynomial in the number of agents nn. Beyond the more general form of it, the advantage of Theorem 3.15 is that it explicitly depends on the parameter δ\delta which measures the convergence speed of consensus. Indeed, Corollary 3.17 uses the bound δ1τ/(2n2)\delta\leq 1-\tau/(2n^{2}) from , which may be very conservative since it applies to a rather general class of graphs. Moreover, any further advances in “consensus theory” deriving improved convergence bounds on consensus would immediately translate into improvements via Corollary 3.17 [cf. (47)], where better bounds would immediately arise with smaller values of δ\delta.

Suppose Assumptions 2, 4, and 5 hold. Also assume that the graphs are time-invariant, undirected and connected (i.e., B=1B=1). Let each W(k)\mathbf{W}(k) be a lazy Metropolis matrix, that is,

Then, with the agents choosing the step-size α(271)\alpha(\frac{2}{71}) [cf. (46) with τ=271\tau=\frac{2}{71} and B=1B=1], the sequence {x(k)}\{{\bf x}(k)\} generated by DIGing converges to the unique optimal solution x{\bf x}^{*} at a global R-linear (geometric) rate of O(λk)O(\lambda^{k}), where λ=11161312n4.5κˉ1.5\lambda=1-\frac{1}{161312n^{4.5}\bar{\kappa}^{1.5}}. In particular, the number of iterations needed to reach ε\varepsilon-accuracy is O(n4.5κˉ1.5ln1ε)O\left(n^{4.5}\bar{\kappa}^{1.5}\ln{\frac{1}{\varepsilon}}\right).

We omit the proof for Corollary 3.19 since it is essentially the same as the proof for Corollary 3.17 in addition to which we further used δ=1171n2\delta=1-\frac{1}{71n^{2}} from Lemma 2.2 of reference .

In certain cases, the bound δ1τ/(2n2)\delta\leq 1-\tau/(2n^{2}) can be conservative. In the most dramatic case, for the (fixed) complete graph this bound tells us that 1δ1-\delta is bounded below by 1/n31/n^{3}, whereas it is not too hard to see that for the complete graph 1δ1-\delta is actually bounded away from zero by a constant.

In general, however, this bound cannot be improved, in the sense that there are graphs for which it is essentially tight. For example, on the (fixed) line or ring graph, the bound gives that 1δ1-\delta is lower bounded by 1/n21/n^{2}, which is the correct scaling up to a constant, as can be seen by observing that it can take two random walks at least Ω(n2)\Omega(n^{2}) steps to intersect on these graphs, and thus the spectral gap has to be at least that much (for more details on making such arguments rigorous, we refer the reader to the introductory chapter of ).

In general, it is not possible to give a nonconservative bound on δ\delta in terms of combinatorial features of the underlying graphs as well as the number of nodes, even for fixed matrices. Such bounds have been explored at great length within Markov chain literature, see e.g., the seminal paper or the monograph , resulting in many different techniques, each giving accurate bounds on some graphs, and others.

Nevertheless, for many sequences when the graph is fixed, the quantity δ\delta can be bounded accurately. The key tool is a connection between δ\delta and the average hitting time (Theorem 11.11 of ) which in turn can be bounded in terms of graph resistance (see ). Using this connection, the following scalings may be obtained; we omit the details.

On the complete graph, δ1Ω(1)\delta\leq 1-\Omega(1).

On the line and ring graphs, we have δ1Ω(1/n2)\delta\leq 1-\Omega(1/n^{2}).

On the 2D grid and on the complete binary tree, δ1Ω(1/(nlogn))\delta\leq 1-\Omega(1/(n\log n)).

On any regular graph, δ1Ω(1/n2)\delta\leq 1-\Omega(1/n^{2}) as a consequence of the hitting time bounds of .

On the star graph and on the two-star graph (defined to be two stars on n/2n/2 nodes with a connection between the centers), δ1Ω(1/n2)\delta\leq 1-\Omega(1/n^{2}).

Our analytical framework also applies to Aug-DGM of when its step-size matrix D\mathbf{D} is set to αI\alpha I (see Subsection 2.2.3). We also find that when the graph is well connected, the ATC strategy employed in Aug-DGM can improve the linear rate. But due to space limit, we omit any detailed discussion on this aspect. In the following design of a variant of DIGing for directed graphs, we also partially employed the ATC strategy to accelerate the convergence. Numerical experiments demonstrating the faster convergence of the ATC variant of DIGing (abbreviated as DIGing-ATC) are also provided in Section 6.

Allowing the agents to use uncoordinated (different) step-sizes (all satisfying an upper bound) is possible; see our very recent subsequent work as well as the related papers . To keep the analysis in this paper concise, we do assume all agents use the same step-size α\alpha.

The step size selection, as instructed by the theorems and corollaries in this paper, requires the agents to agree on a few parameters through preprocessing. It suffices for agents to know a common upper bound on the number of nodes in the network nn and on the connectivity constant BB. The other parameters used in step-size selection are (the lower bound of) μˉ\bar{\mu} and (the upper bound of) LL. Since minima (and maxima) can be easily computed in a time-varying BB-connected network in O(nB)O(nB) rounds by an algorithm wherein each node repeatedly takes minima (or maxima) of in-neighbors, the amount of pre-processing in computing μˉ,L\bar{\mu},L is essentially negligible compared to the worst-case running time of the protocol.

It is possible to improve on this and obtain a geometric rate in tt with probability one. We sketch how this can be done next. The constant BB should be chosen so that the graph obtained by taking all the edges that appear over BB steps is connected with probability at least 1/21/2. This allows us to choose BB to be a constant independent of nn and tt. We then apply the arguments of this paper to the quantities E[q(k)],E[z(k)],E[y(k)],E[[x(k)]E[||\mathbf{q}(k)||],E[||\mathbf{z}(k)||],E[||\overline{\mathbf{y}}(k)||],E[[\overline{\mathbf{x}}(k)||] and use the small gain theorem to obtain that all of these quantities converge to zero at a geometric rate. Markov’s inequality followed by the Borel-Cantelli lemma then yields that asymptotically x(k)\mathbf{x}(k) converges to x{\bf x}^{*} at a geometric rate with probability one .

Distributed Optimization over Directed Graphs

We now focus our attention to directed graphs. For such graphs we want to design a distributed algorithm that can work with mixing matrices that need not be doubly stochastic. To do so, we employ the idea of push-sum protocol which relaxes the requirement of doubly stochastic mixing matrices to column stochastic matrices to achieve average consensus. We then introduce our algorithm that uses push-sum protocol for tracking the gradient average in time. The resulting algorithm is termed Push-DIGing (Algorithm 2), which we analyze later on in Section 5.

Instead, if every agent knows its out-degree, it is possible for the agents to construct a column stochastic matrix C\mathbf{C} and perform the following recursions with initialization u(0)=x(0)\mathbf{u}(0)=\mathbf{x}(0) and v(0)=1\mathbf{v}(0)=\mathbf{1} to achieve the average (push-sum protocol ):

Intuitively, noticing that u(k+1)=u(k)\overline{u}(k+1)=\overline{u}(k) (u(k)\mathbf{u}(k) is sum preserving), rows of u(k)\mathbf{u}(k) are heading towards scaled averages with uneven scaling ratios across the vertices caused by the non-double stochasticity of C\mathbf{C} (the ratios are actually the elements of a right eigenvector of C\mathbf{C} corresponding to the eigenvalue 11), while V(k)\mathbf{V}(k) is recording the ratios. By applying the recorded ratio inverse (V(k))1(\mathbf{V}(k))^{-1} on u(k)\mathbf{u}(k), the algorithm recovers the unscaled average of the rows in x(k)\mathbf{x}(k).

2 The Push-DIGing Algorithm

Next we formally state Push-DIGing in Algorithm 2.

Convergence Analysis for Push-DIGing

We make the following two assumptions for the setup of time-varying directed graphs.

Assumption 6 has been used in distributed optimization over time-varying directed graphs . Similar to the case of undirected graphs, we may have other options for the weights Cij(k)C_{ij}(k). In the existing literature on push-sum consensus protocol, the best understood are the matrices relying on the out-degree information (as in Assumption 7), which we use in establishing the bound on convergence rate of push-sum (see Lemma 5.1 and its proof). Generalizations of Assumption 7 may be of their own interest for the push-sum consensus algorithm, but that is beyond the scope of this paper.

A little algebra shows that the recursion relation of Push-DIGing is equivalent to

where R~(k)=(V(k+1))1C(k)(V(k))\widetilde{\mathbf{R}}(k)=(\mathbf{V}(k+1))^{-1}\mathbf{C}(k)(\mathbf{V}(k)) and h(k)=(V(k))1y(k)\mathbf{h}(k)=(\mathbf{V}(k))^{-1}\mathbf{y}(k). We note here that, under Assumptions 6 and 7, it can be seen that each matrix V(k)\mathbf{V}(k) is invertible, and that

where BB_{\ominus} is the graph connectivity constant defined in Assumption 6. The preceding relation follows from Corollary 2(b) of (here we borrow the notation for the special norm defined in (6)). Also, we note that R~(k)\widetilde{\mathbf{R}}(k) is actually a row stochastic matrix, i.e., every row of R~(k)\widetilde{\mathbf{R}}(k) sums to 11 (see Lemma 4 of).

In what follows, we will use following notation

for any k=0,1,k=0,1,\ldots and b=0,1,b=0,1,\ldots with the convention that Cb(k)=I\mathbf{C}_{b}(k)=I for any needed k<0k<0 and C0(k)=I\mathbf{C}_{0}(k)=I for any kk. The same notation rule applies to R~(k)\widetilde{\mathbf{R}}(k). In the sequel, we will give an upper bound on a norm of (I11/n)R~B(k)(I-\mathbf{1}\mathbf{1}^{\top}/n)\widetilde{\mathbf{R}}_{B}(k), as provided in the next lemma. This lemma comes from the properties of push-sum protocol and can be obtained from references .

Under Assumptions 6 and 7, let BB be an integer satisfying BBB\geq B_{\ominus} and such that

Then, for any k=B1,B,k=B-1,B,\ldots and any matrix b\mathbf{b} with appropriate dimensions, if a=R~B(k)b\mathbf{a}=\widetilde{\mathbf{R}}_{B}(k)\mathbf{b}, i.e., a=(V(k+1))1CB(k)V(k+1B)b\mathbf{a}=(\mathbf{V}(k+1))^{-1}\mathbf{C}_{B}(k)\mathbf{V}(k+1-B)\mathbf{b}, we have a\Lδb\L\|\mathbf{a}\|_{\emph{{{\L}}}}\leq\delta\|\mathbf{b}\|_{\emph{{{\L}}}}.

The convergence analysis of Push-DIGing will be based on analogous recursions illustrated in (48). Similar to the proof in Section 3, we will follow the proof sketch of the small gain theorem around the cycle:

In consensus-based algorithms for optimization over directed graphs, it is difficult to construct monotonically decreasing Lyapunov functions for convergence analysis due to the presence of asymmetric operators in the iterations (arising from the asymmetric weight matrices). Also, to deal with time-varying graphs, one has to resort to time-varying or ergodic metrics which are not easy to construct when the consensus protocol is combined with an optimization algorithm. In this situation, conventional approaches that heavily rely on every step contraction for proving Q-linear convergence are usually inapplicable. However, by defining a special metric and utilizing the small gain theorem, we manage to conveniently analyze the introduced algorithms and establish their linear convergence rates, but without relying on the monotonic decay of any Lyapunov function associated with the recursion.

Noticing that that Lemma 3.6 is a simple consequence of Assumption 4, so it also holds for Push-DIGing. For the sake of reference convenience, we restate it as follows without proof.

Under Assumption 4, we have that for all K=0,1,K=0,1,\ldots and any λ(0,1)\lambda\in(0,1),

The next two lemmata are provided by doing almost identical arguments as those for Lemmata 3.8 and 3.10: indeed, by noticing the similarity between the equivalent recursion of Push-DIGing (48) and the recursion of DIGing (see Algorithm 1), similar bounds to those in Lemmata 3.8 and 3.10 should be obtained by an application of Lemma 5.1, which shows how multiplication by a row stochastic matrix R~(k)\widetilde{\mathbf{R}}(k) shrinks the distance to the consensus subspace.

Let Assumptions 6 and 7 hold, and let λ\lambda be such that δ<λB<1\delta<\lambda^{B}<1, where BB is the constant provided in Lemma 5.1. Then, we have

where Q1Q_{1} and V1max1\|\mathbf{V}^{-1}\|_{\max}^{1} are the constants defined by (50) and (49), respectively.

The equivalent recursion of Push-DIGing involving h\mathbf{h} and z\mathbf{z} is

From (54), using Lemma 5.1, for all kB1k\geq B-1, it follows that

where Q1Q_{1} and V1max1\|\mathbf{V}^{-1}\|_{\max}^{1} are the constants defined in (50) and (49), respectively.

Noticing the similarity between (55) and (14), by the same argument as we have applied in the proof of Lemma 3.8 starting from (15), we can obtain (53).

Let Assumptions 6 and 7 hold, and let λ\lambda be such that δ<λB<1\delta<\lambda^{B}<1, where BB is the constant provided in Lemma 5.1. Then, we have

for all K=0,1,K=0,1,\ldots, where Q1Q_{1} is the constant as introduced in Lemma 5.5 (see (50)).

The equivalent recursions of Push-DIGing involving x\mathbf{x} and h\mathbf{h} is

Noticing the similarity between (57) and (54), with almost identical argument as that illustrated in the proof of Lemma 5.5, we can get (56).

Similar to the proof of Lemma 3.13, Lemma 3.11 also serves as a key ingredient in establishing the last arrow of the proof sketch for Push-DIGing. We state it in the form of a lemma as follows.

Let Assumptions 4, 5, and 7 hold. Also, assume that

where β>0\beta>0 and η>0\eta>0 are some tunable parameters. Then, we have that for all K=0,1,K=0,1,\ldots,

First, by the same argument as in the proof of Lemma 3.13, we have

Applying Lemma 3.11 to the recursion relation of u\overline{u}, namely (59), we obtain

Let us look into the summation in the last term of (60). Since u(k)=V(k)x(k)\mathbf{u}(k)=\mathbf{V}(k)\mathbf{x}(k), it follows that

Finally, we can bound the last term in (62) as follows

2 Linear Convergence of Push-DIGing

Next, we provide a convergence rate estimate for the Push-DIGing algorithm.

Suppose Assumptions 4, 5, 6, and 7 hold. Let BB be a large enough integer constant such that

Also define the constant J2J_{2} as follows:

Then, for any step-size α(0,1.5(1δ)2μˉJ2]\alpha\in\left(0,\frac{1.5(1-\delta)^{2}}{\bar{\mu}J_{2}}\right], the sequence {x(k)}\{\mathbf{x}(k)\} generated by Push-DIGing converges to the unique optimal solution x\mathbf{x}^{*} at a global R-linear (geometric) rate O(λk)O(\lambda^{k}), where λ\lambda is given by

The proof is similar to that of Theorem 3.15. Specifically, we collect all the gains as follows:

To apply the small gain theorem (Theorem 3.2), we need γ1γ2γ3γ4<1\gamma_{1}\gamma_{2}\gamma_{3}\gamma_{4}<1, that is,

The other restrictions on α\alpha, β\beta, η\eta and λ\lambda are the same as those in (38), (39), (40), and (41). Choosing β=2L/μ^\beta=2L/\hat{\mu} and η=1\eta=1, and further using 0.5λ<10.5\leq\lambda<1 and (1λB1)/(1λ)B1(1-\lambda^{B-1})/(1-\lambda)\leq B-1, relations (5.11) and (41) together imply that

Next, similar to the proof of Lemma 3.4 [cf. (45)], it remains to show that there exists some λ(δB,1)\lambda\in(\sqrt[B]{\delta},1) such that

where J2=3Q1V1max1κˉB(δ+Q1(B1))(1+n)(1+4nκˉ)J_{2}=3Q_{1}\|\mathbf{V}^{-1}\|_{\max}^{1}\bar{\kappa}B(\delta+Q_{1}(B-1))(1+\sqrt{n})\left(1+4\sqrt{n}\sqrt{\bar{\kappa}}\right). Noticing the similarity between (65) and (45), by the same argument as in the proof of Theorem 3.15 (starting from (45)), we can obtain the statement of the theorem.

Numerical Experiments

We plot push-gradient method , DIGing, DIGing-ATC (Aug-DGM with step-size matrix αI\alpha I), and Push-DIGing in Fig. 1 and Fig. 2. In the experiments, we observe that DIGing and its variants all have R-linear rates while push-gradient method only has sublinear rate even if the objective is smooth and strongly convex. Since our analyses in the above theorems and corollaries have all been worst-case analyses, the bounds on step-sizes we have given are often conservative. We therefore choose to use hand-optimized step-sizes for the numerical experiments.

We first discuss simulation performed over an undirected communication graph. In this case, we have also tested the EXTRA algorithm from reference . Interestingly, we find the convergence curve of the EXTRA algorithm almost identicalSince the two curves are almost identical, EXTRA is not plotted in the figure. to that of the DIGing algorithm when the same step size is chosen for both algorithms. As commented in Remark 3.21, we have also conducted the numerical experiments for DIGing-ATC and indeed find it faster, in terms of number of iterations, than DIGing (see the left sub-figure of Fig. 2). However, note that the per-iteration communication cost of DIGing-ATC is higher.

For an undirected graph, Push-DIGing reduces to DIGing with partial ATC strategy which only needs one round of communication per iteration. We also plot Push-DIGing for undirected graphs in the left of Fig. 2 and find that Push-DIGing can be as fast as DIGing-ATC in terms of number of iteration but at actually at a half communication cost per iteration. We speculate that this might be because after using the ATC strategy once (in Push-DIGing), the bound of step size/convergence rate has reached the limit of centralized gradient descent and no more improvement can be obtained with ATC structure.

Although DIGing and DIGing-ATC are designed for undirected graphs, we also test them over directed graphs using an asymmetric doubly stochastic matrix. The results are shown in the right sub-figure of Fig. 1. Note that constructing an asymmetric doubly stochastic matrix over a directed graph will require a graph balancing algorithm with considerable overhead. We note that the convergence curve of the DEXTRA/ExtraPush algorithm from reference is not plotted in Fig. 1 since we did not find a stable step size for the algorithm in this case. This might be because the individual functions are not strongly convex.

Finally, for time-varying directed graphs, only Push-DIGing is plotted in the right sub-figure of Fig. 2 as this is the only algorithm over time-varying directed graphs with geometric convergence. We do not plot DIGing/DIGing-ATC since they need real time graph balancing which is not practically implementable when the graph varies.

Conclusion

In this paper, we considered a class of protocols for distributed optimization based on the idea of “distributed inexact gradient” and “gradient tracking”. Under strong convexity, we studied the convergence rates of the algorithms over time-varying directed/undirected graphs. Using the small gain theorem, we showed that our protocols converge at some global R-linear rates for strongly convex functions, and were able to obtain explicit bounds on the rates.

An open question is to obtain improved estimates on the convergence rates of our method, especially as far as scaling with the number of agents, nn, goes. Furthermore, extensions to more complex optimization models containing local constraints, couplings among agents in the objectives would be of considerable interest.

Acknowledgement

We thank César A. Uribe and Thinh T. Doan for their helpful discussions.

References