Inexact Spectral Projected Gradient Methods on Convex Sets
Ernesto G. Birgin
∗
Jos´e Mario Mart´ınez
†
Marcos Raydan
‡
March 26, 2003
Abstract
A new method is introduced for large scale convex constrained optimization. The
general model algorithm involves, at each iteration, the approximate minimization
of a convex quadratic on the feasible set of the original problem and global conver-
gence is obta ined by means of nonmonotone line sear ches. A specific algorithm, the
Inexact Spectral Projected Gradient method (ISPG), is implemented using inexact
projections computed by Dykstra’s alternating projection method and generates inte-
rior iterates. The ISPG method is a generalization of the Spectral Projected Gradient
method (SPG), but can be used when projections are difficult to c ompute. Numerical
results for constra ined least-squares rectangular matrix problems are presented.
Key words: Convex constrained optimization, projected gradient, nonmonotone line
search, spectral gradient, Dykstra’s algorithm.
AMS Subject Classification: 49M07, 49M10, 65K, 90C06, 90C20.
1 Introduction
We consider the problem
Minimize f(x) subject to x ∈ Ω, (1)
where Ω is a closed convex set in IR
n
. Th roughout this paper we assume that f is defin ed
and has continuous partial derivatives on an open set that contains Ω.
The Spectral Projected Gradient (SPG) m ethod [6, 7] was recently proposed for solving
(1), especially for large-scale problems since the storage requirements are minimal. This
∗
Department of Computer Science, Institute of Mathematics and Statistics, University of S˜ao Paulo, Rua
do Mat˜ao 1010 Cidade Universit´aria, 05508-090 S˜ao Paulo, SP - Brazil (egbirgin@ime.usp.br). Sponsored
by FAPESP (Grants 01/04597-4 and 02/00094-0), CNPq (Grant 300151/00-4) and Pronex.
†
Departamento de Matem´atica Aplicada, IMECC-UNICAMP, CP 6065, 13081-970 Campinas SP, Brazil
(martinez@ime.unicamp.br). Sponsored by FAPESP (Grant 01/04597-4), CNPq and FAEP-UNICAMP.
‡
Departamento de Computaci´on, Facultad de Ciencias, Universidad Central de Venezuela, Ap. 47002,
Caracas 1041-A, Venezuela (mraydan@reacciun.ve). Sponsored by the Center of Scientific Computing at
UCV.
1
method has proved to be effective for very large-scale convex programming problems.
In [7] a family of location problems was described with a variable number of variables
and constraints. The SP G method was able to solve problems of this family with up to
96254 variables and up to 578648 constraints in very few seconds of compu ter time. The
computer code that implements S P G and produ ces the mentioned results is pu blished [7]
and available. More recently, in [5] an active-set method which uses SPG to leave the faces
was introduced, and bound-constrained problems with up to 10
7
variables were solved.
The SPG method is related to the practical version of Bertsekas [3] of the classical
gradient projected method of Goldstein, Levitin and Polyak [21, 25]. However, some
critical differences make this method much more efficient than its gradient-projection
predecessors. The main point is that the first trial step at each iteration is taken using
the spectral steplength (also known as the Barzilai-Borwein choice) introduced in [2] and
later analyzed in [9, 19, 27] among others. The spectral step is a Rayleigh quotient related
with an average Hessian matrix. For a review containing the more recent advances on this
special choice of steplength see [20]. The second improvement over traditional gradient
projection methods is that a nonmonotone search must be used [10, 22]. This feature seems
to be essential to preserve the n ice and nonmonotone behaviour of the iterates produced
by single spectral gradient steps.
The reported efficiency of the SPG method in very large problems motivated us to
introduce the inexact-projection version of the method. I n fact, the main drawback of the
SPG method is that it requires the exact projection of an arbitrary point of IR
n
onto Ω
at every iteration.
Projecting onto Ω is a difficult problem unless Ω is an easy set (i.e. it is easy to
project onto it) as a box, an affine subspace, a ball, etc. However, for many important
applications, Ω is not an easy set and the projection can only be achieved inexactly. For
example, if Ω is the intersection of a finite collection of closed and convex easy sets, cycles
of alternating projection methods could be used. This sequence of cycles could be stopped
prematurely leading to an inexact iterative scheme. In this work we are mainly concerned
with extending the machinery developed in [6, 7] for the more general case in which the
projection onto Ω can only be achieved inexactly.
In Section 2 we define a general model algorithm and prove global convergence. In
Section 3 we introduce the ISPG method and we describe the use of Dykstra’s alternat-
ing projection method for obtaining inexact projections onto closed and convex s ets. In
Section 4 we present numerical experiments and in Section 5 we draw some conclusions.
2 A general model algorithm and its global convergence
We say that a point x ∈ Ω is stationary, for problem (1), if
g(x)
T
d ≥ 0 (2)
for all d ∈ IR
n
such that x + d ∈ Ω.
In this work k · k denotes the 2-norm of vectors and matrices, although in some cases it
can be replaced by an arbitrary norm. We also denote g(x) = ∇f (x) and IN = {0, 1, 2, . . .}.
2
Let B be the set of n × n positive definite matrices such that kBk ≤ L and kB
−1
k ≤ L.
Therefore, B is a compact set of IR
n×n
. In the spectral gradient approach, the matrices
will be diagonal. However, the algorithm and theorem that we present below are quite
general. The matrices B
k
may be thought as defining a sequence of different metrics in IR
n
according to which we perform projections. For this reason, we give the name “Inexact
Variable Metric” to the method introduced below.
Algorithm 2.1: Inexact Variable Metric Met hod
Assume η ∈ (0, 1], γ ∈ (0, 1), 0 < σ
1
< σ
2
< 1, M a positive integer. Let x
0
∈ Ω be an
arbitrary initial point. We denote g
k
= g(x
k
) for all k ∈ IN. Given x
k
∈ Ω, B
k
∈ B, the
steps of the k−th iteration of the algorithm are:
Step 1. Compute the search direction
Consider the subproblem
Minimize Q
k
(d) subject to x
k
+ d ∈ Ω, (3)
where
Q
k
(d) =
1
2
d
T
B
k
d + g
T
k
d.
Let
¯
d
k
be the minimizer of (3). (This minimizer exists and is un ique by the strict convexity
of the subproblem (3), but we will see later that we do not need to compute it.)
Let d
k
be such that x
k
+ d
k
∈ Ω and
Q
k
(d
k
) ≤ η Q
k
(
¯
d
k
). (4)
If d
k
= 0, stop the execution of the algorithm declaring that x
k
is a stationary point.
Step 2. Compute the steplength
Set α ← 1 and f
max
= max{f(x
k−j+1
) | 1 ≤ j ≤ min{k + 1, M}}.
If
f(x
k
+ αd
k
) ≤ f
max
+ γαg
T
k
d
k
, (5)
set α
k
= α, x
k+1
= x
k
+α
k
d
k
and finish the iteration. O th erwise, choose α
new
∈ [σ
1
α, σ
2
α],
set α ← α
new
and repeat test (5).
Remark. I n the definition of Algorithm 2.1 the possibility η = 1 corresponds to th e case
in which the subproblem (3) is solved exactly.
Lemma 2.1. The algorithm is well defined.
3
Proof. Since Q
k
is strictly convex and the domain of (3) is convex, the problem (3) has a
unique solution
¯
d
k
. If
¯
d
k
= 0 then Q
k
(
¯
d
k
) = 0. Since d
k
is a feasible point of (3), and, by
(4), Q
k
(d
k
) ≤ 0, it turns out that d
k
=
¯
d
k
. Therefore, d
k
= 0 and the algorithm stops.
If
¯
d
k
6= 0, th en , since Q
k
(0) = 0 and the solution of (3) is unique, it follows that
Q
k
(
¯
d
k
) < 0. Then, by (4), Q
k
(d
k
) < 0. Since Q
k
is convex and Q
k
(0) = 0, it follows that
d
k
is a descent direction for Q
k
, therefore, g
T
k
d
k
< 0. So, for α > 0 small enough,
f(x
k
+ αd
k
) ≤ f(x
k
) + γαg
T
k
d
k
.
Therefore, the condition (5) must be satisfied if α is small enough. This completes the
proof. 2
Theorem 2.1. Assume that the level set {x ∈ Ω | f (x) ≤ f(x
0
)} is bounded. Then, ei-
ther the algorithm stops at some stationary point x
k
, or every limit point of the generated
sequence is stationary.
The proof of Theorem 2.1 is based on the following lemmas.
Lemma 2.2. Assume that the sequence generated by Algorithm 2.1 stops at x
k
. Then, x
k
is stationary.
Proof. If the algorithm stops at some x
k
, we have that d
k
= 0. Therefore, Q
k
(d
k
) = 0.
Then, by (4), Q
k
(
¯
d
k
) = 0. So,
¯
d
k
= 0. Therefore, for all d ∈ IR
n
such that x
k
+ d ∈ Ω we
have g
T
k
d ≥ 0. Thus, x
k
is a stationary point. 2
For the remaining results of this section we assume that the algorithm does not stop.
So, infinitely many iterates {x
k
}
k∈IN
are generated and, by (5), f(x
k
) ≤ f(x
0
) for all
k ∈ IN. Thus, under the hypothesis of T heorem 2.1, the s equ en ce {x
k
}
k∈IN
is bounded.
Lemma 2.3. Assume that {x
k
}
k∈IN
is a sequence generated by Algorithm 2.1. Define,
for all j = 1, 2, 3, . . .,
V
j
= max{f (x
jM−M+1
), f (x
jM−M+2
) . . . , f(x
jM
)},
and ν(j) ∈ {jM − M + 1, jM − M + 2, . . . , jM} such that
f(x
ν(j)
) = V
j
.
Then,
V
j+1
≤ V
j
+ γα
ν(j+1)−1
g
T
ν(j+1)−1
d
ν(j+1)−1
. (6)
for all j = 1, 2, 3, . . ..
Proof. We will prove by ind uction on ℓ that for all ℓ = 1, 2, . . . , M and f or all j = 1, 2, 3, . . .,
f(x
jM+ℓ
) ≤ V
j
+ γα
jM+ℓ−1
g
T
jM+ℓ−1
d
jM+ℓ−1
< V
j
. (7)
4
By (5) we have that, for all j ∈ IN,
f(x
jM+1
) ≤ V
j
+ γα
jM
g
T
jM
d
jM
< V
j
,
so (7) holds f or ℓ = 1.
Assume, as the inductive hypothesis, that
f(x
jM+ℓ
′
) ≤ V
j
+ γα
jM+ℓ
′
−1
g
T
jM+ℓ
′
−1
d
jM+ℓ
′
−1
< V
j
(8)
for ℓ
′
= 1, . . . , ℓ.
Now, by (5), and the definition of V
j
, we have that
f(x
jM+ℓ+1
) ≤ max
1≤t≤M
{f(x
jM+ℓ+1−t
} + γα
jM+ℓ
g
T
jM+ℓ
d
jM+ℓ
= max{f (x
(j−1)M +ℓ+1
), . . . , f(x
jM+ℓ
)} + γα
jM+ℓ
g
T
jM+ℓ
d
jM+ℓ
≤ max{V
j
, f(x
jM+1
), . . . , f(x
jM+ℓ
)} + γα
jM+ℓ
g
T
jM+ℓ
d
jM+ℓ
.
But, by the inductive hypothesis,
max{f(x
jM+1
), . . . , f(x
jM+ℓ
)} < V
j
,
so,
f(x
jM+ℓ+1
) ≤ V
j
+ γα
jM+ℓ
g
T
jM+ℓ
d
jM+ℓ
< V
j
.
Therefore, the inductive proof is complete and, so, (7) is proved. Since ν(j + 1) = jM + ℓ
for some ℓ ∈ {1, . . . , M}, this implies the d esired result. 2
From now on, we define
K = {ν(1) − 1, ν(2) − 1, ν(3) − 1, . . .},
where {ν(j)} is the sequence of indices defined in Lemma 2.3. Clearly,
ν(j) < ν(j + 1) ≤ ν(j) + 2M (9)
for all j = 1, 2, 3, . . ..
Lemma 2.4.
lim
k∈K
α
k
Q
k
(
¯
d
k
) = 0.
Proof. By (6), since f is continuous and bounded below,
lim
k∈K
α
k
g
T
k
d
k
= 0. (10)
But, by (4),
0 > Q
k
(d
k
) =
1
2
d
T
k
B
k
d
k
+ g
T
k
d
k
≥ g
T
k
d
k
∀ k ∈ IN.
5
