ACHIEVING GEOMETRIC CONVERGENCE FOR DISTRIBUTED
OPTIMIZATION OVER TIME-VARYING GRAPHS∗

arXiv:1607.03218v3 [math.OC] 20 Mar 2017

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI†
Abstract. This paper considers the problem of distributed optimization over time-varying graphs. For the case of
undirected graphs, we introduce a distributed algorithm, referred to as DIGing, based on a combination of a distributed
inexact gradient method and a gradient tracking technique. The DIGing algorithm uses doubly stochastic mixing
matrices and employs fixed step-sizes and, yet, drives all the agents’ iterates to a global and consensual minimizer.
When the graphs are directed, in which case the implementation of doubly stochastic mixing matrices is unrealistic,
we construct an algorithm that incorporates the push-sum protocol into the DIGing structure, thus obtaining PushDIGing algorithm. The Push-DIGing uses column stochastic matrices and fixed step-sizes, but it still converges to a
global and consensual minimizer. Under the strong convexity assumption, we prove that the algorithms converge at
R-linear (geometric) rates as long as the step-sizes do not exceed some upper bounds. We establish explicit estimates
for the convergence rates. When the graph is undirected it shows that DIGing scales polynomially in the number of
agents. We also provide some numerical experiments to demonstrate the efficacy of the proposed algorithms and to
validate our theoretical findings.
Key words. distributed optimization, time-varying graphs, linear convergence, small gain theorem, inexact gradient

1. Introduction. This paper focuses on the following distributed convex optimization problem:
n

(1)

minp f (x) =

x∈R

1X
fi (x),
n i=1

where each function fi : Rp → R is held privately by agent i to encode the agent’s objective function.
We assume that the agents are connected through a communication network which can be timevarying. The agents want to collaboratively solve the problem, while each agent can only receive/send
the information from/to its immediate neighbors (to be specified precisely soon). Problems of the
form (1) that require distributed computing have appeared in various domains including information
processing and decision making in sensor networks, networked vehicle/UAV coordination/control, as
well as distributed estimation and learning. Some examples include distributed averaging [7, 50, 72],
distributed spectrum sensing [1], formation control [49, 58], power system control [21, 55], statistical
inference and learning [20, 45, 53]. In general, distributed optimization framework fits the scenarios
where the data is collected and/or stored in a network of agents and having a fusion center is either
inapplicable or unaffordable. In such scenarios, data processing and computing is to be performed in
a distributed but collaborative manner by the agents within the network.
We assume that the functions fi 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 1 .
∗ Initially submitted to the editors on July 13, 2016. Parts of the results have been submitted to the 55th IEEE
Conference on Decision and Control and the 4th IEEE Global Conference on Signal and Information Processing for
possible presentations.
Funding: The work has been partially supported by NSF grant CNS 15-44953, Office of Naval Research under
grant no. N00014-12-1-0998, and Air Force under grant number AF FA95501510394.
† Dr. Angelia Nedić and Dr. Wei Shi are with the School of Electrical, Computer, and Energy Engineering, Arizona
State University, Tempe, AZ 85287, USA. Dr. Alex Olshevsky is with the Department of Electrical Engineering, Boston
University, Boston, MA 02215, USA. (Angelia.Nedich@asu.edu, Wilbur.Shi@asu.edu, alexols@bu.edu)
1 Suppose that a sequence {x(k)} converges to x∗ in some norm k · k. We say that the convergence is: (i) Q-linear
∗
k
≤ λ for all k; (ii) R-linear if there exists λ ∈ (0, 1) and some positive
if there exists λ ∈ (0, 1) such that kx(k+1)−x
kx(k)−x∗ k

constant C such that kx(k) − x∗ k ≤ Cλk for all k. Both of these rates are geometric, and they are often referred to as
1

2

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

1.1. Literature Review. The research on distributed optimization dates back to 1980s [2, 67].
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 [6, 8]. 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 [37–39,54] and incremental proximal methods [3, 68], while a more recent work includes incremental aggregated gradient
methods [23] and its proximal gradient variants [4]. 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 [47], while its stochastic variant was
studied in [56] and its asynchronous variant in [36] 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 [19, 40, 78]. 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 [65] and decentralized alternating direction method of multipliers (ADMM) [5, 32]. Specifically, the decentralized ADMM can employ a fixed step-size and it has
nice provable rates [69]. Under the strong convexity assumption, the decentralized ADMM has been
shown to have linear convergence for time-invariant undirected graphs [63]. Building on (augmented)
Lagrangian, a few improvements have been made via proximal-gradient [10], stochastic gradient [25],
and second-order approximation [34, 35]. In particular, ADMM over a random undirected network
has been shown to have O(1/k) rate for convex functions [25]. However, de-synchronization and extensions of these methods to time-varying undirected graphs are more involved [25, 69], 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 [11] employs multiconsensus inner loop and Nesterov’s acceleration method, which gives a proximal-gradient algorithm
with a rate at least O(1/k). By utilizing multi-consensus inner loop, adapt-then-combine (ATC)
strategy, and Nesterov’s acceleration, the algorithm proposed in [26] is shown to have O ln(k)/k 2
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‘ [59]. 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 [60–62]
use a difference structure to cancel the steady state error in decentralized gradient descent [47, 75],
thereby developing the algorithm EXTRA and its proximal-gradient variant. EXTRA converges at an
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 [40] is the first to propose a different
distributed approach, namely a subgradient-push algorithm that combines the distributed subgradient
global rates to be distinguished from the case when the given relations are valid for some sufficiently large indices k. The
difference between these two types of geometric rate is in that Q-linear rate implies monotonic decrease of kx(k) − x∗ k,
while R-linear rate does not.

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

3

method [47] with the push-sum protocol [27]. While the subgradient-push eliminates the requirement
of graph balancing [22], it suffers from a slow sublinear2 convergence rate even for strongly convex
smooth functions due to its employment of diminishing step-size [41]. On the other hand, noticing
that EXTRA has satisfactory convergence rates for undirected graphs, references [70, 76] combine
EXTRA with the push-sum protocol [27] 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 [70]. Finally, the paper [64] proposed an algorithm with diminishing
step-size for nonconvex optimization over directed graphs based on the push-sum method [27] 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 [73,74] and [15–17]. To be specific, references
[73,74] utilize an Adapt-then-Combine variation [59] of the dynamic average consensus approach [77],
and thereby develop an Aug-DGM algorithm which is capable of employing uncoordinated stepsizes for multi-agent optimization. Independently, a scheme based on difference of gradients was
proposed for the more general class of non-convex functions in [15–17] where a large class of distributed
algorithms is developed. Finally, we note that reference [52], 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.
1.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., [70,76], 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 [62], extending the proof
of [62] 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.
1.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
λk
2 When an algorithm has convergence rate of O(θ(k)), we say that the rate is sublinear if lim
k→+∞ θ(k) = 0 for any
p
constant λ ∈ (0, 1). A typical sublinear rates include O(1/k ) with p > 0.

4

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

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.
1.4. Notation. Throughout the paper, the variable x ∈ Rp in the original problem (1) is viewed
as a column vector. We let agent i hold a local copy of the variable x of the problem in (1), which
is denoted by xi ∈ Rp ; its value at iteration/time P
k is denoted by xi (k). We introduce an aggregate
objective function of the local variables: f (x) , ni=1 fi (xi ), where its argument and gradient are
defined as




— (∇f1 (x1 ))⊤ —
— x⊤
—
1
 — (∇f2 (x2 ))⊤ — 
 — x⊤
— 
2




x,
 ∈ Rn×p ,
 ∈ Rn×p and ∇f (x) , 
..
..




.
.
— x⊤
n

—

— (∇fn (xn ))⊤

—

respectively. Each row i of x and ∇f (x) is associated with agent i. We say that x is consensual if all
of its rows are identical, i.e., x1 = x2 = · · · = xn . The analysis and results of this paper hold for all
p ≥ 1. The reader can assume p = 1 for convenience (so x and ∇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.
We let 1 denote a column vector with all entries equal to one (its size is to be understood from
the context). For any matrix v ∈ Rn×p , we denote its average across the rows (a row corresponds
to an agent) as v = n1 v⊤ 1 ∈ Rp , and its consensus violation as v̌ = v − 1v ⊤ = v − n1 11⊤ v =

I − n1 11⊤ v , Lv, where L = I − n1 11⊤ is a symmetric matrix. For any n × p matrices a and b,
their inner product is denoted as ha, bi = Trace(a⊤ b). For a given matrix a, the Frobenius norm is
given by kakF , the (entry-wise) max norm is given by kakmax , while the spectral norm is given by
kak2 . kak2 equals to the largest
p singular value σmax {a}. We also use kakL to denote its L weighted
(semi)-norm, that is, kakL = ha, Lai. Note that since L = L⊤ L, we always have kakL = kLakF .
For any n × p matrices a, b, c with c = a + b, in view of the definition of the consensus violation, it
holds that
kčkF = kckL = ka + bkL ≤ kakL + kbkL = kǎkF + kb̌kF ,
and we will use this property in the analysis without specific explanation. For any matrix A ∈ Rm×n ,
null{A} , {x ∈ Rn Ax = 0} is the null space of A and span{A} , {y ∈ Rm y = Ax, x ∈ Rn } is the
linear span of all the columns of A.

2. 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.
2.1. The DIGing Algorithm. 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).
Proposition 1 ( [62]). Assume null{I − W} = span{1} where W ∈ Rn×n . If x∗ satisfies
conditions: (i) x∗ = Wx∗ (consensus), (ii) 1⊤ ∇f (x∗ ) = 0 (optimality), then the rows of x∗ are the
same as each other and the transpose of each row is an optimal solution of the 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:
x(k + 1) = Wx(k) − α∇f (x(k)),

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

5

where W is a doubly stochastic mixing matrix and α > 0 is a fixed step-size. The mixing part
“Wx(k)” is necessary for reaching consensus while DGD exhibits undesirable behavior due to its use
of the gradient direction, “−α∇f (x(k))”. To see this,
P let us break the update into steps per agent:
for every agent i, we have xi (k + 1) = Wii xi (k) + j∈Ni Wij xj (k) − α∇fi (xi (k)), where Ni is the
set of the neighbors of agent i in the given graph. Thus, each agent is updating using only the
gradient of its local objective function fi . Suppose now that the values xi (k) have reached consensus
and that xiP
(k) = x∗ for all i and some solution x∗ of the problem (1). Then, the mixing part gives
Wii xi (k) + j∈Ni Wij xj (k) = x∗ for all i. However, the gradient-based term gives −α∇fi (xi (k)) for
all i, which need not be zero in general, thus resulting in xi (k + 1) that will
Pn move away from the
solution x∗ (recall that a solution to the problem (1) is at a point x where j=1 ∇fj (x) = 0 ∀i and
not necessarily a point where ∇fi (x) = 0 ∀i).
Conceptually, one (non-distributed) scheme that bypasses this limitation is the update
(2)

1
x(k + 1) = Wx(k) − α 11⊤ ∇f (x(k))
n

which can be implemented if every agent has access to the average of all the gradients ∇fj (xj (k)),
j = 1, . . . , n (evaluated at each agent’s local copy). One can verify that if (2) converges, its limit point
x(∞) satisfies the optimality conditions as given in Proposition 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, n1 11⊤ ∇f (x(k)), we introduce a
variable y(k) that is updated as follows:
(3)

y(k + 1) = Wy(k) + ∇f (x(k + 1)) − ∇f (x(k)),

with initialization y(0) = ∇f (x(0)) and where each row i of y(k) ∈ Rn×p is associated with agent
i. A similar technique has been introduced in [77] for dynamically tracking the average state of a
multi-agent system, and for tracking some network-wide aggregate quantities in [28, 57]. If x(k + 1)
converges to some point x(∞) and the underlying graph is connected, then it can be seen that the
sequence y(k) generated by the gradient tracking procedure (3) will converge to the point y(∞) given
by
1
y(∞) = 11⊤ ∇f (x(∞)),
n
which is exactly what we need in view of (2). Replacing n1 11⊤ ∇f (x(k)) in (2) by its dynamic
approximation y(k) is exactly what we use to construct the DIGing algorithm. Furthermore, to
accommodate time-varying graphs, the static weight matrix W is replaced by a time varying matrix
W(k), thus resulting in the DIGing algorithm, as given below.
Algorithm 1: DIGing
Choose step-size α > 0 and pick any x(0) ∈ Rn×p ;
Initialize y(0) = ∇f (x(0));
for k = 0, 1, . . . do
x(k + 1) = W(k)x(k) − αy(k);
y(k + 1) = W(k)y(k) + ∇f (x(k + 1)) − ∇f (x(k));
end
Looking at an individual agent i, the initialization of DIGing uses an arbitrary xi (0) ∈ Rp and
sets yi (0) = ∇fi (xi (0)) for all i = 1, . . . , n. At each iteration k, every agent i maintains two vectors,
namely, xi (k), yi (k) ∈ Rp , which are updated as follows:
P
xi (k + 1) = Wii (k)xi (k) + j∈Ni (k) Wij (k)xj (k) − αyi (k),
P
yi (k + 1) = Wii (k)yi (k) + j∈Ni (k) Wij (k)yj (k) + ∇fi (xi (k + 1)) − ∇fi (xi (k)),

where Ni (k) is the set of all neighbors of agent i at time k. At every iteration k, each agent i sends its
current solution estimate xi (k) and average gradient estimate yi (k) to all of its neighbors Ni (k) while

6

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

receiving all of its neighbors’ solution estimates xj (k) and average gradient estimates yj (k), ∀j ∈ Ni (k).
Then, each agent i updates its vector xi (k+1) by mixing its own xi (k) and the neighbors’ copies xj (k),
for j ∈ Ni (k), with specific weights, and adjusting along the direction of −yi (k). Also, each agent i
updates its direction yi (k+1) by mixing its own yi (k) and the neighbors’ directions yj (k), for j ∈ Ni (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)).
2.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.
2.2.1. Connection with EXTRA [62]. When W(k) = W (time-invariant case) and W is
symmetric, the DIGing algorithm shares some similarity with EXTRA. If we eliminate the variables
y(k) in the recursion of DIGing, we will obtain
x(k + 2) = (I + (2W − I))x(k + 1) − W2 x(k) − α[∇f (x(k + 1)) − ∇f (x(k))].
In this form, (2W − I) and W2 will be the two mixing matrices in EXTRA. As long as we have
(2W − I) 4 W2 4 (I + (2W − I)) /2 and W2 ≻ 0,
the convergence properties of DIGing will follow from the results in [62]. It can be seen that, when
W ≻ 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).
2.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 where W is a
symmetric doubly stochastic matrix, and consider the augmented Lagrangian function
(4)

Lα (x, r) = f (x) +

1
1
hr, (I − W)xi +
kxk2I−W2 .
α
2α

If we apply the basic gradient method with a step-size α in Gauss-Seidel-like order for computing the
saddle point of the augmented Lagrangian function (4), we will have
x-update: x(k + 1) =
r-update: r(k + 1) =

x(k) − α∇f (x(k)) − (I − W)r(k) − (I − W2 )x(k),
r(k) + (I − W)x(k + 1).

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

Lα (x, r) = f (x) +

1
1
hr, (I − W)xi +
kxk22I−2W ,
α
2α

and apply the basic gradient method with a step-size α 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, [35], [33], and [24].
Remark 1 (Primal-dual and symmetry). The above discussion assumes a time-invariant matrix W which is symmetric. Even in the time-invariant graph case, but for asymmetric 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).
This suggests the possibility of its extension to the case of directed graphs, which will be addressed in
Section 4.

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

7

2.2.3. Connection with Aug-DGM [74]. A very recent reference that comes to our attention
is [74] 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 [59] and the dynamic average consensus
for gradient tracking of [77]. 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
x-update:
y-update:

x(k + 1) = W (x(k) − Dy(k)) ,
y(k + 1) = W (y(k) + ∇f (x(k + 1)) − ∇f (x(k))) ,

where W is a doubly stochastic matrix and D is a diagonal step-size matrix. When D is chosen as
αI, it turns into an ATC variant of DIGing. With a general (positive) diagonal matrix 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.
3. Convergence Analysis for DIGing over Undirected Graphs. In this section we establish the linear convergence of DIGing over time-varying undirected graphs. Let us formally describe
the assumptions we make on the graphs and on the mixing matrices W(k) that are compatible with
the graphs. Consider a time-varying undirected graph sequence {G un (0), G un (1), . . .}. Every graph instance G un (k) consists of a time-invariant set of agents V = {1, 2, . . . , n} and a set of time-varying edges
E(k). The unordered pair of vertices (j, i) ∈ E(k) if and only if agents j and i can exchange
information

at time (iteration) k. The set of neighbors of agent i at time k is defined as Ni (k) = j (j, i) ∈ E(k) .
In the sequel, we use the following notation:
n
o
[
[ [
Gbun (k) , V, E(k) E(k + 1) · · · E(k + b − 1) for any k = 0, 1, . . . and any b = 1, 2, . . .,
Wb (k) , W(k)W(k − 1) · · · W(k − b + 1) for any k = 0, 1, . . . and any b = 0, 1, . . .,
with the convention that Wb (k) = I for any needed k < 0 and W0 (k) = I for any k.
Next, we give the basic assumption that we impose on the weight matrices.
Assumption 1 (Mixing matrix sequence {W(k)}). For any k = 0, 1, . . ., the mixing matrix
W(k) = [Wij (k)] ∈ Rn×n satisfies the following relations:
(i) (Decentralized property) If i 6= j and the edge (j, i) ∈
/ E(k), then Wij (k) = 0;
(ii) (Double stochasticity) W(k)1 = 1, 1⊤ W(k) = 1⊤ ;
(iii) (Joint spectrum property) There exists a positive integer B such that


1
for all k = 0, 1, . . . .
sup δ(k) < 1 where δ(k) = σmax WB (k) − 11⊤
n
k≥B−1
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 [62]).
In particular, the following two assumptions taken together imply Assumption 1 [43].
Assumption 2 (B̃-connected graph sequence). The time-varying undirected graph sequence
{G un (k)} is B̃-connected. Specifically, there exists some positive integer B̃ such that the undirected
un
graph GB̃
(tB̃) is connected for all t = 0, 1, . . ..
Assumption 2 is typical for many results in multi-agent coordination and distributed optimization
[42]. It is much weaker than the assumption of every G un (k) being connected.
Assumption 3 (Mixing matrix sequence {W(k)}). For any k = 0, 1, . . ., the mixing matrix
W(k) = [Wij (k)] ∈ Rn×n satisfies
(i) (Double stochasticity) W(k)1 = 1, 1⊤ W(k) = 1⊤ ;

8

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

(ii) (Positive diagonal) For all i, Wii (k) > 0;
(iii) (Edge utilization) If (j, i) ∈ E(k), then Wij (k) > 0; otherwise Wij (k) = 0;
(iv) (Non-vanishing weights) There exists some τ > 0 such that if Wij (k) > 0, then Wij (k) ≥ τ ;
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:

 1/ (1 + max{di (k), dj (k)}) , if (j, i) ∈ E(k),
0,
if (j, i) ∈
/ E(k)) and j 6= i,
Wij (k) =
 1−P
W
(k),
if
j
=
i,
il
l∈Ni (k)

where di (k) = |Ni (k)| is the degree of agent i at time k. In this case, Assumption 3 will be satisfied
with the choice of τ = 1/n.
Now let us consider the graph sequence {G un (k)} under Assumption 2. If the constant B mentioned in Assumption 1 is chosen as B ≥ 2B̃ − 1, then for anySk = 0, 1, . . ., the union of a B-length
k+B−1
E(b), is always a super set of
consecutive clip of the edge set sequence from time index k, b=k
S⌈k/B̃⌉B̃+B̃−1
un
E(b) and the graph GB̃ (⌈k/B̃⌉B̃) by assumption is connected, thus it is very easy to
b=⌈k/B̃⌉B̃
un
see that the graph GB
(k) is connected. Here, the notation ⌈·⌉ is used to denote the ceiling function
which rounds a real number to the least succeeding integer. With such a choice of B, Assumptions 2
and 3 together imply Assumption 1. Regarding such a relationship, through out Section 3, we choose
B = 2B̃ − 1
in the context wherever Assumption 2 is used.
The following lemma provides an important relation for later use.
Lemma 2 (B-step consensus contraction). Under Assumption 1, for any k = B − 1, B, . . .,
and any matrix b with appropriate dimensions, if a = WB (k)b, then we have kakL ≤ δ(k)kbkL ,
where δ(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) is doubly stochastic:
kakL

=
=
≤

k(I − n1 11⊤ )WB (k)bkF
k(WB(k) − n1 11⊤ )(I − n1 11⊤ )bkF
σmax WB (k) − n1 11⊤ kbkL .

To make our arguments more concise, we will use δ = supk≥B−1 {δ(k)} in our analysis of the algorithm.
An explicit expression of δ in terms of n can be found in [43] 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.
Assumption 4 (Smoothness). For every agent i, its objective fi : Rp → R is differentiable and
has Lipschitz continuous gradients, i.e., there exists a Lipschitz constant Li ∈ (0, +∞) such that
k∇fi (x) − ∇fi (y)kF ≤ Li kx − ykF for any x, y ∈ Rp .
When Assumption 4 holds, we will also say that each ∇fi is Li -Lipschitz (continuous). In the
forthcoming
Pnanalysis, we will use L , maxi {Li }, which is the Lipschitz constant of ∇f (x), and
L̄ , (1/n) i=1 Li which is the Lipschitz constant of ∇f (x).
Assumption 5 (Strong convexity). For every agent i, its objective fi : Rp → R satisfies
fi (x) ≥ fi (y) + h∇fi (y), x − yi +
where µi ∈ [0, +∞) and at least one µi is nonzero.

µi
kx − yk2F for any x, y ∈ Rp ,
2

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

9

When µi >P0, we will say that fi is µi -strongly convex. In the analysis we will use µ̂ , maxi {µi }
n
and µ̄ , (1/n) i=1 µi . Assumption 5 implies the µ̄-strong convexity of f (x). Under this assumption,
the optimal solution to problem (1) is guaranteed to exist and to be unique since µ̄ > 0. We note
that all the convergence results in our analysis are achieved under Assumption 5. We will also use
κ̄ , L/µ̄.
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 [13]. We will give an
intuition of the whole analytical approach shortly, after stating the small gain theorem at first.
3.1. The Small Gain Theorem. Let us adopt the notation si for the infinite sequence si =
s (0), si (1), si (2), . . . where si (k) ∈ Rn×p , ∀i. Furthermore, let us define
i

ksi kλ,K
,
F

(6)

1 i
ks (k)kF
k=0,...,K λk
max

1 i
ks (k)kF ,
k
k≥0 λ

and ksi kλF , sup

where the parameter λ ∈ (0, 1) will serve as the linear rate parameter later in our analysis. While
ksi kλ,K
is always finite, ksi kλF may be infinite. If n = p = 1, i.e., each si (k) is a scalar, we will just
F
i λ,K
write |s |
and |si |λ for these quantities. Intuitively, ksi kλ,K
is a weighted “ergodic norm” of si .
F
1
Noticing that the weight λk is exponentially growing with respect to k, if we can show that ksi kλF is
bounded, then it would imply that ksi (k)kF → 0 geometrically fast. This ergodic definition enables
us to give analysis to those algorithms which do not converge Q-linearly. Next we will state the small
gain theorem which gives a sufficient condition to for the boundedness of ksi kλF . The theorem is a
basic result in control systems and a detailed discussion about its result can be found in [13]. For the
sake of completeness, we include the proof.
Theorem 3 (The small gain theorem). Suppose that s1 , . . . , sm are sequences such that for
all positive integers K and for each i = 1, . . . , m, we have an arrow si → s(i mod m)+1 , that is
ks(i mod m)+1 kλ,K
≤ γi ksi kλ,K
+ ωi ,
F
F

(7)

where the constants (gains) γ1 , . . . , γm are nonnegative and satisfy γ1 γ2 · · · γm < 1. Then
ks1 kλF

≤

1
1−γ1 γ2 ···γm (ω1 γ2 γ3 · · · γm + ω2 γ3 γ4 · · · γm + · · · + ωm−1 γm + ωm ) .

Proof. By iterating inequality (7) for i from m down to 1, we obtain
ks1 kλ,K
F

≤

γm γm−1 · · · γ1 ks1 kλ,K
+ γm γm−1 · · · γ2 ω1
F
+γm γm−1 · · · γ3 ω2 + · · · + γm ωm−1 + ωm .

Thus,
(8)

ks1 kλ,K
F

≤

1
1−γ1 γ2 ···γm (ω1 γ2 γ3 · · · γm + ω2 γ3 γ4 · · · γm + · · · + ωm−1 γm + ωm ) .

Since (8) holds for all K and its right-hand side does not depend on K, taking K → ∞ implies the
desired relation.
Clearly, the small gain theorem involves a cycle s1 → s2 → · · · → sm → s1 . Due to this cyclic
structure [cf. (7)], similar bounds hold for ksi kλF , ∀i.

Lemma 4 (Bounded norm ⇒ R-linear rate). For any matrix sequence si , if ksi kλF is bounded,
then ksi (k)kF converges at a global R-linear (geometric) rate O(λk ).

Proof. If ksi kλF ≤ U where U is some nonnegative constant, then by the definition we have
supk≥0 λ1k ksi (k)kF ≤ U , thus λ1k ksi (k)kF ≤ U, ∀k. The conclusion follows immediately from
ksi (k)kF ≤ U λk , ∀k.

10

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

3.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 x∗ , 1(x∗ )⊤ where x∗ is the optimal
solution of problem (1). Also, define
q(k) , x(k) − x∗ for any k = 0, 1, . . . ,
which is the optimality residual of the iterates x(k) (at the k-th iteration). Moreover, let us adopt
the notation
z(k) , ∇f (x(k)) − ∇f (x(k − 1)) for any k = 1, 2, . . . ,

and with the convention that z(0) , 0.
We will apply the small gain theorem with the k · kλ,K
metric and a right choice of λ < 1 around
F
the following circle of arrows:
(9)

Algorithm 1: q → z → y̌ → x̌ → q,

where, recall, q is the difference between local copies and the global optimizer, z is the successive
difference of gradients, y̌ is the consensus violation of the estimation of gradient average across agents,
and x̌ is the consensus violation of local copies (see Subsection 1.4 for the definition of operator “ˇ”).
Intuitively: as long q is small, the successive difference of the gradients 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 is small, the structure of DIGing implies that y is close to consensual; as long as
y is close to consensual, then by the structure of DIGing so is x; and, finally, as long as x is close to
consensual, DIGing is very similar to gradient descent and drives the distance to the optimal point q
to zero and thus completes the cycle.
After the establishment of each arrow, we will apply the small gain theorem to conclude that every
corresponding quantities under the metric k · kλ,K
is bounded and hence conclude that all quantities
F
in the “circle of arrows” decay at an R-linear rate O(λk ).
Note that to apply the small gain theorem, we would need to have gains (γi ) that multiply to
less than one. This is achieved by choosing an appropriate step-size α. 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̌), and the dependence of the corresponding gain in that arrow is linear in α. Thus we should
be able to apply the small gain theorem after choosing small enough α.
3.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 5 is a simple consequence of Assumption 4 (namely, it is
a consequence of the fact that the gradient of f is L-Lipschitz).
Lemma 5 (Algorithm 1: The ﬁrst arrow q → z). Under Assumption 4, we have that for all
K = 0, 1, . . . and any λ ∈ (0, 1),


1
λ,K
kqkλ,K
kzkF ≤ L 1 +
F .
λ
Proof. By Assumption 4, ∇f (x) is L-Lipschitz and we have
(10)

k∇f (x(k + 1)) − ∇f (x(k))kF

≤ Lkx(k + 1) − x(k)kF
= Lk(x(k + 1) − x∗ ) − (x(k) − x∗ )kF
≤ Lkx(k + 1) − x∗ kF + Lkx(k) − x∗ kF .

By the definition of z and q, it follows from (10) that
(11)

λ−(k+1) kz(k + 1)kF

≤

−k
kq(k)kF .
Lλ−(k+1) kq(k + 1)kF + L
λλ

Taking maxk=0,1,...,K−1 {·} on both sides of (11) gives


λ,K−1
≤ L 1 + λ1 kqkλ,K
kzkλ,K
≤ Lkqkλ,K
+L
F .
F
F
λ kqkF

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

11

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 2: indeed, a glance at the structure of DIGing
implies that some (semi)-norm of x can be bounded in terms of some (semi)-norm of y, while some
(semi)-norm of y can be bounded in terms of some (semi)-norm of z. This is a fairly straightforward
application of Lemma 2, which shows how multiplication by W(k) shrinks the distance toward the
consensus subspace.
Lemma 6 (Algorithm 1: The second arrow z → y̌). Let Assumption 1 hold, and let δ =
supk≥B−1 {δ(k)}, where δ(k) is as given in Assumption 1(iii). Also, let λ be such that δ < λB < 1.
Then, we have for all K = 0, 1, . . .,
ky̌kλ,K
≤
F

(12)

B
λ(1 − λB )
λB X 1−t
λ,K
kzkF + B
λ ky̌(t − 1)kF .
(λB − δ)(1 − λ)
λ − δ t=1

Proof. The equivalent relation in DIGing involving y and z is
(13)

y(k + 1) = W(k)y(k) + z(k + 1).

From (13), using Lemma 2, for all k ≥ B − 1, it follows that
ky̌(k + 1)kF
(14)

= ky(k + 1)kL
≤ kWB (k)y(k + 1 − B)kL + kWB−1 (k)z(k + 2 − B)kL
+ · · · + kW1 (k)z(k)kL + kW0 (k)z(k + 1)kL
B
P
≤ δky̌(k + 1 − B)kF +
kz(k + 2 − t)kF ,
t=1

and therefore, for all k = B − 1, B, . . .,
(15)

λ−(k+1) ky̌(k + 1)kF ≤ λδB λ−(k+1−B) ky̌(k + 1 − B)kF +

B
P

t=1

1
−(k+2−t)
kz(k + 2 − t)kF .
λt−1 λ

To utilize the norm k · kλ,K
F , we need to take maxk=0,...,K , which in turn requires a relation for
λ−(k+1) ky̌(k + 1)kF with k < B − 1. To obtain such a relation, we complement the initial relation for
(15), i.e.,
λ−(k+1) ky̌(k + 1)kF ≤ λ−(k+1) ky̌(k + 1)kF

(16)

for k = −1, . . . , B − 2. Taking the maximum over k = −1, 0, . . . , B − 2 on both sides of (16) and the
maximum over k = B − 1, . . . , K in (15), and then by combining the obtained relations, we obtain
ky̌kλ,K
F

≤
≤

Hence,
ky̌kλ,K
F

B
B
P
P
λ,K−B
λ,K+1−t
1
δ
+
+
λ1−t ky̌(t − 1)kF
λB ky̌kF
λt−1 kzkF
t=1
t=1
B
B
P
P
λ,K
λ,K
δ
1
ky̌k
+
kzk
+
λ1−t ky̌(t − 1)kF .
t−1
B
F
F
λ
λ
t=1
t=1
λB

≤
=

which is exactly (12).

B
P
1
λt−1
t=1
B
λ −δ

B

kzkλ,K
+ λBλ −δ
F

B
P

λ1−t ky̌(t − 1)kF
t=1
B
P
B
λ,K
λ(1−λB )
λ1−t ky̌(t − 1)kF ,
+ λBλ −δ
(λB −δ)(1−λ) kzkF
t=1

12

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

Lemma 7 (Algorithm 1: The third arrow y̌ → x̌). Let Assumption 1 hold, and let δ =
supk≥B−1 {δ(k)}, where δ(k) is as given in Assumption 1(iii). Furthermore, let λ be such that δ <
λB < 1. Then, we have for all K = 0, 1, . . .,
kx̌kλ,K
≤
F

B
α(1 − λB )
λB X 1−t
λ,K
ky̌k
+
λ kx̌(t − 1)kF .
F
(λB − δ)(1 − λ)
λB − δ t=1

The relation in DIGing involving x and y is given by
(17)

x(k + 1) = W(k)x(k) − αy(k).

Noticing the similarity between (17) and (13), we omit the proof of Lemma 7 since it is almost identical
to that of Lemma 6.
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.
3.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,
n

min g(x) =

x∈Rd

1X
gi (x),
n i=1

where all gi ’s satisfy Assumptions 4 and 5 with fi being replaced by gi . Let us consider the inexact
gradient descent (IGD) on the function g:
n

(18)

pk+1 = pk − θ

1X
∇gi (ski ),
n i=1

where θ 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 pk instead of p(k) to denote the value of p at
iteration k. In particular, we use (p)a instead of pa to denote the a-th power of p when it may cause
confusion. Let p∗ be the global minimum of g, and define
rk , kpk − p∗ kF for any k = 0, 1, . . . .
The main lemma of this subsection is stated next; it is basically obtained by following the ideas
in [14].
Lemma 8 (The error bound on the IGD). Suppose that
s
θµ̄β
1
(19)
,
1−
≤ λ < 1 and θ ≤
β+1
(1 + η)L̄
where β > 0 and η > 0. Then under Assumptions 4 and 5 with fi ’s replaced by gi ’s, the tuple sequence
{rk , pk ; sk1 , sk2 , . . . , skn } generated by the inexact gradient method (18) obeys
 n
q
P
√
L(1+η)
µ̂
kp − si kλ,K
for any K = 0, 1, . . . .
|r|λ,K ≤ 2r0 + (λ n)−1
+
β
F
µ̄η
µ̄
i=1

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

13

Proof. By assumptions, for each i ∈ {1, 2, . . . , n} and k = 0, 1, . . ., we have
gi (p∗ ) ≥ gi (ski ) + h∇gi (ski ), p∗ − ski i +

(20)

µi ∗
kp − ski k2F .
2

β
kpk − p∗ k2F − βkpk − ski k2F where β > 0 is a
Through using the basic inequality kski − p∗ k2F ≥ β+1
tunable parameter, it follows from (20) that

gi (p∗ ) ≥

gi (ski) + h∇gi (ski ), pk − ski i + h∇gi (ski), p∗ − pk i

+ µ2i

β
k
∗ 2
k
k 2
β+1 kp − p kF − βkp − si kF

and therefore
h∇gi (ski ), p∗ − pk i ≤

(21)

gi (p∗ ) − gi (ski ) − h∇gi (ski ), pk − ski i
µi β
− 2(β+1)
kpk − p∗ k2F + µ2i β kski − pk k2F .

Averaging (21) over i through 1 to n gives

≤

(22)

1 Pn
k
∗
k
i=1 h∇gi (si ),
n

p − p i
1 Pn
∗
k
g(p ) − n i=1 gi (si ) + h∇gi (ski ), pk − ski i − µ2i β kski − pk k2F
µ̄β
− 2(β+1)
kpk − p∗ k2F .

On the other hand, we also have that for any vector ∆,

gi ski + ∆ + pk − ski
gi (ski ) + h∇gi (ski ), ∆ + pk − ski i + L2i k∆ + pk − ski k2F
gi (ski ) + h∇gi (ski ), ∆i + h∇gi (ski ), pk − ski i
k∆k2F + Li (1+η)
kpk − ski k2F
+ Li (1+η)
2
2η

gi (pk + ∆) =
≤
≤

where η > 0 is some tunable parameter, and therefore
−h∇gi (ski ), ∆i

(23)

≤ −gi (pk + ∆) + gi (ski ) + h∇gi (ski ), pk − ski i
kpk − ski k2F + Li (1+η)
k∆k2F .
+ Li (1+η)
2η
2

Averaging (23) over i through 1 to n gives
(24)

−h n1

Pn

k
i=1 ∇gi (si ), ∆i

≤

L̄(1+η)
2
−g(pk + ∆)

 + 2 k∆kF
P
n
k
k 2
+ n1 i=1 gi (ski ) + h∇fi (ski ), pk − ski i + Li (1+η)
kp
−
s
k
i F .
2η

Having (22) and (24) at hand, we are ready to show how rk+1 is related to rk . First, plugging
a = pk+1 − p∗ and b = pk − pk+1 into the basic equality kak2F = ka + bk2F − 2ha, bi − kbk2F yields
(rk+1 )2

=
=

(25)
=

(rk )2 − 2hpk+1 − p∗ , pk − pk+1 i − kpk+1 − pk k2F
(substituting (18))
n
P
(rk )2 − 2hpk+1 − p∗ , θ n1
∇gi (ski )i − kpk+1 − pk k2F
i=1
n
P
1
k
∇gi (si ), p∗ − pk i
(r ) + 2θh n
i=1
n
P
∇gi (ski ), pk+1 − pk i − kpk+1 − pk k2F .
−2θh n1
i=1
k 2

14

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

Next, in (25), we substitute (22) for the second term, and we substitute (24) with ∆ = pk+1 − pk for
the third term. Thus, we obtain that
Pn
k
k
k
(rk+1 )2 ≤ (rk )2 + 2θ g(p∗ )− n1 i=1 gi (ski ) + h∇g
 i (si ), p − si i
µ̄β
− µ2i β kski − pk k2F − 2(β+1)
kpk − p∗ k2F
P
n
+2θ −g(pk+1 ) + n1  i=1 gi (ski ) + h∇gi (ski ),pk − ski i
Li (1+η)
+ 2η kpk − ski k2F + L̄(1+η)
kpk+1 − pk k2F − kpk+1 − pk k2F
2
θ
µ̄β
(26)
= (rk )2 + 2θ(g(p∗ ) − g(pk+1 )) − β+1
kpk − p∗ k2F


P
n
θL
(1+η)
i
+ θµi β kpk − ski k2F − (1 − θL̄(1 + η))kpk+1 − pk k2F
+ n1 i=1
η


θ µ̄β
≤
1 − β+1
(rk )2 − 2θ(g(pk+1 ) − g(p∗ ))
 P

n
+ θµ̂β n1 i=1 kpk − ski k2F − (1 − θL̄(1 + η))kpk+1 − pk k2F .
+ θL(1+η)
η
Define ǫk = n1
we have
(27)

Pn

i=1 kp

k

1
such that 1 − θL̄(1 + η) in (26) is nonnegative,
− ski k2F . By choosing θ ≤ (1+η)
L̄





θ µ̄β
(rk+1 )2 ≤ 1 − β+1
(rk )2 − 2θ(g(pk+1 ) − g(p∗ )) + θL(1+η)
+
θ
µ̂β
ǫk .
η

Pn
Let us look into the last two terms of (27). Noticing that µ̄ = (1/n) i=1 µi is a strong convexity
constant of g(p), there are two possibilities that could happen at time k. Possibility A is that


L(1 + η) µ̂
k+1 2
+ β ǫk ,
(r
) ≥
µ̄η
µ̄
while possibility B is the opposite, namely that


L(1 + η) µ̂
k+1 2
(r
) <
+ β ǫk .
µ̄η
µ̄
If possibility A occurs, we have
2θ(g(pk+1 ) − g(p∗ )) ≥ θµ̄kpk+1 − p∗ k2F = θµ̄(rk+1 )2 ≥
which together with (27) implies
(r



θL(1+η)
+ θµ̂β
η



ǫk



θµ̄β
(rk )2 .
) ≤ 1−
β+1

k+1 2

Considering both possibilities A and B, it follows that
n


 o
θ µ̄β
µ̂
(28)
(rk+1 )2 ≤ max 1 − β+1
(rk )2 , L(1+η)
+
β
ǫk .
µ̄η
µ̄

Recursively using the inequality (28) we can see that





k+1
t
L(1+η)
µ̂
θ µ̄β
θ µ̄β
0 2
k−t
k+1 2
+ µ̄ β max
(r ) ,
1 − β+1 ǫ
.
(29)
(r
) ≤ max
1 − β+1
µ̄η
t=0,...,k

Taking square root on both sides of (29) gives us
q


q
q
k+1
t √
L(1+η)
µ̂
θ µ̄β
θ µ̄β
k+1
0
k−t
(30)
r
≤
.
+ µ̄ β max
r +
1 − β+1
1 − β+1
ǫ
µ̄η
t=0,...,k

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

15





θ µ̄β
≤ 1, then from (30) we get
Choose λ that satisfies (19) so to have c , (λ)−2 1 − β+1

(λ)−(k+1) rk+1
(31)


q
n
o
√ k+1 0
√ t√
L(1+η)
µ̂
≤ ( c)
+
β
max (λ)−(k−t) ( c) ǫk−t
r + (λ)−1
µ̄η
µ̄
t=0,...,k

q
n
√ o
L(1+η)
µ̂
0
−1
+ µ̄ β max (λ)−t ǫt .
≤ r + (λ)
µ̄η
t=0,...,k

Further observing that
v
u n
n
√
u1 X
1 X k
k
ǫ =t
kpk − ski k2F ≤ √
kp − ski kF ,
n i=1
n i=1
and combining it with (31), it follows that
 n

P
√ −1 q L(1+η) µ̂
−(k+1) k+1
0
+ µ̄ β
max {(λ)−t kpt − sti kF } .
(λ)
r
≤ r + (λ n)
(32)
µ̄η
i=1 t=0,...,k

Taking maxk=0,1,...,K−1 {·} on both sides of (32) gives
 n
q
P
√
L(1+η)
µ̂
+
β
kp − si kλ,K−1
|r|λ,K ≤ 2r0 + (λ n)−1
F
µ̄η
µ̄
i=1
 n
q
P
√
L(1+η)
≤ 2r0 + (λ n)−1
+ µ̄µ̂ β
kp − si kλ,K
F .
µ̄η
i=1

3.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 8, as a key ingredient.
Lemma 9 (Algorithm 1: The last arrow x̌ → q). Let Assumptions 1, 4, and 5 hold. In
addition, assume that the stepsize α > 0 and the parameter λ are such that
s
αµ̄β
1
,
1−
≤ λ < 1 and α ≤
β+1
(1 + η)L̄
where β > 0 and η > 0 are some tunable parameters. Then, we have


√ q
√
L(1+η)
n
µ̂
(33)
kqkλ,K
≤
1
+
+
β
kx̌kλ,K
+ 2 nkx(0) − x∗ kF for any K = 0, 1, . . . .
F
F
λ
µ̄η
µ̄
Proof. First, let us consider the evolution of x(k). Noticing that
y(k + 1) − n1 (∇f (x(k + 1)))⊤ 1

= y(k) − n1 (∇f (x(k)))⊤ 1
..
.
= y(0) − n1 (∇f (x(0)))⊤ 1
= 0

holds for all k, we then have that

(34)

x(k + 1) =
=
=

x(k) − αy(k)
x(k) − α n1 (∇f (x(k)))⊤ 1
n
P
∇fi (xi (k)).
x(k) − α n1
i=1

16

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

Applying Lemma 8 to the recursion relation of x, namely (34), we obtain
 n

P
√ −1 q L(1+η) µ̂
∗ λ,K
∗
kx − x kF
+ µ̄ β
kx − xi kλ,K
≤ 2kx(0) − x kF + (λ n)
F
µ̄η
i=1 s

q
2

n
P
√
√
L(1+η)
λ,K
µ̂
(35)
+
β
n
−
x
k
≤ 2kx(0) − x∗ kF + (λ n)−1
kx
i
F
µ̄η
µ̄
i=1

q
λ,K
L(1+η)
µ̂
+ µ̄ β kx̌kF .
≤ 2kx(0) − x∗ kF + (λ)−1
µ̄η
Since

q(k)

=
=

x(k) − 1(x(k))⊤ + 1(x(k))⊤ − x∗
x̌(k) + 1(x(k) − x∗ )⊤ ,

it follows that
kqkλ,K
≤ kx̌kλ,K
+
F
F

(36)

√
nkx − x∗ kλ,K
F .

Substituting (35) into (36) yields (33).
3.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 (B, n, and δ), objective parameters (µ̄ and κ̄ = L
µ̄ ), and the algorithmic step-size (α).
Theorem 10 (Algorithm 1: Explicit geometric rate over time-varying graphs). Suppose
that Assumptions 1, 4, and 5 hold. Let




√ √ 
1
and J1 = 3κ̄B 2 1 + 4 n κ̄ .
δ = sup
σmax WB (k) − 11⊤
n
k≥B−1
i

2
, the sequence {x(k)} generated by DIGing algorithm conThen, for any step-size α ∈ 0, 1.5(1−δ)
µ̄J1
verges to the matrix x∗ = 1(x∗ )⊤ , where x∗ is the unique optimal solution of problem (1) at a global
R-linear (geometric) rate O(λk ), where the parameter λ is given by

√
2 #
q

1.5
J12 +(1−δ 2 )J1 −δJ1

αµ̄
2B

,
if α ∈ 0,
1 − 1.5 ,

µ̄J1 (J1 +1)2

#
rq
2
√
λ=

J12 +(1−δ 2 )J1 −δJ1
1.5

B
1.5(1−δ)2
αµ̄J1

.
if α ∈
, µ̄J1


1.5 + δ,
µ̄J1 (J1 +1)2
Proof. Let us collect all the relations/arrows at hand [cf. Lemmata 5, 6, 7, and 9]:

i)
ii)
iii)
iv)

kzkλ,K
≤ γ1 kqkλ,K
+ ω1 where γ1 = L 1 + λ1
F
F



B

and ω1 = 0;
B

)
λ
ky̌kλ,K
≤ γ2 kzkλ,K
+ ω2 where γ2 = (λBλ(1−λ
F
F
−δ)(1−λ) and ω2 = λB −δ
B

B

B
P

t=1
B
P

λ1−t ky̌(t − 1)kF ;

α(1−λ )
λ
kx̌kλ,K
≤ γ3 ky̌kλ,K
+ ω3 where γ3 = (λB
λ1−t kx̌(t − 1)kF ;
F
F
−δ)(1−λ) and ω3 = λB −δ
t=1
√ q
√
+ µ̄µ̂ β and ω4 = 2 nkx(0) − x∗ kF .
kqkλ,K
≤ γ4 kx̌kλ,K
+ ω4 where γ4 = 1 + λn L(1+η)
F
F
µ̄η

To apply the small gain theorem (Theorem 3) to get that kqkλF is bounded, we need γ1 γ2 γ3 γ4 < 1,
that is,


√ q
αL(λ+1)(1−λB )2
L(1+η)
n
µ̂
1+ λ
(37)
+ µ̄ β < 1,
µ̄η
(λB −δ)2 (1−λ)2

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

(38)

where β > 0

17

and η > 0,

along with other restrictions on parameters that appear in Lemmata 6, 7, and 9:
δ < λB < 1;

(39)

(40)

s

(41)

and α ≤

1−

αµ̄β
≤ λ < 1;
β+1
1
.
(1 + η)L̄

We next use relations (37)–(41) with a specific values for the parameters β, η and the stepsize α,
which yields the desired result. Specifically, let β = 2L/µ̂ and η = 1 in relation (38). By further using
0.5 ≤ λ < 1 and (1 − λB )/(1 − λ) ≤ B, (37) and (41) together yield
α≤

(42)

(λB − δ)2
√ √ .
2LB 2 1 + 4 n κ̄

On the other hand, since 1.5 ≥ 1 + 1/β, relation (40) implies that
α≥

(43)

1.5(1 − λ2 )
.
µ̄

√
Using (42) and (43), it remains to show that there exists λ ∈ ( B δ, 1) [cf. (39)] such that
#
"
1.5(1 − λ2 )
(λB − δ)2
,
√ √  6= ∅,
µ̄
2LB 2 1 + 4 n κ̄

which is equivalent to
(44)




1.5(1 − λ2 ) 1.5(λB − δ)2
6 ∅,
=
,
µ̄
µ̄J1

√ √ 
where J1 = 3κ̄B 2 1 + 4 n κ̄ . We consider a smaller interval by enlarging the left-bound of the
interval in (44), i.e., we will prove that


1.5(1 − λ2B ) 1.5(λB − δ)2
6= ∅.
,
(45)
µ̄
µ̄J1
√
B
When λ varies from δ to 1, the left-bound of the interval in (45) is monotonically decreasing from
2
1.5(1−δ 2 )
to 0, while its right-bound is monotonically increasing from 0 to 1.5(1−δ)
. In particular, the
µ̄
µ̄J1
relation (45) is valid when λ (as small as the current choice of all parameters can give) is given by
sp
J12 + (1 − δ 2 )J1 + δ
B
.
λ=
J1 + 1
Thus, for




α ∈ 0,

1.5

2 
p
J12 + (1 − δ 2 )J1 − δJ1 
,
µ̄J1 (J1 + 1)2

18

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

q
αµ̄
we can set λ = 2B 1 − 1.5
, while for


we can use λ =

rq
B

 1.5
α∈


2
p
2
J12 + (1 − δ 2 )J1 − δJ1
1.5(1 − δ) 
,
,
µ̄J1 (J1 + 1)2
µ̄J1

αµ̄J1
1.5 + δ. The rest of the statements follow from Theorem 3 and Lemma 4.

Other possible choices of β, η, α, and λ 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.
Corollary 11 (Algorithm 1: Polynomial network scalability). Under Assumptions 2, 3,
4, and 5, if the agents choose the step-size to be
(46)

α(τ ) =

3τ 2
1.5
√ −
µ̄
128B 2 n4.5 L κ̄



τ2
2
128B n4.5 κ̄1.5

2

,

where τ is the smallest nonzero positive element of the nonnegative matrices W(k) for all k [cf.
Assumption 3]. Then, the sequence {x(k)} generated by DIGing
converges to the unique optimal

solution x∗ at a global R-linear (geometric) rate of O (λ(τ ))k where
r
τ2
B
.
λ(τ ) =
1−
128B 2 n4.5 κ̄1.5
Proof. Define ϕ = 1 − λB , then requiring that the interval in (45) be nonempty is equivalent to
showing that the following inequality has a solution:
(1 − (1 − ϕ)2 )J1 ≤ (1 − δ − ϕ)2 .
Therefore, we can show that an achievable ϕ is
√
2J1 +2(1−δ)− (2J1 +2(1−δ))2 −4(1−δ)2 (J1 +1)
ϕ =
2(J1 +1)
(47)

=
≥

2J1 +2(1−δ)+
(1−δ)2
2(J1 +1) .

√

2(1−δ)2

(2J1 +2(1−δ))2 −4(1−δ)2 (J1 +1)

By Assumption 3, from Lemma 9 of reference [43], we have that δ ≤ 1 − 2nτ 2 . Substituting J1 =
√ √ 
3B 2 κ̄ 1 + 4 n κ̄ into (47) gives us
ϕ≥

8n4 (3B 2 κ̄

τ2
τ2
≥
.
√ √ 
2
128B n4.5 κ̄1.5
1 + 4 n κ̄ + 1)

q
2B
2
)
=
Thus the final rate is λ = B 1 − 128B 2τn4.5 κ̄1.5 , and a step-size to reach this rate is α = 1.5(1−λ
µ̄


2
2
2
τ
3τ
√ − 1.5
.
µ̄
128B 2 n4.5 κ̄1.5
128B 2 n4.5 L κ̄

Corollary 11 explicitly shows how the linear convergence rate of the DIGing algorithm depends
on the condition number κ̄, time-varying graph connectivity constant B, and the network size n. To
reach ε-accuracy, the iteration complexity under conditions of Corollary 11 is O τ −2 B 3 n4.5 κ̄1.5 ln 1ε
which is polynomial in the number of agents n. Beyond the more general form of it, the advantage of
Theorem 10 is that it explicitly depends on the parameter δ which measures the convergence speed

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

19

of consensus. Indeed, Corollary 11 uses the bound δ ≤ 1 − τ /(2n2 ) from [43], 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 11 [cf. (47)], where better bounds would immediately arise with
smaller values of δ.
Corollary 12 (Algorithm 1: Iteration complexity under lazy Metropolis mixing).
Suppose Assumptions 2, 4, and 5 hold. Also assume that the graphs are time-invariant, undirected
and connected (i.e., B = 1). Let each W(k) be a lazy Metropolis matrix, that is,

 1/ (2 max{di (k), dj (k)}) ,
0,
Wij (k) =
 1−P
l∈Ni (k) Wil (k),

if (j, i) ∈ E,
if (j, i) ∈
/ E and j 6= i,
if j = i.

2
2
Then, with the agents choosing the step-size α( 71
) [cf. (46) with τ = 71
and B = 1], the sequence
∗
{x(k)} generated by DIGing converges to the unique optimal solution x at a global R-linear (geometric) rate of O(λk ), where λ =1 − 161312n1 4.5 κ̄1.5 . In particular, the number of iterations needed to reach
ε-accuracy is O n4.5 κ̄1.5 ln 1ε .

We omit the proof for Corollary 12 since it is essentially the same as the proof for Corollary 11
1
in addition to which we further used δ = 1 − 71n
2 from Lemma 2.2 of reference [50].
Remark 2 (Conservatism on the scalability of consensus speed δ). In certain cases, the
bound δ ≤ 1 − τ /(2n2 ) can be conservative. In the most dramatic case, for the (fixed) complete graph
this bound tells us that 1 − δ is bounded below by 1/n3 , whereas it is not too hard to see that for the
complete graph 1 − δ 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 − δ is lower
bounded by 1/n2 , 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 ) 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 [30]).
In general, it is not possible to give a nonconservative bound on δ 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 [18]
or the monograph [29], 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 δ can be bounded accurately.
The key tool is a connection between δ and the average hitting time (Theorem 11.11 of [29]) which
in turn can be bounded in terms of graph resistance (see [9]). Using this connection, the following
scalings may be obtained; we omit the details.
• On the complete graph, δ ≤ 1 − Ω(1).
• On the line and ring graphs, we have δ ≤ 1 − Ω(1/n2 ).
• On the 2D grid and on the complete binary tree, δ ≤ 1 − Ω(1/(n log n)).
• On any regular graph, δ ≤ 1 − Ω(1/n2 ) as a consequence of the hitting time bounds of [12].
• On the star graph and on the two-star graph (defined to be two stars on n/2 nodes with a
connection between the centers), δ ≤ 1 − Ω(1/n2 ).
Remark 3 (Adapt-then-Combine strategy and acceleration). Our analytical framework
also applies to Aug-DGM of [74] when its step-size matrix D is set to α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.

20

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

Remark 4 (Parameter selection). Allowing the agents to use uncoordinated (different) stepsizes (all satisfying an upper bound) is possible; see our very recent subsequent work [44] as well as
the related papers [31, 74]. To keep the analysis in this paper concise, we do assume all agents use the
same step-size α.
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 n and on the connectivity constant B. The other
parameters used in step-size selection are (the lower bound of ) µ̄ and (the upper bound of ) L. Since
minima (and maxima) can be easily computed in a time-varying B-connected network in 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 is essentially negligible compared to the worst-case running time
of the protocol.
Remark 5. Random link activation One important class of time-varying graphs comes from
random link activation. Consider a fixed connected graph G un = (V, E); the graph G un (t) at time t
comes from having each link of G un be present independently with some fixed (constant) probability.
Fix t; it is easy to see that under random link activation the sequence G un (1), . . . , G un (t) (that is,
the graph sequence up to time t) will be B-connected with the probability 1 − 1/poly(t) if we choose
B = O(log(nt)). Applying the results of this paper, for each t can get a bound which holds with high
probability and is geometric in t/ log(nt).
It is possible to improve on this and obtain a geometric rate in t with probability one. We sketch
how this can be done next. The constant B should be chosen so that the graph obtained by taking all
the edges that appear over B steps is connected with probability at least 1/2. This allows us to choose
B to be a constant independent of n and t. We then apply the arguments of this paper to the quantities
E[||q(k)||], E[||z(k)||], E[||y(k)||], E[[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) converges to x∗ at a geometric rate with probability one .
4. 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.
4.1. Motivation. Suppose there are n agents that can communicate over a static strongly connected directed graph. Each agent i initially holds row i of x(0) ∈ Rn×p and would like to calculate
the average n1 1⊤ x(0). To do so, one possible decentralized approach is to construct a doubly stochastic matrix W and perform updates x(k + 1) = Wx(k) starting from x(0). However, in a general
directed graph, construction of a doubly stochastic matrix needs a weight balancing procedure and it
is costly [22]. This becomes even less realistic when the graph is time-varying, as the maintenance of
a doubly stochastic matrix sequence needs a real time weight balancing.
Instead, if every agent knows its out-degree, it is possible for the agents to construct a column
stochastic matrix C and perform the following recursions with initialization u(0) = x(0) and v(0) = 1
to achieve the average (push-sum protocol [27]):
u-update:
v-update:
x-update:

u(k + 1) = Cu(k);
v(k + 1) = Cv(k); V(k + 1) = diag {v(k + 1)} ;
x(k + 1) = (V(k + 1))−1 u(k + 1).

Intuitively, noticing that u(k + 1) = u(k) (u(k) is sum preserving), rows of u(k) are heading towards
scaled averages with uneven scaling ratios across the vertices caused by the non-double stochasticity
of C (the ratios are actually the elements of a right eigenvector of C corresponding to the eigenvalue

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

21

1), while V(k) is recording the ratios. By applying the recorded ratio inverse (V(k))−1 on u(k), the
algorithm recovers the unscaled average of the rows in x(k).
4.2. The Push-DIGing Algorithm. Next we formally state Push-DIGing in Algorithm 2.
Algorithm 2: Push-DIGing
Choose step-size α > 0 and pick any x(0) = u(0) ∈ Rn×p ;
Initialize y(0) = ∇f (x(0)), v(0) = 1 ∈ Rn , and V(0) = diag {v(0)};
for k = 0, 1, . . . do
u(k + 1) = C(k) (u(k) − αy(k));
v(k + 1) = C(k)v(k); V(k + 1) = diag {v(k + 1)};
x(k + 1) = (V(k + 1))−1 u(k + 1);
y(k + 1) = C(k)y(k) + ∇f (x(k + 1)) − ∇f (x(k));
end for
Looking at the individual agents, the initialization of Push-DIGing uses an arbitrary xi (0) =
ui (0) ∈ Rp , and sets yi (0) = ∇fi (xi (0)) and vi (0) = 1 for i = 1, . . . , n. Then, at each iteration k,
every agent i performs updates, as follows:
P
ui (k + 1) = Cii (k)(ui (k) − αyi (k)) + j∈N in (k) Cij (k)(uj (k) − αyj (k)),
i
P
vi (k + 1) = Cii (k)vi (k) + j∈N in (k) Cij (k)vj (k),
i
xi (k + 1) = ui (k + 1)/vi (kP
+ 1),
yi (k + 1) = Cii (k)yi (k) + j∈N in (k) Cij (k)yj (k) + ∇fi (xi (k + 1)) − ∇fi (xi (k)),
i

where Niin (k) is the set of agents that can send information to agent i (in-neighbors of agent i) at
time k, while Niout (k) is the set of agents that can receive the information from agent i (in-neighbors
of agent i) at time k. (Formal definition of the sets Niin (k) and Niout (k) is deferred to Section 5).

At every iteration k, each agent i sends its ui (k) − αyi (k), yi (k), and vi (k) all scaled by Cij (k)
to each of its out-neighbors Niout (k), and receives the corresponding messages from its in-neighbors
Niin (k). Then, each agent i updates its own ui (k + 1) by summing its own Cii (k) (ui (k) − αyi (k))
and the received Cij (k) (uj (k) − αyj (k)) from its in-neighbors Niin (k); a similar strategy applies to
the update of vi (k + 1); then xi (k + 1) is given by scaling ui (k + 1) with (vi (k + 1))−1 ; finally each
agent i updates its own yi (k + 1) by summing its own Cii yi (k) and the received Cij (k)yj (k) from
its in-neighbors Niin (k), and accumulating its current local gradient ∇fi (xi (k + 1)) and subtracting
its previous local gradient ∇fi (xi (k)) (in order to filter in only the new information contained in the
most recent gradient). Unlike DIGing for undirected graphs in which each agent scales the received
variables (xi (k) and yi (k)) and then sums them up, in Push-DIGing, the variables (ui (k) − αyi (k),
yi (k), and vi (k)) are scaled before being sent out. This is due to the fact that, over directed graphs,
usually a scaling weight Cij can only be conveniently determined by the out-degree information of
agent j which is not available to agent i.
5. Convergence Analysis for Push-DIGing. In this section we conduct the convergence
analysis for Push-DIGing over time-varying directed graphs. Consider a time-varying graph sequence
{G dir (0), G dir (1), . . .}. Every graph instance G dir (k) consists of a static set of agents V = {1, 2, . . . , n}
−
→
and a set A(k) of time-varying arcs. An arc (j, i) ∈ A(k) indicates that agent j can send information
to agent i nat time (iteration)
out-neighbors
of agent i at time k are defined as
o k. The set of nin- and
o
−
→
−
→
in
out
Ni (k) = j (j, i) ∈ A(k) and Ni (k) = j (i, j) ∈ A(k) , respectively.
We make the following two assumptions for the setup of time-varying directed graphs.
Assumption 6 (B̃⊖ -strongly connected graph sequence). The time-varying directed graph
sequence G dir (k) is B̃⊖ -strongly connected. Specifically, there exists an integer B̃⊖ > 0 such that for
any t = 0, 1, . . ., the directed graph


B̃⊖ −1
 (t+1)[

dir
(t
B̃
)
,
V,
GB̃
A(ℓ)
⊖
⊖


ℓ=tB̃⊖

22

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

is strongly connected.
Let us denote
B⊖ = 2B̃⊖ − 1.

dir
Note that Assumption 6 implies that GB
(k) is strongly connected for all k = 0, 1, . . ..
⊖
Assumption 7 (Mixing matrix sequence {C(k)}). For any k = 0, 1, . . ., the mixing matrix
C(k) = [Cij (k)] ∈ Rn×n is given by

Cij (k) =

1
−
→
if (j, i) ∈ A(k), and otherwise Cij (k) = 0,
(k)
+
1
dout
j

out
where dout
j (k) = |Nj (k)| is the out-degree of agent j at time k.
Assumption 6 has been used in distributed optimization over time-varying directed graphs [42].
Similar to the case of undirected graphs, we may have other options for the weights Cij (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 13 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

(48)

v(k + 1) = C(k)v(k), V(k + 1) = diag {v(k + 1)} ,
e
x(k + 1) = R(k)
(x(k) − αh(k)) ,
e
h(k + 1) = R(k)h(k)
+ (V(k + 1))−1 (∇f (x(k + 1)) − ∇f (x(k))) ,

e
where R(k)
= (V(k + 1))−1 C(k)(V(k)) and h(k) = (V(k))−1 y(k). We note here that, under Assumptions 6 and 7, it can be seen that each matrix V(k) is invertible, and that
(49)

kV−1 k1max , sup k(V(k))−1 kmax ≤ nnB⊖ ,
k≥0

where B⊖ is the graph connectivity constant defined in Assumption 6. The preceding relation follows
from Corollary 2(b) of [42] (here we borrow the notation for the special norm defined in (6)). Also,
e
e
we note that R(k)
is actually a row stochastic matrix, i.e., every row of R(k)
sums to 1 (see Lemma
4 of [41]).
In what follows, we will use following notation
Cb (k) , C(k)C(k − 1) · · · C(k + 1 − b)
for any k = 0, 1, . . . and b = 0, 1, . . . with the convention that Cb (k) = I for any needed k < 0 and
e
C0 (k) = I for any k. The same notation rule applies to R(k).
In the sequel, we will give an upper
⊤
e
bound on a norm of (I − 11 /n)RB (k), as provided in the next lemma. This lemma comes from the
properties of push-sum protocol and can be obtained from references [46, 48].
Lemma 13 (B-step consensus contraction). Under Assumptions 6 and 7, let B be an integer
satisfying B ≥ B⊖ and such that
B−1

δ , Q1 (1 − τ̃ nB⊖ ) nB⊖ < 1,
where Q1 and τ̃ are given by
(50)

Q1 , 2n

1 + τ̃ −nB⊖
1 − τ̃ nB⊖

and

τ̃ ,

1
.
n2+nB⊖

e B (k)b, i.e.,
Then, for any k = B − 1, B, . . . and any matrix b with appropriate dimensions, if a = R
−1
a = (V(k + 1)) CB (k)V(k + 1 − B)b, we have kakL ≤ δkbkL .

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

23

e
= (V(k + 1))−1 C(k)V(k) is a row
Proof. It has been shown in Lemma 4 of [42] that R(k)
stochastic matrix with entries being R̃ij (k) = Cij (k)vj (k)/vi (k + 1), where vi (k) denotes the i-th
e B (k))⊤ /n which is a stochastic vector. Then it
entry of the vector v(k). Denote d̃(k, B) , (1⊤ R
follows that
e B (k)bkF
kakL = k(I − n1 11⊤ )R
e B (k) − 1(d̃(k, B))⊤ )(I − 1 11⊤ )bkF
= k(R
n
e B (k) − 1(d̃(k, B))⊤ k2 kbkL .
≤ kR

e B (k) − 1(d̃(k, B))⊤ k2 . The first step
Now we can focus on establishment of an upper bound for kR
is to analyse the constituent of R̃ij (k) = Cij (k)vj (k)/vi (k + 1). From Corollary 2(b) of reference [41],
−
→
1
we know that vj (k) ≥ nnB
. Also, obviously we have 1/vj (k) ≥ 1/n, and for any (j, i) ∈ A(k),
⊖
Cij (k) ≥ 1/n. Based on those bounds on Cij (k), vj (k), and 1/vi (k + 1), we can have the following
equivalent arc-utilization lower bound
R̃ij (k) ≥ τ̃ ,

1
n2+nB⊖

−
→
if (j, i) ∈ A(k).

Thus, we have
(51)

e B (k) − 1(d̃(k, B))⊤ k2
kR

≤

≤

e B (k) − 1(d̃(k, B))⊤ kmax
nkR
 B−1
−nB⊖
1 − τ̃ nB⊖ nB⊖ ,
2n 1+τ̃
1−τ̃ nB⊖

where the second inequality comes from Lemma 5 of reference [48]. By (9) and (51), we get kakL ≤
 B−1
−nB⊖
Q1 1 − τ̃ nB⊖ nB⊖ kbkL with Q1 , 2n 1+τ̃
.
1−τ̃ nB⊖

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:
(52)

Algorithm 2: q → z → ȟ → x̌ → q.

Remark 6 (Consensus in the parameters). In consensus-based algorithms for optimization over
directed graphs, it is difficult to construct monotonically decreasing Lyapunov functions for convergence
analysis [51] 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 [66] 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.
5.1. The Establishment of Each Arrow. Noticing that that Lemma 5 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.
Lemma 14 (Algorithm 2: The ﬁrst arrow q → z). Under Assumption 4, we have that for
all K = 0, 1, . . . and any λ ∈ (0, 1),


1
kzkλ,K
≤
L
1
+
kqkλ,K
F
F .
λ
The next two lemmata are provided by doing almost identical arguments as those for Lemmata 6
and 7: indeed, by noticing the similarity between the equivalent recursion of Push-DIGing (48) and

24

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

the recursion of DIGing (see Algorithm 1), similar bounds to those in Lemmata 6 and 7 should be
obtained by an application of Lemma 13, which shows how multiplication by a row stochastic matrix
e
R(k)
shrinks the distance to the consensus subspace.
Lemma 15 (Algorithm 2: The second arrow z → ȟ). Let Assumptions 6 and 7 hold, and
let λ be such that δ < λB < 1, where B is the constant provided in Lemma 13. Then, we have
(53)

kȟkλ,K
≤
F

B
P
Q1 kV−1 k1max λ(1−λB )
λ,K
λB
λ1−t kȟ(t − 1)kF
kzk
+
B
B
F
(λ −δ)(1−λ)
λ −δ
t=1

for all K = 0, 1, . . . ,

where Q1 and kV−1 k1max are the constants defined by (50) and (49), respectively.
Proof. The equivalent recursion of Push-DIGing involving h and z is

e
h(k + 1) = R(k)h(k)
+ (V(k + 1))−1 z(k + 1).

(54)

From (54), using Lemma 13, for all k ≥ B − 1, it follows that
kȟ(k + 1)kF
(55)

=
≤
≤
≤

kh(k + 1)kL
e B (k)h(k + 1 − B)
R

e B−1 (k)(V(k + 2 − B))−1 z(k + 2 − B)
+ R
L
e 1 (k)(V(k)) z(k)kL + kR
e 0 (k)(V(k + 1))−1 z(k + 1)kL
+ · · · + kR
B
P
δkȟ(k + 1 − B)kF + Q1
k(V(k + 2 − t))−1 z(k + 2 − t)kF
L
−1

t=1

δkȟ(k + 1 − B)kF + Q1 kV−1 k1max

B
P

t=1

kz(k + 2 − t)kF ,

where Q1 and kV−1 k1max 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 6 starting from (15), we can obtain (53).
Lemma 16 (Algorithm 2: The third arrow ȟ → x̌). Let Assumptions 6 and 7 hold, and let
λ be such that δ < λB < 1, where B is the constant provided in Lemma 13. Then, we have
(56)

kx̌kλ,K
≤
F

α
B
λ −δ



B
1 − λB−1
λB X 1−t
δ + Q1
kȟkλ,K
+
λ kx̌(t − 1)kF
F
1−λ
λB − δ t=1

for all K = 0, 1, . . ., where Q1 is the constant as introduced in Lemma 15 (see (50)).
The equivalent recursions of Push-DIGing involving x and h is
(57)

e
x(k + 1) = R(k)
(x(k) − αh(k)) .

Noticing the similarity between (57) and (54), with almost identical argument as that illustrated in
the proof of Lemma 15, we can get (56).
Similar to the proof of Lemma 9, Lemma 8 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.
Lemma 17 (Algorithm 2: The last arrow x̌ → q). Let Assumptions 4, 5, and 7 hold. Also,
assume that
s
1
αµ̄β
≤ λ < 1 and α ≤
1−
β+1
(1 + η)L̄
where β > 0 and η > 0 are some tunable parameters. Then, we have that for all K = 0, 1, . . .,


√ q
√
√
L(1+η)
λ,K
n
µ̂
+ µ̄ β kx̌kλ,K
+ 2 nkx(0) − x∗ kF .
(58)
kqkF ≤ (1 + n) 1 + λ
F
µ̄η

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

25

Proof. First, by the same argument as in the proof of Lemma 9, we have
u(k + 1) = u(k) − α n1

(59)

n
P

fi (xi (k)).

i=1

Applying Lemma 8 to the recursion relation of u, namely (59), we obtain
 n

P
√ −1 q L(1+η) µ̂
∗
+
β
n)
≤
2ku(0)
−
x
k
+
(λ
ku − x∗ kλ,K
ku − xi kλ,K
(60)
F
F
F .
µ̄η
µ̄
i=1

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

i=1

(61)

ku − xi kλ,K
=
F
≤
≤
=

n
P

i=1

k(u − x) + (x − xi )kλ,K
F

+
nku − xkλ,K
F

n
P

kx − xi kλ,K
F
i=1
√
λ,K
1 ⊤
1 ⊤
nk n x V1 − n x 1kF + nkx̌kλ,K
F
√
k(1 − v)⊤ xkλ,K
+ nkx̌kλ,K
F
F .

Thus, by (60) and (61) we have
ku − x∗ kλ,K
F

Since

q(k)

=
=

kqkλ,K
F

≤

it follows that

(62)

≤

≤

q

L(1+η)
+ µ̄µ̂ β
µ̄η



kx̌kλ,K
2kx(0) − x kF + (λ)
F

q
√ −1
L(1+η)
+ µ̄µ̂ β k(1 − v)⊤ xkλ,K
+(λ n)
F
µ̄η
∗

−1

x(k) − 1(x(k))⊤ + 1(x(k))⊤ − 1(u(k))⊤ + 1(u(k))⊤ − x∗
x̌(k) + n1 1(1 − v(k))⊤ x(k) + 1(u(k) − x∗ )⊤ ,
√
√
⊤
xkλ,K
+ nku − x∗ kλ,K
kx̌kλ,K
+( n)−1 k(1 − v)
F
F
F

q
√
√
L(1+η)
n
λ,K
µ̂
1+ λ
+ µ̄ β
kx̌kF + 2 nkx(0) − x∗ kF
µ̄η

q

√
L(1+η)
µ̂
+ ( n)−1 + λ1
+
β
k(1 − v)⊤ xkλ,K
F .
µ̄η
µ̄

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

k(1 − v(k))⊤ x(k)kF

=
≤
≤

1
⊤
⊤
k(1
√ − v(k)) (I − n 11 )x(k)kF
2
n − nkx̌(k)kF
nkx̌(k)kF .

Substituting (63) into (62) gives (58).
5.2. Linear Convergence of Push-DIGing. Next, we provide a convergence rate estimate for
the Push-DIGing algorithm.
Theorem 18 (Algorithm 2: Explicit geometric rate over time-varying directed graphs).
Suppose Assumptions 4, 5, 6, and 7 hold. Let B be a large enough integer constant such that

δ , Q1 1 −

1
n(2+nB⊖ )nB⊖

B−1
 nB

⊖

< 1.

Also define the constant J2 as follows:
J2 = 3Q1 kV−1 k1max κ̄B(δ + Q1 (B − 1))(1 +

√


√ √ 
n) 1 + 4 n κ̄ .

26

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

i

2
, the sequence {x(k)} generated by Push-DIGing converges
Then, for any step-size α ∈ 0, 1.5(1−δ)
µ̄J2

to the unique optimal solution x∗ at a global R-linear (geometric) rate O(λk ), where λ is given by

√
2 #
q

1.5
J22 +(1−δ 2 )J2 −δJ2

αµ̄
2B

,
if α ∈ 0,
1 − 1.5 ,

µ̄J2 (J2 +1)2

#
rq
2
√
λ=

2
J22 +(1−δ 2 )J2 −δJ2
1.5

B
1.5(1−δ)
αµ̄J
2

.
if α ∈
, µ̄J2


1.5 + δ,
µ̄J2 (J2 +1)2

Proof. The proof is similar to that of Theorem 10. Specifically, we collect all the gains as follows:

i) γ1 = L 1 + λ1 ;


B−1
;
iii) γ3 = λBα−δ δ + Q1 1−λ
1−λ

Q1 kV−1 k1max λ(1−λB )
;
(λB −δ)(1−λ)


√ q
√
µ̂
iv) γ4 = (1 + n) 1 + λn L(1+η)
+
β
.
µ̄η
µ̄

ii) γ2 =

To apply the small gain theorem (Theorem 3), we need γ1 γ2 γ3 γ4 < 1, that is,





√ q
√
Q1 kV−1 k1max αL(1+λ) 1−λB
L(1+η)
n
µ̂
1−λB−1
1
+
δ
+
Q
(1
+
+
β
< 1.
n)
1 1−λ
1−λ
λ
µ̄η
µ̄
(λB −δ)2
The other restrictions on α, β, η and λ are the same as those in (38), (39), (40), and (41). Choosing
β = 2L/µ̂ and η = 1, and further using 0.5 ≤ λ < 1 and (1 − λB−1 )/(1 − λ) ≤ B − 1, relations (12)
and (41) together imply that
(64)

α≤

(λB − δ)2
√
√ √ .
2Q1 kV−1 k1max LB(δ + Q1 (B − 1))(1 + n) 1 + 4 n κ̄

Next,
√ similar to the proof of Lemma 4 [cf. (45)], it remains to show that there exists some
λ ∈ ( B δ, 1) such that


1.5(1 − λ2B ) 1.5(λB − δ)2
(65)
6= ∅,
,
µ̄
µ̄J2
√ √ 
√
where J2 = 3Q1 kV−1 k1max κ̄B(δ + Q1 (B − 1))(1 + n) 1 + 4 n κ̄ . Noticing the similarity between
(65) and (45), by the same argument as in the proof of Theorem 10 (starting from (45)), we can obtain
the statement of the theorem.
6. Numerical Experiments. Consider a decentralized estimation problem: each agent i ∈
{1, . . . , n} has its own observation yi given by yi = Mi x + ei , where yi ∈ Rmi and Mi ∈ Rmi ×p
are known data, x ∈ Rp is unknown, and ei ∈ Rmi is some noise. The goal is to estimate x.
In this experiment we use the Huber loss, which is known to be robust to outliers, and it allows
us to observe both sublinear
and linear convergence.
The corresponding optimization problem is:
o
Pn nPmi
H
(M
x
−
y
)
,
where
Mi,j and yi,j are the j-th row of matrix Mi
minx∈Rp f (x) = n1 i=1
ξ
i,j
i,j
j=1
and vector yi , respectively. The Huber loss function Hξ is defined by

1 2
if |a| ≤ ξ (ℓ1 zone),
2a ,
Hξ (a) =
ξ(|a| − 12 ξ), otherwise (ℓ22 zone).
In all experiments, we set n = 12, mi = 1 for all i, and p = 3. The data Mi , as well as the noise ei ,
∀i, are firstly generated following the standard normal distribution and then re-normalized so that
Li = 1 for all i. The parameter ξ is set to 2. The algorithm starts from xi (0) = 0 for all i. We
scale ei so that xi (0) is located in the ℓ1 zone for all agents i. The optimal solution x∗ is randomly

27

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS
10 0

4

11

2

7

10 -2

9

12

10
1
6

Residual

3
5

8

4

11

10 -4

10 -6

2

7

Gradient-Push, α(k) = 4.3/k 1/2
DIGing, α = 0.39
DIGing-ATC, α = 0.47
Push-DIGing, α = 0.26

10 -8

9

12

10
1
6

10 -10

3
5

0

500

1000

1500

2000

2500

k

8

10 0

10 0

10 -2

10 -2

10 -4

10 -4

Residual

Residual

Fig. 1. The plots are showing the underlying directed and undirected graphs for experiments. The plot to the right
kx(k)−x∗ k
shows the residuals kx(0)−x∗ k F for the time-invariant directed graph illustrated on the top-left. Step-sizes have been
F
hand-optimized to give faster convergence and more accurate solution for all algorithms.

10 -6

Gradient-Push, α(k) = 10/k 1/2
DIGing, α = 0.37
DIGing-ATC, α = 0.89
Push-DIGing, α = 1.2

10 -8

10 -10

0

500

1000

1500

2000

2500

10 -6

Gradient-Push, α(k) = 2.6/k 1/2
Push-DIGing, α = 0.12
Gradient-Push, α(k) = 3.1/k 1/2
Push-DIGing, α = 0.14

10 -8

10 -10

0

1000

2000

3000

k

4000

5000

6000

k
∗

kF
Fig. 2. The plots to the left and to the right are showing the residuals kx(k)−x
for a time-varying undirected
kx(0)−x∗ kF
graph sequence and a time-varying directed graph sequence, respectively. Step-sizes have been hand-optimized to give
faster convergence and more accurate solution for all algorithms.

but artificially selected so that x∗ is in the ℓ22 zone and kx∗ − xi (0)k = 300, ∀i. Note that while the
objective is restricted strongly convex, not all individual functions are strongly convex.
We conduct three experiments on (i) Time-invariant directed graphs; (ii) Time-varying undirected
dir
graphs; and (iii) Time-varying directed graphs. For (i), the underlying directed graph GTI
= {V, A}
dir
is illustrated in the top left of Fig. 1. GTI is randomly generated with 24 arcs and 12 vertices, and
guaranteed to be strongly connected. The mixing matrix W is chosen as the fastest (asymmetric)
un
linear consensus matrix [71]. For (ii), the underlying undirected graph GTI
= {V, E} is generated
un
dir
has 23
by simply taking out the direction of every arc of GTI (see the bottom left of Fig. 1). GTI
edges in total. At iteration k, we generate the graph instance of the time-varying undirected graph
un
GTV
(k) = {V, E(k)} through randomly uniformly sampling E(k) from E with 40 percent. The mixing
matrix W(k) is generated through Metropolis weights. For (iii), we generate the graph instance of the
dir
time-varying directed graph GTV
(k) = {V, A(k)} through randomly uniformly sampling A(k) from A
dir
of GTI with 80 percent. In all experiments, the mixing matrix C(k) for Push-DIGing is generated
based on the out-degree information [42].

28

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

We plot push-gradient method [41], DIGing, DIGing-ATC (Aug-DGM with step-size matrix α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 [62]. Interestingly, we find the convergence
curve of the EXTRA algorithm almost identical3 to that of the DIGing algorithm when the same
step size is chosen for both algorithms. As commented in Remark 3, 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 subfigure 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 [70, 76] 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.
7. 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, n, 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 César A. Uribe and Thinh T. Doan for their helpful discussions.
REFERENCES
[1] J. Bazerque and G. Giannakis, Distributed Spectrum Sensing for Cognitive Radio Networks by Exploiting
Sparsity, IEEE Transactions on Signal Processing, 58 (2010), pp. 1847–1862.
[2] D. Bertsekas, Distributed Asynchronous Computation of Fixed Points, Mathematical Programming, 27 (1983),
pp. 107–120.
[3] D. P. Bertsekas, Incremental proximal methods for large scale convex optimization, Mathematical Programming,
129 (2011), pp. 163–195.
3 Since the two curves are almost identical, EXTRA is not plotted in the figure.

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

29

[4] D. P. Bertsekas, Incremental aggregated proximal and augmented lagrangian algorithms, tech. report, Laboratory
for Information and Decision Systems Report LIDS-P-3176, MIT, 2015.
[5] D. P. Bertsekas and J. Tsitsiklis, Parallel and Distributed Computation: Numerical Methods, Athena Scientific,
Nashua, 2nd ed., 1997.
[6] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, Distributed Optimization and Statistical Learning
via the Alternating Direction Method of Multipliers, Foundations and Trends R in Machine Learning, 3 (2011),
pp. 1–122.
[7] K. Cai and H. Ishii, Average Consensus on Arbitrary Strongly Connected Digraphs with Time-Varying Topologies, IEEE Transactions on Automatic Control, 59 (2014), pp. 1066–1071.
[8] V. Cevher, S. Becker, and M. Schmidt, Convex Optimization for Big Data: Scalable, Randomized, and Parallel
Algorithms for Big Data Analytics, IEEE Signal Processing Magazine, 31 (2014), pp. 32–43.
[9] A. K. Chandra, P. Raghavan, W. L. Ruzzo, R. Smolensky, and P. Tiwari, The electrical resistance of a
graph captures its commute and cover times, Computational Complexity, 6 (1996), pp. 312–340.
[10] T. Chang, M. Hong, and X. Wang, Multi-Agent Distributed Optimization via Inexact Consensus ADMM, arXiv
preprint arXiv:1402.6065, (2014).
[11] A. Chen and A. Ozdaglar, A Fast Distributed Proximal-Gradient Method, in the 50th Annual Allerton Conference on Communication, Control, and Computing (Allerton), 2012, pp. 601–608.
[12] D. Coppersmith, U. Feige, and J. Shearer, Random walks on regular and irregular graphs, SIAM Journal on
Discrete Mathematics, 9 (1996), pp. 301–308.
[13] C. Desoer and M. Vidyasagar, Feedback Systems: Input-Output Properties, vol. 55, Siam, 2009.
[14] O. Devolder, F. Glineur, and Y. Nesterov, First-Order Methods with Inexact Oracle: The strongly convex
case, tech. report, UCL, 2013.
[15] P. Di Lorenzo and G. Scutari, Distributed nonconvex optimization over networks, in IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP),, 2015, pp. 229–232.
[16] P. Di Lorenzo and G. Scutari, Distributed nonconvex optimization over time-varying networks, in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2016, pp. 4124–4128.
[17] P. Di Lorenzo and G. Scutari, NEXT: In-Network Nonconvex Optimization, IEEE Transactions on Signal and
Information Processing over Networks, 2 (2016), pp. 120–136.
[18] P. Diaconis and D. Stroock, Geometric bounds for eigenvalues of markov chains, The Annals of Applied
Probability, (1991), pp. 36–61.
[19] J. Duchi, A. Agarwal, and M. Wainwright, Dual Averaging for Distributed Optimization: Convergence Analysis and Network Scaling, IEEE Transactions on Automatic Control, 57 (2012), pp. 592–606.
[20] P. Forero, A. Cano, and G. Giannakis, Consensus-Based Distributed Support Vector Machines, Journal of
Machine Learning Research, 59 (2010), pp. 1663–1707.
[21] L. Gan, U. Topcu, and S. Low, Optimal Decentralized Protocol for Electric Vehicle Charging, IEEE Transactions
on Power Systems, 28 (2013), pp. 940–951.
[22] B. Gharesifard and J. Cortes, Distributed Strategies for Generating Weight-Balanced and Doubly Stochastic
Digraphs, European Journal of Control, 18 (2012), pp. 539–557.
[23] M. Gurbuzbalaban, A. Ozdaglar, and P. Parrilo, On the convergence rate of incremental aggregated gradient
algorithms. available on arxiv at http://arxiv.org/abs/1506.02081, 2015.
[24] M. Hong, Decomposing Linearly Constrained Nonconvex Problems by A Proximal Primal Dual Approach: Algorithms, Convergence, and Applications, arXiv preprint arXiv:1604.00543, (2016).
[25] M. Hong and T. Chang, Stochastic Proximal Gradient Consensus over Random Networks, arXiv preprint
arXiv:1511.08905, (2015).
[26] D. Jakovetic, J. Xavier, and J. Moura, Fast Distributed Gradient Methods, IEEE Transactions on Automatic
Control, 59 (2014), pp. 1131–1146.
[27] D. Kempe, A. Dobra, and J. Gehrke, Gossip-Based Computation of Aggregate Information, in Proceedings of
the 44th Annual IEEE Symposium on Foundations of Computer Science, 2003, pp. 482–491.
[28] J. Koshal, A. Nedić, and U. V. Shanbhag, Distributed algorithms for aggregative games on graphs, Operations
Research, 64 (2016), pp. 680–704.
[29] D. A. Levin, Y. Peres, and E. L. Wilmer, Markov chains and mixing times, American Mathematical Soc.,
2009.
[30] T. Lindvall, Lectures on the coupling method, Dover Publications, Inc., 2002.
[31] Q. Lü and H. Li, Geometrical convergence rate for distributed optimization with time-varying directed graphs
and uncoordinated step-sizes, arXiv preprint arXiv:1611.00990, (2016).
[32] G. Mateos, J. Bazerque, and G. Giannakis, Distributed Sparse Linear Regression, IEEE Transactions on Signal
Processing, 58 (2010), pp. 5262–5276.
[33] A. Mokhtari and A. Ribeiro, Dsa: Decentralized double stochastic averaging gradient algorithm, Journal of
Machine Learning Research, 17 (2016), pp. 1–35.
[34] A. Mokhtari, W. Shi, Q. Ling, and A. Ribeiro, DQM: Decentralized Quadratically Approximated Alternating
Direction Method of Multipliers, arXiv preprint arXiv:1508.02073, (2015).
[35] A. Mokhtari, W. Shi, Q. Ling, and A. Ribeiro, A Decentralized Second-Order Method with Exact Linear
Convergence Rate for Consensus Optimization, arXiv preprint arXiv:1602.00596, (2016).
[36] A. Nedić, Asynchronous Broadcast-Based Convex Optimization over a Network, IEEE Transactions on Automatic
Control, 56 (2011), pp. 1337–1351.

30

ANGELIA NEDIĆ, ALEX OLSHEVSKY, AND WEI SHI

[37] A. Nedić and D. Bertsekas, Incremental Subgradient Methods for Nondifferentiable Optimization, SIAM Journal
on Optimization, 12 (2001), pp. 109–138.
[38] A. Nedić and D. P. Bertsekas, Convergence Rate of Incremental Subgradient Algorithms, Kluwer Academic
Publishers, 2000, pp. 263–304.
[39] A. Nedić, D. P. Bertsekas, and V. Borkar, Distributed asynchronous incremental subgradient methods, in Proceedings of the March 2000 Haifa Workshop “Inherently Parallel Algorithms in Feasibility and Optimization
and Their Applications”, Elsevier, Amsterdam, 2001.
[40] A. Nedić and A. Olshevsky, Distributed Optimization over Time-Varying Directed Graphs, in The 52nd IEEE
Annual Conference on Decision and Control, 2013, pp. 6855–6860.
[41] A. Nedić and A. Olshevsky, Stochastic Gradient-Push for Strongly Convex Functions on Time-Varying Directed
Graphs, arXiv preprint arXiv:1406.2075, (2014).
[42] A. Nedić and A. Olshevsky, Distributed Optimization over Time-Varying Directed Graphs, IEEE Transactions
on Automatic Control, 60 (2015), pp. 601–615.
[43] A. Nedić, A. Olshevsky, A. Ozdaglar, and J. Tsitsiklis, On Distributed Averaging Algorithms and Quantization Effects, IEEE Transactions on Automatic Control, 54 (2009), pp. 2506–2517.
[44] A. Nedić, A. Olshevsky, W. Shi, and C. A. Uribe, Geometrically convergent distributed optimization with
uncoordinated step-sizes, arXiv preprint arXiv:1609.05877, (2016).
[45] A. Nedić, A. Olshevsky, and C. Uribe, Fast Convergence Rates for Distributed Non-Bayesian Learning, arXiv
preprint arXiv:1508.05161, (2015).
[46] A. Nedić and A. Ozdaglar, Cooperative Distributed Multi-agent Optimization, Cambridge University Press,
2008, ch. 1, Convex Optimization in Signal Processing and Communications, pp. 3–49.
[47] A. Nedić and A. Ozdaglar, Distributed Subgradient Methods for Multi-agent Optimization, IEEE Transactions
on Automatic Control, 54 (2009), pp. 48–61.
[48] A. Nedić and A. Ozdaglar, Convergence rate for consensus with delays, Journal of Global Optimization, 47
(2010), pp. 437–456.
[49] A. Olshevsky, Efficient Information Aggregation Strategies for Distributed Control and Signal Processing, PhD
thesis, Massachusetts Institute of Technology, 2010.
[50] A. Olshevsky, Linear Time Average Consensus on Fixed Graphs and Implications for Decentralized Optimization
and Multi-Agent Control, arXiv preprint arXiv:1411.4186v6, (2016).
[51] A. Olshevsky and J. Tsitsiklis, On the Nonexistence of Quadratic Lyapunov Functions for Consensus Algorithms, IEEE Transactions on Automatic Control, 53 (2008), pp. 2642–2645.
[52] G. Qu and N. Li, Harnessing smoothness to accelerate distributed optimization, arXiv preprint arXiv:1605.07112,
(2016).
[53] M. Rabbat and R. Nowak, Distributed Optimization in Sensor Networks, in Proceedings of the 3rd international
symposium on Information processing in sensor networks, ACM, 2004, pp. 20–27.
[54] S. Ram, A. Nedić, and V. Veeravalli, Incremental Stochastic Subgradient Algorithms for Convex Optimization,
SIAM Journal on Optimization, 20 (2009), pp. 691–717.
[55] S. Ram, V. Veeravalli, and A. Nedić, Distributed Non-Autonomous Power Control through Distributed Convex
Optimization, in INFOCOM, 2009, pp. 3001–3005.
[56] S. S. Ram, A. Nedić, and V. Veeravalli, Distributed Stochastic Subgradient Projection Algorithms for Convex
Optimization, Journal of Optimization Theory and Applications, 147 (2010), pp. 516–545.
[57] S. S. Ram, A. Nedić, and V. V. Veeravalli, A new class of distributed optimization algorithms: Application to
regression of distributed data, Optimization Methods and Software, 27 (2012), pp. 71–88.
[58] W. Ren, Consensus Based Formation Control Strategies for Multi-Vehicle Systems, in Proceedings of the American Control Conference, 2006, pp. 4237–4242.
[59] A. Sayed, Diffusion Adaptation over Networks, Academic Press Library in Signal Processing, 3 (2013), pp. 323–
454.

[60] W. Shi, Q. Ling, G. Wu, and W. Yin, A Note: the Non-Ergodic o k1 Rate of PG-EXTRA. http://
home.ustc.edu.cn/ ∼shiwei00/papers/PG EXTRA appendix small o 1 k v0.1.pdf.
[61] W. Shi, Q. Ling, G. Wu, and W. Yin, A Proximal Gradient Algorithm for Decentralized Composite Optimization,
IEEE Transactions on Signal Processing, 63 (2015), pp. 6013–6023.
[62] W. Shi, Q. Ling, G. Wu, and W. Yin, EXTRA: An Exact First-Order Algorithm for Decentralized Consensus
Optimization, SIAM Journal on Optimization, 25 (2015), pp. 944–966.
[63] W. Shi, Q. Ling, K. Yuan, G. Wu, and W. Yin, On the Linear Convergence of the ADMM in Decentralized
Consensus Optimization, IEEE Transactions on Signal Processing, 62 (2014), pp. 1750–1761.
[64] Y. Sun, G. Scutari, and D. Palomar, Distributed nonconvex multiagent optimization over time-varying networks, arXiv preprint arXiv:1607.00249, (2016).
[65] H. Terelius, U. Topcu, and R. Murray, Decentralized Multi-agent Optimization via Dual Decomposition, in
IFAC, 2011.
[66] B. Touri, Product of Random Stochastic Matrices and Distributed Averaging, Springer Science & Business Media,
2012.
[67] J. Tsitsiklis, D. Bertsekas, and M. Athans, Distributed Asynchronous Deterministic and Stochastic Gradient
Optimization Algorithms, IEEE Transactions on Automatic Control, 31 (1986), pp. 803–812.
[68] M. Wang and D. P. Bertsekas, Incremental constraint projection-proximal methods for nonsmooth convex
optimization. Lab. for Information and Decision Systems Report LIDS-P-2907, MIT, July 2013; to appear in

GEOMETRIC CONVERGENCE OVER TIME-VARYING GRAPHS

31

SIAM J. on Optimization, 2013.
[69] E. Wei and A. Ozdaglar, On the O(1/k) Convergence of Asynchronous Distributed Alternating Direction Method
of Multipliers, arXiv preprint arXiv:1307.8254, (2013).
[70] C. Xi and U. Khan, On the Linear Convergence of Distributed Optimization over Directed Graphs, arXiv preprint
arXiv:1510.02149, (2015).
[71] L. Xiao and S. Boyd, Fast Linear Iterations for Distributed Averaging, Systems and Control Letters, 53 (2004),
pp. 65–78.
[72] L. Xiao, S. Boyd, and S. Kim, Distributed Average Consensus with Least-mean-square Deviation, Journal of
Parallel and Distributed Computing, 67 (2007), pp. 33–46.
[73] J. Xu, Augmented Distributed Optimization for Networked Systems, PhD thesis, Nanyang Technological University, 2016.
[74] J. Xu, S. Zhu, Y. Soh, and L. Xie, Augmented Distributed Gradient Methods for Multi-Agent Optimization
Under Uncoordinated Constant Stepsizes, in Proceedings of the 54th IEEE Conference on Decision and Control
(CDC), 2015, pp. 2055–2060.
[75] K. Yuan, Q. Ling, and W. Yin, On the Convergence of Decentralized Gradient Descent, arXiv preprint
arXiv:1310.7063, (2013).
[76] J. Zeng and W. Yin, ExtraPush for Convex Smooth Decentralized Optimization over Directed Networks, arXiv
preprint arXiv:1511.02942, (2015).
[77] M. Zhu and S. Martinez, Discrete-Time Dynamic Average Consensus, Automatica, 46 (2010), pp. 322–329.
[78] M. Zhu and S. Martinez, On Distributed Convex Optimization under Inequality and Equality Constraints, IEEE
Transactions on Automatic Control, 57 (2012), pp. 151–164.

