1

The proximal augmented Lagrangian method
for nonsmooth composite optimization

arXiv:1610.04514v3 [math.OC] 4 Mar 2018

Neil K. Dhingra, Sei Zhen Khong, and Mihailo R. Jovanović

Abstract—We study a class of optimization problems in which the
objective function is given by the sum of a differentiable but possibly nonconvex component and a nondifferentiable convex regularization term.
We introduce an auxiliary variable to separate the objective function
components and utilize the Moreau envelope of the regularization term
to derive the proximal augmented Lagrangian – a continuously differentiable function obtained by constraining the augmented Lagrangian
to the manifold that corresponds to the explicit minimization over
the variable in the nonsmooth term. This function is used to develop
customized algorithms based on the method of multipliers (MM) and
a primal-descent dual-ascent gradient method in order to compute
optimal primal-dual pairs. MM algorithm is applicable to a broader
class of problems than proximal gradient methods and it has stronger
convergence guarantees and a more refined step-size update rules than
the alternating direction method of multipliers. These features make it an
attractive option for solving structured optimal control problems. When
the differentiable component of the objective function is (strongly) convex
and the regularization term is convex, we prove global (exponential)
asymptotic stability of the primal-descent dual-ascent algorithm. Finally,
we solve the edge addition in directed consensus networks and optimal
placement problems to demonstrate the merits and the effectiveness of
our approach.

I. I NTRODUCTION
We study a class of composite optimization problems in which the
objective function is a sum of a differentiable but possibly nonconvex
component and a convex nondifferentiable component. Problems of
this form are encountered in diverse fields including compressive
sensing [1], machine learning [2], statistics [3], image processing [4],
and control [5]. In feedback synthesis, they typically arise when a
traditional performance metric (such as the H2 or H∞ norm) is
augmented with a regularization function to promote certain structural
properties in the optimal controller. For example, the `1 norm and
the nuclear norm are commonly used nonsmooth convex regularizers
that encourage sparse and low-rank optimal solutions, respectively.
The lack of a differentiable objective function precludes the use of
standard descent methods for smooth optimization. Proximal gradient
methods [6] and their accelerated variants [7] generalize gradient
descent, but typically require the nonsmooth term to be separable
over the optimization variable. Furthermore, standard acceleration
techniques are not well-suited for problems with constraint sets that
do not admit an easy projection (e.g., closed-loop stability).
An alternative approach is to split the smooth and nonsmooth
components in the objective function over separate variables which
are coupled via an equality constraint. Such a reformulation facilitates the use of the alternating direction method of multipliers
(ADMM) [8]. This augmented-Lagrangian-based method splits the
optimization problem into subproblems which are either smooth or
easy to solve. It also allows for a broader class of regularizers than
proximal gradient and it is convenient for distributed implementation.
However, there are limited convergence guarantees for nonconvex
problems and parameter tuning greatly affects its convergence rate.
Financial support from the National Science Foundation under Awards
ECCS-1739210 and CNS-1544887 and the Air Force Office of Scientific
Research under Award FA9550-16-1-0009 is gratefully acknowledged.
N. K. Dhingra is with the Department of Electrical and Computer Engineering, University of Minnesota, Minneapolis, MN 55455; S. Z. Khong is
with the Institute for Mathematics and its Applications, Minneapolis, MN
55455; M. R. Jovanović is with the Ming Hsieh Department of Electrical
Engineering, University of Southern California, Los Angeles, CA, 90089. Emails: dhin0008@umn.edu, szkhong@ima.umn.edu, mihailo@usc.edu.

The method of multipliers (MM) is the most widely used algorithm
for solving constrained nonlinear programing problems [9]–[11]. In
contrast to ADMM, it is guaranteed to converge for nonconvex
problems and there are systematic ways to adjust algorithmic parameters. However, MM is not a splitting method and it requires joint
minimization of the augmented Lagrangian with respect to all primal
optimization variables. This subproblem is typically nonsmooth and
as difficult to solve as the original optimization problem.
To make this difficult subproblem tractable, we transform the
augmented Lagrangian into the continuously differentiable proximal
augmented Lagrangain by exploiting the structure of proximal operators associated with nonsmooth regularizers. This new form is
obtained by constraining the augmented Lagrangian to the manifold
that corresponds to the explicit minimization over the variable in the
nonsmooth term. The resulting expression is given in terms of the
Moreau envelope of the nonsmooth regularizer and is continuously
differentiable. This allows us to take advantage of standard optimization tools, including gradient descent and quasi-Newton methods, and
enjoy the convergence guarantees of standard MM.
The proximal augmented Lagrangian also enables Arrow-HurwiczUzawa primal-dual gradient flow dynamics. Such dynamics can be
used to identify saddle points of the Lagrangian [12] and have enjoyed
recent renewed interest in the context of networked optimization
because, in many cases, the gradient can be computed in a distributed
manner [13]. Our approach yields a dynamical system with a continuous right-hand side for a broad class of nonsmooth optimization
problems. This is in contrast to existing techniques which employ
subgradient methods [14] or use discontinuous projected dynamics [15]–[17] to handle inequality constraints. Furthermore, since the
proximal augmented Lagrangian is not strictly convex-concave we
make additional developments relative to [18] to show asymptotic
stability. Finally, inspired by recent advances [19], [20], we employ
the theory of integral quadratic constraints (IQCs) [21] to prove global
exponential stability when the differentiable component of the objective function is strongly convex with a Lipschitz continuous gradient.
The rest of the paper is structured as follows. In Section II,
we formulate the nonsmooth composite optimization problem and
provide a brief background on proximal operators. In Section III, we
exploit the structure of proximal operators to introduce the proximal
augmented Lagrangian. In Section III-B, we provide an efficient
algorithmic implementation of the method of multipliers using the
proximal augmented Lagrangian. In Section IV, we prove global (exponential) asymptotic stability of primal-descent dual-ascent gradient
flow dynamics under a (strong) convexity assumption. In Section V,
we use the problems of edge addition in directed consensus networks
and optimal placement to illustrate the effectiveness of our approach.
We close the paper in Section VI with concluding remarks.
II. P ROBLEM FORMULATION AND BACKGROUND
We consider a composite optimization problem,
minimize
x

f (x) + g (T (x))

(1)

where the optimization variable x belongs to a finite-dimensional
Hilbert space (e.g., Rn or Rm×n ) equipped with an inner product
h·, ·i and associated norm k · k. The function f is continuously
differentiable but possibly nonconvex, the function g is convex but

2

potentially nondifferentiable, and T is a bounded linear operator. We
further assume that g is proper and lower semicontinuous, that (1) is
feasible, and that its minimum is finite.
Problem (1) is often encountered in structured controller design [22]–[24], where f is a measure of closed-loop performance,
e.g., the H2 norm, and the regularization term g is introduced
to promote certain structural properties of T (x). For example, in
wide-area control of power systems, f measures the quality of
synchronization between different generators and g penalizes the
amount of communication between them [25]–[27].
P
In particular, for z := T (x) ∈ Rm , the `1 norm, kzk1 :=
|zi |,
is a commonly used convex proxy for promoting
sparsity
of
z.
For
P
z ∈ Rm×n , the nuclear norm, kzk∗ :=
σi (z), can be used to
obtain low-rank solutions to (1), where σi (z) is the ith singular value.
The indicator function, IC (z) := {0, z ∈ C; ∞, z 6∈ C} associated
with the convex set C is the proper regularizer for enforcing z ∈ C.
Regularization of T (x) instead of x is important in the situations
where the desired structure has a simple characterization in the codomain of T . Such problems arise in spatially-invariant systems,
where it is convenient to perform standard control design in the
spatial frequency domain [28] but necessary to promote structure
in the physical space, and in consensus/synchronization networks,
where the objective function is expressed in terms of the deviation
of node values from the network average but it is desired to impose
structure on the network edge weights [23], [24].
Problem (1) is difficult to solve directly because f is, in general,
a nonconvex function and g is typically not differentiable. Since the
existing approaches and our method utilize proximal operators, we
first provide a brief overview; for additional information, see [6].
The proximal operator of the function g is given by
proxµg (v) := argmin g(x) +
x

1
2µ

kx − vk2

(2a)

and the associated optimal value specifies its Moreau envelope,
1
2µ

kproxµg (v) − vk2

(2b)

where µ > 0. The Moreau envelope is a continuously differentiable
function, even when g is not, and its gradient [6] is given by

∇Mµg (v) = µ1 v − proxµg (v) .
(2c)
P
For example, when g is the `1 norm, g(z) = kzk1 =
|zi |, the
proximal operator is determined by soft-thresholding, proxµg (vi ) =
Sµ (vi ) := sign(vi ) max (|vi | − µ, 0), the associated Moreau enve1
lope is the Huber function, Mµg (vi ) = { 2µ
vi2 , |vi | ≤ µ; |vi | −
µ
,
|v
|
≥
µ},
and
the
gradient
of
this
Moreau
envelope is the
i
2
saturation function, ∇Mµg (vi ) = sign(vi ) min (|vi |/µ, 1).
B. Existing algorithms
1) Proximal gradient: The proximal gradient method generalizes
standard gradient descent to certain classes of nonsmooth optimization problems. This method can be used to solve (1) when g(T ) has
an easily computable proximal operator. When T = I, the proximal
gradient method for problem (1) with step-size αl is given by,
xl+1 = proxαl g (xl − αl ∇f (xl )).
When g = 0, the proximal gradient method simplifies to standard
gradient descent, and when g is indicator function of a convex set, it
simplifies to projected gradient descent. The proximal gradient algorithm applied to the `1 -regularized least-squares problem (LASSO)
minimize
x

1
2

minimize

f (x) + g(z)

subject to

T (x) − z = 0.

x, z

(4)

The augmented Lagrangian is obtained by adding a quadratic penalty
on the violation of the linear constraint to the regular Lagrangian
associated with (4),
Lµ (x, z; y) = f (x) + g(z) + hy, T (x) − zi +

1
2µ

kT (x) − zk2

where y is the Lagrange multiplier and µ is a positive parameter.
ADMM solves (4) via an iteration which involves minimization
of Lµ (x, z; y) separately over x and z and a gradient ascent update
(with step-size 1/µ) of y [8],
xk+1 = argmin Lµ (x, z k ; y k )

(5a)

x

z k+1 = argmin Lµ (xk+1 , z; y k )

A. Background on proximal operators

Mµg (v) := g(proxµg (v)) +

αl AT (Axl − b)). This method is effective only when the proximal
operator of g(T ) is easy to evaluate. Except in special cases, e.g,
when T is diagonal, efficient computation of proxµg(T ) does not
necessarily follow from an efficiently computable proxµg . This
makes the use of proximal gradient method challenging for many
applications and its convergence can be slow. Acceleration techniques
improve the convergence rate [7], [29], but they do not perform well
in the face of constraints such as closed-loop stability.
2) Augmented Lagrangian methods: A common approach for
dealing with a nondiagonal linear operator T in (1) is to introduce
an additional optimization variable z

2

kAx − bk + γ kxk1

(3)

where γ is a positive regularization parameter, yields the Iterative Soft-Thresholding Algorithm (ISTA) [7], xl+1 = Sγαl (xl −

y

k+1

= y

k

z

+

1
µ

k+1

(T (x

)−z

k+1

(5b)
).

(5c)

ADMM is appealing because, even when T is nondiagonal, the
z-minimization step amounts to evaluating proxµg , and the xminimization step amounts to solving a smooth (but possibly nonconvex) optimization problem. Although it was recently shown that
ADMM is guaranteed to converge to a stationary point of (4) for
some classes of nonconvex problems [30], its rate of convergence is
strongly influenced by the choice of µ.
The method of multipliers (MM) is the most widely used algorithm
for solving constrained nonconvex optimization problems [9], [31]
and it guarantees convergence to a local minimum. In contrast
to ADMM, each MM iteration requires joint minimization of the
augmented Lagrangian with respect to the primal variables x and z,
(xk+1 , z k+1 ) = argmin Lµ (x, z; y k )
y

k+1

= y

x, z
k

+

1
µ

k+1

(T (x

)−z

(6a)
k+1

).

(6b)

It is possible to refine MM to allow for inexact solutions to the
(x, z)-minimization subproblem and adaptive updates of the penalty
parameter µ. However, until now, MM has not been a feasible choice
for solving (4) because the nonconvex and nondifferentiable (x, z)minimization subproblem is as difficult as the original problem (1).
III. T HE PROXIMAL AUGMENTED L AGRANGIAN
We next derive the proximal augmented Lagrangian, a continuously
differentiable function resulting from explicit minimization of the
augmented Lagrangian over the auxiliary variable z. This brings the
(x, z)-minimization problem (6a) into a form that is continuously
differentiable with respect to both x and y and facilitates the use of a
wide suite of standard optimization tools for solving (1). In particular,
as described below, our approach enables the method of multipliers
and the Arrow-Hurwicz-Uzawa gradient flow dynamics method.
A. Derivation of the proximal augmented Lagrangian
The first main result of the paper is provided in Theorem 1.
We use the proximal operator of the function g to eliminate the
auxiliary optimization variable z from the augmented Lagrangian and
transform (6a) into a tractable continuously differentiable problem.

3

Theorem 1: For a proper, lower semicontinuous, and convex
function function g, minimization of the augmented Lagrangian
Lµ (x, z; y) associated with problem (4) over (x, z) is equivalent to
minimization of the proximal augmented Lagrangian
Lµ (x; y) := f (x) + Mµg (T (x) + µy) −

µ
2

2

kyk

(7)

over x. Moreover, if f is continuously differentiable Lµ (x; y) is
continuously differentiable over x and y, and if f has a Lipschitz
continuous gradient ∇Lµ (x; y) is Lipschitz continuous.
Proof: Through the completion of squares, the augmented Lagrangian Lµ associated with (4) can be equivalently written as
Lµ (x, z; y) = f (x) + g(z) +

1
2µ

kz − (T (x) + µy)k2 −

µ
2

kyk2 .

Minimization with respect to z yields an explicit expression,
zµ? (x, y) = proxµg (T (x) + µy)

(8)

zµ?

and substitution of
into the augmented Lagrangian provides (7),
i.e., Lµ (x; y) = Lµ (x, zµ? (x, y); y). Continuous differentiability
of Lµ (x; y) follows from continuous differentiability of Mµg and
Lipschitz continuity of ∇Lµ (x; y) follows from Lipschitz continuity
of proxµg and boundedness of the linear operator T ; see (2c).
Expression (7), that we refer to as the proximal augmented Lagrangian, characterizes Lµ (x, z; y) on the manifold corresponding to
explicit minimization over the auxilary variable z. Theorem 1 allows
joint minimization of the augmented Lagrangian with respect to x and
z, which is in general a nondifferentiable problem, to be achieved by
minimizing differentiable function (7) over x. It thus facilitates the
use of the method of multipliers in Section III-B and the ArrowHurwicz-Uzawa gradient flow dynamics in Section IV.
Remark 1: The proximal augmented Lagrangian can be derived
even in the presence of a more general linear constraint,
minimize
x1 , x2

subject to

f (x1 ) + g(x2 )

(9a)

T1 (x1 ) + T2 (x2 ) = 0.

Introduction of an additional auxiliary variable z in the nonsmooth
part of the objective function g, can be used to bring this two-block
optimization problem into the following form,
minimize
x1 , x2 , z

subject to

f (x1 ) + g(z)
T1 (x1 ) + T2 (x2 ) = 0, x2 − z = 0.

(9b)

Via an analogous procedure to that described in Theorem 1, explicit
minimization with respect to z can be employed to eliminate it from
the augmented Lagrangian and obtain a continuously differentiable
function of both primal (x1 , x2 ) and dual (y1 , y2 ) variables,
Lµ (x1 , x2 ; y1 , y2 ) = f (x1 ) +

1
2µ

kT1 (x1 ) + T2 (x2 ) + µy1 k2 +

Mµg (x2 + µy2 ) −

µ
2

ky1 k2 −

µ
2

ky2 k2 .

Here, y1 and y2 are the Lagrange multipliers associated with the
respective linear constraints in (9b) and, for simplicity, we use
single parameter µ in the augmented Lagrangian. This approach has
numerous advantages over standard ADMM; e.g., it can be readily
extended to multi-block optimization problems for which ADMM is
not guaranteed to converge in general [32]. These extensions are outside of the scope of the present study and will be reported elsewhere.
B. MM using the proximal augmented Lagrangian
Theorem 1 allows us to solve nondifferentiable subproblem (6a)
by minimizing the continuously differentiable proximal augmented
Lagrangian Lµ (x; y k ) over x. Relative to ADMM, our customized
MM algorithm guarantees convergence to a local minimum and offers
systematic update rules for the parameter µ. Relative to proximal
gradient, we can solve (1) with a general bounded linear operator T
and can incorporate second order information about f .

Using reformulated expression (7) for the augmented Lagrangian,
MM minimizes Lµ (x; y k ) over the primal variable x and updates the
dual variable y using gradient ascent with step-size 1/µ,
xk+1 = argmin Lµ (x; y k )
y

k+1

= y

k

x

+

1
µ

k+1

∇y Lµ (x

(MMa)
k

;y )

(MMb)

where ∇y Lµ (xk+1 ; y k ) = T (xk+1 ) − proxµg (T (xk+1 ) + µy k )
denotes the primal residual, i.e., the difference between T (xk+1 )
and zµ? (xk+1 , y k ).
In contrast to ADMM, our approach does not attempt to avoid
the lack of differentiability of g by fixing z to minimize over x.
By constraining Lµ (x, z; y) to the manifold resulting from explicit
minimization over z, we guarantee continuous differentiability of
the proximal augmented Lagrangian Lµ (x; y). MM is a gradient
ascent algorithm on the Lagrange dual of a version of (4), with
the same constraint, in which the objective function is replaced by
1
f (x)+g(z)+ 2µ
kT (x)−zk2 ; see [8, Section 2.3] and [33]. Since its
closed-form expression is typically unavailable, MM uses the (x, z)minimization subproblem (6a) to evaluate this dual computationally
and then take a gradient ascent step (6b) in y. ADMM avoids
this issue by solving simpler, separate subproblems over x and z.
However, the x and z minimization steps (5a) and (5b) do not
solve (6a) and thus unlike the y-update (6b) in MM, the y-update (5c)
in ADMM is not a gradient ascent step on the “strengthened dual”.
MM thus offers stronger convergence results [8], [9] and may lead
to fewer y-update steps.
Remark 2: The proximal augmented Lagrangian enables MM
because the x-minimization subproblem in MM (MMa) is not more
difficult than in ADMM (5a). For LASSO problem (3), the zupdate in ADMM (5b) is given by soft-thresholding, z k+1 =
Sγµ (xk+1 + µy k ), and the x-update (5a) requires minimization of
the quadratic function [8]. In contrast, the x-update (MMa) in MM
requires minimization of (1/2) kAx − bk2 +Mµk g (x + µk y k ), where
Mµk g (v) is the Moreau envelope associated with the `1 norm; i.e.,
the Huber function. Although in this case the solution to (5a) can
be characterized explicitly by a matrix inversion, this is not true in
general. The computational cost associated with solving either (5a)
or (MMa) using first-order methods scales at the same rate.
1) Algorithm: The procedure outlined in [11, Algorithm 17.4]
allows minimization subproblem (MMa) to be inexact, provides a
method for adaptively adjusting µk , and describes a more refined
update of the Lagrange multiplier y. We incorporate these refinements
into our proximal augmented Lagrangian algorithm for solving (4).
In Algorithm 1, η and ω are convergence tolerances, and µmin is a
minimum value of the parameter µ. Because of the equivalence established in Theorem 1, convergence to a local minimum follows from
the convergence results for the standard method of multipliers [11].
2) Minimization of Lµ (x; y) over x: MM alternates between
minimization of Lµ (x; y) with respect to x (for fixed values of µ and
y) and the update of µ and y. Since Lµ (x; y) is once continuously
differentiable, many techniques can be used to find a solution to
subproblem (MMa). We next summarize three of them.
Gradient descent: The gradient with respect to x is given by,
∇x Lµ (x; y) = ∇f (x) +

1
µ

T † (T (x) + µy − proxµg (T (x) + µy))

where T † is the adjoint of T , hz, T (x)i = T † (z), x . Backtracking
conditions such as the Armijo rule can be used to select a step-size.
Proximal gradient: Gradient descent does not exploit the structure of the Moreau envelope of the function g; in some cases, it
may be advantageous to use proximal operator associated with the
Moreau envelope to solve (MMa). In particular, when T = I, (2a)
and (2c) imply that proxαMµg (v) = x∗ where x∗ satisfies,

x∗ = µ +1 α α proxµg (x∗ ) + µ v . If g is separable and has an

4

Algorithm 1 MM using the proximal augmented Lagrangian.
0

0

input: Initial point x and Lagrange multiplier y
initialize: µ0 = 10−1 , µmin = 10−5 , ω0 = µ0 , and η0 = µ0.1
0
for k = 0, 1, 2, . . .
Solve (MMa) such that
if k∇y Lµk (x

k+1

k∇x Lµ (xk+1 , y k )k ≤ ωk

k

; y )k ≤ ηk

if k∇y Lµk (xk+1 ; y k )k ≤ η and k∇x Lµ (xk+1 , y k )k ≤ ω
stop with approximate solution xk+1
else:
y k+1 = y k +

µk+1 = µk

ηk+1 =

= ωk µk+1

1
∇y Lµk (xk+1 ; y k ),
µk
ηk µ0.9
ωk+1
k+1 ,

else:
y k+1 = y k ,
ηk+1 = µ0.1
k+1 ,

µk+1 = max{µk /5, µmin }
ωk+1 = µk+1

easily computable proximal operator, its Moreau envelope also has
an easily computable proximal operator. In [34], proximal gradient
methods were used for subproblem (MMa) to solve a sparse feedback
synthesis problem introduced in [5]. Computational savings were
shown relative to standard proximal gradient method and ADMM.
Quasi-Newton method: Although proxµg is typically not differentiable, it is Lipschitz continuous and therefore differentiable
almost everywhere [35]. To improve computational efficiency, we
employ the limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) method [11, Algorithm 7.4] which estimates the Hessian
∇xx Lµ (x; y k ) using first-order information and is guaranteed to converge for convex functions with Lipschitz continuous gradients [36].
Remark 3: For regularization functions that do not admit simply
computable proximal operators, proxµg has to be evaluated numerically by solving (2a). If this is expensive, the primal-descent
dual-ascent algorithm of Section IV offers an appealing alternative
because it requires one evaluation of proxµg per iteration. When the
regularization function g is nonconvex, the proximal operator may not
be single-valued and the Moreau envelope may not be continuously
differentiable. In spite of this, the convergence of proximal algorithms
has been established for nonconvex, proper, lower semicontinuous
regularizers that obey the Kurdyka-Łojasiewicz inequality [37].
IV. A RROW-H URWICZ -U ZAWA GRADIENT FLOW
We now consider an alternative approach to solving (1). Instead of
minimizing over the primal variable and performing gradient ascent
in the dual, we simultaneously update the primal and dual variables
to find the saddle point of the augmented Lagrangian. The continuous
differentiability of Lµ (x; y) established in Theorem 1 enables the use
of Arrow-Hurwicz-Uzawa gradient flow dynamics [12],




ẋ
−∇x Lµ (x; y)
=
.
(GF)
ẏ
∇y Lµ (x; y)
In Section IV-A, we show that the continuous-time gradient flow
dynamics (GF) globally converge to the set of saddle points of
the proximal augmented Lagrangian Lµ (x; y) for a convex f . In
Section IV-B, we employ the theory of IQCs to establish global exponential stability for a strongly convex f with a Lipschitz continuous
gradient and estimate convergence rates. Finally, in Section IV-C we
identify classes of problems for which dynamics (GF) are convenient
for distributed implementation and compare/contrast our framework
to the existing approaches.
A. Global asymptotic stability for convex f
We first characterize the optimal primal-dual pairs of optimization
problem (4) with the Lagrangian, f (x) + g(z) + hy, T (x) − zi . The

associated first-order optimality conditions are given by,
0 = ∇f (x? ) + T † (y ? )
0 ∈ ∂g(z ? ) − y ?
0 = T (x? ) − z ?

(10a)
(10b)
(10c)

where ∂g is the subgradient of g. Clearly, these are equivalent to the
optimality condition for (1), i.e., 0 ∈ ∇f (x? ) + T † (∂g(T (x? ))).
To show convergence, we introduce a lemma about proximal
operators which follows almost directly from their definition. Even
though we state the result for x ∈ Rn , g: Rm → R, and T (x) = T x
where T ∈ Rm×n is a given matrix, the proof for x in a Hilbert space
and a bounded linear operator T follows from similar arguments.
Lemma 2: Let g: Rm → R be a proper, lower semicontinuous,
convex function and let proxµg : Rm → Rm be the corresponding
proximal operator. Then, for any a, b ∈ Rm , we can write
proxµg (a) − proxµg (b) = D (a − b)

(11)

where D is a symmetric matrix satisfying I  D  0.
Proof: Let c̃ := a − b, p̃ := proxµg (a) − proxµg (b), and
D := {p̃ p̃T /(p̃T c̃), p̃ 6= 0; 0, otherwise}. Then, by construction,
D = DT  0 and (11) can be written as p̃ = Dc̃. Since proxµg
is firmly nonexpansive [6], p̃T c̃ ≥ kp̃k2 or, equivalently, c̃TD(I −
D)c̃ ≥ 0 for every c̃ ∈ Rm . Positive semi-definiteness of I − D thus
follows from D  0 and commutativity of D and I − D.
Theorem 3: Let f be a continuously differentiable convex function,
and let g be a proper, lower semicontinuous, convex function. Then,
for the primal-descent dual-ascent gradient flow dynamics (GF),



 
ẋ
− ∇f (x) + T T ∇Mµg (T x + µy)
=
(GF1)
ẏ
µ (∇Mµg (T x + µy) − y)
the set of optimal primal-dual pairs (x? , y ? ) of (4) is globally
asymptotically stable, and x? is an optimal solution of (1).
Proof: We introduce a change of variables x̃ := x−x? , ỹ := y−
?
y and a Lyapunov function candidate, V (x̃, ỹ) = 21 hx̃, x̃i+ 21 hỹ, ỹi ,
where (x? , z ? ; y ? ) is an optimal solution to (4) that satisfies (10). The
dynamics in the (x̃, ỹ)-coordinates are given by,




x̃˙
−(∇f (x) − ∇f (x? ) + (1/µ) T T m̃)
=
− µ ỹ
m̃
ỹ˙
(12)
where
m̃ := µ (∇Mµg (T x + µy) − ∇Mµg (T x? + µy ? )).

(13a)

Using expression (11) to construct D such that,
D(T x̃ + µỹ) = proxµg (T x + µy) − proxµg (T x? + µy ? ) (13b)
and definition (2c) of ∇Mµg , we can write
m̃ = (I − D)(T x̃ + µỹ).

(13c)

Thus, the derivative of V along the solutions of (12) is given by
V̇ (x̃, ỹ) = − hx̃, ∇f (x) − ∇f (x? )i
− (1/µ) hT x̃, (I − D) T x̃i − µ hỹ, Dỹi .
Since f is convex, hx̃, ∇f (x) − ∇f (x? )i ≥ 0 and Lemma 2 imply
V̇ ≤ 0.
When ∇f (x) = ∇f (x? ) and when matrices D and T T (I − D)T
have nontrivial kernels, it is possible that V̇ = 0 for a nonzero ỹ ∈
ker{D} and x̃ such that T x̃ ∈ ker{(I − D)}. From (13c), these
conditions imply m̃ = µỹ and (12) simplifies to, x̃˙ = −T T ỹ, ỹ˙ = 0.
Thus, the only invariant set for dynamics (12) is ∇f (x? ) = ∇f (x),
ỹ ∈ ker{T T } ∩ ker{D}, and T x̃ ∈ ker{I − D}. Global asymptotic
stability of these points follows from LaSalle’s invariance principle.
To complete the proof, we show that any x and y that lie in this
invariant set also satisfy the optimality conditions (10) of problem (4)

5

with z = zµ? (x, y) and thus that x solves problem (1).
Since x? and y ? are optimal, (10a) implies ∇f (x? ) + T T y ? = 0.
For any x and y that lie in the invariant set of (12), we can replace
∇f (x? ) with ∇f (x) and add T T ỹ = 0 to obtain
T

?

T

which implies that x and y also satisfy (10a). Furthermore, (I −
D)T x̃ = 0 and Dỹ = 0 can be combined to yield T x̃ − D(T x̃ +
µỹ) = 0, which, by (13b), leads to
T x̃ − (proxµg (T x + µy) − proxµg (T x? + µy ? )) = 0.
By optimality condition (10b), T x? = z ? , and by definition (8),
z ? = zµ? (x? , y ? ) = proxµg (T x? + µy ? ). It thereby follows that
T x − proxµg (T x + µy) = 0 and, thus, T x = proxµg (T x + µy) =
zµ? (x, y) = z which implies that x and z satisfy (10c).
Finally, the optimality condition of the minimization problem (2a)
that defines proxµg (v) is ∂g(z) + µ1 (z − v) 3 0. Setting v = T x +
µy from the expression (8) that characterizes the zµ? -manifold and
T x = z by (10c) yields the final optimality condition (10b).
B. Global exponential stability for strongly convex f
We express (GF), or equivalently (GF1), as a linear system G
connected in feedback with nonlinearities that correspond to the
gradients of f and of the Moreau envelope of g; see Fig. 1. These
nonlinearities can be conservatively characterized by IQCs. Exponential stability of G connected in feedback with any nonlinearity
that satisfies these IQCs implies exponential convergence of (GF) to
the primal-dual optimal solution of (4). In what follows, we adjust
the tools of [19], [20] to our setup and establish global exponential
stability by evaluating the feasibility of an LMI. We assume that the
function f is mf -strongly convex with an Lf -Lipschitz continuous
gradient. Characterizing additional structural restrictions on f and g
with IQCs may lead to tighter bounds on the rate of convergence.
As illustrated in Fig. 1, (GF1) can be expressed as a linear system
G connected via feedback to a nonlinear block ∆,
ẇ = A w + B u, ξ = C w, u = ∆(ξ)






−I − µ1 T T
−mf I
I 0
, B=
A=
, C=
−µI
T µI
0
I
where w := [ xT y T ]T , ξ := [ ξ1T ξ2T ]T , and u := [ uT1 uT2 ]T .
Nonlinearity ∆ maps the system outputs ξ1 = x and ξ2 = T x + µy
to the inputs u1 and u2 via u1 = ∆1 (ξ1 ) := ∇f (ξ1 ) − mf ξ1 and
u2 = ∆2 (ξ2 ) := µ∇Mµg (ξ2 ) = ξ2 − proxµg (ξ2 ).
When the mapping ui = ∆i (ξi ) is the Li -Lipschitz continuous
gradient of a convex function, it satisfies the IQC [19, Lemma 6]

T 


ξi − ξ0
0
L̂i I
ξi − ξ0
≥ 0
(14)
ui − u0
ui − u0
L̂i I −2I
where L̂i ≥ Li , ξ0 is some reference point, and u0 = ∆i (ξ0 ).
Since f is mf -strongly convex, the mapping ∆1 (ξ1 ) is the gradient
of the convex function f (ξ1 ) − (mf /2)kξ1 k2 . Lipschitz continuity
of ∇f with parameter Lf implies Lipschitz continuity of ∆1 (ξ1 )
with parameter L1 := Lf − mf ; thus, ∆1 satisfies (14) with L̂1 ≥
L1 . Similarly, ∆2 (ξ2 ) is the scaled gradient of the convex Moreau
envelope and is Lipschitz continuous with parameter 1; thus, ∆2 also
satisfies (14) with L̂2 ≥ 1. These two IQCs can be combined into
η := [ ξ T uT ]T .

∇f − mf I

ξ2 = T x + µy

(15)

For a linear system G connected in feedback with nonlinearities
that satisfy IQC (15), [20, Theorem 3] establishes ρ-exponential
convergence, i.e., kw(t) − w? k ≤ τ e−ρt kw(0) − w? k for some
τ, ρ > 0, by verifying the existence of a matrix P  0 such that,

 T
 

 T
Aρ P + P Aρ P B
C 0
C
0
+
Π
 0, (16)
0 I
0
I
BT P
0

u2
u1

ξ1 = x

∇f (x) + T (y + ỹ) = ∇f (x) + T y = 0

(η − η0 )T Π (η − η0 ) ≥ 0,

∆
µ∇Mµg

G

Fig. 1: Block diagram of gradient flow dynamics (GF1) where G is
a linear system in feedback with nonlinearities that satisfy (15).
where Aρ := A + ρI. In Theorem 4, we determine a scalar condition
that ensures global exponential stability when T T T is full rank.
Theorem 4: Let f be strongly convex with parameter mf , let its
gradient be Lipschitz continuous with parameter Lf , let g be proper,
lower semicontinuous, and convex, and let T T T be full rank. Then,
if µ ≥ Lf − mf , there is a ρ > 0 such that the dynamics (GF)
converge ρ-exponentially to the optimal point of (4).
Proof: Since any function that is Lipschitz continuous with
parameter L is also Lipschitz continuous with parameter L̂ ≥ L,
we establish the result for µ = L̂1 := Lf − mf and L̂2 = 1.
We utilize [20, Theorem 3] to show ρ-exponential convergence
by verifying matrix inequality (16) through a series of equivalent
expressions (17). We first apply the KYP Lemma [38, Theorem 1]
to (16) to obtain an equivalent frequency domain characterization

∗ 

Gρ (jω)
Gρ (jω)
Π
 0, ∀ ω ∈ R
(17a)
I
I
where Gρ (jω) = C(jωI − Aρ )−1 B. Evaluating the left-hand side
of (17a) for L = µ and dividing by −2 yields the matrix inequality


µm̂ + m̂2 + ω 2
m̂
I
TT
2
2
2
2


m̂ + ω
m̂ + ω
 0 (17b)

m̂/µ
ω 2 − ρµ̂ 
T
∗
TT + 2
I
2
2
2
m̂ + ω
µ̂ + ω
where m̂ := mf − ρ > 0 and µ̂ := µ − ρ > 0 so that Aρ is Hurwitz,
i.e., the system Gρ is stable. Since the (1, 1) block in (17b) is positive
definite for all ω, the matrix in (17b) is positive definite if and only
if the corresponding Schur complement is positive definite,
m̂/µ
ω 2 − ρµ̂
TTT + 2
I  0.
µm̂ + m̂2 + ω 2
µ̂ + ω 2

(17c)

We exploit the symmetry of T T T to diagonalize (17c) via a unitary coordinate transformation. This yields m scalar inequalities
parametrized by the eigenvalues λi of T T T . Multiplying the left-hand
side of these inequalities by the positive quantity (µ̂2 + ω 2 )(µm̂ +
m̂2 + ω 2 ) yields a set of equivalent, quadratic in ω 2 , conditions,
i
i
ω 4 + ( m̂λ
+ m̂2 + µm̂ − ρµ̂)ω 2 + m̂µ̂( µ̂λ
− ρ(µ + m̂)) > 0.
µ
µ
(17d)
Condition (17d) is satisfied for all ω ∈ R if there are no ω 2 ≥ 0
for which the left-hand side is nonpositive. When ρ = 0, both the
constant term and the coefficient of ω 2 are strictly positive, which
implies that the roots of (17d) as a function of ω 2 are either not real
or lie in the domain ω 2 < 0, which cannot occur for ω ∈ R. Finally,
continuity of (17d) with respect to ρ implies the existence a positive
ρ that satisfies (17d) for all ω ∈ R.
Remark 4: Each eigenvalue λi of a full rank matrix T T T is
positive and hence to estimate the exponential convergence rate it
suffices to check (17d) only for the smallest λi . A sufficient condition
for (17d) to hold for each ω ∈ R is positivity of the constant term and
the coefficient multiplying ω 2 . For ρ < min (mf , µ) these can be,
respectively, expressed as the following quadratic inequalities in ρ,

ρ2 − γ ρ + λmin > 0
2ρ2 − (γ + µ + mf ) ρ + γ mf > 0

6

where γ := µ + mf + λmin
. The solutions to these provide the
µ
following estimates of the exponential convergence rate: (i) ρ < ρ1
when mf ≥ µ; and (ii) ρ < min (ρ1 , ρ2 ) when mf < µ, where
p
ρ1 = 21 (γ − γ 2 − 4λmin )
p
ρ2 = 41 (γ + µ + mf − (γ + µ + mf )2 − 8 γ mf ).
Our explicit analytical expressions can be used to determine the
optimal value of µ ≥ Lf − mf to maximize the above decay rates.
Remark 5: A similar convergence rate result can be obtained by
applying [19, Theorem 4] to a discrete-time implementation of the
primal-descent dual-ascent dynamics that results from a forward Euler
discretization of (GF). The proof will be reported elsewhere.
C. Distributed implementation
Gradient flow dynamics (GF) are convenient for distributed implementation. If the state vector x corresponds to the concatenated states
of individual agents, xi , the sparsity pattern of T and the structure of
the gradient map ∇f : Rn → Rn dictate the communicationP
topology
required to form ∇Lµ in (GF). For example, if f (x) =
fi (xi )
is separable over the agents, then ∇fi (xi ) can be formed locally.
If in addition T T is an incidence matrix of an undirected network
with the graph Laplacian T T T , each agent need only share its state
xi with its neighbors and maintain dual variables yi that correspond
to its edges. A distributed implementation is also natural when the
mapping ∇f : Rn → Rn is sparse.
Our approach provides several advantages over existing distributed optimization algorithms. Even for problems (1) with nondifferentiable regularizers g, a formulation based on the proximal
augmented Lagrangian yields gradient flow dynamics (GF) with a
continuous right-hand side. This is in contrast to existing approaches
which employ subgradient methods [14] or use discontinuous projected dynamics [15]–[18]. Note that although the augmented La1
grangian Lµ (x, y; z) contains a quadratic term 2µ
kT (x) − zk2 , it
is not jointly strongly convex in x and z and the resulting proximal
augmented Lagrangian (7) is not strictly convex-concave in x and y.
Furthermore, when T is not diagonal, a distributed proximal gradient
cannot be implemented because the proximal operator of g(T x)
may not be separable. Finally, ADMM has been used for distributed
implementation in the situations where f is separable and T is an
incidence matrix. Relative to such a scheme, our method does not
require solving an x-minimization subproblem in each iteration and
provides a guaranteed rate of convergence.
Remark 6: Special instances of our framework have strong connections with the existing methods for distributed optimization on
graphs; e.g., [13], [14], [39].
P The networked optimization problem
of minimizing fP
(x̄) =
fi (x̄) over a single variable x̄ can be
reformulated as
fi (xi ) + g(T x) where the components fi of the
objective function are distributed over independent agents xi , x is the
aggregate state, T T is the incidence matrix of a strongly connected
and balanced graph, and g is the indicator function associated
with the set T x = 0. The g(T x) term ensures that at feasible
points, xi = xj = x̄ for all i and j. It is easy to show that
∇Mµg (v) = (1/µ) v and that the dynamics (GF) are given by,
ẋ = − ∇f (x) − (1/µ) L x − ỹ
ỹ˙ = β L x

(18)

where β > 0, L := T T T is the graph Laplacian of a connected
undirected network, and ỹ := T T y belongs to the orthogonal
complement of the vector of all ones. The only difference relative
to [13, Eq. (20)] and [39, Eq. (11)] is that −ỹ appears instead of
−Lỹ in equation (18) for the dynamics of the primal variable x.
Remark 7: Forward Euler discretization of (18) is given by
xk+1 = (I − (α/µ)L) xk − α ∇f (xk ) − α ỹ k
ỹ k+1 = ỹ k + α β L xk

(19)

where α is the step-size, and the EXTRA algorithm [40, Equation
(2.13)], which has received significant recent attention,
xk+1 = W xk − α ∇f (xk ) +

k−1
1 X
(W − I) xi
2 i=0

(20)

can be clearly recovered from (19) by setting β = 1/(2αµ) and
taking W = I − (α/µ)L in (20). This demonstrates power and
generality of the proposed approach for distributed optimization.
V. E XAMPLES
We solve the problems of edge addition in directed consensus
networks and optimal placement to illustrate the effectiveness of the
proximal augmented Lagrangian method.
A. Edge addition in directed consensus networks
A consensus network with N P
nodes converges to the average of the
initial node values ψ̄ = (1/N ) i ψi (0) if and only if it is strongly
connected and balanced [41]. Unlike for undirected networks [23],
[24], the problem of edge addition in directed consensus networks is
not known to be convex. The steady-state variance of the deviations
from average is given by the square of the H2 norm of,


Q1/2
ψ̇ = −(Lp + Lx ) ψ + d, ξ =
ψ
−R1/2 Lx
where d is a disturbance, Lp is a weighted directed graph Laplacian
of a plant network, Q := I − (1/N )11T penalizes the deviation
from average, and R  0 is the control weight. The objective is
to optimize the H2 norm (from d to ξ) by adding a few additional
edges, specified by the graph Laplacian Lx of a controller network.
To ensure convergence of ψ to the average of the initial node values, we require that the closed-loop graph Laplacian, L = Lp +Lx , is
balanced. This condition amounts to the linear constraint, 1T L = 0.
We express
P the directed graph
PLaplacian of the controller network as,
Lx =
i 6= j Lij zij =:
l Ll zl where zij ≥ 0 is the added
edge weight that connects node j to node i, Lij := ei eTi − ei eTj ,
ei is the ith basis vector in Rn , and the integer l indexes the
edges such that zl = zij and Ll = Lij . For simplicity, we
assume that the plant network Lp is balanced and connected. Thus,
enforcing thatP
L is balanced amounts to enforcing the linear constraint
1T Lx = 1T ( Ll zl ) =: (Ez)T = 0 on z, where E is the incidence
matrix [41] of the edges that may be added. Any vector of edge
weights z that satisfies this constraint can be written as z = T x where
the columns of T span the nullspace of the matrix E and provide
a basis for the space of balanced graphs, i.e., the cycle space [41].
Each feasible controller Laplacian can thus be written as,

X
X X
X
Lx =
Ll [T x]l =
Ll
(T ek ) xk =:
L̂k xk
l

l

l

k

k

(21a)
where the matrices L̂k are given by L̂k = l Ll [T ek ]l .
Since the mode corresponding to 1 is marginally stable, unobservable, and uncontrollable, we introduce a change of coordinates to the
deviations from average φ = V T ψ where V T 1 = 0 and discard the
average mode ψ̄ = 1T ψ. The energy of the deviations from average
is given by the the H2 norm squared of the reduced system,
D
E
f (x) = V T (Q + LTx RLx ) V, X , Â X + X ÂT + B̂ B̂ T = 0
(21b)
where X is the controllability gramian of the reduced system with
Â := −V T (Lp + Lx )V and B̂ := V T .
To balance the closed-loop H2 performance with the number of
added edges, we introduce a regularized optimization problem
P

xγ = argmin f (x) + γ 1T T x + I+ (T x).
x

(22)

Here, the regularization parameter γ > 0 specifies the emphasis on
sparsity relative to the closed-loop performance f , and I+ is the
indicator function associated with the nonnegative orthant Rm
+ . When
the desired level of sparsity for the vector of the added edge weights
zγ = T xγ has been attained, optimal weights for the identified set
of edges are obtained by solving,
minimize
x

f (x) + IZγ (T x) + I+ (T x)

(23)

where Zγ is the set of vectors with the same sparsity pattern as zγ
and IZγ is the indicator function associated with this set.
1) Implementation: We next provide implementation details for
solving (22) and (23). The proof of next lemma is omitted for brevity.
Lemma 5: Let a graph Laplacian of a directed plant network Lp
be balanced and connected and let Â, B̂, Lx , and V be as defined
in (21a)–(21b). The gradient of f (x) defined in (21b) is given by,
D
E
∇f (x) = 2 vec (R Lx V − V P ) XV T , L̂k
where X and P are the controllability and observability gramians
determined by (21b) and ÂT P + P Â + V T (Q + LTx RLx )V = 0.
The proximal operator associated with the regularization function
gs (z) := γ1T z + I+ (z) in (22) is proxµgs (vi ) =
vi −γµ},
Pmax{0,
2
{v
/(2µ),
vi ≤
the Moreau envelope is given by Mµgs (v) =
i
i
γµ; γ (vi − γµ/2), vi > γµ}, and ∇Mµgs (v) = max {v/µ, γ}.
The proximal operator of the regularization function in (23), gp (z) :=
IZγ (z)+I+ (z), is a projection onto the intersection of the set Zγ and
the nonnegative orthant, proxµgp (v) = PE (v), the Moreau envelope
2
1
is the distance to E := Zγ ∩ Rm
+ , Mµgp (v) = 2µ kv − PE (v)k
and ∇Mµgp (v) is determined by a vector pointing from E to v,
∇Mµgp (v) = µ1 (v − PE (v)).
2) Computational experiments: We solve (22) and (23) using
Algorithm 1, where L-BFGS is employed in the x-minimization
subproblem (MMa). For the plant network shown in Fig. 2a, Fig. 2b
illustrates the tradeoff between the number of added edges and the
closed-loop H2 norm. The added edges are identified by computing
the γ-parameterized homotopy path for problem (22), and the optimal
edge weights are obtained by solving (23). The red dashed lines in
Fig. 2a show the optimal set of 2 added edges. These are obtained for
γ = 3.5 and they yield 23.91% performance loss relative to the setup
in which all edges in the controller graph are used. We note that the
same set of edges is obtained by conducting an exhaustive search.
This suggests that the proposed convex regularizers may offer a good
proxy for solving difficult combinatorial optimization problems.
We also consider simple directed cycle graphs with N = 5 to
50 nodes and m = N 2 − N potential added edges. We solve (22)
for γ = 0.01, 0.1, 0.2, and R = I using the proximal augmented
Lagrangian MM algorithm (PAL), ADMM, and ADMM with an
adaptive heuristic for updating µ [8] (ADMM µ). The x-update in
each algorithm is obtained using L-BFGS. Since gs (T x) and gp (T x)
are not separable in x, proximal gradient cannot be used here.
Figure 3a shows the time required to solve problem (22) in terms
of the total number of potential added edges; Fig. 3b demonstrates
that PAL requires fewer outer iterations; and Fig. 3c illustrates that
the average computation time per outer iteration is roughly equivalent
for all three methods. Even with an adaptive update of µ, ADMM
requires more outer iterations which increases overall solve time relative to the proximal augmented Lagrangian method. Thus, compared
to ADMM, PAL provides computational advantage by reducing the
number of outer iterations (indexed by k in Algorithm 1 and in (5)).
B. Optimal placement problem
To illustrate the utility of our primal-descent dual-ascent approach,
we consider an example in which mobile agents aim to minimize their
Euclidean distances relative to a set of targets {bi } while staying

x1

x4

x2

x5

x6

x7

x3

performance loss (percent)

7

number of added edges
(a) Directed consensus network

(b) Performance-sparsity tradeoff

Fig. 2: (a) A balanced plant graph with 7 nodes and 10 directed edges
(solid black lines). A sparse set of 2 added edges (dashed red lines)
is identified by solving (22) with γ = 3.5 and R = I. (b) Tradeoff
between performance and sparsity resulting from the solution to (22)(23) for the network shown in Fig. 2a. Performance loss is measured
relative to the optimal centralized controller (i.e., all edges are used).
within a desired distance from their neighbors in a network with the
incidence matrix T ,
P
2
minimize
(24)
i (xi − bi ) + I[−1,1] (T x).
x
Here, T x is a vector of inter-agent distances which must be kept
within an interval [−1, 1]. Applying primal-descent dual-ascent update rules to (24) achieves path planning for first-order agents ẋ = u
with u = −∇x Lµ (x; y). The proximal operator is projection onto a
box, proxµI[−1,1] (z) = max(min(z, 1), −1), the Moreau envelope
P 2
1
is the distance squared to that set, MµI[−1,1] (z) = 2µ
S1 (zi ), and
∇MµI[−1,1] (z) = µ1 S1 (z). To update its state, each agent xi needs
information from its neighbors in a network with a Laplacian T T T .
Methods based on the subderivative are not applicable because
the indicator function is not subdifferentiable. Proximal methods are
hindered because the proximal operatorP
of I[−1,1] (T x) is difficult
to compute due to T . Since f (x) =
(xi − bi )2 is separable,
a distributed ADMM implementation can be applied; however, it
may require large discrete jumps in agent positions, which could
be unsuitable for vehicles. Moreover, when f is not separable a
distributed implementation of the x-minimization step (5a) in ADMM
would not be possible.
Figure 4 shows an implementation for a problem with 5 agents
whose set of targets changes position at time 5. The primal-descent
dual-ascent dynamics (GF) are simulated in M ATLAB using ode45.
VI. C ONCLUDING REMARKS
For a class of nonsmooth composite optimization problems that
arise in structured optimal control, we have introduced continuously
differentiable proximal augmented Lagrangian function. This function
is obtained by collapsing the associated augmented Lagrangian onto
the manifold resulting from explicit minimization over the variable in
the nonsmooth part of the objective function. Our approach facilitates
development of customized algorithms based on the method of
multipliers and the primal-descent dual-ascent method.
MM based on the proximal augmented Lagrangian is applicable to
a broader class of problems than proximal gradient methods, and it
has more robust convergence guarantees, more rigorous parameter
update rules, and better practical performance than ADMM. The
primal-descent dual-ascent gradient dynamics we propose are suitable
for distributed implementation and have a continuous right-hand
side. When the differentiable component of the objective function is
(strongly) convex, we establish (exponential) asymptotic convergence.
Finally, we illustrate the efficacy of our algorithms using the edge
addition and optimal placement problems. Future work will focus on
developing second-order updates for the primal and dual variables
and on providing an extension to nonconvex regularizers.

8

m

m

m

(a) Total solve time (s)

(b) Number of outer iterations

(c) Solve time (s) per outer iteration

position

Fig. 3: (a) Total time; (b) number of outer iterations; and (c) average time per outer iteration required to solve (22) with γ = 0.01, 0.1, 0.2
for a cycle graph with N = 5 to 50 nodes and m = 20 to 2450 potential added edges using PAL (–×–), ADMM (- -◦- -), and ADMM
with the adaptive µ-update heuristic [8] (· · ·  · · · ). PAL requires fewer outer iterations and thus a smaller total solve time.

time
Fig. 4: Set of 5 distributed agents tracking targets (black ◦) whose
optimal positions are determined by the solution to (24) (red ×).
R EFERENCES
[1] D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52,
no. 4, pp. 1289–1306, 2006.
[2] T. Hastie, R. Tibshirani, and J. Friedman, “Unsupervised learning,” in
The Elements of Statistical Learning. Springer, 2009, pp. 485–585.
[3] P. J. Bickel and E. Levina, “Regularized estimation of large covariance
matrices,” Ann. Stat., pp. 199–227, 2008.
[4] T. Goldstein and S. Osher, “The split Bregman method for `1 -regularized
problems,” SIAM J. Imaging Sci., vol. 2, no. 2, pp. 323–343, 2009.
[5] F. Lin, M. Fardad, and M. R. Jovanović, “Design of optimal sparse
feedback gains via the alternating direction method of multipliers,” IEEE
Trans. Automat. Control, vol. 58, no. 9, pp. 2426–2431, 2013.
[6] N. Parikh and S. Boyd, “Proximal algorithms,” Foundations and Trends
in Optimization, vol. 1, no. 3, pp. 123–231, 2013.
[7] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding
algorithm for linear inverse problems,” SIAM J. Imaging Sci., vol. 2,
no. 1, pp. 183–202, 2009.
[8] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed
optimization and statistical learning via the alternating direction method
of multipliers,” Found. Trends Mach. Learning, vol. 3, no. 1, pp. 1–124,
2011.
[9] D. P. Bertsekas, Constrained optimization and Lagrange multiplier
methods. New York: Academic Press, 1982.
[10] ——, Nonlinear programming. Athena Scientific, 1999.
[11] J. Nocedal and S. J. Wright, Numerical Optimization. Springer, 2006.
[12] K. J. Arrow, L. Hurwicz, and H. Uzawa, Studies in linear and non-linear
programming. Stanford University Press, 1958.
[13] J. Wang and N. Elia, “A control perspective for centralized and
distributed convex optimization,” in Proceedings of the 50th IEEE
Conference on Decision and Control and the 10th European Control
Conference, 2011, pp. 3800–3805.
[14] A. Nedic and A. Ozdaglar, “Distributed subgradient methods for multiagent optimization,” IEEE Trans. Automat. Control, vol. 54, no. 1, pp.
48–61, 2009.
[15] D. Feijer and F. Paganini, “Stability of primal–dual gradient dynamics
and applications to network optimization,” Automatica, vol. 46, no. 12,
pp. 1974–1981, 2010.
[16] A. Cherukuri, E. Mallada, and J. Cortés, “Asymptotic convergence of
constrained primal–dual dynamics,” Syst. Control Lett., vol. 87, pp. 10–
15, 2016.
[17] A. Cherukuri, E. Mallada, S. Low, and J. Cortes, “The role of convexity
on saddle-point dynamics: Lyapunov function and robustness,” IEEE
Trans. Automat. Control, 2018, doi:10.1109/TAC.2017.2778689.
[18] A. Cherukuri, B. Gharesifard, J., and Cortés, “Saddle-point dynamics:
conditions for asymptotic stability of saddle points,” SIAM J. Control
Optim., vol. 55, no. 1, pp. 486–511, 2017.
[19] L. Lessard, B. Recht, and A. Packard, “Analysis and design of optimization algorithms via integral quadratic constraints,” SIAM J. Optimiz.,
vol. 26, no. 1, pp. 57–95, 2016.

[20] B. Hu and P. Seiler, “Exponential decay rate conditions for uncertain linear systems using integral quadratic constraints,” IEEE Trans. Automat.
Control, vol. 61, no. 11, pp. 3631–3637, 2016.
[21] A. Megretski and A. Rantzer, “System analysis via integral quadratic
constraints,” IEEE Trans. on Automat. Control, vol. 42, no. 6, pp. 819–
830, 1997.
[22] M. R. Jovanović and N. K. Dhingra, “Controller architectures: tradeoffs
between performance and structure,” Eur. J. Control, vol. 30, pp. 76–91,
July 2016.
[23] X. Wu and M. R. Jovanović, “Sparsity-promoting optimal control of
systems with symmetries, consensus and synchronization networks,”
Syst. Control Lett., vol. 103, pp. 1–8, May 2017.
[24] S. Hassan-Moghaddam and M. R. Jovanović, “Topology design for
stochastically-forced consensus networks,” IEEE Trans. Control Netw.
Syst., 2017, doi:10.1109/TCNS.2017.2674962.
[25] F. Dörfler, M. R. Jovanović, M. Chertkov, and F. Bullo, “Sparse and
optimal wide-area damping control in power networks,” in Proceedings
of the 2013 American Control Conference, 2013, pp. 4295–4300.
[26] ——, “Sparsity-promoting optimal wide-area control of power networks,” IEEE Trans. Power Syst., vol. 29, no. 5, pp. 2281–2291,
September 2014.
[27] X. Wu, F. Dörfler, and M. R. Jovanović, “Input-output analysis and decentralized optimal control of inter-area oscillations in power systems,”
IEEE Trans. Power Syst., vol. 31, no. 3, pp. 2434–2444, May 2016.
[28] B. Bamieh, F. Paganini, and M. A. Dahleh, “Distributed control of
spatially-invariant systems,” IEEE Trans. Automat. Control, vol. 47,
no. 7, pp. 1091–1107, 2002.
[29] H. Li and Z. Lin, “Accelerated proximal gradient methods for nonconvex
programming,” in Adv. Neural Inf. Process Syst., 2015, pp. 379–387.
[30] M. Hong, Z.-Q. Luo, and M. Razaviyayn, “Convergence analysis of
alternating direction method of multipliers for a family of nonconvex
problems,” SIAM J. Optimiz., vol. 26, no. 1, pp. 337–364, 2016.
[31] A. R. Conn, N. I. M. Gould, and P. L. Toint, “A globally convergent augmented Lagrangian algorithm for optimization with general constraints
and simple bounds,” SIAM J. Numer. Anal., vol. 28, pp. 545–572, 1991.
[32] C. Chen, B. He, Y. Ye, and X. Yuan, “The direct extension of ADMM
for multi-block convex minimization problems is not necessarily convergent,” Math. Program., vol. 155, no. 1-2, pp. 57–79, 2016.
[33] R. T. Rockafellar, “Monotone operators and the proximal point algorithm,” SIAM J. Control Optim., vol. 14, no. 5, pp. 877–898, 1976.
[34] N. K. Dhingra and M. R. Jovanović, “A method of multipliers algorithm
for sparsity-promoting optimal control,” in Proceedings of the 2016
American Control Conference, Boston, MA, 2016, pp. 1942–1947.
[35] R. T. Rockafellar and R. J.-B. Wets, Variational analysis. Springer
Science & Business Media, 2009, vol. 317.
[36] J.-F. Bonnans, J. C. Gilbert, C. Lemaréchal, and C. A. Sagastizábal,
Numerical optimization: theoretical and practical aspects. Springer
Science & Business Media, 2013.
[37] J. Bolte, S. Sabach, and M. Teboulle, “Proximal alternating linearized
minimization for nonconvex and nonsmooth problems,” Math. Program.,
vol. 146, no. 1-2, pp. 459–494, 2014.
[38] A. Rantzer, “On the Kalman-Yakubovich-Popov lemma,” Syst. Control
Lett., vol. 28, no. 1, pp. 7–10, 1996.
[39] B. Gharesifard and J. Cortés, “Distributed continuous-time convex optimization on weight-balanced digraphs,” IEEE Trans. Automat. Control,
vol. 59, no. 3, pp. 781–786, 2014.
[40] W. Shi, Q. Ling, G. Wu, and W. Yin, “EXTRA: An exact first-order
algorithm for decentralized consensus optimization,” SIAM J. Optimiz.,
vol. 25, no. 2, pp. 944–966, 2015.
[41] M. Mesbahi and M. Egerstedt, Graph theoretic methods in multiagent
networks. Princeton University Press, 2010.

