A Unification and Generalization of Exact
Distributed First Order Methods

arXiv:1709.01317v2 [cs.IT] 25 Dec 2017

Dušan Jakovetić

Abstract—Recently, there has been significant progress in
the development of distributed first order methods. (At least)
two different types of methods, designed from very different
perspectives, have been proposed that achieve both exact and
linear convergence when a constant step size is used – a favorable
feature that was not achievable by most prior methods. In this
paper, we unify, generalize, and improve convergence speed of
these exact distributed first order methods. We first carry out
a novel unifying analysis that sheds light on how the different
existing methods compare. The analysis reveals that a major
difference between the methods is on how a past dual gradient of
an associated augmented Lagrangian dual function is weighted.
We then capitalize on the insights from the analysis to derive
a novel method – with a tuned past gradient weighting –
that improves upon the existing methods. We establish for the
proposed generalized method global R-linear convergence rate
under strongly convex costs with Lipschitz continuous gradients.
Index Terms—Distributed optimization, Consensus optimization, Exact distributed first order methods, R-linear convergence
rate.

I. I NTRODUCTION
Context and motivation. Distributed optimization methods
for solving convex optimization problems over networks,
e.g., [1]–[8], have gained a significant renewed and growing
interest over the last decade, motivated by various applications,
ranging from inference problems in sensor networks, e.g., [9],
[10], to distributed learning, e.g., [11], to distributed control
problems, e.g., [12].
Recently, there have been significant advances in the context
of distributed first order methods. A distinctive feature of
the novel methods, proposed and analyzed in [13]–[26], is
that, unlike most of prior distributed first order algorithms,
e.g., [1], [27], [5], [7], they converge to the exact solution
even when a constant (non-diminishing) step size is used.
This property allows the methods to achieve global linear
convergence rates, when the nodes’ local costs are strongly
convex and have Lipschitz continuous gradients. Among the
exact distributed first order methods, (at least) two different
types of methods have been proposed, each designed through
a very different methodology. The method in [13], dubbed
Extra, see also [25], [23], modifies the update of the standard
distributed gradient method, e.g., [1], by introducing two
different sets of weighting coefficients (weight matrices) for
any two consecutive iterations of the algorithm, as opposed
to a single weight matrix with the standard method [1].
The second type of methods, [14], [19], [18], [20], [21],
replaces the nodes’ local gradient with a “tracked” value
D. Jakovetić is with the Department of Mathematics and Informatics,
Faculty of Sciences, University of Novi Sad, Novi Sad, Serbia. The research is
supported by Ministry of Education, Science and Technological Development,
Republic of Serbia, grant no. 174030. Author’s e-mail: djakovet@uns.ac.rs.

of the network-wide average of the nodes’ local gradients.
The method in [14] has been proved to be convergent under
time-varying networks, e.g., [19], node-varying step sizes,
e.g.,[18], and Nesterov acceleration [22], while such studies
have not been provided to date for the method in [13]. Both
references [13] and [14] establish global linear convergence
rates of the exact distributed first order methods that they
study. Further, reference [28] shows that Extra is equivalent
to a primal-dual gradient-like method (see, e.g., [29], [30],
[31], for primal-dual (sub)gradient methods) applied on the
augmented Lagrangian dual problem of a reformulation of the
original problem of interest. Reference [19] demonstrates that
the method therein, equivalent to the method in [14], can be
put in the form of Extra, with a specific choice of the two
Extra’s weight matrices, and it also provides a primal-dual
interpretation of the exact method therein. (A more detailed
review of works [13]–[26] is provided further ahead.)
Contributions. The main contributions of this paper are
to unify, generalize, and improve convergence speed of exact
distributed first order methods. First, with both the methods
in [13] and [14], we provide an augmented Lagrangian-based
analysis. that gives for both methods a characterization of the
(joint) evolution of primal and dual errors along iterations.
This characterization reveals the effect of the difference between the two methods on their performance and sheds light
on understanding of how the two methods compare under
various problem model scenarios. We further provide a generalized method, parameterized with an additional, easy-to-tune
network-wide weight matrix, that determines the weighting
of the term that corresponds to a past dual gradient of the
associated augmented Lagrangian function. The generalized
method subsumes the two known methods in [13] and [14]
upon setting a weighting parameter matrix to appropriate
specific values. We describe how to set the weighting parameter matrix to improve upon both existing methods. The
proposed parameter matrix tuning is based on an insight from
our analysis on how the negative effect of the primal error
on the dual error can be reduced (compensated) through the
introduced parameter matrix.
The paper carries out convergence rate analysis of the
proposed generalized method when nodes’ local costs are
strongly convex, have Lipschitz continuous gradient, and the
underlying network is static. With the proposed generalized
method, we establish a global R-linear convergence rate,
when the algorithm’s step size is appropriately set. Numerical
examples confirm the insights from our analysis on the mutual
comparison of the methods in [13] and [14], as well as the
improvements of the proposed generalized method over the
two existing ones.

2

Brief literature review. Distributed computation and optimization has been studied for a long time, e.g., [32]. More
recently, reference [1] proposes a distributed first order (subgradient) method for unconstrained problems, allowing for
possibly non-differentiable, convex nodes’ local costs and carries out its convergence and convergence rate analysis for deterministically time varying networks. A distributed projected
subgradient method for constrained problems and possibly
non-differentiable local costs has been proposed and analyzed
in [27]. References [7], [17], [13], [14], [15], [16], [18], [19],
[20], [21], [22] study unconstrained distributed optimization
under more structured local costs. In [7], distributed gradient
methods with an acceleration based on the Nesterov (centralized) gradient method [33] have been proposed and analyzed
under differentiable local costs with Lipschitz continuous and
bounded gradients. The methods in [1], [27], [7], [34] converge
to the exact solution only when a diminishing step-size is
used; when a constant step size is used they converge to a
solution neighborhood. References [35], [36] propose different
types of primal-dual methods, prove their convergence to
the exact solution for a wide class of problems assuming
diminishing step-sizes, and are not concerned with establishing
the methods’ convergence rates. The authors of [37] use
insights from control theory to propose a gradient-like algorithm for which they prove exact convergence under certain
conditions, while the paper is not concerned with analyzing the
method’s convergence rate. Reference [17] proposes several
variants of distributed AL methods that converge to the exact
solution under a constant step size and establishes their linear
convergence rates for twice continuously differentiable costs
with bounded Hessian.
References [13], [14], [15], [16], [18], [19], [20], [21], [22]
develop and/or analyze different variants of exact distributed
first order methods under various assumptions on the nodes’
local costs, algorithm step sizes, and the underlying network.
The papers [13], [14] propose two different exact distributed
first order methods and analyze their convergence rates, as
already discussed above. References [15], [16] develop exact
methods based on diffusion algorithms (see, e.g., [3]) and
through a primal-dual type method on an associated AL
function. Under twice differentiable local costs, each node’s
cost having Lipschitz continuous gradient, and at least one
node’s cost being strongly convex, the papers show linear
convergence rates for the methods therein, allowing for different step sizes across nodes and for a wider range of step
sizes and admissible weight matrices with respect to [13]. The
papers [18], [19], [20], [21] consider strongly convex local
costs with Lipschitz continuous gradients. Under this setting,
reference [18] establishes global linear convergence of the
method studied therein when nodes utilize uncoordinated step
sizes. The paper [19] establishes global linear convergence for
time-varying networks. The authors of [21] prove global linear
rates under both time varying networks and uncoordinated step
sizes, while the issue of uncoordinated step-sizes is previously
considered in [20]. Finally, the paper [22] proposes an accelerated exact distributed first order method based on the Nesterov
acceleration [33] and establishes its convergence rates for local
costs with Lipschitz continuous gradients, both in the presence

and in absence of the strong convexity assumption. Under
strongly convex local costs that have Lipschitz continuous
gradients, the authors of [26] develop optimal distributed
methods, where the optimality is in terms of the number of
oracle calls of a therein appropriately defined oracle. However,
their method is of a different type than [18], [19], [20], [21],
[22] and is different from the method proposed here. Namely,
the method in [26] requires evaluation of Fenchel conjugates of
the nodes’ local costs at each iteration, and hence it in general
has a much larger computational cost per iteration than the
methods in [18], [19], [20], [21], [22] and the method proposed
in this paper. To the best of our knowledge, communication
and computation optimality of distributed first order methods
that do not involve Fenchel conjugates has not been studied
to date.
Paper organization. The next paragraph introduces notation. Section II explains the model that we assume and reviews
existing distributed first order methods. Section III presents the
proposed algorithm and relates it with the existing methods.
Section IV establishes global R-linear convergence rate of the
proposed method, while Section V provides simulation examples. Finally, we conclude in Section VI. Certain auxiliary
proofs are provided in the Appendix.
Notation. We denote by: R the set of real numbers; Rd
the d-dimensional Euclidean real coordinate space; Aij the
entry in the i-th row and j-th column of a matrix A; A>
the transpose of a matrix A; ⊗ the Kronecker product of
matrices; I, 0, and 1, respectively, the identity matrix, the
zero matrix, and the column vector with unit entries; J the
N × N matrix J := (1/N )11> ; A  0 (A  0) means
that the symmetric matrix A is positive definite (respectively,
positive semi-definite); k·k = k·k2 the Euclidean (respectively,
spectral) norm of its vector (respectively, matrix) argument;
λi (·) the i-th largest eigenvalue; Diag (a) the diagonal matrix
with the diagonal equal to the vector a; | · | the cardinality
of a set; ∇h(w) and ∇2 h(w), respectively, the gradient and
Hessian evaluated at w of a function h : Rd → R, d ≥ 1.
Finally, for two positive sequences ηn and χn , we have:
ηn = O(χn ) if lim supn→∞ χηnn < ∞; and ηn = Ω(χn ) if
lim inf n→∞ χηnn > 0.
II. M ODEL AND PRELIMINARIES
Subsection II-A describes the optimization and network
models that we assume. Subsection II-B reviews the standard
(inexact) distributed gradient method in [1], as well as the
exact methods in [13] and [14], and it reviews the known
equivalence (see [28]) of the method in [13] and a primaldual gradient-like method.
A. Optimization and network models
We consider distributed optimization where N nodes in a
connected network solve the following problem:
minx∈Rd f (x) :=

N
X

fi (x).

(1)

i=1

Here, fi : Rd → R is a convex function known only by node i.
Throughout the paper, we impose the following assumption on
the fi ’s.

3

Assumption 1 Each function fi : Rd → R, i = 1, ..., N , is
strongly convex with strong convexity parameter µ, and it
has Lipschitz continuous gradient with Lipschitz constant L,
where L ≥ µ > 0. That is, for all i = 1, ..., N , there holds:
fi (y) ≥ fi (x) + ∇fi (x)> (y − x)
µ
+ kx − yk2 , x, y ∈ Rd
2
k∇fi (x) − ∇fi (y)k ≤ L kx − yk, x, y ∈ Rd .
For a specific result (precisely, Lemma 3 ahead), we additionally assume the following.
Assumption 2 Each function fi : Rd → R, i = 1, ..., N , is
twice continuously differentiable.
Under Assumptions 1 and 2, there holds for every x ∈ Rd
that:
µ I  ∇2 fi (x)  L I.
Further, under Assumption 1, problem (1) is solvable and has
the unique solution x? ∈ Rd .
Nodes i = 1, ..., N constitute an undirected network G =
(V, E), where V is the set of nodes and E is the set of edges.
The presence of edge {i, j} ∈ E means that the nodes i and
j can directly exchange messages through a communication
link. Further, let Ωi be the set of all neighbors of a node i
(including i).
Assumption 3 The network G = (V, E) is connected, undirected and simple (no self-loops nor multiple links).
We associate with network G a N × N symmetric, (doubly)
stochastic weight matrix W . Further, we let: Wij = Wji > 0,
if {i, j} ∈ E,
/ E, i 6= j; and
P i 6= j; Wij = Wji = 0, if {i, j} ∈
Wii = 1 − j6=i Wij > 0, for all i = 1, ..., N . Denote by λi ,
i = 1, 2, ..., N , the eigenvalues of W , ordered in a descending
order; it can be shown that they obey 1 = λ1 > λ2 ≥ ... ≥
λN > −1.
Throughout the paper, we consider several iterative distributed methods to solve (1). An algorithm’s iterations are
(k)
indexed by k = 0, 1, 2, ... Further, we denote by xi ∈ Rd the
estimate of the solution to (1) available to node i at iteration k.
To avoid notational clutter, we will keep the same notation
(k)
xi across different methods, while it is clear from context
which method
is in question. For compact
notation, we use


(k)

(k)

(k)

>

∈ RN d to denote
x(k) = (x1 )> , (x2 )> , ..., (xN )>
the vector that stacks the solution estimates by all nodes at
iteration k. We use analogous notation for certain auxiliary
(k)
sequences that the algorithms maintain (e.g., see ahead si
and s(k) with (14).) Again, to simplify notation, with all
methods to be considered we will assume equal initialization
(0)
(0)
(0)
across all nodes, i.e., we let x1 = x2 = ... = xN ; e.g.,
(0)
nodes can set xi = 0, for all i.
For future reference, we introduce the following quantities.
We let the (N d) × (N d) matrix W = W ⊗ I, where I is the
d × d identity matrix. That is, W is a N × N block matrix
with d × d blocks, such that the block at the (i, j)-th position

equals Wij I.1 Note that the (i, j)-th and (j, i)-th block of W
equal to zero if {i, j} ∈
/ E, i 6= j. In other words, we say that
W respects the sparsity pattern of the underlying graph G.
Further, recall that J = N1 11> is the ideal consensus matrix,2
f = W −J = (W −J)⊗I the
and let J = J ⊗I. Denote by W
matrix that describes how far is W from the ideal matrix J .
f = max{λ2 , −λN } =: σ ∈ [0, 1).
It can be shown that kWk
•
?
Next, let x = 1 ⊗ x , where we recall that 1 is an all-ones
vector (here of size N × 1), and x? is the (d × 1) solution to
(1). In other words, x• concatenates N repetitions of x? on
top of each other; the i-th repetition corresponds to node i in
the network. Our goal is that, with a distributed method, we
(k)
have x(k) → x• , or, equivalently, xi → x? , for all nodes
i = 1, ..., N . Further, we define function F : RN d → R, by
F (x) := F (x1 , ..., xN ) =

N
X

fi (xi ).

(2)

i=1

Note that, under Assumption 1, F is strongly convex with
strong convexity parameter µ, and it has Lipschitz continuous
gradient with Lipschitz constant L.
B. Review of distributed first order methods
We review three existing distributed first order methods
that are of relevance to our studies – the standard distributed
gradient method in [1], the Extra method in [13], and the
method in [14].
Standard distributed gradient method in [1]. The method
(k)
updates the solution estimate xi at each node i by weightaveraging node i’s solution estimate with the estimates of its
immediate neighbors j ∈ Ωi \ {i}, and then by taking a step
in the negative local gradient’s direction. This corresponds to
the following update rule at arbitrary node i:
X
(k+1)
(k)
(k)
xi
=
Wij xj − α ∇fi (xi ), k = 0, 1, ... (3)
j∈Ωi

In matrix format, the network-wide update is as follows:3
x(k+1) = W x(k) − α ∇F (x(k) ), k = 0, 1, ...

(4)

Here, constant α > 0 is the algorithm’s step size. A drawback
of algorithm (4) is that (when constant step size α is used)
it does not converge to the exact solution x• := 1 ⊗ x? ,
but only to a point in a O(α) solution neighborhood; see,
e.g., [39]. The algorithm can be made convergent to x• by
1 We will frequently work with Kronecker products of types A ⊗ B, a ⊗ B,
and a ⊗ b, where A is an N × N matrix, a – N × 1 vector, B – d × d matrix,
and b – d×1 vector. In other words, the Kronecker products throughout always
appear with the first argument of dimension either N × N or N × 1, and the
second argument either d × d or d × 1. As the capital letters denote matrices
and the lower case letters denote vectors, the dimensions of the Kronecker
products arguments will be clear.
2 The consensus algorithm, e.g., [38], computes the global average of nodes’
local quantities through a linear system iteration with a stochastic system
matrix W . When W = J, consensus converges in a single iteration, hence
we name J the ideal consensus matrix. Clearly, J is not realizable over a
generic graph as it is not sparse, but it is usually desirable to have W as
close as possible to J in an appropriate sense, given the constraints on the
network sparsity.
3 To save space, we will present all subsequent algorithms in compact forms
only. From a compact form, it is straightforward to recover local update forms
at any node i.

4

taking an appropriately set diminishing step size (e.g., square
summable, non-summable), but the convergence rate of the
resulting method is sublinear.
Extra [13]. The method modifies the update rule (1) and
achieves exact (global linear) convergence with a fixed step
size α. Extra works as follows. It uses the following update
rule:4
x(1)
x

(k+1)

= W x(0) − α ∇F (x(0) )

(5)

= 2 W x(k) − α ∇F (x(k) ) − W x(k−1)

(6)

+ α ∇F (x

(k−1)

), k = 1, 2, ...

Reference [28] demonstrates that algorithm (5)–(6) is a
primal-dual gradient-like method; see, e.g., [29], [30], [31], for
primal-dual (sub)gradient methods. Denote by L := I − W,
and consider the following constrained problem:
minimize F (x)

1 1/2
L x = 0.
α

subject to

(7)

It can be shown, e.g., [28], that L1/2 x = 0, if and only if x1 =
x2 = ... = xN , where xi ∈ Rd is the i-th consecutive block
of x = ( (x1 )> , ..., (xN )> )> . Therefore, (7) is equivalent to
(1). Introduce the augmented Lagrangian function A : RN d ×
RN d → R (with the penalty parameter equal to α) associated
with (7):
1 >
1
x Lx,
(8)
A(x, µ) = F (x) + µ> L1/2 x +
α
2α
and the corresponding dual problem:
maximizeµ∈RN d

inf A(x, µ).

x∈RN d

(9)

Consider the following primal-dual method to solve (9):5
x(k+1)
(k+1)

µ

=
=

x(k) − α ∇x A(x(k) , µ(k) )
(k)

µ

+ α ∇µ A(x

(k+1)
(0)

(k)

,µ

(10)

=

= x(k) − α




1
L x(k) + ∇F (x(k) ) + u(k) (12)
α

1
L x(k+1) , k = 0, 1, ...
(13)
α
It turns out that (12)–(13) with a proper initialization is
equivalent to (5)–(6).7
u(k+1)

= u(k) +

Lemma 1 ([28]) The sequence of iterates {x(k) } generated by
(0)
(0)
algorithm (5)–(6), with initialization xi = xj , for all i, j,
is the same as the sequence of iterates generated by (12)–(13),
(0)
with the same initialization of xi , i = 1, ..., N , and with u(k)
initialized to zero.
Under Assumptions 1 and 3 an appropriately chosen step size
α, one has that x(k) → x• , and u(k) → −∇F (x• ), at an
R-linear rate [13].
The exact method in [14]. The authors of [14], see
also [18], [19], [20], [21], [22], consider a distributed first
order method that, besides solution estimate x(k) ∈ RN d , also
maintains an auxiliary variable s(k) ∈ RN d . Here, at each node
(k)
i, quantity si ∈ Rd serves to approximate the network-wide
PN
(k)
gradient average N1 i=1 ∇fi (xi ); then, the gradient con(k)
tribution −α∇F (x ) with the standard distributed gradient
method (4) is replaced with s(k) . More precisely, the update
rule for k = 0, 1, ... is as follows:
x(k+1)

=

W x(k) − α s(k)

(k+1)

=

W s(k) + ∇F (x(k+1) ) − ∇F (x(k) ), (15)

s

(14)

with s(0) = ∇F (x(0) ). The method also achieves a global
R-linear convergence under appropriately chosen step-size α,
when Assumptions 1 and 3 are in force.

Nd

III. T HE PROPOSED METHOD
Subsection III-A presents the method that we propose and
explains how to recover the existing methods in [13], [14] by
a particular setting of a proposed method’s parameter matrix.
Subsection III-B gives further insights into the proposed
method and explains how to tune the parameter matrix. Finally,
Subsection III-C provides a primal-dual interpretation of the
proposed method and the methods in [13], [14].

µ(k) + L1/2 x(k+1) , k = 0, 1, ...

Introducing the new variable u(k) :=
4 The

x(k+1)

), k = 0, 1, ...(11)

with step size α > 0, arbitrary x ∈ R , and µ(0) = 0.
Here, ∇x and ∇µ denote the partial derivatives with respect
to x and µ, respectively. Evaluating the partial derivatives, we
arrive at the following method:

1
L x(k)
x(k+1) = x(k) − α
α

L1/2 (k)
+ ∇F (x(k) ) +
µ
α
µ(k+1)

at the following method:6

1 1/2
αL

µ(k) , one arrives

Extra method utilizes two different weight matrices in general. When
these matrices are tuned as it is suggested in [13], Extra can be represented
as in (5)–(6). The formulas here and in [13] appear slightly different due to
f = 1 (I + W ) in [13] corresponds to W here. It is
a different notation: W
2
required in [13] for the analysis of (5)–(6) that (in our notation) (2W − I)
be doubly stochastic and that W be, besides being doubly stochastic, also
positive definite. The latter assumptions are not imposed here when analyzing
the proposed method (16)–(17) (Section IV).
5 More precisely, under appropriate conditions, with (10) one has that

x(k) , µ(k) converges to a saddle point of function A(x, µ), which then
implies that x(k) solves (7) and µ(k) solves (9).

A. The proposed method and its relation with existing algorithms
We now describe the generalized exact first order method
that we propose. The method subsumes the known methods [13] and [14] upon a specific choice of the tuning
parameters, as explained below. The algorithm maintains over
iterations k the primal variable (solution estimate) x(k) ∈ RN d
6 See also [17] for similar methods with multiple primal updates per each
dual update.
7 Reference [28] shows the equivalence under a more general initialization;
we keep the one as in Lemma 1 for simplicity.

5

>
(k)
(k)
(ex )> , (eu )> , for k = 0, 1, ..., satisfies the following
recursion:
"
# 
 " (k) #
(k+1)
ex
W − α Hk
−α I
ex
(k+1) = (W − I) (H − B) W − J
(k) . (18)
k
eu
eu

and the dual variable u(k) ∈ RN d , initialized with the zero
vector. The update rule is for k = 0, 1, ... given as follows:


x(k+1) = W x(k) − α ∇F (x(k) ) + u(k)
(16)


u(k+1) = u(k) − L ∇F (x(k) ) + u(k) − B x(k) .(17)



Here, quantities W, L, and α are the same as before. We
note that one can use two different (doubly stochastic) weight
matrices W1 and W2 in (16) and (17) above; more precisely,
one can replace W with W1 in (16) and L with I − W2 in
(17). However, throughout we keep a single weight matrix W
for simplicity of the analysis and presentation. Quantity B
is a (N d) × (N d) symmetric matrix that respects the blocksparsity pattern of the underlying graph G and satisfies the
property that, for any y ∈ Rd , there exists some c ∈ R, such
that B (1 ⊗ y) = c (1 ⊗ y). Specifically, we will consider
the following choices (that clearly obey the latter conditions):
1) B = b I, where b ≥ 0 is a scalar parameter; and 2)
B = b0 W, for b0 ≥ 0. These choices are easy to implement
and incur no additional communication overhead while lead
to efficient algorithms (see also Section V). That is, with
the above choices of matrix B, the proposed method (16)–
(17) requires communicating 2 vectors of size d per node,
per iteration. We note that, for a generic choice of matrix B,
method (16)–(17) incurs one additional communication (of a
d-dimensional vector) per node, per iteration, i.e., it requires
communicating 3 vectors of size d per node, per iteration.
The next lemma, proved in the Appendix, explains how
to recover the existing exact distributed first order methods
from (16)–(17).8

Lemma 3 expresses the primal dual error e(k) as a recursion,
guided by a time varying matrix:


W − α Hk
−α I
Mk :=
,
(19)
(W − I) (Hk − B) W − J

Lemma 2 Consider algorithm (16)–(17). Then, the following
holds.
(a) Algorithm (16)–(17) with B = 0 is equivalent to
method (14), proposed in [14].
(b) Algorithm (16)–(17) with B = α1 W is equivalent to
method (5)–(6), proposed in [13].
B. Further insights into the proposed method and parameter
tuning
We give further intuition and insights into the proposed
method, the existing algorithms in [14] and [13], and we also
describe how to improve upon existing methods by appropri(k)
(k)
ately setting matrix B. Denote by ex := x(k) −x• and eu :=
u(k) + ∇F (x• ) the primal and dual errors, respectively. Also,

>
(k)
(k)
let the primal-dual error vector e(k) := (ex )> , (eu )> ,

R1
and Hk := t=0 ∇2 F x• + t (x(k) − x• ) dt. We have the
following Lemma.9
Lemma 3 Let Assumptions 1–3 hold, and consider algorithm (16)–(17). Then, the primal-dual error vector e(k) :=
8 The equivalence claimed in Lemma 2 is in the sense that the methods
generate the same sequence of iterates {x(k) } under the same initialization
x(0) , the appropriate initialization of methods’ auxiliary variables, and the
same step size α.
9 In Section IV, we show that, under Assumptions 1 and 3, e(k) converges
to zero R-linearly.

that we refer to as the error dynamics matrix.10 Equation (19)
shows the partitioning of the error dynamics matrix as a
2 × 2 block matrix with blocks of size (N d) × (N d). The
(1, 1)-th block (W − α Hk ) describes how the current primal
error affects the next primal error, the (1, 2)-th block (−α I)
describes how the current dual error affects the next primal
error, and so on. Specifically, with the methods in [14]
and [13], the blocks on the positions (1, 1), (1, 2), and (2, 2)
are of the same structure for both methods, and are of the
same structure as in (19). For clarity of explanation, assume
that the fi ’s are strongly convex quadratic functions, so that
Hk = H = ∇2 F (x) = const, for any x ∈ RN d , with
µ I  H  L I, and of course the same H appearing in
the error dynamics matrices for each of the three methods.
Then, the blocks (1, 1), (1, 2), and (2, 2) match completely
for the three methods. However, the error dynamics matrices
for different methods differ in the (2, 1)-th block. Specifically,
with [14], the (2, 1)-th block of the corresponding error
dynamics matrix equals (W − I) H. On the
 other hand, with
[13] the block equals (W − I) H − α1 W . Intuitively, one
may expect that if [Mk ]2,1 is smaller (as measured by an
appropriate matrix norm), then the algorithm’s convergence
is likely to be faster.11 This provides an intuition on the
comparison between [14] and [13]. Namely, if H is small
relative to (H − α1 W) (for instance, when α is very small
and so α1 W has a very large norm), then we expect that the
method in [14] is faster than the method in [13], and vice
versa. This intuition is confirmed in Section V by numerical
examples.
We can go one step further and seek to tune matrix B such
that [Mk ]2,1 = (W − I)(H − B) is smallest in an appropriate
sense. We consider separately the cases B = b I and B =
b0 W. For the former, a possible worst case-type approach is as
follows: choose parameter b that solves the following problem:


min b≥0 sup kH − b Ik ,
(20)
H∈H

where H is the set of all (N d) × (N d) symmetric block
diagonal matrices H with arbitrary d × d diagonal blocks that
in addition obey the following condition: µ I  H  L I. It is
easy to show (see the Appendix for details) that the solution
µ+L
to (20) is b? = µ+L
2 . Note that the choice B = 2 I in (16)–
(17) does not match either [14] or [13] method and thus
10 In Section IV, we show that, under appropriate choice of parameters α
and B, the primal-dual error e(k) converges to zero R-linearly.
11 This intuition is corroborated more formally in Theorem 4 and Remark 4.

6

represents a novel algorithm. For the latter choice B = b0 W,
the analogous problem:


min b0 ≥0 sup kH − b0 Wk
(21)
H∈H

is more challenging. A sub-optimal choice for λN > 0, as
L+µ
shown in the Appendix, is b0 = 1+λ
, where we recall that λN
N
is the smallest eigenvalue of matrix W . Extensive simulations
(see Also Section V) show that the simple choice b0 = L
works well in practice.
Note that, ideally, one would like to minimize with respect
to b the following quantity: supH: µ IHL I kMk k, i.e., one
wants to take into account the full matrix Mk . This problem
is challenging in general and is hence replaced here by
method (20) (or, similarly, (21)), i.e., by considering the
(2, 1)-th block of Mk only. Note that a smaller norm of
the (2, 1)-th block of Mk might not necessarily imply a
smaller norm of the full matrix Mk . However, extensive
numerical experiments on quadratic and logistic losses (see
also Section 5) demonstrate that tuning method (20) yields
fast algorithms while at the same time is very cheap.
C. Primal-dual interpretations
We now provide a primal-dual interpretation of the proposed
method. The interpretation builds on construction (8)–(11)
from [28]. It is worth noting that [19] provides a primal-dual
interpretation of the method therein, equivalent to the method
in [14]. The construction in [19] starts from a different problem
reformulation than (8) and utilizes a different quadratic penalty
term for the Lagrangian function.
To start, we write the methods (5)–(6) and (14)–(15) in
another equivalent form (see the Appendix as to why this
equivalence also holds.) Namely, (5)–(6) can be equivalently
represented as follows (this is essentially a re-write of (12)–
(13)) but is useful to present it here):


x(k+1) = W x(k) − α ∇F (x(k) ) + u(k)
(22)
1
L x(k+1) , k = 0, 1, ...
(23)
α
Similarly, (14) can be, for k = 0, 1, ..., equivalently represented as:


x(k+1) = W x(k) − α ∇F (x(k) ) + u(k)
(24)
u(k+1)

u(k) +

=

1
1
L x(k+1) − W L x(k) . (25)
α
α
We can re-interpret (24) as a primal-dual gradient-like method
for solving (9). Namely, due to the fact that matrices W and
L commute, it is easy to see that (24)–(25) corresponds to the
following method:
u(k+1)

x(k+1)
(k+1)

µ

u(k) +

=

=
=

x(k) − α ∇x A(x(k) , µ(k) )
(k)

µ

+ α ∇µ A(x

(k+1)

(k)

,µ

(26)
)

(27)

−α W ∇µ A(x(k) , µ(k) ).
Hence, (24) is a primal-dual gradient-like method that modifies the dual update step to also incorporate the (weighted)
previous dual gradient term. Clearly, the proposed generalized

method (16)–(17) also incorporates the (weighted) previous
dual gradient term. It is shown in the Appendix that, assuming
that matrices B and L commute (which is the case for the two
specific choices B = b I and B = b0 W considered here), we
have that (16)–(17) is equivalent to the following primal-dual
method:
x(k+1)
(k+1)

µ

= x(k) − α ∇x A(x(k) , µ(k) )
(k)

= µ

(k+1)

+ α ∇µ A(x

,µ

−α (W − α B) ∇µ A(x

(28)

(k)

(k)

)

(29)
(k)

,µ

).

IV. C ONVERGENCE RATE ANALYSIS
This Section presents the results on global R-linear convergence of the proposed method (16)–(17), and provides
the needed intermediate results and proofs. The Section is
organized as follows. First, Subsection IV-A states the result
(Theorem 4) on global R-linear convergence of the proposed
method (16)–(17) and discusses the implications of the result.
Subsection IV-B sets up the analysis and gives preliminary
Lemmas. Finally, Subsection IV-C proves a series of intermediate Lemmas, followed by the proof of Theorem 4.
A. Global R-linear convergence rate: Statement of the result
We show that, under appropriate choice of step size α, the
proposed method (16)–(17) converges to the exact solution at
a global R-linear rate. Recall quantity σ = max{λ2 , −λN } ∈
[0, 1), with λi the i-th largest eigenvalue of W .
Theorem 4 Consider algorithm (16)–(17) with B = b I,
and nlet Assumptions o1 and 3 hold. Further, let α <
1/2
µ (1−σ)2 µ
0
2
2
.
min (1−σ)
19 L2 , 192 L0 L , where L = L + b − 2 b µ
Then, the sequence of iterates x(k) generated by algo•
?
rithm (16)–(17) converges to
 x = 1⊗x R-linearly, i.e., there
(k)
•
k
holds kx − x k = O r , with r ∈ (0, 1). The convergence
factor r is at most max{1 − α2µ , 1+σ
2 } +  , where  > 0 is
arbitrarily small.
Several Remarks on Theorem 4 are now in order.
Remark 1 When B = b I is replaced with a generic symmetric
matrix B that respects the sparsity pattern of graph G and
satisfies the property that, for any y ∈ Rd , there exists some
c ∈ R, such that B (1 ⊗ y) = c (1 ⊗ y), Theorem 4 continues
to hold with L0 replaced with constant (L + kBk); see the
Appendix for the proof.
Remark 2 The maximal admissible step size for which Theorem 4 guarantees R-linear convergence can be very small
for poorly conditioned problems (L/µ large) and/or weekly
connected networks (σ close to one). However, extensive
simulations show that the proposed method converges (Rlinearly) in practice with large step sizes, e.g., α = 31L . Similar
theoretical and practical admissible values for step size α are
reported in earlier works for the methods studied therein [14],
[13].

7

Remark 3 Theorem 4 generalizes existing R-linear convergence rate results of exact distributed first order methods, e.g.,
[14], to a wider class of algorithms. Recall that for b = 0
we recover the method in [14]. Note that for b = 0 we
have L0 = L. Abstracting universal constants, the convergence
factor obtained in [14] (see Theorem 1 and Lemma 2 therein)
and the one obtained here are the same and equal


(1 − σ)2 µ2
.
(30)
1−Ω
L2
This bound is obtained here by setting the maximal
step

 size
(1−σ)2 µ
α permitted by Theorem 4, which is α = Ω
.
L2
Remark 4 Adapting the results in [18] to our setting by letting
the nodes’ step-sizes and their strong convexity parameters be
equal, the results from [18] yield the convergence factor:


(1 − σ)2 µ3/2
√
.
1−Ω
N L3/2
Comparing this with bound (30) obtained here, we can see that
our bound exhibits better scaling with respect to the number
of nodes N ; on the other hand, the bounds in [18], [19] exhibit
better scaling in terms of condition number (L/µ).
Remark 5 It is interesting to observe what happens in Theorem 4 when b = µ – the case that corresponds to a novel
1/2
exact method. Then, one has L0 = L2 − µ2
, which is
arbitrarily small when L becomes close to µ. Hence, for wellconditioned problems (L close to µ), the maximal

admissible
step size in Theorem 4 becomes α = Ω (1−σ)µ
, and the
2
 L 2
µ
corresponding convergence factor is 1 − Ω (1−σ)
. Hence,
L2
according to the upper bounds derived here and in [14], the
proposed method with
 b =2 µ
improvesconvergence
 factor
(1−σ) µ2
(1−σ) µ2
over [14] from 1−Ω
to 1−Ω
for wellL2
L2
conditioned problems (L sufficiently close to µ). Simulations
confirm large gains in convergence speed of the new method
over [14] (see Section V).
B. Setting up analysis and preliminary Lemmas
The proof of Theorem 4 is based on the small gain Theorem [40]. It is worth noting that analytical approaches using
small gain theorem have been developed in [18] and [19] prior
to our work and hold for a more general network and stepsize models than what we study here. Namely, reference [19]
carries out analysis under undirected and directed time varying
graphs, while reference [18] allows for node-dependent stepsizes. However, the analysis here is different. First, it is carried
out for a different algorithm in general (the proposed method
(16)–(17) that involves matrix B); as such, it reveals the novel
insight that the negative coupling from the primal to the dual
error can be partially compensated by tuning matrix B; see
also Remark 5. Second, when the results here are specialized
to the method in [18], [19], the obtained convergence factors
are different. (See Remark 4 for details.)
Besides the small gain Theorem, this Subsection also reviews some known results on the convergence of inexact
(centralized) gradient methods that will be used subsequently.

We start by reviewing the setting of the small gain theorem [40]. Denote by a := a(0) , a(1) , ..., a(k) , ... an infinite
sequence of vectors, a(k) ∈ Rp , k = 0, 1, ... For a fixed δ ∈
(0, 1), define the following quantities:


1 (k)
ka
k
(31)
kakδ,K :=
max
k=0,...,K
δk


1 (k)
kakδ := sup
ka
k
.
(32)
δk
k≥0
0

Clearly, we have that kakδ,K ≤ kakδ,K ≤ kakδ , for
K 0 ≥ K ≥ 0. It is also clear that, if kakδ is finite, then the sequence a(k) converges to zero R-linearly. Indeed, provided that
kakδ ≤ Ca < ∞, there holds: ka(k) k ≤ Ca δ k , k = 0, 1, ...
Hence, if for some δ ∈ (0, 1) we have that kakδ is finite,
then the sequence {a(k) } converges to zero R-linearly with
convergence factor at most δ. We state the small gain theorem
for two infinite sequences a and b as this suffices for our
analysis; for the more general Theorem involving an arbitrary
(finite) number of infinite sequences, see, e.g., [40], [18].
Theorem 5 Consider two infinite sequences a = a(0) , a(1) , ...,
and b = b(0) , b(1) , ..., with a(k) , b(k) ∈ Rp , k = 0, 1, ...
Suppose that for some δ ∈ (0, 1), and for all K = 0, 1, ...,
there holds:
kakδ,K
kbk

δ,K

≤ γ1 kbkδ,K + ω1
≤ γ2 kak

δ,K

+ ω2 ,

(33)
(34)

where γ1 , γ2 ≥ 0 and γ1 γ2 < 1. Then, there holds:
kakδ ≤

1
(ω2 γ1 + ω1 ) .
1 − γ1 γ2

(35)

We will frequently use the following simple Lemma.
Lemma 6 Consider two infinite sequences a = a(0) , a(1) , ...,
and b = b(0) , b(1) , ..., with a(k) , b(k) ∈ Rp . Suppose that, for
all k = 0, 1, ..., there holds:
ka(k+1) k ≤ c1 ka(k) k + c2 kb(k) k,

(36)

where ci ≥ 0, i = 1, 2. Then, for all K = 0, 1, ..., for any
δ ∈ (0, 1), we have:
c1
c2
kakδ,K ≤
kakδ,K +
kbkδ,K + ka(0) k.
(37)
δ
δ
1
Proof: Divide inequality (36) by δk+1
, δ ∈ (0, 1). The
resulting inequality implies, for all K = 1, 2, ... that:






c1
1 (k)
(k+1)
ka
k
ka
k
≤
max
k=0,...,K−1
δ k+1
δ k=0,...,K−1 δ k


c2
c1
1 (k)
c2
max
kb k = kakδ,K−1 + kbkδ,K−1
δ k=0,...,K−1 δ k
δ
δ
c1
c2
≤ kakδ,K + kbkδ,K .
(38)
δ
δ
max

+

1

Note that (38) implies that, for all K = 0, 1, ...


1
(k+1)
kakδ,K =
max
ka
k
k=−1,0,...,K−1 δ k+1
c1
c2
≤
kakδ,K + kbkδ,K + ka(0) k,
δ
δ
which is precisely what we wanted to show.

8

The use of the Lemma will be to bound kakδ,K by kbkδ,K .
Namely, whenever c3 := cδ1 < 1, (36) implies the following
bound:
c2 /δ
1
kakδ,K ≤
kbkδ,K +
ka(0) k.
(39)
1 − c3
1 − c3
We will also need the following Lemma from [41] on the
convergence of inexact (centralized) gradient methods.
Lemma 7 Consider unconstrained minimization of function
φ : Rp → R, where φ is assumed to be strongly convex
with strong convexity parameter m, and it also has Lipschitz
continuous gradient with Lipschitz constant M , M ≥ m > 0.
Consider the following inexact gradient method with step size
1
:
γ≤M


y (k+1) = y (k) − γ ∇φ(y (k) ) + (k) , k = 0, 1, ...,
with arbitrary initialization y (0) ∈ Rp and (k) ∈ Rp . Then,
for all k = 0, 1, ..., there holds:
ky (k+1) − y ? k ≤ (1 − γ m) ky (k) − y ? k + γ k(k) k,

Proof: Consider (17). Multiplying the equality from the
left by N1 (1 ⊗ I)> , using (1 ⊗ I)> L = (1> ⊗ I)((I − W ) ⊗ I)
= (1> (I − W )) ⊗ I = 0, we obtain that:
u(k+1) = u(k) , k = 0, 1, ...

(41)

Recall that u(0) = 0, by assumption. Thus, the result.
Lemma 9 Consider algorithm
(16)–(17), with α ≤ 1/L. Then,

for any δ ∈ 1 − α2µ , 1 , there holds:
4L
2
kex kδ,K ≤ √
ke(0) k.
ke
xkδ,K +
αµ x
Nµ

(42)

Proof: Consider (16). Multiplying the equation from the
left by N1 ( 1 ⊗ I )> , using the fact that N1 ( 1 ⊗ I )> W =
1
1
>
> (k)
= u(k) = 0,
N ( 1 ⊗ I ) , and the fact that N ( 1 ⊗ I ) u
for all k (by Lemma 8), we obtain:
x(k+1) = x(k) −

(40)

N
N
α X
α X
(k)
∇fi (xi ) = x(k) −
∇fi (x(k) )
N i=1
N i=1

N

α X
(k)
∇fi (xi ) − ∇fi (x(k) )
N i=1

α 
∇f (x(k) ) + (k) , where
= x(k) −
N
N 

X
(k)
∇fi (xi ) − ∇fi (x(k) ) .
=

−

where y ? = arg miny∈Rp φ(y).
Quantity (k) in (40) is an inexactness measure that says how
far is the employed search direction from the exact gradient
at the iterate y (k) .

(k)

i=1

(k)

C. Intermediate Lemmas and proof of Theorem 4
We now carry out convergence proof of Theorem 4 through
a sequence of intermediate Lemmas. We first split theprimal
(k)
error as follows: ex = x(k) − x• = x(k) − 1 ⊗ x(k) +1⊗

(k)
(k)
?
(k)
e(k) = x(k) −1⊗x(k)
x − x =: x
e +1⊗ex . Quantity x
(k)
says how mutually different are the solution estimates xi ’s at
(k)
different nodes; quantity ex says how far is the global average
PN (k)
1
(k)
x = N i=1 xi from the solution x? . We decompose the
(k)
dual error eu = u(k) + ∇F (x• ) in the following way:
e(k)
e(k) + 1 ⊗ e(k)
u .
u =u
P
(k)
(k)
(k)
N
Here, eu = N1 (1 ⊗ I)> eu = N1 i=1 [eu ]i , and u
e(k) =
(k)
(k)
(k)
(k)
(I − J )eu = eu − 1 ⊗ eu . Note that eu = u(k) +
PN
PN (k)
1
(k)
?
, where u(k) = N1 i=1 ui .
i=1 ∇fi (x ) = u
N
We will be interested in deriving bounds on kex kδ,K =
(k)
maxk=0,...,K δ1k kex k, for some δ ∈ (0, 1), and on the
analogous quantities that correspond to other primal and dual
errors that we defined above. Specifically, the proof path is as
(k)
follows. First, Lemma 8 shows that eu = 0, for all k, which
simplifies further analysis. Our goal is to apply Theorem 5
e,
with the following identification of infinite sequences: a → x
e . Then, Lemmas 9-11 are devoted to deriving
and b → u
a bound like in (33), while Lemma 12 devises a bound
that corresponds to (34). As shown below, this sequence of
Lemmas will be sufficient to complete the proof of Theorem 4.
(k)

Lemma 8 With algorithm (16)–(17), there holds: eu
u(k) = 0, for all k = 0, 1, ...

=

Now, applying Lemma 7, with ex = x(k) − x? , we obtain:
α (k)
ke(k+1)
k ≤ (1 − α µ) ke(k)
k k.
(43)
x
x k+
N
We next upper bound
α (k)
k k
N

≤

α
N

k(k) k as follows:

N
α X
(k)
∇fi (xi ) − ∇fi (x(k) )
N i=1

≤

N
α X
(k)
L xi − x(k)
N i=1

(44)

αL
≤√
x
e(k) .
(45)
N
Inequality (44) is by the Lipschitz continuity of the ∇fi ’s,
while (45) is by noting that
N
X

(k)

kxi

i=1

− x(k) k =

N
X

(k)

ke
xi k ≤

√

N ke
x(k) k.

i=1

Substituting the last bound in (43), we obtain:
α L (k)
ke(k+1)
k ≤ (1 − α µ) ke(k)
ke
x k.
x
x k+ √
N
Now, applying Lemma 6, we obtain:

(46)

1
αL 1
(1 − α µ) kex kδ,K + √
ke
xkδ,K + ke(0)
x k.
δ
Nδ
(47)
From (47), using α ≤ 1/L, one can verify that, for all δ ≥
1 − α2µ , there holds:

α µ
kex kδ,K ≤
1−
kex kδ,K
2
kex kδ,K ≤

9

αL
2 √ ke
xkδ,K + ke(0)
x k.
N
The last bound yields the desired result.

Lemma 10 Let α <
holds:
ke
xk

δ,K

≤

1−σ
3L ,

and δ ∈ 1 −

αµ
2 ,1



(48)

. Then, there

√
12 α L N
kex kδ,K
5 1−σ
12 α
2
+
ke
ukδ,K +
ke
x(0) k.
5 1−σ
1−σ

Proof: Consider (16). Subtracting x• = 1 ⊗ x? from both
sides of the equation, and noting that W x• = (W ⊗I)(1⊗x? )
= (W 1) ⊗ (I x? ) = 1 ⊗ x? = x• , we obtain:


(k)
(k)
(k)
.
(49)
e(k+1)
=
W
e
−
α
∇F
(x
)
+
u
x
x
Next, note that


∇F (x(k) ) + u(k) = ∇F (x(k) ) − ∇F (x• )


+ ∇F (x• ) + u(k)
(50)


= ∇F (x(k) ) − ∇F (x• ) + u
e(k) .
(51)
(k)

e(k) (by
In (51), we use that ∇F (x• ) + u(k) = eu = u
Lemma 8). Substituting (51) into (49), we get:


(k)
e(k+1)
= W e(k)
) − ∇F (x• ) − α u
e(k) (52)
x
x − α ∇F (x
We next multiply (52) from the left by I − J = (I − J) ⊗ I,
(k)
e(k) . Noting that
and we use that [ (I − J) ⊗ I ]ex = x
[ (I − J) ⊗ I ] W

=

[ (I − J) ⊗ I ] [ W ⊗ I ]


Next, note that for δ ∈ 1 − α2µ , 1 , and α <
holds:
1
σ+1
(σ + α L) <
.
δ
2

2

5 (1−σ)µ
12 α L
Next, note that, for α < 96
L2 , we have that (1−σ) µ <
1/2. Substituting the latter bound in (59) and manipulating the
terms, the desired result follows.

1+σ
Lemma 12 Let δ
∈
, 1 , and recall L0
=
2
1/2
2
2
L + b − 2bµ
. Then, the following holds:

fx
W
e(k) − α (I − J ) ∇F (x(k) )

−

∇F (x• )) − α u
e(k) .
(k)

•

(53)
(k)
ex

(k)

(k)
+1⊗ex ,

Next, use the decomposition x −x =
=x
e
and Lipschitz continuity of ∇F , to note that:
√
k∇F (x(k) ) − ∇F (x• )k ≤ L ke
x(k) k + L N ke(k)
x k.

(54)

Using the latter bound and taking the 2-norm in (53), while
using its sub-additive and sub-multiplicative properties, we
get:
√
u(k) k.
ke
x(k+1) k ≤ (σ + α L) ke
x(k) k + α L N ke(k)
x k + α ke
(55)
Now, similarly to Lemma 6, it is easy to see that the last
equation implies:
1
ke
xk
≤ (σ + α L) ke
xkδ,K
δ
√
1
α
+ α L N kex kδ,K + ke
ukδ,K + ke
x(0) k. (56)
δ
δ
δ,K

≤

40 L0 L
2
ke
xkδ,K +
ke
u(0) k
µ (1 − σ)
1−σ
√
16 L0 N
kex (0) k.
+
1 − σ αµ

Proof: Consider (17). Adding ∇F (x• ) to both sides of
(k)
e(k) ,
the equality, and recalling that eu = u(k) + ∇F (x• ) = u
we obtain:



=

(57)

Proof: We combine Lemmas 9 and 10, to obtain:
√
48 α L2
24
NL
δ,K
δ,K
ke
xk
≤
ke
xk
+
ke(0) k
5 (1 − σ) µ
5 (1 − σ) µ x
12 α
2
+
ke
ukδ,K +
ke
x(0) k.
(59)
5 1−σ
1−σ

= [ (I − J)W ] ⊗ I
f ⊗ I = W,
f
=W

x
e(k+1)

there

1
Also, as α < 1−σ
3L < 3L , we have that 1/δ ≤ 6/5. Substituting
the last two bounds in (56), we obtain:
√
1+σ
6
ke
xkδ,K + α L N kex kδ,K
ke
xkδ,K ≤
2
5
6α
δ,K
+
ke
uk
+ ke
x(0) k.
5
Rearranging terms in the last equality, the desired result
follows.

αµ
5 (1−σ)µ
Lemma 11 Let α < 96
L2 , and δ ∈ 1 − 2 , 1 . Then,
there holds:
4
24 α
ke
ukδ,K +
ke
x(0) k (58)
ke
xkδ,K ≤
5 1−
σ
1
−
σ
√
48 L N
+
ke(0) k.
5 µ (1 − σ) x

kukδ,K

and (I − J )e
u(k) = u
e(k) , we get:

1−σ
3L ,

u
e(k+1)

= u
e(k) − L ( ∇F (x(k) ) − ∇F (x• ) + u(k)
+ ∇F (x• ) − b (x(k) − x• ) )
(60)
= Wu
e(k) − L ( ∇F (x(k) ) − ∇F (x• ) )

+ b L (x(k) − x• )
fu
= W
e(k) − L ( ∇F (x(k) ) − ∇F (x• ) )
+ b L (x(k) − x• ).

(61)

Here, (60) uses the fact that L x• = 0, while (61) holds
fu
because W
e(k) = (W − J ) u
e(k) = W u
e(k) , as J u
e(k) =
J (I − J )eu (k) = (J − J )eu (k) = 0. We nextupper bound
the term ∇F (x(k) ) − ∇F (x• ) − b (x(k) − x• ) as follows:
k∇F (x(k) ) − ∇F (x• ) − b (x(k) − x• )k2
= k∇F (x(k) ) − ∇F (x• )k2 + b2 kx(k) − x• k2

> 

−2 b ∇F (x(k) ) − ∇F (x• )
x(k) − x•

10

≤ L2 kx(k) − x• k2 + b2 kx(k) − x• k2

max{ 1+σ
2 , 1−

αµ
2 }, 1



, as desired.

− 2 b µ kx(k) − x• k2 ,
where the last inequality holds by the Lipschitz continuity of
∇F , and by the strong monotonicity of ∇F :

V. S IMULATIONS

This Section provides a simulation example on learning a
linear
classifier via minimization of the `2 -regularized logistic
(∇F (x) − ∇F (y)) (x − y) ≥ µ kx−yk2 , for all x, y ∈ Rd . loss. Simulations confirm the insights gained through the
theoretical analysis and demonstrate that the proposed generHence, we have that:
alized method improves convergence speed over the existing
k∇F (x(k) ) − ∇F (x• ) − b (x(k) − x• )k ≤ L0 kx(k) − x• k.
methods [13], [14].
The simulation setup is as follows. We consider distributed
Next, taking the norm in (61), using kLk ≤ 2, exploiting its
learning
of a linear classifier via the `2 -regularized logistic
sub-additive and sub-multiplicative properties, and using the
loss,
e.g.,
[11]. Each node i has J data samples {aij , bij }Jj=1 .
Lipschitz continuity of ∇F , we obtain:
Here, aij ∈ Rd−1 , d ≥ 2, is a feature vector, and bij ∈
(k+1)
(k)
0
(k)
ke
u
k ≤ σ ke
u k + 2 L kex k.
(62) {−1, +1} is its class label. The goal is to learn a vector x =
>
d−1
(x>
, and xP
0 ∈ R, d ≥ 1, such that the total
(k)
(k)
1 , x0 ) , x1 ∈ R
Decomposing ex = x
e(k) + 1 ⊗ ex , and applying Lemma 6,
N
`2 -regularized surrogate loss i=1 fi (x) is minimized, where,
we obtain:
for i = 1, ..., N , we have:
2 L0
σ
δ,K
δ,K
δ,K
 1
ke
uk
+
ke
xk
ke
uk
≤
+ Rkxk2 .
fi (x) = ln 1 + exp −bij ( a>
δ √
δ
ij x1 + x0 )
2
2 L0 N
kex kδ,K + ke
+
u(0) k.
Here, R is a positive regularization parameter. We can take the
δ
strong convexity constant
µ = R, while a Lipschitz constant
Pas
N PJ
Next, applying Lemma 9, we get:
1
>
L
can
be
taken
as
k
i=1
j=1 cij cij k + R, where cij =
4N
 0 

>
>
σ
2L
4L
(bij aij , bij ) .
ke
ukδ,K ≤
ke
ukδ,K +
1+
ke
xkδ,K
With all experiments, we test the algorithms on a connected
δ
δ
µ
√
network
with N = 30 nodes and 123 links, generated as
0
4 NL
ke(0)
u(0) k.
+
a realization of the random
geometric graph model with
x k + ke
p
αµδ
communication radius ln(N )/N .

From now on, assume that δ ∈ 1+σ
We generate data and set the algorithm parameters as
2 , 1 . Then, it is easy to
.
Also,
note that 1δ ≤ 2. Thus, follows. Each node i has J = 2 data points whose dimension
see that there holds: σδ < 1+σ
2
we have:
is d − 1 = 5. The aij ’s are generated independently over
0
i and j; each entry of aij is drawn independently from the
1
+
σ
20
L
L
ke
ukδ,K ≤
ke
ukδ,K +
ke
xkδ,K
standard normal distribution. We generate the “true” vector
2
µ
√
x? = ((x?1 )> , x?0 )> by drawing its entries independently from
0
8L N
(0) δ,K
(0)
kex k
+ ke
u k.
+
standard normal distribution. The class
labels are generated

αµ
as bij = sign (x?1 )> aij + x?0 + ij , where ij ’s are drawn
After rearranging expressions, the desired result follows.
independently from normal distribution with zero mean and
We are now ready to prove Theorem 4.
variance 0.4. We set the regularization parameter as R = 0.03.
(0)
With all algorithms, we initialize xi to zero for all
Proof of Theorem 4: We apply Theorem 5 with the
e and b → u
e , by utilizing Lemma 11 i = 1, ..., N . The auxiliary variables for each algorithm
identification a → x
and Lemma 12. Assume the algorithm parameters obey the are initialized as described in Subsection II-B. Further, the
1
weight matrix is as follows: Wij = 2 (max{deg(i),deg(j)}+1)
,
conditions of Lemmas 11 and 12, namely that:


for
i
=
6
j,
{i,
j}
∈
E;
W
=
0,
for
i
=
6
j,
{i,
j}
∈
/
E;
and
ij
5 µ (1 − σ)
1+σ
αµ
P
α<
and δ ∈ max{
, 1−
}, 1 . Wii = 1 − j6=i Wij , for i = 1, ..., N . Here, deg(i) is the
2
96
L
2
2
number of neighbors of node i (excluding i).
Then, by the two Lemmas, the product of gains equals:
As an error metric, we use the following quantity:
6α 40 L0 L
?
PN kx(k)
γ1 γ2 = 1−σ
1
i −x k
µ (1−σ) We need that γ1 γ2 < 1 in order for (35)
, x? 6= 0, that we refer to as the relative
i=1
N
kx? k
µ (1−σ)2 µ(1−σ)
?
to hold. Therefore, when α < min{ 192 L0 L , 19 L2 }, we have error. Quantity x is obtained numerically beforehand by a
that ke
xkδ ≤ C < ∞, for a constant C ∈ (0, ∞). Note also centralized Nesterov gradient method [33].
that, by Lemma 9, we have:
We compare four methods: the method in [14], that we
refer
to here as “harnessing”; the Extra method in [13]; the
4L
2
kex kδ ≤ √
ke
xkδ +
ke(0)
(63) proposed method (16)–(17) with B = L+µ I – that we refer to
x k.
2
αµ
Nµ
as the modified “harnessing”; and the proposed method (16)–


(k)
(k)
As ex = x
e(k) +1⊗ex , we have: kex kδ ≤ 1 + √4NLµ ke
xkδ (17) with B = L W – that we refer to as the modified Extra.


Figure 1 plots the relative error versus number of iterations k
(0)
(0)
+ α2µ kex k = 1 + √4NLµ C+ α4µ kex k=: C 0 < +∞. for the four methods, for different values of step sizes: the
(k)
Therefore, kex k ≤ C 0 δ k , for all k, for any δ ∈ top Figure: α = 1/(3L); middle: α = 1/(9L); and bottom:
>

11

α = 1/(15L). First, on the top Figure, we can see that the
proposed modifications yield improvements in the convergence
speed over the respective original methods in [14] and [13].
While the improvement is not very large for [13], it is quite
significant for the method in [14].
As the step size decreases (the middle and bottom Figures),
we can see that the gain of the proposed method is reduced,
and the four methods tend to behave mutually very similarly.
Next, while for the large step size (the top Figure) Extra [13]
performs better than the method in [14], for the small step
size (bottom Figure) the performance of the two methods
is reversed, as predicted by our theoretical considerations.
(Though the difference between the methods is quite small
for the small step size.)
Figure 2 repeats the experiment for a 100-node, 561-link,
connected network (generated also as
p an instance of a random
geometric graph model with radius ln(N )/N ), and for stepsizes α = 1/(6L) (top Figure); α = 1/(18L) (middle); and
α = 1/(54L) (bottom). We can see that a similar behavior of
the four methods can be observed here, as well.
VI. C ONCLUSION
We considered exact first order methods for distributed optimization problems where N nodes collaboratively minimize
the aggregate sum of their local convex costs. Specifically, we
unified, generalized, and improved convergence speed of the
existing methods, e.g., [13], [14]. While it was known that
the method in [13] is equivalent to a primal-dual gradient-like
method, we show here that this is true also with [14], where the
corresponding primal-dual update rule incorporates a weighted
past dual gradient term, in addition to the current dual gradient term. We then generalize the method by proposing an
optimized, easy-to-tune weighting of the past dual gradient
term, and show both theoretically and by simulation that the
modification yields significant improvements in convergence
speed. We establish for the proposed exact method global
R-linear convergence rate, assuming strongly convex local
costs with Lipschitz continuous gradients and static networks.
Possible future work directions include extensions to uncoordinated step-sizes across nodes and time varying, random,
and directed networks. Specifically, with directed networks,
it would be interesting to extend the methods to work with
singly stochastic matrices, i.e., to relax the requirement on
the double stochasticity, as doubly stochastic methods become
undesirable when dealing with time-varying or directed graph
topologies.
.
A PPENDIX
A. Proof of Lemma 2
We first prove part (b), i.e., we set B = α1 W. We must
show that (16)–(17) is equivalent to (5)–(6). It suffices to show
that (16)–(17) is equivalent to (12)–(13), due to Lemma 1.
First, note that, due to identity L = I − W, we have that (16)
and (12) are the same. Next, from (16), we have that:
1
1
∇F (x(k) ) + u(k) = − x(k+1) + W x(k) .
α
α

Fig. 1: Relative error versus number of iterations k, for three
different values of step-size α: Top: α = 31L ; middle: α = 91L ;
and bottom: α = 151L .

12

Substituting this into (17), and using B = α1 W, we recover (13). Hence, the equivalence of (16)–(17) and (5)–(6)
with B = α1 W. Next, we prove part (a). That is, we consider
(14) and (16)–(17) with B = 0. The key is to note that s(k) ,
k = 0, 1, ..., can be written as: s(k) = u(k) +∇F (x(k) ), where
u(k) is defined by the following recursion:
u(k+1)

=

W u(k) + (W − I) ∇F (x(k) )

=

u(k) − L (u(k) + ∇F (x(k) )),

(64)

for k = 0, 1, ..., and u(0) = 0. Substituting (64) into (16),
yields the desired equivalence.
B. Proof of Lemma 3
Consider (16). Subtracting x• from both sides of the equality, noting that W x• = x• , and adding and subtracting
∇F (x• ) to the term in the parenthesis on the right hand side
of the equality, we obtain:
e(k+1)
x

=



(k)
W e(k)
) − ∇F (x• ) + ∇F (x• ) + u(k)
x − α ∇F (x

=

(k)
(k)
W e(k)
x − α Hk ex − α eu ,

(65)

by the
 definition of Hk , Hk =
R 1where2 (65)• follows
(k)
•
∇
F
x
+
t
(x
−
x
)
dt, and by the definition of
t=0
(k)
eu . Next, consider (17). Using identity L = I − W, the
equality can be equivalently written as:


u(k+1) = W u(k) − L ∇F (x(k) ) − B x(k) .
(66)
Next, add ∇F (x• ) to both sides of the equality, and express
the quantity on the right hand side as ∇F (x• ) = W ∇F (x• )+
L ∇F (x• ). We obtain:

e(k+1)
= W (u(k) + ∇F (x• )) − L ∇F (x(k) )
u
−

∇F (x• )) + L B x(k)

(k)
(k)
= W e(k)
.
u − L Hk ex + L B x

(67)

Next, by the property B x• = c x• , for some c ∈ R, it follows
that: L B x• = c L x• = c [ (I − W ) ⊗ I ] [1 ⊗ x? ] = c [ (I −
W )1 ] ⊗ [ I ⊗ x? ] = 0. Hence, we can subtract L B x• = 0
from the right hand side of (67) to obtain the following:
e(k+1)
u

(k)
(k)
= W e(k)
u − L Hk ex + L B ex
(k)
= W e(k)
u + L (B − Hk ) ex .

(68)

Finally, note from (68) that:
J e(k+1)
= J e(k)
u
u , k = 0, 1, ...,
because J L = [ J ⊗ I ] [ (I − W ) ⊗ I ] = [ J(I − W ) ] ⊗ I =
(0)
[J − J] ⊗ I = 0. Note that J eu = J (0 + ∇F (x• )) =
P
(k)
N
1
?
i=1 ∇fi (x ) ) = 0. Thus, we conclude that J eu = 0,
N 1(
for all k. Applying the latter fact to (68), we obtain:
(k)
e(k+1)
= (W − J ) e(k)
u
u + L (B − Hk ) ex .

(69)

The relations (65) and (69) yield the claim of the Lemma.
Fig. 2: Relative error versus number of iterations k, for three
different values of step-size α: Top: α = 61L ; middle: α =
1
1
18 L ; and bottom: α = 54 L .

C. Derivation of the solution to (20) and of an approximate
solution to (21)
We first consider (20). Note that, for any H ∈ H,
we have that: kb I − Hk = maxi=1,...,N d |b − hi | =

13

max {|b − hN d |, |h1 − b|} ≤ max {|b − µ|, |L − b|} . Here,
hi denotes the i-th largest eigenvalue of H. In the last
inequality above, we used the fact that µ I  H  L I, for all
H ∈ H. Therefore, we have that:
max kb I − Hk = max{|b − µ|, |L − b|}.
H∈H

(70)

The maximum in (70) is attained, e.g., for H =
Diag(L, µ, ..., µ), where µ is repeated (N d − 1) times. The
quantity (70) is clearly minimized over b ≥ 0 at b = b? =
L+µ
2 .
Now, consider (21), and assume that λN > 0. Note that
the maximal eigenvalue of W = W ⊗ I equals one, and the
minimal eigenvalue of W equals λN . We have:
kb0 W − Hk ≤ max {|b0 − µ|, |L − b0 λN |} .

(71)

We choose a sub-optimal b0 that minimizes the upper bound
in (71) on the desired function supH∈H kb0 W − Hk. It is easy
L+µ
.
to see that the corresponding value is b0 = 1+λ
N
D. Proof of equivalence of (16)–(17) and (28)–(29)
Consider (16). From the equation,
we have: ∇F (x(k) ) +

1
(k+1)
(k)
u = −α x
−Wx
. Substituting the latter relation
in (17), and using the fact that matrices L and W commute,
as well as that L and B commute, (17) leads to (25). Now,
consider (28). Proceeding by the same steps as in deriving
(12) from (10), it is straightforward to verify that (28) leads
to (24). Thus, the equivalence between (16)–(17) and (28)–
(29). The respective equivalence for the method in [14] follows
by setting B = 0 in (28)–(29).
(k)

E. Proof of Theorem 4 for generic matrices B
We show here that, when B = b I is replaced with a
generic symmetric matrix B that respects the sparsity pattern
of graph G and obeys that for any y ∈ Rd , there exists some
c ∈ R, such that B (1 ⊗ y) = c (1 ⊗ y), Theorem 4 continues
to hold with L0 replaced with constant (L + kBk). Namely,
it is easy to see that Lemma 8 continues to hold unchanged.
Further, Lemmas 9–11 pertain to the update equation (16) that
does not depend on B, and hence they also hold unchanged.
The only modifications occur with Lemma 12. Namely, (61)
fu
becomes: u
e(k+1) = W
e(k) − L ( ∇F (x(k) ) − ∇F (x• ) )
(k)
•
+L B (x − x ). Using Lipschitz continuity of ∇F and
the fact that kLk ≤ 2, the latter equality implies: ke
u(k+1) k
≤ σ ke
u(k) k + 2 L kx(k) − x• k +2kBk kx(k) − x• k. The proof
of the modified Lemma 12 then proceeds in the same way
as the remaining part of the proof of Lemma 12. Finally, the
proof of the modified Theorem 4 then also proceeds in the
same way as the proof of Theorem 4.
R EFERENCES
[1] A. Nedic and A. Ozdaglar, “Distributed subgradient methods for multiagent optimization,” IEEE Transactions on Automatic Control, vol. 54,
no. 1, pp. 48–61, January 2009.
[2] C. Lopes and A. H. Sayed, “Adaptive estimation algorithms over
distributed networks,” in 21st IEICE Signal Processing Symposium,
Kyoto, Japan, Nov. 2006.
[3] F. Cattivelli and A. H. Sayed, “Diffusion LMS strategies for distributed
estimation,” IEEE Trans. Sig. Process., vol. 58, no. 3, pp. 1035–1048,
March 2010.

[4] A. H. Sayed, S.-Y. Tu, J. Chen, X. Zhao, and Z. Towfic, “Diffusion
strategies for adaptation and learning over networks,” IEEE Sig. Process.
Mag., vol. 30, no. 3, pp. 155–171, May 2013.
[5] J. Duchi, A. Agarwal, and M. Wainwright, “Dual averaging for distributed optimization: Convergence and network scaling,” IEEE Trans.
Aut. Contr., vol. 57, no. 3, pp. 592–606, March 2012.
[6] I. Lobel and A. Ozdaglar, “Convergence analysis of distributed subgradient methods over random networks,” in 46th Annual Allerton
Conference onCommunication, Control, and Computing, Monticello,
Illinois, September 2008, pp. 353 – 360.
[7] D. Jakovetic, J. Xavier, and J. M. F. Moura, “Fast distributed gradient
methods,” IEEE Trans. Autom. Contr., vol. 59, no. 5, pp. 1131–1146,
May 2014.
[8] ——, “Convergence rates of distributed Nesterov-like gradient methods
on random networks,” IEEE Transactions on Signal Processing, vol. 62,
no. 4, pp. 868–882, February 2014.
[9] S. Kar, J. M. F. Moura, and K. Ramanan, “Distributed parameter estimation in sensor networks: Nonlinear observation models and imperfect
communication,” IEEE Transactions on Information Theory, vol. 58,
no. 6, pp. 3575–3605, June 2012.
[10] M. Rabbat and R. Nowak, “Distributed optimization in sensor networks,”
in IPSN 2004, 3rd International Symposium on Information Processing
in Sensor Networks, Berkeley, California, USA, April 2004, pp. 20 –
27.
[11] 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 in Machine Learning, Michael
Jordan, Editor in Chief, vol. 3, no. 1, pp. 1–122, 2011.
[12] F. Bullo, J. Cortes, and S. Martinez, Distributed control of robotic
networks: A mathematical approach to motion coordination algorithms.
Princeton University Press, 209.
[13] W. Shi, Q. Ling, G. Wu, and W. Yin, “EXTRA: An exact first-order
algorithm for decentralized consensus optimization,” SIAM J. Optim.,
vol. 25, no. 2, pp. 944–966, 2015.
[14] G. Qu and N. Li, “Harnessing smoothness to accelerate distributed
optimization,” to appear in IEEE Transactions on Control of Network
Systems, 2017, DOI: 10.1109/TCNS.2017.2698261.
[15] K. Yuan, B. Ying, X. Zhao, and A. H. Sayed, “Exact diffusion for
distributed optimization and learning — Part I: Algorithm development,”
2017, arxiv preprint, arXiv:1702.05122.
[16] ——, “Exact diffusion for distributed optimization and learning — Part
II: Convergence analysis,” 2017, arxiv preprint, arXiv:1702.05142.
[17] D. Jakovetic, J. Xavier, and J. M. F. Moura, “Linear convergence rate of
a class of distributed augmented Lagrangian algorithms,” IEEE Trans.
Autom. Contr., vol. 60, no. 4, pp. 922–936, April 2015.
[18] A. Nedic, A. Olshevsky, W. Shi, and C. A. Uribe, “Geometrically convergent distributed optimization with uncoordinated step-sizes,” 2016,
arXiv preprint arXiv:1609.05877.
[19] A. Nedic, A. Olshevsky, and W. Shi, “Achieving geometric convergence
for distributed optimization over time-varying graphs,” 2016, arXiv
preprint arXiv:1607.03218.
[20] J. Xu, S. Zhu, Y. Soh, and L. Xie, “Augmented distributed gradient
methods for multi-agent optimization under uncoordinated constant
stepsizes,” in 54th IEEE Conference on Decision and Control (CDC),
2015, pp. 2055–2060.
[21] Q. Lu and H. Li, “Geometrical convergence rate for distributed optimization with time-varying directed graphs and uncoordinated step-sizes,”
2016, arXiv preprint arXiv:1611.00990.
[22] G. Qu and N. Li, “Accelerated distributed Nesterov gradient descent,”
2017, arxiv preprint arXiv:1705.07176.
[23] C. Xi and U. A. Khan, “Dextra: A fast algorithm for optimization over
directed graphs,” IEEE Transactions on Automatic Control, 2017, to
appear, DOI: 10.1109/TAC.2017.2672698.
[24] J. Zeng and W. Yin, “Extrapush for convex smooth decentralized
optimization over directed networks,” Journal of Computational Mathematics, vol. 35, no. 4, pp. 381–394, 2017.
[25] W. Shi, Q. Ling, G. Wu, and W. Yin, “A proximal gradient algorithm
for decentralized composite optimization,” Journal of Machine Learning
Research, vol. 63, no. 22, pp. 6013–6023, 2015.
[26] K. Scaman, F. Bach, S. Bubeck, Y. T. Lee, and L. Massoulie, “Optimal
algorithms for smooth and strongly convex distributed optimization in
networks,” 2017, arxiv preprint, arXiv:1702.08704.
[27] A. Nedic, A. Ozdaglar, and A. Parrilo, “Constrained consensus and
optimization in multi-agent networks,” IEEE Transactions on Automatic
Control, vol. 55, no. 4, pp. 922–938, April 2010.

14

[28] A. Mokhtari and A. Ribeiro, “DSA: Decentralized double stochastic
averaging gradient algorithm,” Journal of Machine Learning Research,
vol. 17, pp. 1–35, 2016.
[29] A. Nedic and A. Ozdaglar, “Subgradient methods for saddle point
problems,” Journal of Optimization Theory and Applications, vol. 145,
no. 1, pp. 205–228, July 2009.
[30] M. Kallio and C. H. Rosa, “Large-scale convex optimization via saddlepoint computation,” Oper. Res., pp. 93–101, 1999.
[31] H. Uzawa, “Iterative methods in concave programming,” 1958, in Arrow,
K., Hurwicz, L., Uzawa, H. (eds.) Studies in Linear and Nonlinear
Programming, pp. 154-165. Stanford University Press, Stanford.
[32] J. Tsitsiklis, D. Bertsekas, and M. Athans, “Distributed asynchronous
deterministic and stochastic gradient optimization algorithms,” IEEE
Trans. Autom. Contr., vol. 31, no. 9, pp. 803–812, Sep. 1986.
[33] Y. E. Nesterov, “A method for solving the convex programming problem
with convergence rate O(1/k2 ),” Dokl. Akad. Nauk SSSR, vol. 269, pp.
543–547, 1983, (in Russian).
[34] D. Jakovetic, D. Bajovic, N. Krejic, and N. K. Jerinkic, “Distributed
gradient methods with variable number of working nodes,” IEEE Transactions on Signal Processing, vol. 64, no. 15, pp. 4080–4095, August
2016.
[35] M. Zhu and S. Martinez, “On distributed convex optimization under
inequality and equality constraints,” IEEE Transactions on Automatic
Control, vol. 57, no. 1, pp. 151–164, Jan. 2012.
[36] T.-H. Chang, A. Nedic, and A. Scaglione, “Distributed constrained optimization by consensus-based primal-dual perturbation method,” IEEE
Transactions on Automatic Control, vol. 59, no. 6, pp. 1524–1538, June
2014.
[37] J. Wang and N. Elia, “Control approach to distributed optimization,”
in 48th Annual Allerton Conference onCommunication, Control, and
Computing, Monticello, IL, Oct. 2010.
[38] A. Dimakis, S. Kar, J. M. F. Moura, M. Rabbat, and A. Scaglione,
“Gossip algorithms for distributed signal processing,” Proceedings of
the IEEE, vol. 98, no. 11, pp. 1847–1864, 2010.
[39] K. Yuan, Q. Ling, and W. Yin, “On the convergence of decentralized
gradient descent,” SIAM J. Optim., vol. 26, no. 3, pp. 1835–1854, 2016.
[40] C. Desoer and M. Vidyasagar, Feedback Systems: Input-Output Properties. SIAM, 2009.
[41] M. Schmidt, N. L. Roux, and F. Bach, “Convergence rates of inexact
proximal-gradient methods for convex optimization,” in Advances in
Neural Information Processing Systems 24, 2011, pp. 1458–1466.

