arXiv:1711.08020v1 [math.OC] 21 Nov 2017

First-order methods for constrained convex programming based on
linearized augmented Lagrangian function∗
Yangyang Xu†

Abstract
First-order methods have been popularly used for solving large-scale problems. However,
many existing works only consider unconstrained problems or those with simple constraint.
In this paper, we develop two first-order methods for constrained convex programs, for which
the constraint set is represented by affine equations and smooth nonlinear inequalities. Both
methods are based on the classic augmented Lagrangian function. They update the multipliers
in the same way as the augmented Lagrangian method (ALM) but employ different primal
variable updates. The first method, at each iteration, performs a single proximal gradient step
to the primal variable, and the second method is a block update version of the first one.
For the first method, we establish its global iterate convergence as well as global sublinear
and local linear convergence, and for the second method, we show a global sublinear convergence
result in expectation. Numerical experiments are carried out on the basis pursuit denoising and
a convex quadratically constrained quadratic program to show the empirical performance of the
proposed methods. Their numerical behaviors closely match the established theoretical results.
Keywords: augmented Lagrangian method (ALM), nonlinearly constrained programming,
first-order method, global convergence, iteration complexity
Mathematics Subject Classification: 90C06, 90C25, 90C30, 68W40.

1

Introduction

Recent years have witnessed the surge of first-order methods partly due to the increasingly big
data involved in modern applications. Compared to second or higher-order methods, first-order
ones only require gradient information and generally have much lower per-iteration complexity.
However, many existing works on first-order methods are about problems without constraint or
with easy-to-project constraint and/or with affine constraint.
∗
†

This work is partly supported by NSF grant DMS-1719549.
xuy21@rpi.edu. Department of Mathematical Sciences, Rensselaer Polytechnic Institute, Troy, New York.

1

In this paper, we consider the generally constrained convex programming
min f0 (x) ≡ g(x) + h(x), s.t. Ax = b, fj (x) ≤ 0, j = 1, . . . , m,
x

(1)

where g and fj for j = 1, . . . , m are convex and Lipschitz differentiable functions, and h is a proper
closed convex (possibly nondifferentiable) function. For practical efficiency of our algorithms, we
will assume h to be simple in the sense that its proximal mapping is easy to compute. However,
our convergence results do not require this assumption.
Applications that can be formulated into (1) appear in many areas including operations research,
statistics, machine learning, engineering, just to name a few. Towards finding a solution to (1), we
design algorithms that only need zeroth and first-order information of g and fj , j = 1, . . . , m, and
the proximal mapping of h.

1.1

Augmented Lagrangian method

Our algorithms are based on augmented Lagrangian function of (1). In the literature, there are
several different augmented Lagrangian functions (see [1] for example), and we use the classic one.
Let

 uv + β u2 , if βu + v ≥ 0,
2
ψβ (u, v) =
2
v
 − ,
if βu + v < 0,
2β
and
Ψβ (x, z) =

m
X

ψβ (fj (x), zj ).

j=1

Then the classic augmented Lagrangian function of (1) is
Lβ (x, y, z) = g(x) + h(x) + y > (Ax − b) +

β
kAx − bk2 + Ψβ (x, z),
2

(2)

where y and z are Lagrangian multipliers, and β > 0 is the penalty parameter.
The augmented Lagrangian method (ALM) for (1), at each iteration, renews x-variable by minimizing Lβ with respect to x while y and z are fixed and then perform an augmented dual gradient
ascent update to the multipliers y and z, namely,
xk+1 ∈ arg min Lβ (x, y k , z k ),

(3a)

y k+1 = y k + ρy ∇y Lβ (xk+1 , y k , z k ),

(3b)

z k+1 = z k + ρz ∇z Lβ (xk+1 , y k , z k ).

(3c)

x

In general, it is difficult to solve the x-subproblem exactly or to a high accuracy. In one recent
work [19], we show that if (3a) is solved to a certain error tolerance, a global sublinear convergence
of the inexact ALM can be established. In this work, we propose to perform one single proximal
gradient update to (3a), and a sublinear convergence can still be shown.
2

1.2

Related work

ALM has been popularly used to solve constrained optimization problems; see books [2,3]. However,
most works on first-order methods in the ALM framework consider affinely constrained problems,
and only a few study the methods for generally constrained problems in the form of (1). We review
these works below.
For smooth affinely constrained convex programs, [10] analyzes the iteration complexity of an
inexact ALM, where each primal subproblem is approximately solved by Nesterov’s optimal first7
order method [14]. It shows that to reach an ε-optimal solution (see Definition 1.1 below), O(ε− 4 )
gradient evaluations are sufficient. In addition, it shows that O(ε−1 | log ε|) gradient evaluations can
guarantee an ε-optimal solution by an inexact proximal ALM. Although the number of gradient
evaluations is not explicitly given, [11, 12] also consider inexact ALM. They specify the accuracy
that each primal subproblem need be solved to and estimate the outer iteration complexity of the
inexact ALM. Within the ALM framework, [17] perform a single proximal gradient update to primal
variable at each iteration and establish O(ε−1 ) complexity result to have an ε-optimal solution for
affinely constrained composite convex programs. This linearized ALM also appears as a special case
of the methods in [5–7, 9, 20], which perform Gauss-Seidel or randomized block coordinate update
to the primal variable in the ALM framework.
Towards finding solutions of general saddle-point problems, [13] gives a subgradient method. If
√
both primal and dual constraint sets are compact, the method has O(1/ k) convergence rate in
terms of primal-dual gap, where k is the number of iterations. It also discusses how to apply
the subgradient method to convex optimization problems with nonlinear inequality constraint.
On smooth constrained convex problems, [21] proposes a primal-dual type first-order method (see
(69) in section 6.2). Assuming compactness of the constraint set, it establishes O(ε−1 ) iteration
complexity result to produce an ε-optimal solution. Recently, [19] studies an inexact ALM for
(1) and proposes to use Nesterov’s optimal first-order method to approximately solve each x3
subproblem. When the constraint set is bounded, it shows that nearly O(ε− 2 ) gradient evaluations
suffice to obtain an ε-optimal solution, and for the smooth case, the result can be improved to
O(ε−1 | log ε|). Compared to these works, our iteration complexity results will be better under
weaker assumptions.

1.3

Contributions

This paper mainly makes the following contributions.
• We propose a first-order method, named LALM, for solving composite convex problems with
both affine equality and smooth nonlinear inequality constraints. The method is based on
proximal linearization of the classic augmented Lagrangian function. Under mild assumptions,
we show global iterate sequence convergence of LALM to a primal-dual optimal solution.

3

• Also, we analyze the iteration complexity of the proposed method. We show that to reach
an ε-optimal solution, O(ε−1 ) gradient evaluations are sufficient. In addition, we establish its
local linear convergence by assuming the existence of a non-degenerate primal-dual solution
and positive definiteness of Hessian of the augmented Lagrangian function near the nondegenerate primal-dual solution.
• Furthermore, as the problem has the so-called coordinate friendly structure, we propose a
block update version of LALM. At each iteration, the method renews a single block coordinate
while keeping all the other coordinates unchanged and then immediately performs an update
to dual variables. We show that in expectation, an ε-optimal solution can be obtained by
O(ε−1 ) gradient evaluations.
• We implement LALM and its block update version and apply them to the basis pursuit
denoising problem and a quadratically constrained quadratic program. On both problems,
we notice better performance of the block-LALM in terms of iteration number. In addition,
when the iterate is far away from optimality, sublinear convergence is observed, and while
the iterate approaches to optimality, both methods converge linearly.

1.4

Notation and organization

We focus on finite-dimensional Euclidean space, but our analysis can be directly extended to a
Hilbert space. We use [m] as the set {1, 2, . . . , m}, and [a]+ = max(0, a) denotes the positive part
of a real number a. We use I √
as the identity matrix. Given a symmetric positive definite (SPD)
matrix P , we define kxkP = x> P x, and if P = I, we simply write it as kxk. Also, given a
P
nonnegative vector ` = [`1 , . . . , `n ] ∈ Rn , we define kxk2` = ni=1 `i kxi k2 if x is partitioned into
˜ (x) as its subgradient and ∂f (x)
n blocks (x1 , . . . , xn ). For any convex function f (x), we use ∇f
˜ (x)
the subdifferential of f at x, i.e., the set of all subgradients at x. When f is differentiable, ∇f
coincides with the gradient of f , and we simply write it to ∇f (x). The indicator function of a set
X is defined as ιX (x) = 0 if x ∈ X and +∞ otherwise. Bγ (x) represents a ball with radius γ and
center x. Eik denotes the expectation about ik conditioned on all previous history.
For ease of notation, we use w as the triple (x, y, z) and denote the smooth part of Lβ as
Fβ (w) = Lβ (w) − h(x).
In addition, we define
Φ(x̄; x, y, z) = f0 (x̄) − f0 (x) + y > (Ax̄ − b) +

m
X

zj fj (x̄).

(4)

j=1

Definition 1.1 (ε-optimal solution) Let f0∗ be the optimal value of (1). We call x̄ an ε-optimal
solution to (1) if
m
X
|f0 (x̄) − f0∗ | ≤ ε, kAx̄ − bk +
[fj (x̄)]+ ≤ ε.
j=1

4

Organization. The rest of the paper is organized as follows. Section 2 gives several technical
results that will be used to prove our main theorems. We propose a linearized ALM for (1) in
section 3 and a block linearized ALM in section 4. Convergence results are also given. In section
5, we discuss a few applications and how the proposed methods can be applied. Numerical results
are given in section 6, and finally section 7 concludes the paper.

2

Technical assumptions and preliminary results

A point w = (x, y, z) satisfies the Karush-Kuhn-Tucker (KKT) conditions for (1) if
0 ∈ ∇g(x) + ∂h(x) + A> y +

m
X

zj ∇fj (x),

(5a)

Ax = b,

(5b)

zj ≥ 0, fj (x) ≤ 0, zj fj (x) = 0, ∀j ∈ [m].

(5c)

j=1

If w satisfies the above conditions, we call it a KKT point. For convex programs, the conditions
in (5) are sufficient for x to be an optimal solution of (1). If a certain qualification condition (e.g.,
the Slater condition) holds, they are also necessary.

2.1

Technical assumptions

Throughout the paper, we assume the existence of a KKT point.
Assumption 1 There exists a point w∗ = (x∗ , y ∗ , z ∗ ) satisfying the KKT conditions in (5).
Under the above assumption, it follows from the convexity of f0 that
Φ(x; w∗ ) ≥ 0, ∀x,

(6)

where Φ is defined in (4).
In addition, we make the following assumption, which holds if dom(h) is bounded.
Assumption 2 There are constants Lg , L1 , . . . , Lm and B1 , . . . , Bm such that
k∇g(x̂) − ∇g(x̃)k ≤ Lg kx̂ − x̃k, ∀x̂, x̃ ∈ dom(h),
k∇fj (x̂) − ∇fj (x̃)k ≤ Lj kx̂ − x̃k, ∀x̂, x̃ ∈ dom(h), ∀j ∈ [m],
k∇fj (x)k ≤ Bj , ∀x ∈ dom(h), ∀j ∈ [m].

(7)
(8)
(9)

From the mid-point theorem, the boundedness of ∇fj implies the Lipschitz continuity of fj , i.e.,
|fj (x̂) − fj (x̃)| ≤ Bj kx̂ − x̃k, ∀x̂, x̃ ∈ dom(h), ∀j ∈ [m].
5

(10)

2.2

Preparatory lemmas

In this subsection, we give several lemmas that will be used multiple times in our convergence
analysis. First we show the Lipschitz continuity of ∇x Ψ(w) with respect to x.
Lemma 2.1 Under Assumption 2, we have
k∇x Ψ(x̂, z) − ∇x Ψ(x, z)k ≤ LΨ (x, z)kx̂ − xk, ∀x̂, x, z,
where
LΨ (x, z) =

m 
X


 
βBj2 + Lj βfj (x) + zj +

(11)

(12)

j=1

Proof. First we notice that

∂
∂u ψβ (u, v)

= [βu + v]+ , and thus for any v,

∂
∂
ψβ (û, v) −
ψβ (ũ, v) ≤ β|û − ũ|, ∀û, ũ.
∂u
∂u
Let hj (x, zj ) = ψβ (fj (x), zj ), j = 1, . . . , m. Then
k∇x hj (x̂, zj ) − ∇x hj (x, zj )k
∂
∂
=
ψβ (fj (x̂), zj )∇fj (x̂) −
ψβ (fj (x), zj )∇fj (x)
∂u
∂u
∂
∂
≤
ψβ (fj (x̂), zj )∇fj (x̂) −
ψβ (fj (x), zj )∇fj (x̂)
∂u
∂u
∂
∂
+
ψβ (fj (x), zj )∇fj (x̂) −
ψβ (fj (x), zj )∇fj (x)
∂u
∂u
∂
≤β|fj (x̂) − fj (x)| · k∇fj (x̂)k +
ψβ (fj (x), zj ) · k∇fj (x̂) − ∇fj (x)k
∂u


≤βBj2 kx̂ − xk + Lj βfj (x) + zj + · kx̂ − xk.

(13)

Hence,
k∇x Ψβ (x̂, z) − ∇x Ψβ (x, z)k ≤

m
X

k∇x hj (x̂, zj ) − ∇x hj (x, zj )k ≤ LΨ (x, z)kx̂ − xk,

j=1

which completes the proof.



Remark 2.1 Note that the Lipschitz constant LΨ (x, z) in (11) depends on the point (x, z) and
is not a universal constant. We will set its value at the iterate of the algorithm. Together with
the next lemma, the inequality in (11) implies that a sufficient progress can be obtained after each
x-update.
Lemma 2.2 For a continuously differentiable function φ(u) and a given v, if k∇φ(u) − ∇φ(v)k ≤
Lφ (v)ku − vk, ∀u, then
φ(u) ≤ φ(v) + h∇φ(v), u − vi +
6

Lφ (v)
ku − vk2 .
2

The following result is easy to show (c.f., [4, Prop. 2.3]). It will be used for establishing iterate
convergence of the proposed algorithm.
Lemma 2.3 Let {P k } be a sequence of SPD matrices, and there are SPD matrices P and P such
that P  P k  P . Let W be a nonempty set. If the sequence {wk } satisfies
kwk+1 − wk2P k+1 ≤ kwk − wk2P k , ∀w ∈ W,
and {wk } has a cluster point w̄ in W, then wk converges to w̄.
The result below will be used to establish convergence rate of our algorithms. It is similar to a
deterministic result in [19] and can be shown in the same way. We omit its proof.
Lemma 2.4 Assume (x∗ , y ∗ , z ∗ ) is a KKT point of (1). Let x̄ be a stochastic point such that for
any y and any z ≥ 0,
E[Φ(x̄; x∗ , y, z)] ≤ α + c1 kyk2 + c2 kzk2 ,
(14)
where α and c1 , c2 are nonnegative constants independent of y and z. Then


m
X
− α + 4c1 ky ∗ k2 + 4c2
(zj∗ )2  ≤ E[f0 (x̄) − f0 (x∗ )] ≤ α,

(15)

j=1

EkAx̄ − bk +

m
X

m
X
2
2
E[fj (x̄)]+ ≤ α + c1 1 + ky k + c2
1 + zj∗ .
∗

j=1

3

(16)

j=1

Linearized augmented Lagrangian method

In this section, we propose a linearized augmented Lagrangian method (LALM). Different from
the step in (3a), it updates x-variable by a single proximal gradient descent of the augmented
Lagrangian function. The method is summarized in Algorithm 1, where δ ≥ 0 is a constant and
LkF = Lg + βkAk2 + LΨ (xk , z k )
with LΨ defined in (11).
Note that the setting of η k is for simplicity of our analysis. Practically, one can choose it by starting
from η k−1 and then backtracking such that
Fβ (xk+1 , y k , z k ) ≤ Fβ (wk ) + ∇x Fβ (wk ), xk+1 − xk +

η k k+1
kx
− xk k2 ,
2

(17)

and all our convergence results can still be shown. When η k ≥ LkF , the above inequality always
holds from Lemma 2.2.

7

Algorithm 1: Linearized augmented Lagrangian method (LALM) for (1)
1
2
3
4

Initialization: choose x0 , y 0 , z 0 and β, ρy , ρz , δ ≥ 0; set η −1 = 0
for k = 0, 1, . . . do
Let η k = max(η k−1 , LkF + δ)
Perform the updates
xk+1 = arg min h(x) + ∇x Fβ (wk ), x +
x

ηk
kx − xk k2 ,
2

y k+1 = y k + ρy (Axk+1 − b), ,
zjk+1 = zjk + ρz · max −

3.1

zjk
β

(18a)
(18b)

!
, fj (xk+1 ) , j = 1, . . . , m.

(18c)

Global convergence analysis

To show the convergence results of Algorithm 1, we need the following two lemmas, which can be
found in [19].
Lemma 3.1 Let y and z be updated by (18b) and (18c) respectively. Then for any k, it holds

1  k+1
ky
− yk2 − ky k − yk2 + ky k+1 − y k k2 − hy k+1 − y, rk+1 i = 0, (19)
2ρy
m
zjk
 X

1  k+1
kz
− zk2 − kz k − zk2 + kz k+1 − z k k2 −
(zjk+1 − zj ) · max − , fj (xk+1 ) = 0, (20)
2ρz
β
j=1

where rk = Axk − b.
Lemma 3.2 For any z ≥ 0, we have
m
m
X
X
zjk

β − 2ρz k+1
k+1
k 2
k+1 k
k+1
kz
−
z
k
≤
Ψ
(x
,
z
)
−
z
f
(x
)
−
(z
−
z
)
·
max
−
, fj (xk+1 ) .
j
j
j
β
j
2
2ρz
β
j=1

j=1

(21)
Using the above two lemmas, we establish a fundamental result on Algorithm 1.
Theorem 3.1 (One-iteration progress of LALM) Let {wk } be the sequence generated from
Algorithm 1. Then for any x such that Ax = b and fj (x) ≤ 0, ∀j ∈ [m], any y, and any z ≥ 0, it

8

holds that
1
1
η k k+1
kx
− xk2 +
ky k+1 − yk2 +
kz k+1 − zk2
2
2ρy
2ρz
β − ρy k+1 2 β − ρz k+1
1
+
kr k +
kz
− z k k2 + kxk+1 − xk k2(ηk −Lg −Lk )I−βA> A
2
Ψ
2
2ρz
2
k
η
1
β
1
≤ kxk − xk2 − krk k2 +
ky k − yk2 +
kz k − zk2 ,
2
2
2ρy
2ρz
Φ(xk+1 ; w) +

(22)

where Φ is defined in (4), and LkΨ = LΨ (xk , z k ) with LΨ defined in (12).
Proof. From the update in (18a), it follows that
0 ∈ ∂h(xk+1 ) + ∇g(xk ) + A> y k + βA> rk + ∇x Ψ(xk , z k ) + η k (xk+1 − xk ).
By the convexity of h, we have
D
E
k+1
˜
xk+1 − x, ∇h(x
) ≥ h(xk+1 ) − h(x),

(23)

(24)

From the convexity of g and Ψ(·, z), and also Lemmas 2.1 and 2.2, we have
D
E
xk+1 − x, ∇g(xk ) + ∇x Ψ(xk , z k )
D
E D
E
= xk+1 − xk , ∇g(xk ) + ∇x Ψ(xk , z k ) + xk − x, ∇g(xk ) + ∇x Ψ(xk , z k )
≥ g(xk+1 ) + Ψ(xk+1 , z k ) − g(xk ) − Ψ(xk , z k ) −

Lg + LkΨ k+1
kx
− xk k2
2

+ g(xk ) + Ψ(xk , z k ) − g(x) − Ψ(x, z k )
= g(xk+1 ) + Ψ(xk+1 , z k ) − g(x) − Ψ(x, z k ) −

Lg + LkΨ k+1
kx
− xk k2 .
2

(25)

For x such that Ax = b, it holds that
hxk+1 − x, A> y k + βA> rk i
= hrk+1 , y k+1 − ρy rk+1 + βrk i
= y > rk+1 + hy k+1 − y, rk+1 i + (β − ρy )krk+1 k2 + βhrk+1 , rk − rk+1 i
i
β h k+1 2
= y > rk+1 + hy k+1 − y, rk+1 i + (β − ρy )krk+1 k2 −
kr k − krk k2 + krk+1 − rk k2 . (26)
2
Adding (24), (25), (26), and the following equation
i
η k h k+1
kx
− xk2 − kxk − xk2 + kxk+1 − xk k2 ,
xk+1 − x, η k (xk+1 − xk ) =
2
we have from (23) that
f0 (xk+1 ) − f0 (x) + Ψ(xk+1 , z k ) − Ψ(x, z k ) + y > rk+1 + hy k+1 − y, rk+1 i
i
η k − Lg − LkΨ k+1
ηk
β h k+1 2
kr k + krk+1 − rk k2 + (β − ρy )krk+1 k2 + kxk+1 − xk2 +
kx
− xk k2
−
2
2
2
ηk
β
≤ kxk − xk2 − krk k2 .
(27)
2
2
9

The desired result is obtained by noting Ψ(x, z k ) ≤ 0, adding (19), (20), and (21) to the above
inequality, and rearranging terms.

The next lemma shows the upper boundedness of η k .
Lemma 3.3 Let {wk } be the sequence generated from Algorithm 1 with zj0 ≥ 0, ∀j ∈ [m]. If
ρy , ρz ∈ (0, β], then η k ≤ η̄, ∀k ≥ 0, where η̄ is a constant satisfying
!
m
m
0 − y∗k
0 − z∗k
X
X
ky
kz
η̄ ≥ δ + Lg + βkAk2 + β
Bj2 + β
+ p
Bj Lj kx0 − x∗ k + p
ρy η 0
ρz η 0
j=1
j=1
v
(28)
!
s
r
uX
 0
um 2 √
ρ
η̄
η̄
z
0
∗
0
∗
∗
∗
+t
ρz η̄kx − x k +
kz − z k .
Lj
ky − y k + kz k + max 1,
ρy η 0
η0
j=1

Proof. Since zj0 ≥ 0, ∀j ∈ [m], we have zjk ≥ 0, ∀j ∈ [m] from the update of z and the condition
ρz ∈ (0, β]. Note fj (xk ) ≤ fj (x∗ ) + Bj kxk − x∗ k ≤ Bj kxk − x∗ k. It follows from the increasing
monotonicity of [a]+ that
m
X

Lj [βfj (xk ) + zjk ]+ ≤

j=1

m
X


βBj Lj kxk − x∗ k + Lj zjk ,

j=1

and thus
LkΨ ≤ β

m
X
j=1

≤β

m
X
j=1

Bj2 +

m
X

βBj Lj kxk − x∗ k + Lj zjk



j=1

Bj2 + β

m
X

v
uX
um 2 ∗
k
∗
Lj (kz k + kz k − z ∗ k).
Bj Lj kx − x k + t

(29)

j=1

j=1

We next show the desired result by induction. First, the result for k = 0 directly follows from (28)
and (29). Assume η k ≤ η̄, ∀k ≤ K − 1. Then letting w = w∗ in (22), we have from (6) and by
dropping nonnegative terms on the left hand side that
η k k+1
1
1
kx
− x∗ k2 +
ky k+1 − y ∗ k2 +
kz k+1 − z ∗ k2
2
2ρy
2ρz
≤

ηk k
1
1
kx − x∗ k2 +
ky k − y ∗ k2 +
kz k − z ∗ k2 .
2
2ρy
2ρz

(30)

Since η k+1 ≥ η k , dividing by η k on both sides of the above inequality yields
1 k+1
1
1
kx
− x∗ k2 +
ky k+1 − y ∗ k2 +
kz k+1 − z ∗ k2
k+1
2
2ρy η
2ρz η k+1
1
1
1
≤ kxk − x∗ k2 +
ky k − y ∗ k2 +
kz k − z ∗ k2 .
k
2
2ρy η
2ρz η k
10

(31)

Repeatedly using (31) and also from (30), we have
1 K
1
1
ky K − y ∗ k2 +
kz K − z ∗ k2
kx − x∗ k2 +
K−1
2
2ρy η
2ρz η K−1
1
1
1
≤ kx0 − x∗ k2 +
ky 0 − y ∗ k2 +
kz 0 − z ∗ k2 .
0
2
2ρy η
2ρz η 0
The above inequality together with η K−1 ≤ η̄ implies
ky 0 − y ∗ k kz 0 − z ∗ k
kxK − x∗ k ≤ kx0 − x∗ k + p
+ p
ρy η 0
ρz η 0
and
kz

K

∗

−z k≤

√

s
0

∗

ρz η̄kx − x k +

ρz η̄ 0
ky − y ∗ k +
ρy η 0

r

η̄ 0
kz − z ∗ k.
η0

Hence, η̄ ≥ LK
F + δ from (28), (29), and the above two inequalities. This completes the proof.



We are now ready to show our main convergence and rate results.
Theorem 3.2 (Iterate convergence of LALM) Under Assumptions 1 and 2, let {wk } be the
sequence generated from Algorithm 1 with any x0 , y 0 , and zj0 ≥ 0, ∀j ∈ [m]. If ρy , ρz ∈ (0, β) and
δ > 0, then wk converges to a KKT point w̄ = (x̄, ȳ, z̄) of (1).
Proof. Letting w = w∗ in (22) and dividing by η k , we have from (6) that
1
1
1 k+1
kx
− x∗ k2 +
ky k+1 − y ∗ k2 +
kz k+1 − z ∗ k2
k
2
2ρy η
2ρz η k
β − ρy k+1 2
1
1
+
kr k + 2 k (β − ρz )kz k+1 − z k k2 + k kxk+1 − xk k2(ηk −Lg −Lk )I−βA> A
Ψ
2η k
2ρz η
2η
1
1
1
≤ kxk − x∗ k2 +
ky k − y ∗ k2 +
kz k − z ∗ k2 .
(32)
2
2ρy η k
2ρz η k
Summing up (32) over k and noting η k+1 ≥ η k , we have from the condition δ > 0, ρy , ρz ∈ (0, β)
and Lemma 3.3 that
lim xk+1 − xk = 0,

k→∞

lim y k+1 − y k = lim ρy rk+1 = 0,

k→∞

k→∞

lim z k+1 − z k = 0.

k→∞

(33)

In addition, it follows from (32) that {wk } is bounded and must have a cluster point w̄. Hence,
there is a subsequence {wk }k∈K convergent to w̄. Since {η k } is increasing and bounded, it must
converge to a number η ∞ .
Below we show that w̄ is a KKT point. First, we have Ax̄ − b = 0 from (33), i.e, x̄ satisfies (5b).
Secondly, from the update of z and ρz < β, it follows zjk ≥ 0, ∀j ∈ [m], ∀k, and thus z̄j ≥ 0, ∀j ∈ [m].
If fj (xk+1 ) > 0, then fj (xk+1 ) = ρ1z (zjk+1 − zjk ) → 0 that indicates [fj (xk+1 )]+ → 0. Hence,
11

z̄

fj (x̄) ≤ 0, ∀j ∈ [m] follows from the continuity of fj ’s. For any j ∈ [m], if z̄j > 0, then zjk > 2j ,

zk
as k ∈ K is sufficiently large. It follows from max − βj , fj (xk+1 ) → 0 that fj (xk+1 ) → 0 as
K 3 k → ∞. Hence, fj (x̄) = 0. Therefore, (x̄, z̄) satisfies (5c).
Thirdly, from the optimality of xk+1 , it holds that
∇x F (wk ), xk+1 + h(xk+1 ) +

η k k+1
ηk
kx
− xk k2 ≤ ∇x F (wk ), x + h(x) + kx − xk k2 , ∀x.
2
2

Taking limit infimum over k ∈ K on both sides of the above equation, we have from the lower
semicontinuity of h and continuity of g and fj ’s that
∇x F (w̄), x̄ + h(x̄) ≤ ∇x F (w̄), x + h(x) +
namely,
x̄ = arg min ∇x F (w̄), x + h(x) +
x

η∞
kx − x̄k2 , ∀x,
2

η∞
kx − x̄k2 .
2

Therefore, w̄ satisfies (5a) from the optimality condition of the above minimization problem, and
thus w̄ is a KKT point of (1).
Hence, (31) holds with w∗ replaced by w̄, and thus wk converges to w̄ from Lemma 2.3.



Theorem 3.3 (Sublinear convergence rate of LALM) Under Assumptions 1 and 2, let {wk }
be the sequence generated from Algorithm 1 with any x0 , y 0 = 0 and z 0 = 0. If ρy , ρz ∈ (0, β], then

|f0 (x̄k+1 ) − f0 (x∗ )| ≤

kAx̄k+1 − bk +

m
X
[fj (x̄k+1 )]+ ≤
j=1

η∞



η 0 (k + 1)



η0
2

kx0 − x∗ k2 +



2ky ∗ k2
ρy

m
X
2(zj∗ )2

+

j=1

ρz

m


,
(34a)


∗ 2
0
X (1 + zj∗ )2
η∞
 η kx0 − x∗ k2 + (1 + ky k) +
,
η 0 (k + 1) 2
2ρy
2ρz
j=1

(34b)
where η ∞ ≤ η̄ is the limit of {η k } and x̄k+1 =

xt+1
1
t=0 Pk

Pk

.

t=0 η t

Proof. Letting x = x∗ in (22), dividing by η k , and summing it up give
!


k
k
X
X
1
1
1
1
1 η0 0
∗ 2
2
2
k+1 ∗
t+1 ∗
Φ(x̄ ; x , y, z) ≤
Φ(x ; x , y, z) ≤ 0
kx − x k +
kyk +
kzk .
ηt
ηt
η
2
2ρy
2ρz
t=0

Since
c1 =

Pk

t=0

1
k+1
t=0 η t ≥ η ∞ , the desired results
∞
η∞
, and c2 = 2ρz ηη0 (k+1) .
2ρy η 0 (k+1)

are obtained from Lemma 2.4 with α =

η∞
0
2(k+1) kx

− x∗ k2 ,


12

Theorem 3.3 implies that to reach an ε-optimal solution of (1), it is sufficient to evaluate the
gradients of g and fj , j ∈ [m] and proximal mapping of h for K times, where



2 

2

 ∞
m
∗
∗
∗
∗
0
X

η
max(1 + zj , 2zj )
max(1 + ky k, 2ky k)
η
0
∗ 2

.


+
K=
kx − x k +
0
εη
2
2ρy
2ρz
j=1

3.2

Local linear convergence of LALM for constrained smooth problems

In this subsection, we assume that h(x) = ιX (x) for a closed convex set X and g, f1 , . . . , fm are twice
continuously differentiable. We show local linear convergence of Algorithm 1 under the following
assumption.
Assumption 3 There is a KKT point w∗ and a subset J ⊂ [m] such that x∗ ∈ int(X ), and
1. fj (x∗ ) = 0, zj∗ > 0, ∀j ∈ J and fj (x∗ ) < 0, zj∗ = 0, ∀j 6∈ J;


P
2. x> ∇2 g(x∗ ) + j∈J zj∗ ∇2 fj (x∗ ) x > 0 for any nonzero vector x ∈ Null(D> ), where


D = A> , ∇f1 (x∗ ), . . . , ∇fm (x∗ )
is column full-rank.
When item 1 holds in the above assumption, we have

X
∇2x Lβ (w∗ ) = ∇2 g(x∗ ) + βA> A +
zj∗ ∇2 fj (x∗ ) + β∇fj (x∗ )[∇fj (x∗ )]> ,
j∈J

and thus if in addition item 2 holds, then ∇2x Lβ (w∗ ) is positive definite. We denote µ > 0 as its
smallest eigenvalue.
From the continuity of ∇2x Lβ , we have the following result.
Proposition 3.1 There is γ > 0 such that if max(kx − x∗ k, ky − y ∗ k, kz − z ∗ k) ≤ γ, then
x ∈ int(X ); ∇2x Lβ (w) 

µ
I; βfj (x) + zj > 0, ∀j ∈ J; βfj (x) + zj < 0, ∀j 6∈ J.
2

(35)

Hence, for any x ∈ Bγ (x∗ ),
f0 (x) − f0 (x∗ ) + hy ∗ , Ax − bi +

β
µ
kAx − bk2 + Ψ(x, z ∗ ) − Ψ(x∗ , z ∗ ) ≥ kx − x∗ k2 .
2
4

The next lemma can be easily verified from the definition of Ψ. We omit its proof.

13

(36)

Lemma 3.4 If xk+1 ∈ Bγ (x∗ ) and z k ∈ Bγ (z ∗ ), then
X
(zjk − zj∗ )fj (xk+1 ) = Ψ(xk+1 , z k ) − Ψ(xk+1 , z ∗ ) − Ψ(x∗ , z k ) + Ψ(x∗ , z ∗ ).

(37)

j∈J

From the update rule of z, we have following result.
Lemma 3.5 If xk+1 ∈ Bγ (x∗ ) and z k ∈ Bγ (z ∗ ), then
m
X

(zjk+1

−

zj∗ )

j=1

X
zjk

1
(zjk − zj∗ )fj (xk+1 ).
· max − , fj (xk+1 ) ≤ kz k+1 − z k k2 +
β
ρz

(38)

j∈J

Proof. When xk+1 ∈ Bγ (x∗ ) and z k ∈ Bγ (z ∗ ), it follows from (35) that
βfj (xk+1 ) + zjk > 0, ∀j ∈ J; βfj (xk+1 ) + zjk < 0, ∀j 6∈ J.
Hence,
m
X

(zjk+1 − zj∗ ) · max −

j=1

=

X

zjk

, fj (xk+1
β

(zjk + ρz fj (xk+1 ) − zj∗ )fj (xk+1 ) +

X

(1 −

j6∈J

j∈J

≤

X

(zjk + ρz fj (xk+1 ) − zj∗ )fj (xk+1 ) + ρz

zjk 
ρz k 
)zj −
β
β

X zjk 2
β
j6∈J

j∈J

X
1
(zjk − zj∗ )fj (xk+1 ),
= kz k+1 − z k k2 +
ρz
j∈J

which completes the proof.
Let η = Lg + βkAk2 +

Pm

2
j=1 βBj .


Then we have the next theorem.

Theorem 3.4 Let {wk } be the sequence generated from Algorithm 1 with w0 satisfying the following
condition:
1
1
1
1 
kx0 − x∗ k2 +
ky 0 − y ∗ k2 +
kz 0 − z ∗ k2 ≤ γ 2 · min 1,
,
,
(39)
ρy η
ρz η
ρy η̄ ρz η̄
where γ is given in Proposition 3.1. For any θ ∈ (0, 1), if 0 < ρy ≤ β and 0 < ρz ≤ β(1 − θ), then
for any k, it holds that
θµ k+1
ηk
β
1
1
kx
− x∗ k2 + kxk+1 − x∗ k2 − krk+1 k2 +
ky k+1 − y ∗ k2 +
kz k+1 − z ∗ k2
4
2
2
2ρy
2ρz


1 k+1
θ  ρy
k 2
+ kx
− x k(ηk −Lg −Lk )I−βA> A + β 1 −
−
krk+1 k2
Ψ
2
2
2
β
1
1
ηk
≤ kxk − x∗ k2 − krk k2 +
ky k − y ∗ k2 +
kz k − z ∗ k2 .
2
2
2ρy
2ρz
14

(40)

Proof. We first note max(kxk − x∗ k, ky k − y ∗ k, kz k − z ∗ k) ≤ γ, ∀k ≥ 0 from (31), (39), and the
following inequality
1
1  k
(kx − x∗ k2 + ky k − y ∗ k2 + kz k − z ∗ k2 )
min 1,
,
ρy η̄ ρz η̄
1
1
ky k − y ∗ k2 +
kz k − z ∗ k2
≤ kxk − x∗ k2 +
k
ρy η
ρz η k
1
1
≤ kx0 − x∗ k2 +
ky 0 − y ∗ k2 +
kz 0 − z ∗ k2 .
ρy η
ρz η
Adding (19) with y = y ∗ , (20) with z = z ∗ , θ times of (37) and (38), and 1 − θ times of (21) to
(27) with (x, y) = (x∗ , y ∗ ), we have by rearranging terms that


β k+1 2
k+1
∗
∗ k+1
k+1 ∗
∗ ∗
θ f0 (x ) − f0 (x ) + hy , r i + kr k + Ψ(x , z ) − Ψ(x , z )
2


m
X
+ (1 − θ) f0 (xk+1 ) − f0 (x∗ ) + hy ∗ , rk+1 i +
zj∗ fj (xk+1 ) − Ψ(x∗ , z k )
j=1

ηk

β k+1 2
1
1
kr k +
ky k+1 − y ∗ k2 +
kz k+1 − z ∗ k2
2
2
2ρy
2ρz




1 k+1
θ  ρy
1
β(1 − θ)
k 2
k+1 2
+ kx
− x k(ηk −Lg −Lk )I−βA> A + β 1 −
−
kr k +
− 1 kz k+1 − z k k2
Ψ
2
2
2
2ρz
ρz
ηk
β
1
1
≤ kxk − xk2 − krk k2 +
ky k − y ∗ k2 +
kz k − z ∗ k2 .
2
2
2ρy
2ρz
+

kxk+1 − x∗ k2 −

From (6), (36), the above inequality, and Ψ(x∗ , z k ) ≤ 0, the desired result follows.



In addition, we can bound ky k − y ∗ k2 and kz k − z ∗ k2 by x-terms.
Lemma 3.6 Let ν > 0 be the smallest eigenvalue of D> D. Under the assumption of Theorem 3.4,
we have

ν ky k − y ∗ k2 + kz k − z ∗ k2


X
(41)
≤ 4L2g + 8|J|
(β 2 Bj4 + |zjk |2 L2j ) kxk − x∗ k2 + 4β 2 kAk2 krk k2 + 4η̄ 2 kxk+1 − xk k2 .
j∈J

Proof. Note that
∇x Ψ(xk , z k ) =

m
X
X
[βfj (xk ) + zjk ]+ ∇fj (xk ) =
(βfj (xk ) + zjk )∇fj (xk ).
j=1

j∈J

Hence, from the update of x and the fact xk+1 ∈ int(X ), it follows
X
∇g(xk ) + A> y k + βA> rk +
(βfj (xk ) + zjk )∇fj (xk ) + η k (xk+1 − xk ) = 0.
j∈J

15

In addition, since x∗ ∈ int(X ), it holds that
∇g(x∗ ) + A> y ∗ +

X

zj∗ ∇fj (x∗ ) = 0.

j∈J

From the above two equations, it follows that
2

D[y k ; z k ] − D[y ∗ ; z ∗ ]

= ∇g(xk ) − ∇g(x∗ ) + βA> rk +

X

(βfj (xk ) + zjk )∇fj (xk ) −

j∈J

X

zjk ∇fj (x∗ ) + η k (xk+1 − xk )

2

j∈J




X

≤4 k∇g(xk ) − ∇g(x∗ )k2 + kβA> rk k2 +
(βfj (xk ) + zjk )∇fj (xk ) − zjk ∇fj (x∗ )

2

+ kη k (xk+1 − xk )k2 

j∈J

≤4L2g kxk − x∗ k2 + 4β 2 kAk2 krk k2 + 4η̄ 2 kxk+1 − xk k2 + 4

X


(βfj (xk ) + zjk )∇fj (xk ) − zjk ∇fj (x∗ )

2

. (42)

j∈J

Note fj (x∗ ) = 0, ∀j ∈ J. Hence,
X

(βfj (xk ) + zjk )∇fj (xk ) − zjk ∇fj (x∗ )

2

j∈J

≤|J|

X

=|J|

X

(βfj (xk ) + zjk )∇fj (xk ) − zjk ∇fj (x∗ )

2

j∈J


β fj (xk ) − fj (x∗ ) ∇fj (xk ) + zjk ∇fj (xk ) − zjk ∇fj (x∗ )

2

j∈J

≤2|J|

X
(β 2 Bj4 + |zjk |2 L2j )kxk − x∗ k2 .
j∈J

Plugging in the above inequality into (42) and noting ν [y k ; z k ]−[y ∗ ; z ∗ ]
we obtain the desired result.

2

2

≤ D[y k ; z k ]−D[y ∗ ; z ∗ ] ,


If necessary, taking a smaller γ, we can assume
X
j∈J

Lj (βfj (x̂) + ẑj ) −

X

Lj (βfj (x̃) + z̃j ) ≤

j∈J

µ
, ∀x̂, x̃ ∈ Bγ (x∗ ), ∀ẑ, z̃ ∈ Bγ (z ∗ ).
8

(43)

Then we have the local linear convergence of Algorithm 1 as follows.
Theorem 3.5 (Local linear convergence) Under Assumptions 2 and 3, let {wk } be the sequence generated from Algorithm 1 with w0 satisfying (39), ρy = ρz = β2 , and δ > 0. Let


X
C = L2g + 2|J| 
β 2 Bj4 + 2L2max (|z ∗ |2 + γ 2 ) ,
j∈J

where Lmax = maxj Lj . For any α > 0 such that


µ δ
1
α < min
, 2,
,
8C η̄ βkAk2
16

(44)

it holds φ(xk+1 , y k+1 , z k+1 ) ≤ σ · φ(xk , y k , z k ), where
φ(xk , y k , z k ) =


µ
1
ηk  k
kx − x∗ k2 +
ky k − y ∗ k2 + kz k − z ∗ k2
+
16
2
β

and


σ = max

αβν
αC + η̄
,1 −
µ
8
8 + η̄


< 1.

Proof. Adding α8 of (41) to (40) with θ = 12 , and noting αη̄ 2 I  (η k − Lg − LkΨ )I − βA> A, ∀k and
αβ 2 kAk2 ≤ β gives
µ η k  k+1
kx
− x∗ k2 +
+
8
2
αC
ηk  k
≤
+
kx − x∗ k2 +
2
2
Let
2

ηmax = δ + Lg + βkAk +


1
ky k+1 − y ∗ k2 + kz k+1 − z ∗ k2
β

1
αν  k
−
ky − y ∗ k2 + kz k − z ∗ k2 .
β
8

m
X

βBj2 + max ∗

X

x∈Bγ (x )
z∈Bγ (z ∗ ) j∈J

j=1

(45)

Lj (βfj (x) + zj ).

From the setting of η k , we have that
η k = max(L0F , L1F , . . . , LkF ) + δ ≤ ηmax , ∀k.
In addition, (43) indicates η k+1 − η k ≤ ηmax − (LkF + δ) ≤ µ8 . Hence,
µ
η k+1
µ ηk
+
≤ + .
16
2
8
2
k

η
Since αC
2 + 2 ≤ σ
definition of φ.

µ
16

+

ηk 
2

and

1
β

−

αν
8

(46)

≤ βσ , we have the desired result from (45), (46), and the


Remark 3.1 In Theorem 3.5, the setting of ρy = ρz = β2 is for simplicity of the analysis. The local
linear convergence can be obtained for any ρy , ρz ∈ (0, β). Therefore, from Theorems 3.2 and 3.5,
the algorithm may eventually converge linearly. This phenomenon is observed from our numerical
experiments; see Figures 1 and 2.

4

Block linearized augmented Lagrangian method

In this section, we assume that in (1), x can be partitioned into n disjoint blocks and the nondifferentiable part h(x) is separable, i.e.,
x = (x1 , x2 , . . . , xn ),

h(x) =

n
X

hi (xi ).

i=1

Correspondingly, A can be written as the block matrix format [A1 , . . . , An ].
17

4.1

Algorithm

Towards a solution of the block structured problem, we propose a block linearized augmented
Lagrangian method (BLALM). At each iteration, it randomly picks one block primal variable to
update and then immediately renews the multipliers. The method is summarized in Algorithm 2.
Algorithm 2: Block linearized augmented Lagrangian method for (1)
1
2
3

Initialization: choose x0 , y 0 , z 0 and β, ρy , ρz , η = [η1 , · · · , ηn ]; let r0 = Ax0 − b
for k = 0, 1, . . . do
Pick ik ∈ [n] uniformly at random and perform the updates
(
xk+1
i

=

arg min hi (xi ) + h∇xi Fβ (wk ), xi i +
xi

ηi
2 kxi

− xki k2 ,

xki ,

if i = ik

(47a)

if i 6= ik

k
rk+1 = rk + Aik (xk+1
ik − xik ),

y k+1 = y k + ρy rk+1 ,

(47b)
!

zjk+1 = zjk + ρz · max −

zjk
, fj (xk+1 ) , j = 1, . . . , m.
β

(47c)

To make Algorithm 2 efficient, we require (1) to have the so-called coordinate friendly structure [16].
Roughly speaking, computing all n block partial gradients ∇xi Fβ has nearly the same complexity
as a full gradient evaluation. In addition, f (xk+1 ) can be easily calculated from xk , f (xk ) and the
change of xik .
We let `ki be the Lipschitz constant of ∇xi g(x)+∇xi Ψ(x, z k ) with respect to xi for every i = 1, . . . , n
and `k = [`k1 , · · · , `kn ]. In general, `ki can be significantly smaller than the Lipschitz constant of
∇g(x) + ∇x Ψ(x, z), and thus a larger stepsize can be made if a single block is updated instead of
all blocks.

4.2

Convergence analysis

To show the convergence results of Algorithm 2, we first establish a fundamental result that is
similar to Theorem 3.1.
Theorem 4.1 (One-iteration result of BLALM) Let {wk } be the sequence generated from Al-

18

gorithm 2. Then for any x such that Ax = b and fj (x) ≤ 0, ∀j ∈ [m], it holds
h
i
Eik f0 (xk+1 ) − f0 (x) + hy k+1 , rk+1 i + (β − ρy )krk+1 k2 + Ψβ (xk+1 , z k )
h
i β


1
+ Eik kxk+1 − xk2η − kxk − xk2η + kxk+1 − xk k2η−`k − Eik krk+1 k2 − krk k2 + kxk+1 − xk k2A> A
2
2

1 
k
k k
k 2
k k
≤ 1−
f0 (x ) − f0 (x) + hy , r i + βkr k + Ψβ (x , z ) .
(48)
n
Proof. From the update of xik , we have
k+1
k
k
> k
k
k k
0 ∈ ∂hik (xk+1
ik ) + ∇xik g(x ) + Aik (y + βr ) + ∇xik Ψβ (x , z ) + ηik (xik − xik ).

(49)

Note that for any x,
k+1
k+1
˜
Eik xk+1
ik − xik , ∇hik (xik ) ≥ Eik [hik (xik ) − hik (xik )]
k
k
= Eik [hik (xk+1
ik ) − hik (xik ) + hik (xik ) − hik (xik )]
1
= Eik [h(xk+1 − h(xk )] + [h(xk ) − h(x)],
n

(50)

and
k
k k
Eik xk+1
ik − xik , ∇xik g(x ) + ∇xik Ψβ (x , z )
k k
k
k
k
= Eik xk+1
ik − xik + xik − xik , ∇xik g(x ) + ∇xik Ψβ (x , z )


1
≥ Eik g(xk+1 ) + Ψβ (xk+1 , z k ) − g(x) − Ψβ (xk , z k ) − kxk+1 − xk k2`k
(51)
2

1
+ g(xk ) + Ψβ (xk , z k ) − g(x) − Ψβ (x, z k ) .
n
In addition, for any x such that Ax = b, we have from [18, Lemma 3.2] that


1 k k
> k
k
(hy , r i + βkrk k2 ) + Eik hy k , rk+1 i + βkrk+1 k2
Eik xk+1
ik − xik , Aik (y + βr ) = − 1 −
n


β
− Eik krk+1 k2 − krk k2 + kxk+1 − xk k2A> A .
(52)
2
Furthermore,
 k+1

1
k+1
k
Eik hxk+1
− xk2η − kxk − xk2η + kxk+1 − xk k2η .
(53)
ik − xik , ηik (xik − xik )i = Eik kx
2
Adding (50) through (53), we have from (49) that
h
i
Eik f0 (xk+1 ) − f0 (x) + hy k , rk+1 i + βkrk+1 k2 + Ψβ (xk+1 , z k )
h
i β


1
+ Eik kxk+1 − xk2η − kxk − xk2η + kxk+1 − xk k2η−`k − Eik krk+1 k2 − krk k2 + kxk+1 − xk k2A> A
2
2
 1
1 
k
k k
k 2
k k
≤ 1−
f0 (x ) − f0 (x) + hy , r i + βkr k + Ψβ (x , z ) + Ψβ (x, z k ).
n
n

Since y k+1 = y k + ρy rk+1 and Ψβ (x, z k ) ≤ 0, (48) is obtained from the above inequality.
We also need the next lemma.
19



Lemma 4.1 For any ρz ≤ β,
−

1 k+1
kz
− z k k2 ≤ Ψβ (xk+1 , z k ) − Ψβ (xk+1 , z k+1 )
ρz

(54)

Proof. Since ρz ≤ β, we have zjk ≥ 0, ∀k. Let
J1k = {j ∈ [m] : βfj (xk+1 ) + zjk ≥ 0, βfj (xk+1 ) + zjk+1 ≥ 0},

(55a)

J2k = {j ∈ [m] : βfj (xk+1 ) + zjk ≥ 0, βfj (xk+1 ) + zjk+1 < 0},

(55b)

J3k = {j ∈ [m] : βfj (xk+1 ) + zjk < 0}.

(55c)

For any j ∈ J1k ∪ J2k , zjk+1 = zjk + ρz fj (xk+1 ), and for any j ∈ J3k , zjk+1 = 1 −
βfj

(xk+1 )

+

zjk+1

ρz  k
β zj

and

< 0. Hence,

Ψβ (xk+1 , z k ) − Ψβ (xk+1 , z k+1 )
 X k 2
X 
(zj )
β
k+1 2
k
k+1
zj fj (x ) + [fj (x )] −
=
2
2β
k
k
k
j∈J1 ∪J2

j∈J3


X
β
k+1
k+1 2
k+1
−
zj fj (x ) + [fj (x )] +
2
k
j∈J1

=−

X

ρz [fj (x

k+1 2

)] +

j∈J1k

−

X

zjk fj (xk+1 )

j∈J2k

X
j∈J2k ∪J3k

(zjk+1 )2
2β

(zjk+1 )2
β
k+1 2
+ [fj (x )] +
2
2β

!
(56)

X (zjk )2 − (zjk+1 )2
j∈J3k

2β

For j ∈ J2k , we have
(zjk+1 )2
β
[fj (xk+1 )]2 +
2
2β
k
(zj + ρz fj (xk+1 ))2
β
= zjk fj (xk+1 ) + [fj (xk+1 )]2 +
2
2β
k
2
(zj )
ρz  k
β
ρ2 
+ 1+
+ z [fj (xk+1 )]2
=
zj fj (xk+1 ) +
2β
β
2
2β
zjk fj (xk+1 ) +

≥ − ρz [fj (xk+1 )]2 ,

(57)

where the inequality follows from the Young’s inequality. For j ∈ J3k , we have
−

(zjk )2 − (zjk+1 )2
2β



ρz 2 (zjk )2
ρz
=− 1− 1−
≥ − 2 (zjk )2 .
β
2β
β

20

(58)

Plugging (57) and (58) into (56) gives
X

Ψβ (xk+1 , z k ) − Ψβ (xk+1 , z k+1 ) ≥ −

ρz [fj (xk+1 )]2 −

j∈J1k ∪J2k

X ρz
1
(zjk )2 = − kz k+1 − z k k2 ,
2
β
ρz
k

j∈J3

which completes the proof.



The following results are easy to show from the Young’s inequality and the update rule of z.
Lemma 4.2 For any y and z ≥ 0,
hy, rk+1 i ≤ hy k , rk+1 i +
and
0 ≤ Ψβ (x

k+1

k

,z ) −

m
X

β k+1 2
1 k
kr k +
ky − yk2 ,
2
2β

(59)

1 k
kz − zk2
2β

(60)

zj fj (xk+1 ) +

j=1

Proof. The inequality in (59) directly follows from the Young’s inequality.
k = J k ∪ J k and J k = J k , where J k , J k and J k are defined in (55). Then
Let J+
−
1
2
3
1
2
3

Ψβ (x

k+1

k

,z ) −

m
X

zj fj (xk+1 ) +

j=1

1 k
kz − zk2
2β

#
 X"
(zjk )2
β
1 k
1 k
k+1 2
2
k+1
2
=
− zj )fj (x ) + [fj (x )] +
(z − zj ) +
− zj fj (x ) +
(z − zj )
−
2
2β j
2β
2β j
k
k
j∈J+
j∈J−
#
"
X
X 1
(zjk )2
1 k
≥
− zj fj (xk+1 ) +
(zj − zj )2 ≥
(zj )2 ,
−
2β
2β
2β
k
k
X

(zjk

k+1

j∈J−

j∈J−

where the first inequality follows from the Young’s inequality, and the second one holds because
zk

k and z ≥ 0, ∀j. This completes the proof.
fj (xk+1 ) ≤ − βj , ∀j ∈ J−
j



Using the previous establish results, we are now able to show the convergence rate of Algorithm 2.
Theorem 4.2 (Sublinear convergence of BLALM) Under Assumptions 1 and 2, let {wk } be
the sequence from Algorithm 2 with y 0 = 0 and z 0 = 0. Assume `ki is upper bounded by `¯i for any
β
i ∈ [n] and any k. If ρy ∈ (0, βn ], ρz ∈ (0, 2n
], and ηi ≥ `¯i + βkAi k2 , ∀i ∈ [n], then


1
2ky ∗ k2 2kz ∗ k2
k+1
∗
E[f0 (x̄ ) − f0 (x )] ≤
C x0 +
+
, (61a)
nρy
nρz
1 + nk




m
m
∗
2
X
X
2
1 
(1 + ky k)
1
E kAx̄k+1 − bk +
[fj (x̄k+1 )]+  ≤
C x0 +
+
1 + zj∗  , (61b)
k
2nρy
2nρz
1+ n
j=1

j=1

21

where x̄k+1 =

C x0

1
k
1+ n

Pk

t=0 x

t+1 ,

and




m
X
1
β
1
= 1−
f0 (x0 ) − f0 (x∗ ) + kr0 k2 +
[fj (x0 )]2+  + kx0 − x∗ k2η .
n
2
2
j=1

Proof. Since ηi ≥ `¯i + βkAi k2 , ∀i ∈ [n] and xk+1
= xki , ∀i 6= ik , it holds
i
kxk+1 − xk k2η−`k ≥ βkxk+1 − xk k2A> A .
Hence, taking expectation on both sides of (48) with x = x∗ and summing it up give
h
i β
E f0 (xk+1 ) − f0 (x∗ ) + hy k , rk+1 i + Ψβ (xk+1 , z k ) + Ekrk+1 k2
2
k−1
k−1

X
1X 
β
t+1
∗
t+1 t+1
t+1 t
+
E f0 (x ) − f0 (x ) + hy , r i + Ψβ (x , z ) +
− ρy
Ekrt+1 k2
n
n
t=0

t=0

k−1
X


 1
1
E Ψβ (xt+1 , z t ) − Ψβ (xt+1 , z t+1 ) + Ekxk+1 − x∗ k2η
n
2
t=0


1
β 0 2
1
β
0
∗
0 0
0 0
≤ 1−
f0 (x ) − f0 (x ) + hy , r i + kr k + Ψβ (x , z ) + kx0 − x∗ k2η − kr0 k2 , (62)
n
2
2
2
+ 1−

where in the first line, we have used y k+1 = y k − ρy rk+1 . Summing (19), (20), and (21) gives
# k−1
"
k−1
X
X
1
k
2
0
2
t+1
t 2
ky − yk − ky − yk +
ky
−y k −
hy t+1 − y, rt+1 i
2ρy
t=0
t=0
"
#
k−1
k−1
X
1
β − 2ρz X t+1
kz
− z t k2
+
kz k − zk2 − kz 0 − zk2 +
kz t+1 − z t k2 +
2ρz
2ρ2z
t=0
t=0


m
k−1
X
X
Ψβ (xt+1 , z t ) −
zj fj (xt+1 ) .
≤
t=0

j=1

Since y 0 = 0 and z 0 = 0, adding
β
n

≥ ρy ,

β−ρz
2nρ2z

≥

1
1− n
ρz

1
n

of the above inequality to (62), using Lemma 4.1, and noting

from the choice of ρy , ρz , we have

h
i β
1
E f0 (xk+1 ) − f0 (x∗ ) + hy k , rk+1 i + Ψβ (xk+1 , z k ) + Ekrk+1 k2 +
Eky k − yk2
2
2nρy
k−1

1
1
1X
+ Ekxk+1 − x∗ k2η +
Ekz k − zk2 +
E[Φ(xt+1 ; x∗ , y, z)]
2
2nρz
n
t=0

 1
1
β
≤ 1−
f0 (x0 ) − f0 (x∗ ) + βkr0 k2 + Ψβ (x0 , 0) + kx0 − x∗ k2η − kr0 k2
n
2
2
1
1
2
2
+
kyk +
kzk .
2nρy
2nρz
22

(63)

P
β
0 2
Note Ψβ (x0 , 0) = m
j=1 [fj (x )]+ . Since ρy , ρz ≤ n , plugging (59) and (60) into (63) and using the
convexity of fi ’s yield


1
1
1
k+1 ∗
2
2
E[Φ(x̄ ; x , y, z)] ≤
C x0 +
kyk +
kzk .
2nρy
2nρz
1 + nk
Therefore, we complete the proof by Lemma 2.4.



Remark 4.1 If n block updates of x costs roughly the same as one full update to x, then the results
in (61) are comparable to those in (34) by noting their differences in choosing ρy , ρz . One drawback
of Theorem 4.2 is the assumption on the upper bound of `k . From (13), we see that the upper
bound can be pre-calculated if fj (x), ∀j ∈ [m] are affine. However, in general, it is unknown and
dependent on the iterates. Numerically, we can gradually increase ηi by a fixed amount or ratio if
ηi < `ki + βkAi k2 is detected or by backtracking until the following inequality holds:
k
Fβ (xk+1 , y k , z k ) ≤ Fβ (wk ) + ∇xik Fβ (wk ), xk+1
ik − x ik +

ηik k+1
kxik − xkik k2 .
2

(64)

After finitely many increases, `ki + βkAi k2 ≤ ηi , ∀i, will hold in high probability for every k. This
can be explained by the following arguments.
Let ηi = ζ ≥ 1, ∀i ∈ [n]. Since nρz ≤ β2 , then from (6), (63) with (y, z) = (y ∗ , z ∗ ), (59), and (60),
it follows that
q
∗
∗
Cx0
0
∗
√ky k + √kz k ,
Ekxk − x∗ k ≤
ζ + kx − x k +
nρz ζ
nρy ζ
q
q
p
√
β
0
∗
k
∗
Ekz − z k ≤ βCx0 + βζkx − x k + nρy ky ∗ k + nρβ z kz ∗ k.
√
Hence, we have from (29) that E[`ki ] = O( ζ). By the Markov inequality, for every k, `ki ≤ ηi , ∀i ∈
[n] holds in high probability if ηi  ζ, ∀i ∈ [n], where ζ ≥ 1 satisfies


m
m
X
X
p
ky ∗ k
kz ∗ k
0
∗
2
2
Cx0 + kx − x k + √
ζ ≥ Lg + βkAk + β
Bj + β
Bj Lj
+√
nρy
nρz
j=1
j=1
v
s
s
!
uX
q
p
um 2
2β
2β
∗
∗
∗
0
∗
ky k +
kz k .
+t
Lj kz k + 2βCx0 + 2βζkx − x k +
nρy
nρz
j=1

Remark 4.2 If Assumption 3 is satisfied, we can also show a local linear convergence result of
Algorithm 2 following the analysis in section 3.2 and that in [20]. We do not expand details here
but leave it to interested readers.

5

Applications

In this section, we give a few applications that can be formulated in the form of (1) and discuss
how Algorithm 1 and/or Algorithm 2 can be applied.
23

5.1

Basis pursuit denosing

Suppose we observe a noisy measurement b = Aθo + ξ of a signal θo , where A is a measuring matrix,
and ξ is a noise vector. Assume θo can be sparsely represented by a dictionary D. Then we can
recover the signal through solving the so-called basis pursuit denoising (BPDN) problem:
min kxk1 , s.t. kADx − bk2 ≤ δ,
x

(65)

where δ measures the noise level. Upon obtaining a solution x∗ to (65), we let θr = Dx∗ be the
recovered signal. Depending on the application, one can impose certain bounds on x to make
the recovered signal physically meaningful. In this case, all conditions in Assumption 2 holds. In
addition, assuming b ∈ Range(AD), then Slater condition holds, and thus Assumption 1 is satisfied.
Hence, Algorithm 1 is applicable, and the x-subproblem (18a) has closed-form solution by shrinkage
or soft-thresholding. If A and D are stored as matrices, (65) is coordinate friendly, and we can
also use Algorithm 2. However, for certain signal processing problems, evaluating Aθ and/or Dx
may not require explicit form of A or D but can be efficiently realized, such as a partial circulant
A and/or a discrete cosine dictionary D. For this case, Algorithm 2 will not be as efficient as
Algorithm 1 since evaluating coordinate gradient of kADx − bk2 may require full gradient.

5.2

Quadratically constrained quadratic programming

The quadratically constrained quadratic programming (QCQP) can be formulated as
1 >
x Q0 x + c>
0 x + d0
2
1
s.t. x> Qj x + c>
j x + dj ≤ 0, ∀j ∈ [m],
2
Ax = b,

min

x∈Rp

(66)

li ≤ xi ≤ ui , ∀i ∈ [p].
Let X = [l1 , u1 ] × · · · × [lp , up ] and h(x) = ιX (x). Then (66) can be written as (1) by adding h(x)
into the objective. When every Qj is positive semidefinite, the problem is convex, and if all li ’s and
ui ’s are finite, then X is bounded and all conditions in Assumption 2 hold. Hence, we can apply
Algorithm 1 to find a solution of (66), and the solution of x-subproblem (18a) can be explicitly
given by performing projection to a box constraint. In addition, the problem is coordinate friendly
1
since evaluating the partial derivative of 12 x> Qj x + c>
j x + dj about each xi costs roughly p of
computing the full gradient. Furthermore, if we maintain Qj x, then calculating the function value
is negligible compared to the gradient computation. Therefore, we can also apply Algorithm 2 to
the QCQP.

24

5.3

Finite minimax problems

Many applications can be formulated as a finite minimax problem (e.g., see [15] and the references
therein):
min max fj (x),
(67)
x∈X 1≤j≤m

where each fj is a smooth convex function. Although all fj ’s are differentiable, the objective of
(67) is generally not differentiable due to the max operation. Introducing variable t and requiring
max1≤j≤m fj (x) ≤ t, one can express the minimax problem equivalently to
min t, s.t. fj (x) − t ≤ 0, ∀j ∈ [m].

x∈X ,t

(68)

For any x ∈ int(X ), each inequality constraint holds strictly at (x, maxj fj (x) + 1), and thus the
Slater condition holds. Hence, Assumption 1 is satisfied. In addition, if X is bounded, then all
conditions in Assumption 2 also hold. Therefore, we can use Algorithm 1 to find a solution of (68)
and equivalently (67), and every iteration requires performing a projection to X . Depending on
applications, one may also apply Algorithm 2 if the problem is coordinate friendly, for example,
every fj is a quadratic function.

6

Numerical experiments

In this section, we test Algorithms 1 and 2 on BPDN (65) and QCQP (66) to show their numerical performance. The two algorithms are named as LALM and BLALM respectively. For both
algorithms, we choose the parameter η by backtracking. More precisely, at each iteration k, for
LALM, we start from η k = η k−1 and multiply it by 1.5 if (17) fails, and for BLALM, we initialize
ηikk = ηik−1
and multiply it by 1.5 if (64) does not hold. For both tests, we run the compared
k
methods to 105 epochs, where one epoch is equivalent to n block updates. Optimal solutions to
both tested problems are computed by CVX [8] with high precision.

6.1

Basis pursuit denoising

In this test, we show the convergence speed of LALM and BLALM on solving BPDN (65). For
simplicity, we set D = I. The matrix A ∈ R50×100 is randomly generated according to the standard
Gaussian distribution, and the underlying sparse signal xo has 5 nonzero components following the
standard Gaussian distribution. Then we let b = Ax0 + 0.1ξ, where ξ is a unit Gaussian noise
vector. For BLALM, we evenly partition the variable x into 10 blocks. The parameter β is simply
β
set to 1 for both methods, and ρz = β is set for LALM and ρz = 10
for BLALM. Note that for the
latter, the value of ρz is larger than that given by the theorem, and the algorithm still works well.
This may indicate that our analysis is not tight.

25

10 5

Violation of constraint

distance of objective to optimal value

10 5

10 0

10 -5
LALM ergodic
BLALM ergodic
LALM nonergodic
BLALM nonergodic
1/k rate

10 -10

10 -15
0
10

10

1

2

3

10
10
10
Number of Epochs

4

10

10 0

10 -5

10 -10

10 -15
0
10

5

LALM ergodic
BLALM ergodic
LALM nonergodic
BLALM nonergodic
1/k rate

10

1

2

3

10
10
10
Number of Epochs

4

10

5

Figure 1: Convergence behaviors of Algorithm 1 (named LALM) and Algorithm 2 (named BLALM)
on the BPDN problem (65). Left: distance of objective value to the optimal value |f0 (x) − f0 (x∗ )|;
P
k
Right: constraint residual m
j=1 [fj (x)]+ . “ergodic” curves are measured by averaged iterate x̄
and “nonergodic” ones by actual iterate xk . The missing part on each constraint violation curve
corresponds to zero residual.
Figure 1 plots the objective values and constraint residuals produced by both algorithms, where the
curve corresponding to “ergodic” is obtained by using the averaged iterates x̄k and “nonergodic”
by the actual iterate xk . The missing part on each constraint violation curve corresponds to zero
residual. Since LALM and BLALM have similar per-epoch complexity, their comparison in terms
of running time is similar to that in Figure 1. From the figure, we see that BLALM is better than
LALM in terms of both ergodic and nonergodic iterates. The ergodic convergence speed of both
methods is precisely the order of k1 and matches our theorems. However, the nonergodic convergence
is significantly faster, especially as the iterate approaches to optimality. This is possibly because
the iterate enters a region where the algorithm has linear convergence as indicated by the analysis
in section 3.2. For this reason, we use the actual iterate in the next test.

6.2

Quadratically constrained quadratic programming

In this subsection, we test LALM and BLALM on the QCQP problem (66) and compare them
to the recently proposed first-order primal-dual type method by Yu&Neely [21]. They consider a
smooth constrained convex program in the form of (1) without an explicit linear equality constraint.
Their method that we name as PD-YN iteratively performs the updates:



m
X

1
xk+1 = PX xk − ∇f0 (xk ) +
λkj + fj (xk ) ∇fj (xk ) ,
(69a)
η
j=1

λk+1
= max − fj (xk ), λkj + fj (xk ) , ∀j ∈ [m],
(69b)
j
26

10 5

Violation of constraint

distance of objective to optimal value

10 5

10 0

10 -5

LALM
BLALM
PD-YN
1/k rate

10 -10

10 -15
0
10

10

1

2

3

10
10
10
Number of Epochs

4

10

10 0

10 -5

10 -15
0
10

5

LALM
BLALM
PD-YN
1/k rate

10 -10

10

1

2

3

10
10
10
Number of Epochs

4

10

5

Figure 2: Convergence behaviors of Algorithm 1 (named LALM), Algorithm 2 (named BLALM),
and the primal-dual type method (named PD-YN) in [21] on the QCQP problem (66). Left:
distance of objective value to the optimal value |f0 (x) − f0 (x∗ )|; Right: constraint residual
Pm
j=1 [fj (x)]+ . The missing part on each constraint violation curve corresponds to zero residual.
where λ0j = max(0, −fj (x0 )), ∀j ∈ [m], and η is the step size. In the test, we also choose η
adaptively by backtracking such that
φ(xk+1 , z k ) ≤ φ(xk , z k ) + ∇x φ(xk , z k ), xk+1 − xk +

η k k+1
kx
− xk k2 ,
2

P
k
k
k
where φ(x, z) = f0 (x) + m
j=1 zj fj (x) and zj = λj + fj (x ). Although [21] does not show the
convergence of PD-YN with the above adaptive η k , we observe its better performance than that
with a fixed η.
The problem size is set to m = 10 and p = 2000. We randomly generate SPD matrices Qj , j =
0, 1, . . . m. A is set to zero, i.e., there is no linear equality constraint. The vector cj ’s are generated
according to the standard Gaussian distribution, and d0 = 0 and each dj is a negative number for
j ∈ [m]. Also we set li = −10 and ui = 10 for each i ∈ [p]. Hence, the zero vector is an interior
point of X and makes every inequality hold strictly, namely, the Slater condition holds. We set
β = 0.1 for both LALM and BLALM, and for the latter, we evenly partition the variable into 200
β
blocks. The parameter ρz is set to β and 200
for the two algorithms respectively.
Figure 2 plots the results by the three compared algorithms. Both the proposed methods perform
significantly better than PD-YN, and BLALM is the best among the three. In addition, we notice
that LALM and BLALM converge linearly when the iterate approaches to optimality, as indicated
by Theorem 3.5.

27

7

Conclusions

We have presented a first-order method for solving composite convex programming with both equality and smooth nonlinear inequality constraints. The method is derived from proximal linearization
of the classic augmented Lagrangian function. Its global iterate convergence and global sublinear
and local linear convergence results have been established. For the problem that has coordinate
friendly structure, we have also proposed a first-order randomized block update method and shown
its global sublinear convergence in expectation. In addition, we have implemented the two methods
on solving the basis pursuit denoising problem and the convex quadratically constrained quadratic
programming. Global sublinear and local linear convergence are both observed in the numerical
experiments.

References
[1] A. Ben-Tal and M. Zibulevsky. Penalty/barrier multiplier methods for convex programming
problems. SIAM Journal on Optimization, 7(2):347–366, 1997. 2
[2] D. P. Bertsekas. Nonlinear programming. Athena scientific Belmont, 1999. 3
[3] D. P. Bertsekas. Constrained optimization and Lagrange multiplier methods. Academic press,
2014. 3
[4] P. L. Combettes and J.-C. Pesquet. Stochastic quasi-fejér block-coordinate fixed point iterations with random sweeping. SIAM Journal on Optimization, 25(2):1221–1248, 2015. 7
[5] Y. Cui, X. Li, D. Sun, and K.-C. Toh. On the convergence properties of a majorized alternating
direction method of multipliers for linearly constrained convex optimization problems with
coupled objective functions. Journal of Optimization Theory and Applications, 169(3):1013–
1041, 2016. 3
[6] X. Gao, Y. Xu, and S. Zhang. Randomized primal-dual proximal block coordinate updates.
arXiv preprint arXiv:1605.05969, 2016. 3
[7] X. Gao and S.-Z. Zhang. First-order algorithms for convex optimization with nonseparable objective and coupled constraints. Journal of the Operations Research Society of China,
5(2):131–159, 2017. 3
[8] M. Grant, S. Boyd, and Y. Ye. CVX: Matlab software for disciplined convex programming,
2008. 25
[9] M. Hong, T.-H. Chang, X. Wang, M. Razaviyayn, S. Ma, and Z.-Q. Luo. A block successive
upper bound minimization method of multipliers for linearly constrained convex optimization.
arXiv preprint arXiv:1401.7079, 2014. 3
28

[10] G. Lan, D. Renato, and C. Monteiro. Iteration-complexity of first-order augmented lagrangian
methods for convex programming. Mathematical Programming, 155(1-2):511–547, 2016. 3
[11] Y.-F. Liu, X. Liu, and S. Ma. On the non-ergodic convergence rate of an inexact augmented
lagrangian framework for composite convex programming. arXiv preprint arXiv:1603.05738,
2016. 3
[12] V. Nedelcu, I. Necoara, and Q. Tran-Dinh. Computational complexity of inexact gradient
augmented lagrangian methods: application to constrained mpc. SIAM Journal on Control
and Optimization, 52(5):3109–3134, 2014. 3
[13] A. Nedić and A. Ozdaglar. Subgradient methods for saddle-point problems. Journal of optimization theory and applications, 142(1):205–228, 2009. 3
[14] Y. Nesterov. Gradient methods for minimizing composite functions. Mathematical Programming, 140(1):125–161, 2013. 3
[15] E. Pee and J. O. Royset. On solving large-scale finite minimax problems using exponential
smoothing. Journal of optimization theory and applications, 148(2):390–421, 2011. 25
[16] Z. Peng, T. Wu, Y. Xu, M. Yan, and W. Yin. Coordinate-friendly structures, algorithms and
applications. Annals of Mathematical Sciences and Applications, 1(1):57–119, 2016. 18
[17] Y. Xu. Accelerated first-order primal-dual proximal methods for linearly constrained composite
convex programming. SIAM Journal on Optimization, 27(3):1459–1484, 2017. 3
[18] Y. Xu.
Asynchronous parallel primal-dual block update methods.
arXiv:1705.06391, 2017. 19

arXiv preprint

[19] Y. Xu. Global convergence rates of augmented lagrangian methods for constrained convex
programming. preprint, 2017. 2, 3, 7, 8
[20] Y. Xu and S. Zhang. Accelerated primal-dual proximal block coordinate updating methods
for constrained convex optimization. arXiv preprint arXiv:1702.05423, 2017. 3, 23
[21] H. Yu and M. J. Neely. A primal-dual type algorithm with the O(1/t) convergence rate for
large scale constrained convex programs. In Decision and Control (CDC), 2016 IEEE 55th
Conference on, pages 1900–1905. IEEE, 2016. 3, 26, 27

29

