DIPARTIMENTO DI MATEMATICA
Complesso universitario, via Vivaldi - 81100 Caserta
ON MUTUAL IMPACT OF NUMERICAL LINEAR ALGEBRA
AND LARGE
-
SCALE OPTIMIZATION
WITH FOCUS ON INTERIOR POINT METHODS
Marco D’Apuzzo
1
,Valentina De Simone
1
and Daniela di Serafino
1
PREPRINT: n. 1, impaginato nel mese di marzo 2008
CLASSIFICAZIONE AMS: 65F10, 65K05, 90C30, 90C06
1
Dipartimento di Matematica, Seconda Università di Napoli, via Vivaldi 43, I-81100 Caserta, Italy
On mutual impact of numerical linear algebra
and large-scale optimization
with focus on interior point methods
Marco D’Apuzz o, Valentina De Simone, Daniela di Serafino
Department of Mathematics, Second University of Naples, Caserta, Italy
E-mail: {marco.dapuzzo,valentina.desimone,daniela.diserafino}@unina2.it
Abstract
The solution of KKT systems is ubiquitous in optimization methods and often
domi nates the computation time, especially when l arge-scale problems are consid-
ered. Thus, the effective implementation of such methods is highl y dependent on
the availability of effective linear algebra algorithms and software, that are able, in
turn, to take into account specific needs of optimization. In this paper we discuss
the mutual impact of linear algebra and optimization, focusing on interior point
methods and on the iterative solution of the KKT system. Three critical issues
are addressed: preconditioning, termination control for the inner iterations, and
inertia control.
Keywords: large-scale optimization, interior point methods, KKT system, constraint
preconditioners, adaptive stopping criteria, inertia control.
1 Introduction
The strong interplay between numerical linear algebra and optimization has been ev-
ident for a long time. Much progress in numerical linear algebra has been spurred by
the need of solving linear systems with special features in the context o f optimization,
and many optimization codes have benefited, in terms of both efficiency and robust-
ness, from advances in numerical linear algebra, coming out also fro m needs in other
fields of scientific computing. This interplay is clearly recog nized in the textbook [41]
by Gil l et al., where a presentation of numerical optim ization and numerical linear
algebra techniques is provided, highlig hting the relations between the two fields in t he
broader context of scientific computing. A general di s cussion of the r ole of numerical
linear algebra in optimization in the 20th century is in t he essay by O’Leary [65]. She
1
points out that in a ny optimization algorithm the work involved “in generati ng points
approximating an o pt imal point” is often “dom inated by linear algebra, usually in the
form of solution of a linear system or least squares problem and updat ing of matrix
information”. By looking at the connections between the advances in numerical linear
algebra and in optimization, O’Leary comes to the conclusion that t here is a symbiosis
between the two fields and foresees that it will continue in the current century.
In our opinio n, this symbiosis is getting stronger and st ronger, especially in the con-
text of large-scale o ptimization problems, where the solution of linear algebra problems
often dominates the computation time. A clear signal of this trend is also the organi-
zation of events gathering people working i n the two fields (see, e.g., [82, 83, 84]). The
aim of this paper is j ust to discuss the mutual impact of recent developments in numer-
ical linear algebra and optimization, focusing on the so lution of large-scale nonlinear
optimization problems and on iterative linear algebra techniques, where much progress
has been made in the last years (see [2, 5, 48, 74] and the references therein).
Thus, we consider the following general nonlinear optimizati on problem:
minimize
x
f(x)
s. t. c
I
(x) ≥ 0,
c
E
(x) = 0,
(1)
where f : <
n
−→ < is the objective function, c
I
: <
n
−→ <
m
I
and c
E
: <
n
−→ <
m
E
are the inequality and equality constraints, respectively, and m
I
+m
E
≤ n. We assume
that f, c
I
and c
E
are twice continuously differentiable and some constraint qualification
holds, such as the Linear Independence or the Mangasari an-Fromovitz one, so that a
solution of problem (1) satisfies the Karush-Kuhn-Tucker (KKT) conditions (see, e.g.,
[64, Chapter 12]).
General problems of type (1) are often hard t o solve and in the last years many
research efforts have been devoted to improve optimizati on algorithms with the double
goal of success and high performance over a wide range of pro blems. However, no single
approach has resulted uniformly robust and efficient in tackling nonli near optimization
problems. Among the various methods developed for such problems, two approaches
have emerged: Sequential Quadratic Programming (SQP) and Interior Point (IP). Both
approaches are based on the idea of moving toward a (local) solution of problem (1) by
approximately solving a sequence of “simpler” problems; of course, they strongly differ
for the characterist ics of these problems and the strategies used to solve them.
As explained in the next section, SQP and IP methods have a commo n linear algebra
kernel; at each iteration they require the solution of the so-called KKT lin ear system:
H −J
T
−J −D
v
w
=
c
d
, (2)
2
where H ∈ <
n×n
, D ∈ <
m×m
and J ∈ <
m×n
, wi th m ≤ m
I
+ m
E
. The matr ix H
is usually (an approximation of) the Hessian of the Lagrangian of problem (1) at the
current iteration, and hence it is symmetric and possi bly indefinite, J is the Jacobian
of some or all the constraints, and D is diagonal and positive semidefinite, possibly
null. Note that, in large-scal e problems, system (2) is usually sparse. In the following,
the matrix of this system is denoted by K and is called KKT matrix.
The solution of KKT systems is often the most computationally expensive task in
SQP and IP methods. Thus, the effective implement ation of such methods is highly
dependent on the availability of effective linear algebra algorithms and software, that
are able, in turn, to t ake into account specific needs of the optimiza tion solvers. Note
that KKT sys tems arise also in the solution of other optimization problems, such as lea st
squares ones, and in the more general context of saddle-point problems [7], therefore
very l arge interest is devoted to this subject. For these reasons, our discussion on linear
algebra and optimi zation is centred on the KKT system.
The remainder of the paper is organized as follows. In Section 2 we show how the
KKT sys tem arises in SQP and IP methods for solvi ng problem (1). In Section 3 we
report main properties of the KKT matr ix, which must be taken into account in s olving
the related system in the context of optimization. In Section 4 fundamental issues in
solving the KKT system are discussed, focusing on interior p oint methods and iterative
linear algebra solvers; preconditioning, adaptive termination of the inner iterations and
inertia control are addressed. Finally, in Section 5 we report some experiences i n the
application of iterative linear algebra techniques in the context of a potential reduction
method for quadratic programming. Concluding remarks are given in Section 6.
2 Linear algebra in SQP and IP methods
In order to show how KKT sys tems arise in SQP and IP methods, we give a sketch of
both, presenting only their basic ideas applied to the general pro blem (1). A deeper
discussion of these methods is outside the scope of the paper; many details can be found
in the surveys [ 33, 48, 50] and in the references therein.
Henceforth we use the following notations: g(x) = ∇f(x) (gradient of the objective
function), L(x, y) = f(x) − y
T
I
c
I
(x) − y
T
E
c
E
(x) (Lagrangi an function of the probl em),
H(x, y) ≈ ∇
xx
L(x, y) (approximation to) the Hessian of the Lagrangian function w ith
respect to x, J
I
(x) = ∇c
I
(x) and J
E
(x) = ∇c
E
(x) (Jacobian matrices of the inequality
and equality constraints, respectively). Furthermo re, the identity matrix is denoted
by I and the vector of all 1’s by e; for any vector v the diagonal ma trix di ag(v) is
denoted by the corresponding uppercase letter V , and, for any vectors v and w, (v, w)
is a shorthand for (v
T
, w
T
)
T
.
The basic idea of a SQP m ethod is to generate a sequence of approximate (local)
3
solutions of problem (1), by solving, at each iteration, a Quadratic Programming (QP)
problem, such as
minimize
δx
q(δx) ≡ δx
T
g(x) +
1
2
δx
T
H(x, y)δx
s. t. c
I
(x) + J
I
(x)δx ≥ 0,
c
E
(x) + J
E
(x)δx = 0,
(3)
where q(δx) is a quadratic (e.g. Quasi-Newt on) approximation of L(x, y) and δx is a
search direction. A commonly used strategy to solve the SQP subproblem is based on
the active-set approach. It tries to predict the inequality constraints that are active
at the solution and solves an equality constrained optimization problem, hence it is
called Sequential Equality-constrained Quadratic Programming (SEQP). The quadrati c
problem (3) reduces to
minimize
δx
q(δx) ≡ δx
T
g(x) +
1
2
δx
T
H(x, y)δx
s.t. c
A
(x) + J
A
(x)δx = 0,
(4)
where A is an estim ate of the active set at x, c
A
are t he co nstraints corresponding t o
A and J
A
is t he related Jacobi an. An optimal solution of problem (4) satisfies the first
order optimality conditions for this problem, i.e. it is solution of the linear system
H(x, y) −J
A
(x)
T
−J
A
(x) 0
δx
y
A
=
−g(x)
c
A
(x)
, (5)
where y
A
is the vector of Lagrangian multipliers for the quadratic problem (4), that
provides an approximation of the Lagrangian multipliers corresponding to c
A
in the
original problem (1). This system has the form (2 ), with H = H(x, y) and J = J
A
(x).
Systems of this type are obtained also in the case of Sequential Inequality-constrained
Quadratic Programming met hods, where no a priori prediction of the active set is
made [48]. The step δx resulting from the solution of (5) is used to update the current
approximation of the solution; actually, a linesearch or t rust-region approach must b e
applied to obtain useful updates.
A fundamental aspect for the effectiveness of SQP methods is the choice of the
Hessian approximation H(x, y ) at each iterat ion. The SQP methods also have several
criti cal shortcomings, such as the possibility that the subproblem is not convex, the
linearized constraints are inconsist ent and the iterates do not converge. For a discussion
on these issues and the strategies to deal with them the reader is referred to [4 8] and
the references therein. We only note that any variant of the SQP method outlined here
requires the solution of KKT systems.
The key idea of IP met hods is to approach a so lution of problem (1) by approxi-
mately solving a sequence of barrier problems (BPs), depending on a parameter µ > 0.
4
