You are on page 1of 106

Institute of Applied Mathematics

LECTURE NOTES SERIES

1

NUMERICAL OPTIMIZATION
Constrained Optimization
Bülent Karasözen
Gerhard-Wilhelm Weber

Design & Layout: Ömür Ugur

Preface
These are lecture notes offered to the students of the course Numerical Op-
timization at the Institute of Applied Mathematics (IAM) of Middle East
Technical University (METU) in Summer Semester 2003. In these months,
this course was held for the first time at our new and for Turkey pioneering
institute which was founded in Autumn 2002. There has been the institute
teacher’s conviction that optimization theory is an important key technology
in many modern fields of application from science, engineering, operational
research and economy and, forthcoming, even in social sciences.
To be more precise, these lecture notes are prepared on the course’s sec-
ond part which treated the case of constrained continuous optimization from
the numerical viewpoint. Here, we pay attention to both the cases of lin-
ear and nonlinear optimization (or: programming). In future, extensions of
these notes are considered, especially in direction of unconstrained optimiza-
tion. Herewith, our lecture notes are much more a service for the students
than a complete book. They essentially are a selection and a composition of
three textbooks’ elaborations: There are the works “Lineare und Netzwerkop-
timierung. Linear and Network Optimization. Ein bilinguales Lehrbuch” by
H. Hamacher and K. Klamroth (2000) used in the parts about linear pro-
gramming, “Linear and Nonlinear Programming” by S.G. Nash and A. Sofer
(1996) and “Numerical Optimization” by J. Nocedal and S.J. Wright (1999)
used for the parts about foundations and nonlinear programming.
During Summer Semester 2003, these lecture notes were given to the
students in the handwritten form of a manuscript. We express our deep
gratitude to Dr. Ömür Uğur from IAM of METU for having prepared this
LATEX-typed version with so much care and devotion. Indeed, we are looking
forward that in future many of students of IAM will really enjoy these notes
and benefit from them a lot. Moreover, we thank the scientific and institu-
tional partners of IAM very much for the financial support which has made
the lecture notes in the present form possible.
With friendly regards and best wishes,

Bülent Karasözen and Gerhard-Wilhelm Weber,

Ankara, in October 2003
ii
Contents

1 Introduction 1
1.1 Some Preparations . . . . . . . . . . . . . . . . . . . . . . . . 2

3 Linear Programming: Interior-Point Methods 31

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2 Primal-Dual Methods . . . . . . . . . . . . . . . . . . . . . . . 32
3.3 A Practical Primal-Dual Algorithm . . . . . . . . . . . . . . . 40

4 Nonlinear Programming: Feasible-Point Methods 45

4.1 Linear Equality Constraints . . . . . . . . . . . . . . . . . . . 45
4.2 Computing the Lagrange Multipliers λ . . . . . . . . . . . . . 52
4.3 Linear Inequality Constraints . . . . . . . . . . . . . . . . . . 57
4.4 Sequential Quadratic Programming . . . . . . . . . . . . . . . 67
4.5 Reduced-Gradient Methods . . . . . . . . . . . . . . . . . . . 74

5 Nonlinear Programming: Penalty and Barrier Methods 81

5.1 Classical Penalty and Barrier Methods . . . . . . . . . . . . . 82
iv CONTENTS
List of Tables

2.1 Chocolate production. . . . . . . . . . . . . . . . . . . . . . . 10

2.2 Corner points and values. . . . . . . . . . . . . . . . . . . . . 12

4.1 SQP method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.2 Reduced-gradient method. . . . . . . . . . . . . . . . . . . . . 78

5.1 Barrier function minimizers. . . . . . . . . . . . . . . . . . . . 87

vi LIST OF TABLES
List of Figures

2.1 Simplex algorithm, feasible region. . . . . . . . . . . . . . . . . 11

2.2 Simplex algorithm. . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3 Feasible region. . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4 Basis change. . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.1 Central path, projected into x-space, showing a typical neigh-

bourhood N . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2 Iterates of Algorithm IPM2 plotted in (x, s)-space. . . . . . . . 40
3.3 Central path C, and a trajectory H from the current (noncen-
tral) point (x, y, s) to the solution set Ω. . . . . . . . . . . . . 41

4.1 Sequence of movements in an active-set method. . . . . . . . . 61

4.2 Illustration of active-set algorithm. . . . . . . . . . . . . . . . 66
4.3 Zigzagging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

5.1 Effect of barrier term. . . . . . . . . . . . . . . . . . . . . . . 84

5.2 Contours of the logarithmic barrier function. . . . . . . . . . . 90
viii LIST OF FIGURES
Chapter 1

Introduction

In this chapter with its various aspects of consideration, we take into account
constraints additionally to our problem of minimizing an objective function
f . Actually, we become concerned with the problem

 minimize f (x)

subject to
(P)

 hi (x) = 0 ∀i ∈ I := {1, 2, . . . , m} ,

gj (x) ≥ 0 ∀j ∈ J := {1, 2, . . . , s} .

Here, f, hi and gj are supposed to be smooth real-valued function on Rn . By

smooth, we usually think of being one- or two-times continuously differen-
tiable. The first group of constraints, where we demand . . . = 0, are called
equality constraints. The second group of constraints, where we ask . . . ≥ 0,
are called inequality constraints. Denoting the feasible set, where we restrict
the objective function f on, by


M := x ∈ Rn hi (x) = 0 (i ∈ I), gj (x) ≥ 0 (j ∈ J) ,
our constrained optimization problem can be written as follows:

or equivalently,

(P) min f (x).

x∈M
2 CHAPTER 1. INTRODUCTION

Depending on the context, namely on the various assumptions on f, h, g, I

and J we make, we shall sometimes denote the functions, the problem and
its feasible set a bit differently or specific. For example, this will be the case
when f, h, g are linear (to be more precise: linearly affine). In fact, we are
going to distinguish the linear and the nonlinear case, speak about linear
optimization and nonlinear optimization. As in our course, the practical
character and aspect of optimization is more emphasized than the analytical
or topological ones, we also talk about (non-)linear programming. This may
remind us of both the motivation of optimization by concrete applications,
and the numerical solution algorithms. Because of the priority we gave to the
numerical-algorithmical aspect, we make the following introductory Section
not so long.

1.1 Some Preparations

By the following definition we extend corresponding notions from uncon-
strained optimization almost straightforwardly:
Definition 1.1. Let a vector (or, point) x∗ ∈ Rn be given.

(i) We call x∗ a local solution of (P), if x∗ ∈ M and there is a neighbour-

hood N (x∗ ) of x∗ such that

(ii) We call x∗ a strict local solution of (P), x∗ ∈ M and there is a neigh-

bourhood N (x∗ ) of x∗ such that


f (x∗ ) < f (x) ∀x ∈ M ∩ N (x∗ ) \ {x∗ } .

(iii) We call x∗ an isolated local solution of (P), x∗ ∈ M and there is a

neighbourhood N (x∗ ) of x∗ such that x∗ is the only local solution of
(P) in M ∩ N (x∗ ).

Sometimes, we also say (strict or isolated) local minimizer for a (strict or

local) local solution.

Do you understand and geometrically distinguish these three conditions?

We shall deepen our understanding in the exercises and in the following
sections.
1.1. SOME PREPARATIONS 3

Example 1.1. Let us consider the problem

 
minimize f (x) := x1 + x2 
(P)
subject to x21 + x22 − 2 = 0 i.e., h(x) := x21 + x22 − 2 ,
Here, | I | = {1} (h1 = h) and√J = ∅. We see by inspection that the feasible
set M is the circle of radius 2 centered at the origin (just the boundary
of the corresponding disc, not its interior). The solution x∗ is obviously
(−1, −1)T , this is the global solution of (P): From any other point on the
circle, it is easy to find a way to move that stays feasible √
(i.e., remains on the
circle) while decreasing f . For instance, from the point ( 2, 0)T any move in
the clockwise direction around the circle has the desired effect. Would you
please illustrate this and the following graphically? Indeed, we see that at
x∗ , the gradient ∇h(x∗ ) is parallel to the gradient ∇f (x∗ ):
 
1
( = ) ∇f (x∗ ) = λ∗ ∇h(x∗ )
1
where λ∗ = λ∗1 = − 21 .

The condition found in the previous example to be necessary for a (local)
minimizer (or, local minimum) can be expressed as follows:

∇x L(x∗ , λ∗ ) = 0,
where L : Rn × Rm −→ R is defined by L(x, λ) := f (x) − λh(x) and called
the Lagrange function. The value λ∗ is called the Lagrange multiplier (at x∗ ).
This Lagrange multiplier rule is a first-order necessary optimality condition
1
(NOC) which can be extended to cases where I is of some other cardi-
nality m ≤ n. in fact, provided that the Linear Independence Constraint
Qualification (a regularity condition) holds at x∗ , saying that the gradients

∇h1 (x∗ ), ∇h2 (x∗ ), . . . , ∇hm (x∗ ),

regarded as a family (i.e., counted by multiplicity), are linearly independent,
then there exist so-called Lagrange multipliers λ∗1 , λ∗2 , . . . , λ∗m such that

m
X

∇f (x ) = λ∗i ∇h(x∗ ).
i=1
1
In the next sections, we also consider feasibility to be a part of it.
4 CHAPTER 1. INTRODUCTION

The latter equation is just

∇x L(x∗ , λ∗ ) = 0,
where L : Rn × Rm −→ R is the so-called Lagrange function L(x, λ) :=
f (x) − λT h(x), with λ := (λ1 , λ2 , . . . , λm )T and h := (h1 , h2 , . . . , hm )T .
Next, we slightly modify Example 1.1. We replace the equality constraint
by an inequality.
Example 1.2. Let us consider the problem

 
minimize f (x) := x1 + x2 
(P)
subject to 2 − x21 − x22 ≥ 0 i.e., g(x) := 2 − x21 − x22 ,
here, I = ∅ and | J | = 1 (g = g1 ). The feasible set consists of the cir-
cle considered in Example 1.1 and its interior. We note that the gradient
∇g(x) = (−2x1 , −2x2 )T points from the √ circle in direction of the √ interior.
For example, inserting the point x = ( 2, 0) gives ∇g(x) = (−2 2, 0)T ,
T

pointing in the direction of the negative x1 − axis, i.e., inwardly (regarded

from x). By inspection, we see that the solution of (P) is still x∗ = (−1, −1)T ,
and we have now

∇f (x∗ ) = µ∗ ∇g(x),
1
where µ∗ = µ∗1 = 2
≥ 0.
Again, we can write our necessary equation by using a Lagrange function,
namely, L(x, µ) := f (x) − µg(x):

∇x L(x∗ , µ∗ ) = 0, where µ∗ ≥ 0.

Again, our necessary optimality condition can be generalized by admitting

another finite cardinality s of J, and m = | I | possibly to be positive (m ≤ n).
For this purpose, we need a regularity condition (constraint qualification) at
our local minimizer x∗ , once again. In the presence of both equality and
inequality constraints, we can again formulate LICQ. For this purpose, we
introduce the set of active inequality constraints at some feasible point x:


J0 (x) := j ∈ J gj (x) = 0 .
Then, LICQ at x∗ means the linear independence of the vectors (together
regarded as a family)
1.1. SOME PREPARATIONS 5


∇hi (x∗ ) (i ∈ I), ∇gj (x∗ ) j ∈ J0 (x∗ )
(so that the gradients of the active inequalities are taken into account to the
gradients of the equality constraints, additionally). Another famous regular-
ity condition is somewhat weaker than LICQ and called MFCQ (Mangasarian-
Fromovitz Constraint Qualification). The geometrical meaning of MFCQ is
the following: At each point of M we have an inwardly pointing direction (rel-
ative to the surface given by the zero set of h). In the sequel, we refer to LICQ
for simplicity. Now, our first-order necessary optimality conditions (NOC)
are called Karush-Kuhn-Tucker
 conditions: There exist λ∗ = (λ∗1 , . . . , λ∗m )T
and µ∗ = µ∗j j∈J0 (x∗ ) such that
 m
 X X
 ∇f (x∗ ) = λ∗i ∇hi (x∗ ) + µ∗j ∇gj (x∗ ),
 i=1 j∈J0 (x∗ )
 ∗ ∗
µj ≥ 0 ∀j ∈ J0 (x ).
The second ones of the Lagrange multipliers, namely µ∗j (j ∈ J0 (x∗ )), neces-
sarily are nonnegative. Here, we may interpret this as by the existence of (rel-
atively) inwardly pointing direction along which f does not decrease, when
starting from our local minimizer x∗ . Let us summarize our consideration
by referring
 to the Lagrange function L(x, λ, µ) := f (x) − λT h(x) − µT g(x),
where µ∗j j∈J0 (x∗ ) is filled up by additional parameters µj (j 6∈ J0 (x∗ )) (at
x = x∗ , we have µ∗j = 0, j 6∈ J0 (x∗ )).

Theorem 1.1. Suppose that x∗ is a local solution of (P), and that LICQ
holds at x∗ . Then, there uniquely exist Lagrange multiplier vectors λ∗ ∈
Rm , µ∗ ∈ Rs , such that the following conditions are satisfied:

 ∇x L(x∗ , λ∗ , µ∗ ) = 0,

 h(x∗ ) = 0,
(NOC) g(x∗ ) ≥ 0,

 µ∗ ≥ 0,

µ∗j gj (x∗ ) = 0 ∀j ∈ J.

The latter multiplicative conditions are local complementarity conditions. To

which case study do they give rise?

Proof. See Nocedal, Wright (1999); the proof of uniqueness is left to the
6 CHAPTER 1. INTRODUCTION

Often, the compact form of (NOC), where we do not explicitly refer to

the active inequalities j ∈ J0 (x∗ ), is more convenient. In other times, a direct
reference to J0 (x∗ ) is better for us. We note that our first-order necessary
optimality conditions are not the only ones. In the subsequent sections and in
the presentation given in lectures, we shall briefly indicate them sometimes.
These conditions are, more or less, generalization of the condition (taught to
us in the unconstrained case) that the gradient vanishes at a local solution:
∇f (x∗ ) = 0, if I = J = ∅. Provided that all the defining functions are
twice continuously differentiable, then we know the second-order optimality
condition saying that the Hessian matrix has to be positive semi-definite:

ρT ∇2xx f (x∗ )ρ ≥ 0 ∀ρ ∈ Rn .

But how does this condition look like in our case of constrained optimization?
Let us define the tangent space of M at x∗ as follows:


Tx∗ := ρ ∈ Rn ∇T hi (x∗ )ρ = 0 ∀i ∈ I, ∇T gj (x∗ )ρ = 0 ∀j ∈ J0 (x∗ ) .

Then, our second-order condition states that the Hessian of the Lagrange
function L at x∗ is positive semi-definite over Tx∗ . Here, we take the vectors
ρ ∈ Rn for left- or right-multiplication from Tx∗ . This can also be expressed
as the positive semi-definiteness of this Hessian, when firstly left- and right-
multiplied by B T and B, respectively, over Rn . Here, B is any matrix whose
columns constitute a basis of the linear space Tx∗ . Herewith, we can express
our second-order optimality condition so, where n b = n − m,b mb = | J0 (x∗ ) |:

In literature, we also find some refined condition, where in the statement

of positive-definiteness the reference space Tx∗ is replaced by the set Cx•∗
which is a (tangent) cone. Namely, Cx•∗ comes from substituting the gradient
conditions on the active inequalities by

∇T gj (x∗ )ρ ≥ 0 ∀j ∈ J0 (x∗ ) with µ∗j = 0.

1.1. SOME PREPARATIONS 7

Theorem 1.2. Suppose that x∗ is a local solution of (P), and that LICQ
is satisfied at x∗ . Let λ∗ and µ∗ be the corresponding Lagrange multiplier
vectors so that (NOC) is (necessarily) fulfilled.
Then, the Hessian ∇2xx L(x∗ , λ∗ , µ∗ ) is positive semi-definite over Tx∗ , or
(refined version) over Cx•∗ .

Proof. See, e.g., Nash, Sofer (1996) and Nocedal, Wright (1999). 
This theorem plays an important role for detecting that some found or
submitted stationary point x∗ (i.e., x∗ fulfills our first-order necessary opti-
mality conditions) is not a local solution, but a local minimizer or a saddle
point. How can we use Theorem 1.2 for this purpose? A saddle point is in
between of a local minimizer and a local maximizer: there are feasible direc-
tions of both an increasing and a decreasing behaviour of f . But how can we
detect that a stationary point, i.e., a candidate for being a local minimizer
really is such a one? what we know from the unconstrained case gives rise to
assume that we should turn from positive semi-definiteness (over Tx∗ or Cx•∗ )
to positive definiteness (over Tx∗ \ {0} or Cx•∗ \ {0}). In fact, by this sharp-
ening of (NOC)s.o. we obtain a condition which together with our first-order
(NOC) is sufficient for x∗ to become a local solution. In particular, we arrive
at

We state the following theorem on second-order sufficient optimality condi-

tions:

Theorem 1.3. Suppose that for some point x∗ ∈ M there are Lagrange mul-
tiplier vectors λ∗ , µ∗ such that the conditions (NOC) (of first-order) are satis-
fied and µ∗j > 0 (j ∈ J0 (x∗ )). Suppose also that the Hessian ∇2xx L(x∗ , λ∗ , µ∗ )
is positive definite over Tx∗ , or over Cx•∗ (i.e., over Tx∗ \ {0}, or Cx•∗ \ {0}).
Then, x∗ is a strict local solution for (P).

Proof. See, e.g., Nash, Sofer (1996) and Nocedal, Wright (1999). 
In our exercises, we shall consider a few examples on finding local or
global solutions for nonlinear optimization problem, where the utilization of
numerical-algorithmical methods is still not necessary. Let us also conclude
our Section 1.1 with another such example.
Example 1.3. Let us consider the problem (P) from Example 1.2:
8 CHAPTER 1. INTRODUCTION

 
minimize f (x) := x1 + x2 
(P)
subject to 2 − x21 − x22 ≥ 0 i.e., g(x) := 2 − x21 − x22 .

L(x, µ) := (x1 + x2 ) − µ(2 − x1 − x2 ),

and it is easy to show that the Karush-Kuhn-Tucker conditions (NOC) are
satisfied by x∗ = (−1, −1)T , with µ∗ = 21 . The Hessian of L (with respect to
x) is
   
2µ∗ 0 1 0
∇2xx L(x∗ , µ∗ ) = = .
0 2µ∗ 0 1
This matrix is positive definite, i.e., it satisfies ρT ∇2xx L(x∗ , µ∗ )ρ > 0 for all
ρ ∈ Rn \ {0} (in particular, for all elements of Tx∗ , Cx•∗ , except 0). So, it
certainly satisfies the conditions of Theorem 1.3.
We conclude that x∗ = (−1, −1)T is a strict local solution for (P). Let us
remark that x∗ is a even a global solution, since (P) is a convex programming
problem. 

In this course, we do not pay much attention to stability aspects which,

however, can be well analyzed by LICQ and (SOC)s.o . Next, we concentrate
of the linear case of (P).
Chapter 2

Linear Programming:
Foundations and Simplex
Method

In this chapter, based on the book of Hamacher, Klamroth (2000), we con-

sider easy examples which lead to optimization (“programming”) problems:
Linear programs play a central role in modelling of optimization problems.
Example 2.1. The manufacturer “Chocolate & Co.” is in the process of
reorganizing its production. Several chocolate products which have been
produced up to now are taken out off production and are replaced by two
new products. The first product P1 is fine cocoa, and the second one P2 is
dark chocolate. The company uses three production facilities F1, F2, and F3
for the production of the two products. In the central facility F1 the cocoa
beans are cleaned, roasted and cracked. In F2, the fine cocoa and in F3 the
dark chocolate are produced from the preprocessed cocoa beans. Since the
chocolate products of the company are known for their high quality, it can
be assumed that the complete output of the company can be sold on the
market. The profit per sold production unit of P1 (50kg of cocoa) is 3 (¤),
whereas it is 5 per sold production unit of P2 (100kg of chocolate). However,
the capacities of F1–F3 are limited as specified in Table 2.1.
(For example, the row corresponding to F1 implies that the per unit P1 and
P2, 3% and 2%, respectively, of the daily capacity of production facility F1
are needed, and that 18% of the daily capacity of F1 is available for the
production of P1 and P2.)
The problem to be solved is to find out how many units x1 of products
P1 and x2 of product P2 should be produced per day in order to maximize
CHAPTER 2. LINEAR PROGRAMMING: FOUNDATIONS AND
10 SIMPLEX METHOD
P1 P2 available capacity
(in % of the daily capacity)
F1 3 2 18
F2 1 0 4
F3 0 2 12
Table 2.1: Chocolate production.

the profit (while satisfying the capacity constraints.) This problem can be
formulated as a linear program (LP):

 maximize 3x1 + 5x2 =: cT x

 subject to the constraints

3x1 + 2x2 ≤ 18 (I)
(LP)

 x 1 ≤ 4 (II)

 2x2 ≤ 12 (III)

x1 , x2 ≥ 0 (IV), (V)
The function cT x is the (linear) objective function of the LP. The constraints
are partitioned into functional constraints ((I)–(III)) and nonnegativity con-
straints ((IV),(V)). Each x which satisfies all the constraints, is called a fea-
sible solution of the LP and cT x is its objective (function) value. Instead of
maximizing the function, it is often minimized. In this case, the coefficients
ci of the objective function can be interpreted as unit costs. The functional
constraints can also be equations (“=”) or inequalities with “≥”.
Let us continue by introducing a graphical procedure for the solution of
(LP). In a first step, we draw the set of feasible solutions (the feasible region
P) of (LP),
  that is the set of all points (vectors) (x1 , x2 ), in matrix notation:
x1
, satisfying all the constraints; see Figure 2.1.
x2
If Ai x ≤ bi is one of the constraints (e.g., Ai = (3, 2), b1 = 18) we draw the
line corresponding to the equation

Ai x = b i .
This space separates the space R2 into two half-spaces

Ai x ≤ bi , and Ai x ≥ bi
The set P of feasible solutions is a subset of the half-space Ai x ≤ bi . We
obtain P by taking the intersection of all these half-spaces including the half-
11

spaces x1 ≥ 0 and x2 ≥ 0. A set of points in R2 , which is obtained in such a

way is called a convex polyhedron.
In a second step, we draw the objective function z := cT x = 3x1 + 5x2 .
For any given value of z we obtain a line, and for any two different values
of z we obtain two parallel lines for which z is as large as possible. The
following Figure 2.2 shows the feasible region P and, additionally, the line
corresponding to the objective function 3x1 + 5x2 , i.e., the line corresponding
to the objective value z = 15. The arrow perpendicular to this line indicates,
in which direction the line should be shifted in order to increase z.

cT x = z with P in some x ∈ P corresponds

The intersection of any line 
x1
to a feasible solution x = with objective value z. Hence, in order to
x2
maximize the objective function, the line is shifted parallel as far as possible
without violating the condition
CHAPTER 2. LINEAR PROGRAMMING: FOUNDATIONS AND
12 SIMPLEX METHOD


x ∈ R2 cT x = z ∩ P 6= ∅.
 
∗ 2
This just gives the line which passes through the point x = . Hence,
6
z ∗ = cT x∗ = 3 · 2 + 5 · 6 = 36 is the maximum possible profit.
As a consequence, the company “Chocolate & Co.” will produce 2 pro-
duction units of cocoa and 6 production units of the dark chocolate per day.


In this example, we observe an important property of linear programs:

There is always an optimal solution (i.e., a feasible solution for which cT x is
maximal) which is a corner point of P –or there is no optimal solution at all.
This general property of LPs is a fundamental for the method introduced
below as a general solution for LPs: the simplex method :
Imagine to move from corner point to corner point of P (a corner point
is an element of P which is an intersection of two or more of its facets)
using a procedure which guarantees that the objective value is improved in
every step. As soon as we have reached a corner point in which no further
improvement of the objective value is possible, the algorithm stops and the
corner point corresponds to an optimal solution of the LP. In Example 2.1,
such a sequence of corner points is, e.g., as stated in Table 2.2.

corner
 point
 objective value 3x1 + 5x2
0
1. 0
 0 
4
2. 12
 0 
4
3. 27
 3 
2
4. 36 (STOP), no further improvement possible.
6
Table 2.2: Corner points and values.

This is exactly the idea on which the simplex method is based, which was
developed by G. Dantzig in 1947. It yields an optimal solution (whenever it
exists) for only LP. An important step in the development of the preceding
idea for arbitrary LPs is the application of methods from linear algebra.

Definition 2.1. A given LP is said to be in standard form if and only if it

is given as
13

 minimize cT
(LP)st subject to Ax = b

x ≥ 0 (i.e., xi ≥ 0 ∀i ∈ {1, . . . , n}),
where A in an m × n matrix, m ≤ n, and rank(A) := dim(Im T ) = m (T
being the linear transformation represented by A relative to standard bases).

If m ≤ n or rank(A) = m are not the case, then we can reintroduce

m̃ ≤ n, rank(A) =: m̃ by omitting redundant (i.e., linearly on other m̃ ones
depending) constraints. In the following, we show how any given LP can be
transformed into standard form.
Assume that an LP in general form is given:

 minimize c1 x1 + · · · + cn xn

 subject to ai1 x1 + · · · + ain xn = bi (i ∈ {1, . . . , p})

ai1 x1 + · · · + ain xn ≤ bi (i ∈ {p + 1, . . . , q})
(LP)gf

 ai1 x1 + · · · + ain xn ≥ bi (i ∈ {q + 1, . . . , m})

 xj ≥ 0 (j ∈ {1, . . . , r})

xj ≤ 0 (j ∈ {r + 1, . . . , n})

(a) The “≤” constraints can be transformed into equality constraints by

introducing slack variables

xn+i−p := bi − ai1 x1 − · · · − ain xn (i ∈ {p + 1, . . . , q}).

We obtain


ai1 x1 + · · · + ain xn + xn+i−p = bi
(i ∈ {p + 1, . . . , q}).
xn+i−p ≥ 0

(b) Analogously, the “≥” constraints can be transformed into equality con-
straints by introducing surplus variables

xn+i−p := ai1 x1 + · · · + ain xn − bi (i ∈ {q + 1, . . . , m}).

We obtain


ai1 x1 + · · · + ain xn − xn+i−p = bi
(i ∈ {q + 1, . . . , m}).
xn+i−p ≥ 0
CHAPTER 2. LINEAR PROGRAMMING: FOUNDATIONS AND
14 SIMPLEX METHOD

(c) An LP in which the objective function is to be maximized can be trans-

formed into standard form by using the identity

 
max cT x · · · = − min (−c)T x · · · ,
where “· · · ” stand for properties required for x, and by solving an LP
with the coefficients −cj in standard form.

replace it by

xj =: x+ −
j − xj with x+ −
j , xj ≥ 0,

in connection with LPs, we will usually denote as follows (see (LP)st ):

 
A = aij i∈{1,...,m} = aij : m ≤ n, rank(A) = m,
j∈{1,...,n}

a1j
 
Aj =  ...  , Ai = (ai1 . . . ain ),
amj
  
c1 b1
   
c =  ...  , b =  ...  ,
cn bm


P = x ∈ Rn Ax = b, x ≥ 0 .

Basic Solutions: Optimality Test and Basic Exchange


Definition 2.2. A basis of A is a set B := AB(1) , . . . , AB(m) of m linear
independent columns of A. The  index set is given in an arbitrary but fixed
order: B := B(1), . . . , B(m) .
Often, we call (somewhat inexact) the index set B itself a basis of A.

By AB := AB(1) , . . . , AB(m) we denote the regular m × m sub-matrix
of A corresponding to B. The corresponding variables xj are called basic
variables, and they are collected in the vector
15

 
xB(1)
 
xB :=  ...  .
xB(m)
The remaining indices are contained in a set which is again in arbitrary but
fixed order: N := N (1), . . . , N (n − m) . The variables xj with j ∈ N are
called non-basic variables.

We immediately obtain

x is a solution of Ax = b ⇔ 

(xB , x
N ) is asolution of AB xB + AN xN  =b 
xB xB 

(i.e., is a solution of (AB , AN ) = b). 
xN xN
(This can easily be seen by changing the order of the columns Aj suitably.)
By multiplying the latter equation by A−1
B we obtain:

x is a solution of Ax = b ⇔
(xB , xN ) is a solution of xB = A−1 −1
B xB − A B A N xN ; (∗)
(∗) is called the basic representation of x (with respect to B).

Definition 2.3. For any choice of the non-basic variables xN (j) (j ∈ {1, . . . , n − m})
we obtain, according to (∗), uniquely defined values of the basic variable
xB(j) (j ∈ {1, . . . , m}). The basic solution (with respect to B) is the solution
of Ax = b with xN = 0 and, xB = A−1 B b.

Example 2.2. We consider the LP

 maximize x1

subject to − x 1 + x2 ≤ 1
(LP)

 x1 + x 2 ≤ 3

x1 , x2 ≥ 0.

The feasible region is illustrated in Figure 2.3.

We can easily see that x∗ = (3, 0)T is the optimal solution. First, we
transform this LP into standard form by introducing slack variables x3 , x4
and by transforming the maximization into a minimization (the “−” sign
CHAPTER 2. LINEAR PROGRAMMING: FOUNDATIONS AND
16 SIMPLEX METHOD

Figure 2.3: Feasible region.

in front of “min” can be dropped, but we must not forget this “−” when
interpreting the final result). Hence we obtain

 minimize − x1

subject to − x 1 + x2 + x3 = 1
(LP)st

 x1 + x 2 + x 4 = 3

x1 , x2 , x3 , x4 ≥ 0,

 
−1 1 1 0
A= .
1 1 0 1

(i) For B = (B(1), B(2)) = (3, 4) we obtain the basic solution

     −1  
xB(1) x3 1 0 1
xB = = = A−1
B b = b=b= ,
xB(2) x4 0 1 3

   
xN (1) x1
xN = = = 0,
xN (2) x2
hence, x = (0, 0, 1, 3)T .
(ii) For B = (1, 2) we get by a proposition from linear algebra
   
−1 1 1 −1 1
AB = =⇒ A−1
B = ,
1 1 2 1 1
17

therefore,

        
xB(1) x1 1 −1 1 1 1
xB = = = A−1
B b = = ,
xB(2) x2 2 1 1 3 2

   
xN (1) x3
xN = = = 0,
xN (2) x4
hence, x = (1, 2, 0, 0)T .

(iii) For B = (4, 1) we get (exercise)

   −1    
x4 0 −1 1 4
xB = = A−1
B b = = ,
x1 1 1 3 −1
 
x2
xN = = 0,
x3
hence, x = (−1, 0, 0, 4)T .

Evaluation: In the cases (i), (ii), the basic solutions correspond to corner
points of P, namely, (i): x corresponds to the corner point (0, 0), and (ii): x
corresponds to the corner point (1, 2). Since xB ≥ 0 in both cases, (i) and
(ii) yield basic feasible solutions.
In the case (iii) however, xB ≥ 0 is not satisfied. So, (xB , xN ) is a basic
solution which is not feasible. In Figure 2.3, it can be easily seen that the
corresponding point (x1 , x2 ) = (−1, 0) is not contained in P.


Now, we consider a basic feasible solution (xB , xN ) and use the basic
representation with respect to B to derive an optimality criterion. First, we
partition the coefficients of c into

cN := (cN (1) , . . . , cN (n−m) )T .

Then, the objective value cT x can be written as
CHAPTER 2. LINEAR PROGRAMMING: FOUNDATIONS AND
18 SIMPLEX METHOD

cT x = cTB xB + cTN xN 
= cTB A−1B b − A −1
B A N x N + cTN
 xN
T −1 T T −1
= c B A B b + c N − c B A B A N xN .
In the current basic feasible solution we have that xN = 0 and, therefore, its
objective value is

cTB xB = cTB A−1

B b.
T
 solutions the objective value differs from this value by cN −
For all other
T −1
cB AB AN xN . If we increase the value of xN (j) = 0 to xN (j) = δ, δ > 0, then
the objective value is changed by

=:zN (j)
z }| {
T −1
δ cN (j) − cB AB AN (j) .
| {z }
=:Π

Therefore, the objective value of the given basic feasible solution increases if
cN (j) − zN (j) > 0, and it decreases if cN (j) − zN (j) < 0. The values

cN (j) := cN (j) − zN (j) ,

called the reduced or relative costs of the non-basic variable xN (j) , thus con-
tain the information whether it is useful to increase the value of the non-basic
variable xN (j) from 0 to a value δ > 0. In particular, we obtain

Theorem 2.1 (Optimality Condition for Basic Feasible Solutions).

If x is a basic feasible solution with respect to B and if

cN (j) = cN (j) − zN (j) = cN (j) − cTB A−1

B AN (j) ≥ 0 ∀j ∈ {1, . . . , n − m} ,

then x is an optimal solution of the LP

 minimize cT
(LP)st subject to Ax = b

x ≥ 0.

As we have seen in Example 2.2, basic feasible solutions of Ax = b cor-

bx ≤ b. Here, A
respond to extreme points of the feasible region Ab b is derived
from A by deleting those columns of A which correspond to a unit matrix,
19

and xb are the components of x which correspond to the columns of A. b Ac-

cording to the idea mentioned at the beginning of this chapter, we will move
to a new basic feasible solution if the optimality condition is not satisfied for
the current basic feasible solution.

Example 2.3 (Continuation of Example 2.2). We apply our optimality test

to two basic feasible solutions of Example 2.2:

(a) B = (1, 2): In (ii) we computed that

 
1 −1 1
A−1
B = .
2 1 1

Thus, we obtain:


zN (1) = cB(1) , cB(2) A−1 B AN (1)
−1
= (c1 , c2 )AB A3  
1 −1 1 1
= (−1, 0) 2
  1 1 0
1
= 21 (1, −1)
0
1
= 2 (i.e., Π = 21 (1, −1)), N (1) = 3
=⇒ cN (1) = c3 − z3 = 0 − 21 < 0.

So, the optimality condition is violated. Note that, using Theorem 2.1
one cannot conclude that the corresponding basic feasible solution is
not optimal, since our theorem gives only a sufficient optimality condi-
tion.

   
−1 1 0 1
AB = =⇒ A−1
B = .
1 0 1 1

We check the optimality condition by computing c := cN − zN with

T
zN := cTB A−1
B AN , and check whether cN ≥ 0.
CHAPTER 2. LINEAR PROGRAMMING: FOUNDATIONS AND
20 SIMPLEX METHOD

T
zN = cTB A−1
B AN
= (c1 , c3 )A−1
B (A1 , A )
4 
0 1 1 0
= (−1, 0)
 1 1  1 1
1 1
= (−1, 0)
2 1
= (−1, −1)
=⇒ cN = (0, 0)T − (−1, −1)T = (1, 1)T ≥ 0.
Hence, the basic feasible solution
      
x1 0 1 1 3
xB = = A−1
B b = = ,
x3 1 1 3 4
   
x2 0
xN = =
x4 0
corresponding to B is optimal. This basic feasible solution corresponds
to the corner point
   
∗ x1 3
x = = ,
x2 0
which we have already identified in Example 2.2 as optimal using the
graphical-geometrical procedure.

In the following, we show how to obtain a new feasible solution if the

current one does not satisfy the optimality condition.
Suppose that cN (s) = cN (s) − zN (s) < 0 for some s ∈ N . Since

cT x = cTB A−1 T T −1
B b + (cN − CB AB AN )xN

increasing xN (s) by one unit will decrease cT x by cN (s) . Since our goal is to
minimize the objective function cT x, we want to increase xN (s) by as many
units as possible.
How large can xN (s) be chosen? This question is answered by the basic
representation (∗) and the requirement xB ≥ 0. If we keep xN (j) , j 6= s,
equal to 0 and increase xN (s) from 0 to δ ≥ 0, then we obtain for the resulting
solution of Ax = b
21

xB = A−1 −1
B b − AB AN (s) δ.

Denoting the ith component of A−1 −1

B b and AB AN (s) by b̃i and ãiN (s) , respec-
tively, we have to choose δ such that

xB(i) = b̃i − ãiN (s) δ ≥ 0. (⊗)

In order to choose δ as large as possible, we thus compute δ as
( )
b̃i
δ = xN (s) := min ãiN (s) > 0 . (min ratio rule)
ãiN (s)
While computing δ with respect to the min ratio rule, two cases may occur:

Case 1: ãiN (s) ≤ 0 ∀i ∈ {1, . . . , m}.

Then, δ can be chosen arbitrarily large without violating any of the nonneg-
ativity of constraints. As a consequence, cT x can be made arbitrarily small
and the LP is unbounded. Thus we obtain:

Theorem 2.2 (Criterion on Unbounded LPs). If x is a basic feasible

solution with respect to B and if

B AN (s) ≤ 0

for some s ∈ N , then the LP

 minimize cT
(LP)st subject to Ax = b

x ≥ 0

is unbounded.

Case 2: ∃i ∈ {1, . . . , m} : ãiN (s) > 0.

In this case, the minimum in the computation of xN (s) by the min ratio rule
b̃r
is attained. Suppose that we get δ = xN (s) = . (If the index r is not
ãrN (s)
b̃r
uniquely determined, then we choose any of the indices such that = δ.)
ãrN (s)
According to (⊗) and the min ratio rule, the new solution is
CHAPTER 2. LINEAR PROGRAMMING: FOUNDATIONS AND
22 SIMPLEX METHOD

b̃r
xN (s) = , xN (j) = 0 ∀j 6= s,
ãrN (s)

b̃r
xB(i) = b̃i − ãiN (s) xN (s) = b̃i − ãiN (s) ∀i.
ãrN (s)
In particular, xB(r) = 0. It is easy to check that


B 0 = B 0 (1), . . . , B 0 (r − 1), B 0 (r), B 0 (r + 1), . . . , B 0 (m)
with

0 B(i) if i 6= r
B (i) :=
N (s) if i = r
defines a new basis for A (argument: ãrN (s) 6= 0). Thus, the computation of
xN (s) has induced a basis change. The variable xB(r) has left the basis, and
xN (s) has entered the basis. The new index set of non-basic variables is


N 0 = N 0 (1), . . . , N 0 (s − 1), N 0 (s), N 0 (s + 1), . . . , N 0 (n − m)

with

0 N (j) if j 6= s
N (j) :=
B(r) if j = s.
Our basis change is indicated by Figure 2.4.

Figure 2.4: Basis change.

Example 2.4 (Continuation of Example 2.3). In (a), we saw that the opti-
mality condition is violated for B = (1, 2). Since c3 < 0, we want to let
xN (1) = x3 enter into the basis. We compute (cf. Example 2.3):
23

      
ã1N (1) 1 −1 1 1 1 −1
= A−1
B AN (1) = = ,
ã2N (1) 2 1 1 0 2 1

      
b̃1 1 −1 1 1 1
b̃ = = A−1
B b = =
b̃2 2 1 1 3 2

b̃2
=⇒ xN (1) = x3 = = 4.
ã2N (1)
(Note that here the minimization of the min ratio rule is taken over a single-
element set.) Thus, we obtain B 0 = (1, 3) as our new index set of the basis
(x2 has left the basis, x3 has entered the basis). The corresponding basic
feasible solution is

xB 0 (2) = xN (1) = x3 = 4,
1
xB 0 (1) = b̃1 − ã1N (1) xN (1) = 1 − − 4 = 3,
2
1
xN 0 (1) = xB(2) = b̃2 − ã2N (1) xN (1) = 2 − 4 = 0,
2
xN 0 (2) = xN (2) = 0.
This gives the same solution which was obtained directly in Example 2.3 (b)
by using the definition of the basic feasible solution with respect to B 0 =
(1, 3).

It is the idea of the simplex method to move iteratively from the basic
feasible solutions to basic feasible solution until an optimal basic feasible
solution is reached. Nevertheless, it still remains to be shown that there
is always an optimal basic feasible solution and that the procedure is, in-
deed, finite. But first we show how to organize the basic exchange and the
optimality test efficiently.
If we write the objective function as −z + c1 x1 + · · · , +cn xn = 0, then
the objective function and the constraints can be stored in a matrix which
we write in tableau form, called the starting tableau T = tij i∈{0,1,...,m} :
j∈{0,1,...,n,n+1}

1 c 1 . . . cn 0
0 a11 . . . a1n b1 1 cT 0
T = .. .. .. =
. . . 0 A b
0 am1 . . . amn bm
CHAPTER 2. LINEAR PROGRAMMING: FOUNDATIONS AND
24 SIMPLEX METHOD

Here, T represents a system of linear equations with m + 1 equations.

The (n + 1)-st column contains the information about the right-hand sides
of the equations. If B is a basis, then we denote by TB the nonsingular
(m + 1) × (m + 1) matrix

 
1 cTB
TB := .
0 AB

It is easy to verify that

 
1 −cTB A−1
TB−1:= B
,
0 A−1
B 
−1 1 cT − cTB A−1
B A −c T −1
A
B B b
TB T = =: T (B).
0 A−1
B A A−1B b

We call T (B) the simplex tableau associated with the basis B. Since TB−1
is nonsingular, T (B) represents the same system of linear equations as the
starting tableau T . The entries of T (B) can be interpreted as follows:

(i) The first column is always the first standard vector E1T . It emphasizes
the character of the 0th row as an equation. Later on, we will omit this
column.

(ii) For j = B(i) we have:

A−1 T
B Aj = E i (being a column).

Furthermore,

cj − cTB A−1
B Aj = c j − c j = 0 then

T (B) contains in the column corresponding to the ith basic variable

xB(i) the value 0 in row 0 and, then, the ith unit vector EiT with m
components.

(iii) For j = N (i) we have:

A−1
B Aj = (e amj )T .
a1j , ..., e
25

Furthermore,
t0j = cj − cTB A−1
B Aj = c j

(iv) In the bast column A−1

B b is the vector of the basic variables with respect
to B. Consequently, cTB A−1B b is the negative of the objective value of
the current basic feasible solution.

Example 2.5. (Continuation of Example 2.4). If we consider again Exam-

ple 2.2 with B = (1, 2), then

1 -1 0 0 0 0
T = 0 -1 1 1 0 1
0 1101 3
and because of
 
− 21 1
A−1
B = 1
2
and
1
2  1
2
1

 −2 2 
cTB A−1
B = −1 0   = 1
2
− 12
1 1
2 2

we get
 
1 −0.5 0.5
TB−1 :=  0 −0.5 0.5  .
0 0.5 0.5

Hence, the simplex tableau corresponding to B is

1 0 0 -0.5 0.5 1
T (B) = TB−1 T = 0 1 0 -0.5 0.5 1
0 0 1 0.5 0.5 2

Following the interpretation of T (B), the reduced costs c3 = − 21 , c4 =

0.5 of the non-basic variables can be taken from the 0th of T (B). It can
be easily seen that the optimality condition is not satisfied (which we know
already from Example 2.3).
Looking at the last column of the tableau, it can be seen that x1 = 1 , x2 = 2
CHAPTER 2. LINEAR PROGRAMMING: FOUNDATIONS AND
26 SIMPLEX METHOD

are the values of the basic variables in the basic feasible solution, yielding an
objective value of −t0n+1 = −1.


If ∃j ∈ { 1, ..., n } : t0j < 0, then we try to more the non-basic variable

xnj into the basis. Since the entries of the tableau are

t1j = e
a1j , ..., tmj = e
amj ,

and
t1n+1 = eb1 , ..., tmn+1 = ebm ,
we can perform the min ratio test using the simplex tableau in a very simple
way:
( )
ti n+1 e
e
δ = xj := min ti j > 0 .
e
ti j

Thus, an unbounded objective function can be recognized by the fact that

the column corresponding to one of the non - basic variables xj with t0j < 0
contains only entries ≤ 0 (cf. Theorem 2.2).
If δ = trtn+1
rj
, a pivot operation is carried out with the element tr j > 0,
i.e., we transform the j th column of T (B) into the j th standard (column)
vector EjT using only elementary row operations. The resulting tableau is
the simplex tableau T (B 0 ) with respect to the new basis B 0 .
Example 2.6. (Continuation of Example 2.5). Let us come back to Exam-
ple 2.5:
Since t03 = − 12 , we try to move x3 into the basis. The min ratio rule yields

t25 2
δ = x3 = = = 4.
t23 0.5

Hence, we pivot the last tableau (see above) with the element t23 = 0.5.

1 0 0 -0.5 0.5 1 1 0101 3

T (B) = 0 1 0 -0.5 0.5 1 ∼ 0 1 1 0 1 3 = T (B 0 ).
0 0 1 0.5 0.5 2 0 0211 4
27

In T (B 0 ), all reduced cost values (= t0j , j ∈ {1, ..., m} are ≥ 0: The

corresponding basic solution with
x1 = 3, x3 = 4 , x2 = x4 = 0
is optimal.

If t0j ≥ 0, ∀j ∈ {1, ..., n} and tin+1 ≥ 0, ∀i ∈ { 1, ..., m }, then T (B) is
called an optimal (simplex) tableau. We summarize:

Algorithm (Simplex Method for the Solution of LPs of the type min cT x :
Ax = b, x ≥ 0)

(Input) Basic feasible solution (xB , xN ) with respect to basis B.

(1) Compute the simples tableau T (B).
(2) If t0j ≥ 0 ∀j = 1, . . . , n,
(STOP) (xB , xN ) with xB(i) = tin+1 (i = 1, . . . , m) xN = 0
and with the objective value
−t0n+1 is an optimal solution of LP.
(3) Choose j with t0j < 0.
(4) If tij ≤ 0 ∀i = 1, . . . , m (STOP) The LP is unbounded.
(5) Determine r ∈ {1, . . .n, m} with o
trn+1

trn+1
tij
= min tij
tij > 0
and pivot with trj .
Goto (2).

For the well-definiteness of (1) we state:

Theorem 2.3 (Fundamental Theorem of LP s). Suppose that an LP

 minimize cT x
(LP)st subject to Ax = b

x ≥ 0
is given and that P := { x ∈ Rn | Ax = b, x ≥ 0 } 6= ∅ . Then, there exists
a basic feasible solution.

Proof. cf., e.g., Hamacher, Klamroth, (2000). 

CHAPTER 2. LINEAR PROGRAMMING: FOUNDATIONS AND
28 SIMPLEX METHOD

We conclude our first chapter by an example:

Example 2.7. Solve the LP

 maximize 3x1 + x2

 subject to −x1 + x2 ≤ 1
(LP) x1 − x 2 ≤ 1

 x2 ≤ 2

x1 , x2 ≥ 0

by the Simplex Method.

Solution:
First we transform LP into standard form:

 (-) minimize −3x1 − x2

 subject to −x1 + x2 +x3 = 1
(LP)st x1 − x 2 + x4 = 1

 x 2 + x5 = 2

x1 , x2 , x3 , x4 , x5 ≥ 0

(The ”-” in front of ”min” will be omitted in the following.) As a conse-

quence, the objective function value of the current basic feasible solution is
−(−t06 ) = t06 in the following tableaus (and not −t06 ).

A basic feasible starting solution is given by B = (3, 4, 5), i.e., the slack
variables are the basic variables. The first simplex tableau is therefore (with-
out the O th column):

-3 -1 0 0 0 0
-1 1 1 0 0 1
T (B) =
1 -1 0 1 0 1
01001 2

While applying the Simplex Method, we indicate the pivot column (i.e.,
the column corresponding to a non-basic variable xj with t0j < 0) and the
pivot row (i.e., the row r determined by (5)) by representing the pivot ele-
ment in boldface.
29

0 -4 0 3 0 3 00 0 3 4 11
00110 2 00 1 1 0 2
T (B) ∼ ∼ optimal tableau
1 -1 0 1 0 1 10 0 1 1 3
01001 2 01 0 0 1 2

(STOP) (xB , xN ) with

   
x3 2    
x4 0
xB =  x1  =  3  and xN = = ,
x5 0
x2 2

hence,    
x1 2
 x2   3 
   
x := 

 x3  = 
  2 

 x4   0 
x5 0

is optimal with objective value 11.


With our next chapter we remain in the field of LP, but follow a basically
different approach in solving the linear program.
CHAPTER 2. LINEAR PROGRAMMING: FOUNDATIONS AND
30 SIMPLEX METHOD
Chapter 3

Linear Programming:
Interior-Point Methods

3.1 Introduction
In the 80s of 20th century, it was discovered that many large linear pro-
grams could be solved efficiently by formulating them as nonlinear problems
and solving them with various modifications of nonlinear algorithms such as
Newton’s method. One characteristic of these methods was that they re-
quired to satisfy inequality constraints strictly. So, they soon became called
interior-point methods. By the early 90s of last century, one class —so-called
primal-dual methods— had distinguished itself as the most efficient practical
approach and proved to be a strong competitor to the simplex method on
large problems. These methods are the focus of this chapter.
The simplex method can be quite inefficient on certain problems: The
time required to solve a linear program may be exponential in the problem’s
size, as measured by the number of unknowns and the amount of storage
needed for the problem data. In practice, the simplex method is much more
efficient than this bound would suggest, but its poor worst-case complexity
motivated the development of new algorithms with better guaranteed perfor-
mance. Among them is the ellipsoid method proposed by Khachiyan (1979),
which finds a solution in time that is at worst polynomial in the problem size.
Unfortunately, this method approaches its worst-case bound on all problems
as is not competitive with the simplex method. Karmarkar (1984) announced
a projective algorithm which also has polynomial complexity, but it came
with the added inducement of a good practical behaviour. The initial claims
of excellent performance on large linear programs were never fully borne out,
CHAPTER 3. LINEAR PROGRAMMING: INTERIOR-POINT
32 METHODS

but the announcement prompted a great deal of research activity and a wide
array of methods described by such labels as “affine-scaling ”, “logarithmic-
barrier”, “potential-reduction”, “path-following”, “primal-dual”, and “infea-
sible interior point”. Many of the approaches can be motivated and described
independently of the earlier ones of Karmarkar or of so-called log-barrier, etc..
Interior-point methods share common features that distinguish them from
the simplex method. Each interior-point iteration is expensive to compute
and can make a significant progress towards the solution, while the simplex
method usually requires a large number of inexpensive iterations. The sim-
plex method works its way around the boundary of the feasible polytope,
testing a sequence of vertices in turn until it finds the optimal one. Interior-
point methods approach the boundary of the feasible set in the limit. They
may approach the solution either from the interior or from the exterior of the
feasible region, but they never actually lie on the boundary of this region.

3.2 Primal-Dual Methods

Outline: Let our LP be given in standard form:

 minimize cT x
(LP)st subject to Ax = b

x ≥ 0;

cf. Chapter 2. The dual (linear) problem, LDP, is defined by

 maximize bT y
T
(DP) subject to A y + s = c,

s ≥ 0,

where y ∈ Rm (dual variable) and s ∈ Rn (slack variable). “Primal-dual”

solutions of (LP)st , (DP) fulfill the (first-order) necessary optimality condi-
tions from general optimization theory which are called Karush-Kuhn-Tucker
conditions (cf. Section 1.1). We state them here as follows (exercise):
 T

 A y + s = c,

Ax = b,
(NOC)

 xj sj = 0 ∀j ∈ {1, 2, . . . , n} ,

(x, s) ≥ 0.1

1

The condition (x, s) ≥ 0 just means nonnegativity of all xj , sj .

3.2. PRIMAL-DUAL METHODS 33

Primal-dual methods find solutions (x∗ , y ∗ , s∗ ) of (NOC) by applying vari-

ants of Newton’s method to the three groups of equalities, and modifying
the search directions and step lengths so that the inequalities (x, s) ≥ 0
are satisfied strictly at every iteration. The equations in (NOC) are only
mildly nonlinear and so are not difficult to solve by themselves. However, the
problem becomes more difficult when we add the nonnegativity requirement
(x, s) ≥ 0. The nonnegativity condition is the source of all the complications
in the design and analysis of interior-point methods.
To derive primal-dual interior-point methods, we restate (NOC) in a
slightly different form by a mapping F : R2n+m −→ R2n+m :

  

 AT y + s − c
  =! 0,
F (x, y, s) :=  Ax − b
(NOC)

 XS 1

(x, s) ≥ 0,

where

X := diag(x1 , x2 , . . . , xn ),

S := diag(s1 , s2 , . . . , sn ), and

1 := (1, 1, . . . , 1)T .
Primal-dual methods generate iterates (xk , y k , sk ) that satisfy the bounds
(x, s) ≥ 0 strictly, i.e., xk > 0 and sk > 0. By respecting these bounds,
the interior-point methods avoid spurious solutions, i.e., points that sat-
isfy F (x, y, s) = 0 but not (x, s) ≥ 0. Spurious solutions abound and do
not provide useful information about solutions of (LP)st and (DP). So, it
is reasonable to exclude them altogether from the region of search. Many
interior-point methods actually require (xk , y k , sk ) to be strictly feasible, i.e.,
to satisfy the linear equality constraints for the primal and dual problems.
We put


M := (x, y, s) Ax = b, AT y + s = c, (x, s) ≥ 0 ,


M0 := (x, y, s) Ax = b, AT y + s = c, (x, s) > 0 ;
CHAPTER 3. LINEAR PROGRAMMING: INTERIOR-POINT
34 METHODS

so, M0 is some relative interior of M. Herewith, the strict feasibility condi-

tion can be written concisely as

(SF) (xk , y k , sk ) ∈ M0 ∀k ∈ N0 .

(b) a measure of desirability of each point in search space.

The search direction procedure has its origins in Newton’s method, here,
!
applied to nonlinear equations F = 0 in (NOC). Newton’s method forms a
linear model for F around the current point and obtains the search direction
(∆x, ∆y, ∆s) by solving the following system of linear equations (exercise):
 
∆x
DF (x, y, s)  ∆y  = −F (x, y, s).
∆s
If (x, y, s) ∈ M0 (strict feasibility!), then the Newton step equations become

    
0 AT I ∆x 0
(ASD)sf  A 0 0   ∆y  =  0 .
S 0 X ∆s −XS 1

Usually, a full step along this direction is not permissible, since it would
violate the bound (x.s) ≥ 0. To avoid this difficulty, we perform a line
search along the Newton direction so that the new iterate is

(x, y, s) + α(∆x, ∆y, ∆s)

for some line search parameter α ∈ (0, 1]. Unfortunately, we often can take
only a small step along the direction (i.e., α << 1) before violating the
condition (x, s) > 0. Hence, the pure Newton iteration (ASD), which is
known as the affine scaling direction, often does not allow us to make much
progress towards a solution. Primal-dual methods modify the basic Newton
procedure in two ways:
3.2. PRIMAL-DUAL METHODS 35

(i) They bias the search direction towards the interior of the “nonnegative
orthant” (x, s) ≥ 0, so that we can move further along the direction
before one of the components of (x, s) becomes negative.

(ii) They keep the components of (x, s) from moving “too close” to the
boundary of the nonnegative orthant.

The Central Path (embedding our framework): The central path C

consists of strictly feasible points (xτ , yτ , sτ ) ∈ C parametrized by a scalar
τ > 0, where each point solves the following system:

 T

 A y+s = c,

Ax = b,
(NOC)τsf

 xj s j = τ ∀j ∈ {1, 2, . . . , n} ,

(x, s) > 0.

From (NOC)τsf we can define the central path as C := (xτ , yτ , sτ ) τ > 0 .
In fact, it can be shown that (xτ , yτ , sτ ) is uniquely defined for each τ > 0
if and only if M0 6= ∅. A plot of C for a typical problem, projected into the
space of primal variables x, is shown in Figure 3.1:

Figure 3.1: Central path, projected into x-space, showing a typical neigh-
bourhood N .
CHAPTER 3. LINEAR PROGRAMMING: INTERIOR-POINT
36 METHODS

Another way of defining C is to use our mapping F (see (NOC)) and write

 
0
(NOC)τsf F (xτ , yτ , sτ ) =  0  , (xτ , sτ ) > 0.
τ1
The equations (NOC)τsf approximate the ones in (NOC) more and more as
τ → 0. If C converges to anything as τ → 0, then it must converge to a
primal-dual solution of the linear problem. The central path thus guides us
to a solution along a route that steers clear of spurious solutions by keeping all
x- and s- components strictly positive and decreasing the pairwise products
xj sj to 0 at roughly the same rate.
Primal-dual algorithms take Newton steps towards points on C for which
τ > 0, rather that pure Newton steps for F . Since these steps are biased
towards the interior of the nonnegative orthant (x, s) ≥ 0, it usually is pos-
sible to take longer steps along them than along the pure Newton steps for
F , before violating the positivity condition.
To describe the biased search direction, we introduce a centering param-
eter σ ∈ [0, 1] and a duality measure µ defined by the average

n
1X 1
µ := xj sj = xT s.
n j=1 n
By writing τ := σµ and applying Newton’s method to (NOC)τsf , we obtain

    
0 AT I ∆x 0
(ASD)σµ
sf
 A 0 0   ∆y  =  0 .
S 0 X ∆s −XS 1 + σµ1
The step (∆x, ∆y, ∆s) is a Newton step towards the point (xσµ , yσµ , sσµ ) ∈ C,
at which the pairwise products xj sj are equals to σµ.
If σ = 1, then (ASD)τsf defines a centering direction, a Newton step
towards the point (xµ , yµ , sµ ) ∈ C, at which xj sj = µ for all j. Centering
directions are usually biased strongly towards the interior of the nonnegative
orthant and make little progress in reducing the duality measure µ. By
moving closer to C, however, they set the scene for substantial progress on
the next iteration.
If σ = 0, then we get the standard Newton step (ASD)sf . many algo-
rithms use intermediate values of σ ∈ (0, 1) to trade off between the twin
3.2. PRIMAL-DUAL METHODS 37

A Primal-Dual Framework: With our concepts from above, we can

define a general framework for primal-dual algorithms:

Framework IPM1 Primal-Dual

Given (x0 , y 0 , s0 ) ∈ M0
for k = 0, 1, 2, . . .
Solve     
0 AT I ∆xk 0
(ASD)σsfk µk  A 0 0   ∆y k  =  0 ,
k k k k k
S 0 X ∆s −X S 1 + σk µk 1
k T k
where σk ∈ [0, 1] and µk = (x ) s /n;
Set
(xk+1 , y k+1 , sk+1 ) = (xk , y k , sk ) + αk (∆xk , ∆y k , ∆sk ),
choosing αk such that (xk+1 , sk+1 ) > 0.
end(for).

The choices of centering parameter σk and step length αk are crucial to

the performance of the method. Techniques for controlling these parameters
give rise to a wide variety of methods with varying theoretical properties.
So far, we have assumed that the starting point (x0 , y 0 , s0 ) is strictly
feasible and, in particular, that it satisfies the linear equations Ax0 = b,
AT y 0 + s0 = c. All subsequent iterates also respect these constraints, because
σk µ k
of the zero right-hand-side terms in (ASD)sf . However, for most problems,
a strictly feasible starting point is difficult to find! Infeasible-interior-point
methods require only that the components of x0 and s0 be strictly positive.
The search direction needs to be modified so that it improves feasibility as
well as centrality at each iteration, but this requirement entails only a slight
change to the step equation (ASD)σµ sf .

If we define the residuals for the two linear equations as

rb := Ax − b, rc := AT y + s − c,
then the modified step equation is

    
0 AT I ∆x −rc
(ASD)σµ
inf
 A 0 0   ∆y  =  −rb .
S 0 X ∆s −XS 1 + σµ1
CHAPTER 3. LINEAR PROGRAMMING: INTERIOR-POINT
38 METHODS

The search direction is still a Newton step towards the point (xσµ , yσµ , sσµ ) ∈
C. It tries to correct all the infeasibility in the equality constraints in a single
step. If a full step is taken at any iteration (i.e., ∃k : αk = 1), the residuals
rb and rc become 0, and all subsequent iterates remain strictly feasible.

Path-Following Methods: These algorithms explicitly restrict the it-

erates to a neighbourhood of the central path C and follow C to a solution
of our LP. By preventing the iterates from coming too close to the bound-
ary of the nonnegative orthant, they ensure that search directions calculated
from each iterate make at least some minimal amount of progress towards
the solution.
In path-following algorithms, the duality measure µ fills the role (to be
satisfied by any optimization algorithm) of (b) introduced above. The duality
measure µk is forced to 0, µk → 0 (k → ∞), so the iterates (xk , y k , sk ) come
closer and closer to satisfy the Karsh-Kuhn-Tucker conditions (NOC)!
There are two most interesting neighbourhoods of C:


(α) N2 (θ) := (x, y, s) ∈ M0 || XS 1 − µ1 ||2 ≤ θµ :

2-norm-neighbourhood, of some θ ∈ [0, 1), || · ||2 : Euclidean norm,

 

(β) N−∞ (γ) := (x, y, s) ∈ M
0
max xj sj ≥ γµ :
j∈{1,2,...,n}

one-sided ∞- (or max-) norm-neighbourhood, of some γ ∈ (0, 1]. Typical

1
values are θ = 0.5 and γ = 1,000 . If a point lies in N−∞ (γ), each xj sj must be
at least some small multiple γ of their average value µ. This requirement is
quite modest, and we can make N−∞ (γ) encompass most of the feasible region
M by choosing γ close to 0. The neighbourhood N2 (θ) is more restrictive,
since certain points in M0 do not belong to N2 (θ), no matter how close θ is
chosen to its upper bound of 1.
The projection of neighbourhood N onto the x-space for a typical problem
is shown as the region between the dotted lines in Figure 3.1.
By keeping all iterates inside one or another of these neighbourhoods,
path-following methods reduce all the xj sj to 0 at more or less the same
rate. Our following so-called long-step path-following algorithm can make
rapid progress because of its use of wide neighbourhood N−∞ (γ), for γ ≈ 0.
It depends on two parameters σmin and σmax , which are upper and lower
3.2. PRIMAL-DUAL METHODS 39

bounds on the centering parameter σk . As usual, the search direction is

obtained by solving (ASD)σsfk µk , and we choose the step length αk to be as
large as possible, subject to staying inside N−∞ (γ).
Here and later on, we use the notation

(xk (α), y k (α), sk (α)) := (xk , y k , sk ) + α(∆xk , ∆y k , ∆sk ),

µk (α) := n1 xk (α)sk (α)
(cf. our corresponding definitions from above).

Algorithm IPM2 (Long-Step Path-Following)

Given γ, σmin , σmax with γ ∈ (0, 1), 0 < σmin < σmax < 1,
and (x0 , y 0 , s0 ) ∈ N−∞ (γ);
for k = 0, 1, 2, . . .
Choose σk ∈ [σmin , σmax ];
σk µ k
Solve (ASD)sf to obtain (∆xk , ∆y k , ∆sk );
Choose αk as the largest value of α in [0, 1] such that

Set (xk+1 , y k+1 , sk+1 ) = (xk (αk ), y k (αk ), sk (αk ));

end(for).

The typical behaviour of our algorithm is illustrated in Figure 3.2 for the
case of n = 2. The horizontal and vertical axes in this figure represent the
pairwise product x1 s1 and x2 s2 ; so the central path C is the line emanating
from the origin at an angle of 45o . (Here, a point at the origin is a primal-
dual solution if it also satisfies the feasibility conditions from (NOC) different
from xj sj = 0.) In the usual geometry of Figure 3.2, the search directions
(∆xk , ∆y k , ∆sk ) transform to curves rather than lines. As Figure 3.2 shows,
the bound σmin ensures that each search direction stands out by moving
away from the boundary of N−∞ (γ) and into the relative interior of this
neighbourhood. That is, small steps along the search direction improve the
centrality.
Larger values of α take us outside the neighbourhood again, since the error
in approximating the nonlinear system (NOC)τsf by the linear step equations
(ASD)σµsf becomes more pronounced as α increases. Still, we are guaranteed
that a certain minimum step can be taken before we reach the boundary of
N−∞ (γ). A complete analysis of Algorithm IPM2 can be found in Nocedal,
CHAPTER 3. LINEAR PROGRAMMING: INTERIOR-POINT
40 METHODS

Wright (1999). It shows that the primal-dual methods can be understood

without recourse to profound mathematics. This algorithm is fairly efficient
in practice, but with a few more changes it becames the basis of a truly
competitive method.
An infeasible-interior-point variant of Algorithm IPM2 can be constructed
by generalizing the definition of N−∞ (γ) to allow violation of the feasibility
conditions. In this extended neighbourhood, the residual norms || rb || and
|| rc || are bounded by a constant multiple of the duality measure µ. By
squeezing µ to 0, we also force rb and rc to 0, so that the iterates approach
complementarity and feasibility simultaneously.

3.3 A Practical Primal-Dual Algorithm

Most existing interior-point codes are based on a predictor-corrector algo-
rithm proposed by Mehrotra (1992). The two key features of this algorithm
are:

(I) addition of a corrector step to the search direction of Framework IPM1,

3.3. A PRACTICAL PRIMAL-DUAL ALGORITHM 41

so that the algorithm more closely follows a trajectory to the primal-

dual solution set;

(II) adaptive choice of the centring parameter σ.

In (I), we shift the central path C so that it starts at our current iterate
(x, y, s) and ends at the set Ω of solution points. Let us denote this modified
trajectory by H and parametrize it by the parameter τ ∈ [0, 1], so that


x(τ ), yb(τ ), sb(τ )) τ ∈ [0, 1) ,
H = (b
with x
b(0), yb(0), sb(0)) = (x, y, s) and

lim (b
τ →1−

(see Figure 3.3).

Figure 3.3: Central path C, and a trajectory H from the current (noncentral)
point (x, y, s) to the solution set Ω.

Algorithms from Framework IPM1 are first-order methods, in that they

find the tangent to a trajectory like H and perform line search along it.
This tangent is knows as the predictor step. Mehrotra’s algorithm takes the
next step of calculating the curvature of H at the current point, thereby
obtaining a second-order approximation of this trajectory. The curvature is
used to define the corrector step. It can be obtained at a low marginal cost,
since it reuses the matrix factors from the calculation of the predictor step.
CHAPTER 3. LINEAR PROGRAMMING: INTERIOR-POINT
42 METHODS

In (II), we chose the centering parameter σ adaptively, in contrast to algo-

rithms from Framework IPM1, which assign a value to σk prior to calculating
the search direction. At each iteration, Mehrotra’s algorithm first calculates
the affine-scaling direction (the predictor step) and assesses its usefulness as
a search direction. If this direction yields a large reduction in µ without
violating the positivity condition (x, s) > 0, then the algorithm concludes
that little centering is needed. So, it chooses σ close to zero and calculates a
centered search direction with this small value. If the affine-scaling direction
is not so productive the algorithm enforces a larger amount of centering by
choosing a value of σ closer to 1.
The algorithm thus combines three steps to form the research direction:
a predictor step which allows us to determine the centering parameter σk , a
corrector step using second-order information of the path H leading towards
the solution, and centering step in which the chosen value of σk is substituted
in (ASD)σµinf . Hereby, the computation of the centered direction and the
corrector step can be combined, so that adaptive centering does not add
further to the cost of each iteration.
For the computation of the search direction (∆x, ∆y, ∆s) we proceed as
follows:
First, we calculate the predictor step (∆xaf f , ∆y af f , ∆saf f ) by setting
σ = 0 in (ASD)σµ
inf , i.e.,

    
0 AT I ∆xaf f −rc
0,af f  A 0 0   ∆y af f  =  −rb  .
(ASD)inf
S 0 X ∆saf f −XS 1
pri dual
To measure the effectiveness of this direction, we find αaf f and αaf f to the
largest length that can be taken along this direction before violating the
nonnegativity conditions (x, s) ≥ 0, and an upper bound of 1:
( )
pri xj
αaf f := min 1, min − af f ,
j:∆xaf
j
f
<0 ∆xj
( )
dual sj
αaf f := min 1, min − .
j:∆saf
j
f
<0 ∆saf
j
f

Furthermore, we define µaf f to be the value of µ that would be obtained by

a full step to the boundary, i.e.,
3.3. A PRACTICAL PRIMAL-DUAL ALGORITHM 43

1 pri af f
 pri af f

µaf f :=x + αaf f ∆x s + αaf f ∆s ,
n
and set the damping parameter σ to be
 3
µaf f
σ := .
µ
When good progress is made along the predictor direction, then we have
µaf f << µ, so this σ is small; and conversely.
0,af f
The corrector step is obtained by replacing the right-hand-side of (ASD)inf
by (0, 0, −∆X af f , ∆S af f 1), while the centering step requires a right-hand-
side of (0, 0, σµ1). So we can obtain the complete Mehrotra step, which
includes the predictor, corrector and centering step components, by adding
the right-hand-sides of these three components and solving the following sys-
tem:

    
0 AT I ∆x −rc
(⊕)  A 0 0   ∆y  =  −rb .
af f af f
S 0 X ∆s −XS 1 − ∆X ∆S 1 + σµ1

We calculate the maximum steps that can be taken along these directions
before violating the nonnegativity condition (x, s) > 0 by formulae similar
pri dual
to the ones for αaf f and αaf f ; namely,

( )
pri
xkj
αmax := min 1, min − ,
j:∆xj <0 ∆xj
( )
dual
skj
αmax := min 1, min − ,
j:∆sj <0 ∆sj
and then choose the primal and dual step lengths as follows


αkpri := min 1, ηαmax
pri
,


αkdual := min 1, ηαmax
dual
,
where η ∈ [0.9, 1.0) is chosen so that η → 1 near the solution, to accelerate
the asymptotic convergence.
CHAPTER 3. LINEAR PROGRAMMING: INTERIOR-POINT
44 METHODS

usual format.

Given (x0 , y 0 , s0 ) with (x0 , s0 ) > 0;

for k = 0, 1, . . .
Set (x, y, s) = (xk , y k , sk ) and
0,af f
solve (ASD)inf for (∆xaf f , ∆y af f , ∆saf f );
pri dual
Calculate αaf f , αaf f , and µaf f ;
Set centering parameter to σ = (µaf f /µ)3 ;
Solve (⊕) for (∆x, ∆y, ∆s);
Calculate αkpri and αkdual ;
Set
xk+1 = xk + αkpri ∆x,
(y k+1 , sk+1 ) = (y k , sk ) + αkdual (∆y, ∆s);
end(for).

It is important to note that no convergence theory is available for Mehro-

tra’s algorithm in its previous form. In fact, there are examples for which
the algorithm diverges. Simple safeguards could be incorporated into the
method to force it into the convergence framework of existing methods. How-
ever, most programs do not implement these safeguards, because the good
practical performance of Mehrotra’s algorithm makes them unnecessary.
to Nocedal, Wright (1999) and to Nash, Sofer (1996).
Chapter 4

Nonlinear Programming:
Feasible-Point Methods

In this chapter, we examine methods which solve constrained (nonlinear) op-

timization problems by attempting to remain feasible at every iteration. If
all the constraints are linear, maintaining feasibility is straightforward. We
discuss this case first. When nonlinear constraints are present, then more
elaborate procedures are required. We discuss two such approaches: se-
quential quadratic programming and reduced-gradient methods. Both these
approaches generalize the techniques for linear constraints. Although they
are motivated by the idea of maintaining feasibility at every iteration, the do
not always achieve this.

4.1 Linear Equality Constraints

The majority of methods for solving problems with linear equality constraints
are feasible-point methods: They start from a feasible point and move along
feasible descent directions to consecutively better feasible points. There are
two features making this approach particularly attractive:

(i) It is practically advantageous that all iterates are feasible. Even if the
algorithm fails to solve the problem to the desired accuracy, it might
still provide a feasible solution that is usable.

(ii) By restricting movement to feasible directions, the equality-constrained

problem is transformed to an unconstrained one in the null space of the
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
46 METHODS

constraints. This new problem may then be solved using unconstrained

minimization techniques.

Let us write our problem as follows:


minimize f (x)
(NLP)
subject to Ax = b,

where f is twice-continuously differentiable: f ∈ C 2 , and A is an m×n matrix

of full row rank, where m ≤ n. As in the constrained case, the methods we
describe are only guaranteed to find a stationary point of the problem. In
the special case where f is convex, this point will be a global minimizer of f .
Let x be a feasible point, i.e., Ax = b. Since any other feasible point can
be reached from x by moving in a feasible direction, the solution to (NLP)
can be written as x∗ = x + p, where p solves the problem

minimize f (x + p)
subject to Ap = 0,
(Note: Ap = A(x∗ − x) = Ax∗ − Ax = b − b = 0).
Let Z denote an n × (n − m) basis matrix for the null-space (kernel) of
A. Then, p = Zv for some (n − m)-dimensional vector v. This problem is
equivalent to the unconstrained problem

minimize ϕ(v) := f (x + Zv).

So, we have reduced the problem of finding the best n-dimensional vector p
to the unconstrained problem of finding the best (n − m)-dimensional vector
v.
Conceptually, it is possible to minimize the reduced function ϕ using any
of the unconstrained methods. In practice it is not necessary to provide an
explicit expression for ϕ(v). Instead, it is possible to work directly with the
original variable x, using

(∗∗) ∇2 ϕ(v) = Z T ∇2 f (x)Z,

4.1. LINEAR EQUALITY CONSTRAINTS 47

where x = x + Zv. Let us shortly demonstrate this:

1
ϕ(v) = f (x + Zv) = f (x) + v T Z T ∇f (x) + v T Z T ∇2 f (x)Zv + · · · .
2
Ideally, we would like to use the Newton direction; it is obtained by mini-
mizing the quadratic approximation to ϕ(v) obtained from the Taylor series.
Setting the gradient of the quadratic approximation to 0 gives the following
linear system in v:


Z T ∇2 f (x)Z v = −Z T ∇f (x),
called the reduced Newton (or null-space) equation. Its solution is just

−1
v = − Z T ∇2 f (x)Z Z T ∇f (x),
being an estimate of the best. In turn, it provides an estimate of the best p:

−1
p = Zv = −Z Z T ∇2 f (x)Z Z T ∇f (x),
called the reduced Newton direction at x. Now, we can derive the equality-
constrained analog of the classical Newton method. The method sets

xk+1 := xk + pk ,
 −1
where pk = −Z Z T ∇2 f (xk )Z Z T ∇f (xk ) is the reduced Newton direction
at xk . This is just the mathematical formulation of the method, and in
practice, explicit inverses are not normally computed. The method does not
require that the reduced function be formed explicitly.
Example 4.1. We consider the problem:

 minimize f (x) = 12 x21 − 21 x23 + 4x1 x2 + 3x1 x3 − 2x2 x3
(NLP)= subject to x1 − x2 − x3 = −1

(i.e., A := (1, −1, −1), b = (−1)).

As a basis matrix for the null-space of A we choose

 
1 1
Z =  1 0 .
0 1
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
48 METHODS

The reduced gradient at x := (1, 1, 1)T , being the feasible point which we
consider, is
 
  8  
1 1 0  2  = 10 ,
Z T ∇f (x) =
1 0 1 8
0
and the reduced Hessian matrix at x is

  
  1 4 3 1 1  
1 1 0  4 0 −2   1 0  = 9 6
Z T ∇2 f (x)Z = .
1 0 1 6 6
3 −2 −1 0 1

The reduced Newton equation yields

 
− 32
v= ;
− 23
hence, the reduced Newton direction is

− 43
p = Zv =  − 23  .
− 23
Since the objective function is quadratic and the reduced Hessian matrix
is positive definite, a step length of α = 1 leads to the optimum x∗ =
T
− 31 , 13 , 13 . At x∗ , the reduced gradient is Z T ∇f (x∗ ) = Z T (2, −2, 2)T = 0
as expected. The corresponding Lagrange multiplier is λ∗ = 2, because

∇f (x∗ ) = 2AT .


The reduced Newton direction is invariant with respect to the null-space

matrix Z: Any choice of the basis matrix Z will yield the same search di-
rection p. Numerically, however, the choice of Z can have dramatic effect on
the computation.
The classical reduced Newton method has all the properties of the clas-
sical Newton method. In particular, if the reduced Hessian matrix at the
solution is positive definite, and if the starting point is sufficiently close to
4.1. LINEAR EQUALITY CONSTRAINTS 49

the solution, then the iterates will converge quadratically. In the more gen-
eral case, however, the method may diverge or fail. Then, some globalization
strategy should be used.
If xk is not a local solution and if the reduced Hessian matrix is positive
definite, then the reduced Newton direction is a descent direction, since
−1
pT ∇f (xk ) = −∇f (xk )Z Z T ∇2 f (xk )Z Z T ∇f (xk ),
< 0.

If the Hessian matrix is not positive definite, then the search direction
may not be a descent direction, and worse still, it may not be defined. Then,
the modified factorizations can be applied to the reduced Hessian matrix to
provide a descent direction; see Nash, Sofer (1996).
Other compromises on Newton’s method may be made to obtain cheaper
iterations. The simplest of all methods is of course the steepest-descent
method. For the reduced function, this strategy gives the direction

v = −Z T ∇f (xk )
in the reduced space, which yields the reduced steepest-descent direction

p = −ZZ T ∇f (xk )
in the original space. Here, Z may be any null-space matrix for A. However,
the direction will vary with the particular choice of Z, unlike the reduced
Newton direction!
Example 4.2. The reduced gradient at the initial point of Example 4.1 is
Z T ∇f (x) = (10, 8)T ; hence, the reduced steepest-descent direction is p =
−ZZ T ∇f (x) = (−18, −10, −8)T . Had we chosen
 
2 0
b= 1 4 
Z
1 −4
as the null-space matrix for A, the reduced gradient would be ZbT ∇f (x) =
bT (8, 2, 0)T = (18, 8)T , and the reduced steepest-descent direction would be
Z
p = −ZbZ bT ∇f (x) = (−36, −50, 14)T . 
The reduced steepest-descent method has the same properties as that of
its classical unconstrained counterpart. The iterations are cheap, but con-
vergence may be very slow. Quasi-Newton methods are a more sophisticated
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
50 METHODS

compromise. A common approach is to construct an approximation Bk to

the reduced Hessian matrix at xk . Let Z be a basis matrix for the null-space
of A. The search direction is computed as p = Zv, where v is obtained by
solving

Bk v = −Z T ∇f (xk ).

The approximation Bk is updated in much the same way as in the un-

constrained case, except that all quantities are in the reduced space. For
example, the symmetric rank-one update formula becomes

(y k − Bk sk )(y k − Bk sk )T
Bk+1 = Bk + ,
(y k − Bk sk )T sk

where y k = Z T ∇f (xk+1 ) − ∇f (xk ) and sk = Z T (xk+1 − xk ).
In principle, it is possible to solve an equality-constrained problem using
any classical technique from the unconstrained case. In practice, the nu-
merical difficulties encountered when solving equality-constrained problems
are not quite the same as those encountered when solving unconstrained
problems, and it is not always possible to solve a large equality-constrained
problem by simple applying general-purpose software for unconstrained op-
timization. One reason is that the reduced Hessian matrix in a constrained
problem is often different in structure from the Hessian matrix in an uncon-
strained minimization problem:
Example 4.3. Consider the quadratic problem


minimize f (x) = 21 x21 + x22 + 2x23 + 4x24
(NLP)=
subject to x1 + x2 + x3 + x4 = 1.

Taking

 1 1 1

2 2 2
 − 12 − 12 
1
Z=
 −
2 
2
1 1
2
− 12 
− 21 − 12 12

as a basis for the null-space of the constrained matrix A = (1 1 1 1), we

obtain
4.1. LINEAR EQUALITY CONSTRAINTS 51

 
1 0 0 0
 0 2 0 0 
∇2 f (x) = 
 0
, and
0 4 0 
0 0 0 8
 15 3

4 4
− 54
Z T ∇2 f (x)Z =  3
4
15
4
− 94  .
− 45 − 94 15
4

Thus, the reduced matrix Z T ∇2 f (x)Z is dense, even though the matrix
∇2 f (x) is sparse. The special diagonal structure of the Hessian matrix is
destroyed by the reduction! 

To overcome these problems, special implementations that are tailored

to the equality-constrained problem may be needed. For example, if a
conjugate-gradient method is used to solve the reduced Newton equation,
then it is not necessary to form the reduced Hessian matrix explicitly.
Once an optimal solution to the equality-constrained problem is obtained,
the associated vector of Lagrange multipliers is computed. There are several
reasons for this.

(i) The Lagrange multipliers measure the sensitivity of the solution to

changes in the constraints.

(ii) The equality-constrained problem could be one of the sequence of prob-

lems generated by an algorithm for solving a problem with inequality
constraints. In this case, the Lagrange multipliers indicate how to im-
prove the current solution.

The optimality conditions for (NLP)= can be used directly to derive al-
gorithms. The conditions are (see Section 1.1):


∇f (x) − AT λ = 0, and
(NOC)=
b − Ax = 0,

where λ is the vector of Lagrange multipliers. If Newton’s method is applied

to this nonlinear system, then

xk+1 = xk + pk ,
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
52 METHODS

λk+1 = λk + νk ,

where updates pk and νk are the solutions to the Newton equations

   
∇2 f (xk ) −AT p AT λk − ∇f (xk )
.
−A 0 ν Axk − b

Some algorithms for constrained optimization work directly with this lin-
ear system, although some care must be taken to ensure that descent di-
rections are obtained. This linear system is closely related to the reduced
Newton equation derived above.

4.2 Computing the Lagrange Multipliers λ

We still consider the linear equality-constrained problem (NLP)= from Sec-
tion 4.1. Assume, that the linear independence constraint qualification holds,
i.e., the regularity condition saying that the rows of A are linearly indepen-
dent. Consider now the optimality condition

This is a system of n equations in m ≤ n unknowns, so it cannot always be

expected to have a solution. At most feasible points x∗ , this over-determined
system will be inconsistent, but if x∗ is a local solution of the optimization
problem, then the system will have a solution. How can such a solution λ∗
be computed?
A useful tool is a matrix known as the right inverse. We define an n × m
matrix Ar to be a right inverse for the m × n matrix A, if

AAr = Im .

It is easy to see that a matrix A has a right inverse only if it has full
row rank. In this case, and if m = n, then the right inverse is unique, and
Ar = A−1 . If m < n, the right inverse is generally not unique. For example,
the matrices
4.2. COMPUTING THE LAGRANGE MULTIPLIERS λ 53

   
0.75 0 1 0
 −0.25 0   0 0 
  and  
 0 0.5   0 1 
0 0.5 0 0
are both right inverses for the matrix
 
1 −1 0 0
A= .
0 0 1 1

To see how right inverses are of use in solving the system ∇f (x∗ ) = AT λ∗ ,
suppose that a solution to this system exists. If both sides of the equation
are multiplied by ATr , then we obtain

λ∗ = ATr ∇f (x∗ ).
If the system ∇f (x∗ ) = AT λ is consistent, its solution λ∗ = ATr ∇f (x∗ )
is unique, even though the right inverse may not be unique. (Note, that
∇f (x∗ ) = AT λ implies that AAT λ = A∇f (x∗ ), so the unique solution is
λ∗ = (AAT )−1 A∇f (x∗ ). In fact, if A has full row rank, then AAT is positive
definite and, hence, its inverse exists.)
The linear system ∇f (x∗ ) = AT λ is consistent if and only if ∇f (x∗ ) is
a linear combination of the rows of A. Hence, λ∗ = ATr ∇f (x∗ ) will be a
solution to the system if and only if


I − AATr ∇f (x∗ ) = 0.

In practice, we will almost never find a point x∗ that satisfies the optimal-
ity conditions to within some specified tolerance. The point xk will be an es-
timate of the optimal solution. Correspondingly, the vector λk = ATr ∇f (xk )
will only be an estimate of the vector of Lagrange multipliers at the solu-
tion. It is sometimes called a first-order estimate, because for sufficiently
∗ ||
small , if || xk − x∗ || = O() (i.e., || xk −x

is bounded under  → 0) then

|| λk − λ || = O() also.
The rest of this section is an excursion about computing a right inverse
matrix.
Variable Reduction Method: Here, the variables are partitioned into
m basic and n − m nonbasic variables. The matrix A is partitioned into basic
and nonbasic columns correspondingly. Assuming that the first m columns
are basic, then,
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
54 METHODS

A = ( B N ),
where B is an m × n nonsingular matrix, and the n × (n − m) matrix
 
−B −1 N
Z :=
I
is a basis matrix for the null space of A. The matrix
 
B −1
Ar :=
0
is a right-inverse matrix for A that is available with no additional computa-
tion.
Orthogonal Projection Matrix: Let the n × n matrix

P := I − AT (AAT )−1 A
be the orthogonal projection matrix into the null space of A. A right inverse
for A associated with the orthogonal projection is the matrix

Ar := AT (AAT )−1 .
This matrix, which we will denote by A+ , is a special right inverse. It satisfies
the following four conditions:

A+ AA+ = A+ , (A+ A)T = A+ A.

It can be shown that, for any m × n matrix A, there is a unique n × m
matrix A+ that satisfies these conditions. We call A+ be the Penrose-Moore
generalized inverse of A. If A has full row rank, then,

A+ = AT (AAT )−1 ,
and if A has full columns rank, then,

A+ = (AT A)−1 AT .
Formulas for A+ can also be developed when A does not have full row or
column rank; cf. Golub, Van Loan (1989).
4.2. COMPUTING THE LAGRANGE MULTIPLIERS λ 55

Given a point xk , the vector of Lagrange multipliers estimates (A+ )T ∇f (xk )

obtained from the Penrose-Moore generalized inverse has the appealing prop-
erty that it solves the problem

minimize ∇f (xk ) − AT λ .
m
λ∈R 2

For this reason it is called the least-squares Lagrange multiplier estimate at

xk .
As the condition number of AAT is the square of the condition number of
A, the computation of (AAT )−1 is potentially unstable. The QR factorization
provides a stable approach to computing this matrix that is practical for
smaller problems.
Non-Orthogonal Projection: Let D be a positive n × n matrix. Then,
the n × n projection matrix

PD := I − DAT (ADAT )−1 A

is a null space matrix for A. A right inverse for A associated with this
projection is

QR Factorization (known from numerical mathematics): The QR fac-

torization represents AT as a product of an orthogonal matrix Q and an
upper triangular matrix R. Denoting the first m columns of Q by Q1 and
the last n − m columns by Q2 , we have
 
T R1
A = QR = (Q1 Q2 ) ,
0
where R1 is an m × m triangular matrix. The n × (n − m) matrix

Z := Q2
is an orthogonal basis for the null space of A. The matrix


Ar := Q1 R1−T , R1−T := (R1−1 )T
is a right inverse for A available from the QR factorization at little additional
cost. In fact, this matrix need not be formed explicitly: A computation of
the form λk = ATr ∇f (xk ) may be done by first computing y1 = QT1 ∇f (xk )
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
56 METHODS

and then solving the triangular system R1 λ = y1 . It is easy to show that

Ar = AT (AAT )−1 ; hence, this right-inverse is in fact the Penrose-Moore
generalized inverse of A.
Example 4.4. Let us construct two right inverses for
 
1 −1 0 0
A := .
0 0 1 1

If variable reduction is used, with columns 2 and 3 of A being selected as

the basic columns, then
   
−1 0 1 0
B= , N= .
0 1 0 1
From these we determine that
   
1 0 0 0
 1 0   −1 0 
Z= 
 0 −1  and Ar = 
 0
.
1 
0 1 0 0
If a QR factorization of AT is used, then (exercise)
 
− √12 0 − 12 − 12
 √1 0 − 12 − 12 
 2 
Q=  and
 0 − √12 12 − 12 
0 − √12 − 12 12
√  
− 2 0

 0 − 2 
R=
 0
.

0
0 0

 1 
− 2 − 12
 −1 −1 
Z= 2
 1 −1  .
2 
2 2
− 12 12

The right inverse is obtained from the formula

4.3. LINEAR INEQUALITY CONSTRAINTS 57

Ar = Q1 R1−T ,
where Q1 consists of the first two columns of Q:
 
− √12 0
 √1 0 
 2 
Q1 =  
 0 − √12 
0 − √12
and R1 consists of the first two rows of R:
 √ 
− 2 0

R1 = .
0 − 2
Hence,
 1

2
0
 −1 0 
Ar = 
 0
2 
1 .
2
1
0 2

If

∇f (x) = (7, −7, −2, −2)T ,

then, for the above right inverses
 
7
λ= ATr ∇f (x) = .
−2
No matter which right inverse has been used, the same values of the Lagrange
multipliers are obtained. 

4.3 Linear Inequality Constraints

In this section, we discuss methods for solving problems with linear inequal-
ities. The problem is now written in the following form:


minimize f (x)
(NLP)A
>
subject to Ax ≥ b,
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
58 METHODS

where f is still a C 2 function. For simplicity, we assume that the problem

has only inequality constraints. (The extension to the case with equality
constraints additionally is easy, and mentioned later on).
Suppose that the point x∗ is a local solution to this problem. Let A b be

a matrix whose rows are the coefficients of the active constraints at x , and
let Z be a null-space matrix for A.b The first-order optimality condition of
Karush-Kuhn-Tucker state that there exists a vector µ∗ such that


bT µ∗
∇f (x∗ ) = A (⇔ Z T ∇f (x∗ ) = 0),
(NOC)≥
µ∗ ≥ 0 (and Ax∗ ≥ b).

Problems that have inequality constraints are significantly more difficult

to solve than problems in which all constraints are equations. The reason
is that it is not known in advance which inequality constraints are active
at the solution. If we know a priori the correct active set, then we could
ignore the inactive constraints and minimize the objective function with all
active inequalities (i.e., where “=” holds) treated as equalities. In practice,
unfortunately, we do not know what the active set is.
How can we resolve this combinatorial issue? A brute force approach
would be to solve the equality constrained problem for all possible selections
of active constraints, and then choose the best solution. Even for a small
problem, however, the number of such subproblems is enormous, and the
amount of work could be prohibitive!
Active-set methods attempt to overcome this difficulty by moving se-
quentially from one choice of active constraints to another choice that is
guaranteed to produce at least a good solution. The hope is that only a
fraction of the potential subproblems will be considered.
The most commonly used active-set methods are feasible-point methods.
An initial feasible point, if none is provided, can be obtained much as in
linear programming; see Section 1.1 and Nash, Sofer (1996).
At each iteration of the active set method, we select a working set J w of
constraints that are assumed to be active at the optimum. We attempt to
minimize f with all constraints in the working set as equalities. All other
constraints are considered inactive and temporarily ignored.
In general, the working set at the current point x is a subset of the
constraints that are active at x, so that x is a feasible point for the working
set. There may also be constraints that are active at x but that are not
included in the working set; hence, the working set is not necessarily equal
4.3. LINEAR INEQUALITY CONSTRAINTS 59

to the active set.

Example 4.5.

 minimize f (x)

 subject to x 1 + x2 + x3 ≥ 1
(NLP)A
> x1 ≥ 0

 x2 ≥ 0

x3 ≥ 0.

At xa = (0, 21 , 12 )T , the first two inequalities are active. If both constraints

are chosen for the working set, i.e., J w := {1, 2}, then we will attempt to
minimize f on the set


x ∈ R3 x1 + x2 + x3 = 1, x1 = 0 .
If only the first constraint is selected to be in the working set, then, we will
attempt to minimize f on the set


x ∈ R 3 x1 + x 2 + x 3 = 1 .


The problem of minimizing f subject to the constraints defined by a

working set is an equality-constrained problem. Therefore, we can use any
of the techniques describe in Section 4.1 to obtain a feasible direction p. We
could then use this search direction within a line search method, and find an
appropriate step length α. However, in our problem (NLP)A > with inequality
constraints, an acceptable step might lead to a point that is infeasible, i.e.,
a point that violates one or more of the constraints that we have ignored.
Geometrically, as we move from x along p, we may encounter the boundary
of some constraint for some step length α. The value α represents the largest
possible step that may be taken without violating feasibility. The step length
must never exceed α.
It is possible, that the best acceptable step length that does not exceed
α is α itself. This step leads to a point on the boundary of the constraint.
The constraint encountered is now satisfied exactly at the new point, and
it is added to the working set. With the step length determined, and any
necessary adjustments made to the working set, the entire process is repeated.
It is also possible, that the best acceptable step length is 0. This is an
exceptional case, and occurs when α = 0. If this occurs, no step is taken and
a constraint is added to the working set (since α = α = 0).
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
60 METHODS

Suppose, now, that we have found a point x that minimizes f on a given

working set. Then, the first-order optimality conditions for the equality-
constrained problem are satisfied at x, and we can compute the Lagrange
multipliers corresponding to the constraints in the working set. (The La-
grange multipliers for all other constraints are assumed to be 0.) If the
Lagrange multipliers are just all nonnegative, then x is also a solution to the
original inequality-constrained problem and the problem is solved. However,
if some Lagrange multiplier is negative, then x is not an optimal solution.
The negative multiplier indicates that the function can be decreased if we
move away from the corresponding constraint into the interior of the feasible
region. Hence, we can drop this constraint from the working set. Now, we
have a new working set, and the process is repeated. (For simplicity, we al-
low only one constraint to enter or leave the working set at a time, although
alternatives are possible.)
Example 4.6. Figure 4.1 illustrates a possible sequence of movements in an
active-set method for minimizing a convex function f on the box

3
!
Y
[0, 1]3 = [0, 1] ⊆ R3 .
j=1

The directions of movement in this example are arbitrary feasible descent

directions and do not correspond to any specific method.
Let A be the starting point. The initial working set is chosen as the active
constraints at A. It consists of the upper bound constraint on x3 , which is
the single constraint that is active at A. Suppose that the first search ends
with B. If B is still not a minimizer of f on x3 = 1, then another search
is made, starting from B; etc.. Assume, finally, that E is a local minimizer
of f along the constraints x1 = 1 and x2 = 1. if the Lagrange multipliers
corresponding to these constraints are both nonnegative, then, E is a local
minimizer for the original problem, and the algorithm is terminated. 
Now, let us discuss some of the above ideas in further detail. We as-
sume that we have a point x which is feasible for the inequality constrained
problem. Denote the working set at x by J := J w . Furthermore, denote the
coefficient matrix for the constraints in J by A (= Aw ) and the corresponding
right-hand-side vector by b (= bw ). Let Z be a null-space matrix for A. The
equality constrained problem for J reads as follows:


A minimize f (x)
(NLP)>
subject to Ax ≥ b.
4.3. LINEAR INEQUALITY CONSTRAINTS 61

This problem is commonly solved using some feasible direction method.

Thus, if p is in the search direction at x, p satisfies Ap = 0. The step-
size procedure will attempt to find an acceptable step length, while retaining
feasibility with respect to all constraints.
It is easy to compute the maximum feasible step length that can be taken
along a direction p using a min ratio rule (see Chapter 2):

α = maxα x + αp is feasible 
aTi x − bi T
= min ai p < 0, i ∈
6 J .
−aTi p

Now, we outlet a simple active-set method. Assume, that a feasible start-

ing point x0 is given and let J still be the index set of the active constraints
at x0 , now. Let A be the corresponding constraint matrix, and let Z be a
null-space matrix for A. Set k = 0.

1. Optimality Test: If Z T ∇f (xk ) = 0, then:

(a) If no constraints are active, then the current point is a local (un-
constrained) stationary point
—STOP.
(b) Else, compute Lagrange multipliers:
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
62 METHODS

T
µ = Ar ∇f (xk ).
(c) If µ ≥ 0, then, STOP: a local stationary point has been reached.
Otherwise, drop a constraint corresponding to a negative multi-
plier from the active set, update J, A, Z and Ar .

2. Search Direction: Compute a descent direction p that is feasible with

respect to the constraints in J .

f (xk + αp) < f (xk )

and

α ≤ α,
where α is the maximum feasible step along p.

4. Update: Find the new point

xk+1 = xk + αp.
If a new constraint boundary is encountered (α = α), then, add it to
the working set and update J, A, Z and Ar accordingly.
(if more than one constraint boundary is encountered, then, pick one
of the constraints to enter the working set; this is a degenerate case.)
Set k = k + 1 and return to 1.

Example 4.7. Consider the problem

 minimize f (x) := 12 (x1 − 3)2 + (x2 − 2)2

subject to 2x1 − x2 ≥ 0
(NLP)A
>

 −x1 − x2 ≥ −4

x2 ≥ 0.

We use an active-set method to solve this problem. The equality constrained

subproblems will be solved using a reduced Newton method. Since the ob-
jective function is quadratic, we use a step α = 1 whenever it is feasible, i.e.,
whenever α ≥ 1. Otherwise, we shall take the step α = α.
4.3. LINEAR INEQUALITY CONSTRAINTS 63

Furthermore, we use the variable reduction method to compute the null-

space matrix for the active constraints, and for simplicity always use the
left-hand submatrix of A as the basis matrix. This gives:
   
−B −1 N B −1
A = (B N ), Z= , Ar = .
I 0

Let x0 = (0, 0)T be our feasible starting point. Since the constraints j = 1
and j = 3 are active, we let J = {1, 3}, so
   1 1

2 −1 2 2
A= , Ar = .
0 1 0 1
The matrix Z is empty, hence the reduced gradient Z T ∇f (x0 ) vanishes triv-
ially. We therefore compute Lagrange multipliers (for simplicity, suppressing
index k)

   1
   
µ1 T 0 −3 − 32
µ= = Ar ∇f (x0 ) = 2
1 = .
µ3 2
0 −4 − 11
2

Both multipliers are negative, thus, we should drop one of the constraints
from J. We drop j = 3, because its multiplier is more negative. Updating
the working set J := {1} gives
 1
  1

A = (2, 1), Z= 2 , Ar = 2 ,
1 0
so that the new reduced gradient is Z T ∇f (x0 ) = − 11
2
. We now compute the
reduced Newton search direction to obtain
−1 T
p = −Z Z T∇2 f (x0 )Z Z ∇f (x0 )
1 
9 −1

= − 2 4
− 11
2
1
 11 
= 9 .
22
9

The maximum step to the boundary of the constraints is α = 12 11

; hence,
a step length of 1 is feasible (α = 1). The full Newton step is taken to reach
the point
   11
  11

0 9 9
x1 = +1 22 = 22 .
0 9 9
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
64 METHODS

The next iteration begins with the optimality test:

 
T 1 − 16
Z ∇f (x1 ) = ( 1) 9
8 = 0.
2 9

Thus, the reduced gradient vanishes at x1 , as expected. Since a local mini-

mum of f with respect to the working set has been found, we compute the
Lagrange multiplier corresponding to the active constraint:
 
T 1 − 16 8
µ = (µ1 ) = Ar ∇f (x1 ) = ( 0) 8
9 =− .
2 9 9
As this multiplier is negative, we drop the constraint from J. Now, we are
left with J = ∅, which means that the problem is locally unconstrained.
The reduced gradient is simply the gradient itself (Z = I), and the search
direction is simply the unconstrained Newton direction:

2 −1
p = − (∇
 f (x1 ))  ∇f (x
1)
16
1 0 −9
= − 1 8
0
 16 2 9

= 9 .
− 94
Since the largest feasible step to the boundary is
 
1 11 1
α = min , = ,
4 2 4
we use α = 41 , and at the new point
 11
  16
  5

9
1 9 3
x2 = 22 + = ,
9 4 − 49 7
3

the constraint j = 2 is active. We now update J = {2},

   
−1 −1
A = (−1 1), Z= , Ar = .
1 0
Again testing for optimality, we find that the reduced gradient at the new
point
 
T − 34
Z ∇f (x2 ) = (−1 1) 2 = 2.
3
4.3. LINEAR INEQUALITY CONSTRAINTS 65

Since it is not zero, we continue with a search in the Newton direction. This
gives

−1
p = −Z Z T ∇2 f (x2 )Z Z T ∇f (x2 )
−1
= − (3)−1 2
 2  1
= 3 .
− 23

The maximum step to the boundary is α = 72 , hence, the step length is

α = 1. This gives
 5
  2
  7

x3 = 3 +1 3 = 3 .
7
3
− 32 5
3

At the new point, the reduced gradient is

 
T − 23
Z ∇f (x3 ) = (−1 1) = 0.
− 23
Since

µ = (µ2 )
T
= Ar ∇f (x
3 ) 2 
−3
= (−1 0)
− 23
2
= 3 > 0,
the point x3 satisfies the first-order optimality conditions from (NOC)> , and
we terminate. Since the objective function f is strictly convex, the solution
x∗ = ( 73 , 53 )T is a strict global minimizer. The Lagrange multipliers corre-
sponding to the three constraints are µ∗1 = 0, µ∗2 = 32 , µ∗3 = 0.
The progress of the algorithm is shown in Figure 4.2.


One possible modification is to solve the equality constrained subprob-

lems inexactly, whenever there is reason to believe that the working set is
not the optimal active set. The rationale is that faster progress may be made
by obtaining a better working set than by getting a few extra digits of accu-
racy on an incorrect set. This idea requires care in implementation, however.
The reason is that as the solution to a problem becomes less accurate, the
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
66 METHODS

computed Lagrange multipliers also become less accurate. These inaccura-

cies can effect the sign of a computed Lagrange multiplier. Consequently, a
constraint may erroneously be deleted from the working set, thereby out any
potential savings.
Another possible danger is zigzagging. This phenomenon can occur if
the iterated cycle repeatedly between two working sets. This situation is
depicted in Figure 4.3. Zigzagging can not occur if the equality-constrained
problems are solved sufficiently accurately before constraints are dropped
from the working set.

Figure 4.3: Zigzagging

4.4. SEQUENTIAL QUADRATIC PROGRAMMING 67

Finally, we conclude by indicating how the active-set method can be

adopted to solve a problem of the form

 minimize f (x)
A1,2
(NLP) subject to A1 x = b1 ,
{>=} 
A2 x ≥ b 2 ,
where equality constraints arise in addition to the inequality constraints. In
this case, the equality constraints are kept permanently in the working set J
since they must be kept satisfied at every iteration. The Lagrange multipliers
for equality constraints can be positive or negative, so they do not play a role
in the optimality test. The equality constraints also do not play a role in
the selection of the maximum allowable step length α. These are the only
changes that need be made to the active-set method!
At the present stage of our considerations, we could consider linear pro-
gramming (cf. Chapter 2) to a great extent as special case of the methods
considered in the present Chapter 4. We recall that (LP) is just (NLP){=}
>
with f (x) := cT x, A1 := A, b1 := b, A2 := I and b2 := 0. This is (by some
hints) left as an exercise.

4.4 Sequential Quadratic Programming

Sequential quadratic programming is a popular and successful technique for
solving nonlinear constrained problems. The main idea is to obtain a search
direction by solving a quadratic program, i.e., a problem with a quadratic
objective function and linear constraints. This approach is a generalization
of Newton’s method for unconstrained minimization.
Methods for solving our problem


minimize f (x)
(NLP)=
subject to h(x) = 0
can be derived by applying Newton’s method to the corresponding optimality
conditions. Here, f and h := (h1 , h2 , . . . , hm )T are assumed to be of class C 2 .
The Lagrangian for (NLP)= , arising in these Karush-Kuhn-Tucker conditions
(here, Lagrange multiplier rule) is

L(x, λ) := f (x) − λT h(x),

CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
68 METHODS

and this first-order optimality condition is

(NOC)= ∇L(x, λ) = 0.

Then, the formula for Newton’s method is

     
xk+1 xk pk
= + ,
λk+1 λk νk
where pk and νk are obtained as the solution to the linear system
 
2 pk
∇ L(xk , λk ) = −∇L(xk , λk ).
νk
This system has the form (exercise)

 
∇2xx L(xk , λk ) −∇h(xk )    
(⊕)  − (∇h(xk ))T 0  pk = −∇x L(xk , λk ) ,
| {z } νk h(xk )
=:∇T h(xk )

and it represents the first-order optimality conditions for the (auxiliary) op-
timization problem

 1 T
 
k minimize 2
p∇2xx L(x
 k , λ k ) p + p T
∇ x L(x k , λ k )
(QP)
subject to ∇T h(xk ) p + h(xk ) = 0,
with νk being the vector of Lagrange multipliers. This optimization problem
is a quadratic program; i.e., it is the minimization of a quadratic function
subject to linear constraints. The quadratic function is a Taylor series ap-
proximation to the Lagrangian at (xk , λk ), and the constraints are a linear
approximation to h(xk + p) = 0.
In a sequential quadratic programming (SQP) method, at each iteration
a quadratic problem is solved to obtain (pk , νk ). These are used to update
(xk , λk ), and the process repeats at the new point. Each of the quadratic
programs is solved using the technique described in Chapter 2.
Example 4.8. We apply the SQP method to the problem


minimize f (x) := e3x1 +4x2
(NLP)=
subject to h(x) := x21 + x22 − 1 = 0.
4.4. SEQUENTIAL QUADRATIC PROGRAMMING 69

The solution to this problem is x∗ = (− 53 , − 54 )T with λ∗ = − 25 e−5 ≈ −.016845.

We use the initial guess x0 = (−7, −7)T and λ0 = −.01. At this point
   
3x1 +4x2 3 .02234
∇f = e = ,
4 .02979
  
2 9 12
3x1 +4x2 .06702 .08936
∇f = e = ;
12 16 .08936 .11915
= x21 + x22 − 1 = −.02,
     
2x1 1.4 2 2 0
∇h = = , ∇ h= ;
2x2 1.4 0 2
 
.008340
∇x L = ∇f − λ∇h = ,
.015786
 
.08702 .08936
∇2xx L 2
= ∇ f − λ∇ h = 2
.
.08936 .13915

The corresponding quadratic program is

 1 T
 
0 minimize 2
p ∇2xx L p + pT ∇x L
(QP)
subject to ∇T h p + h = 0,

Its solution can be found using

    
∇2xx L −∇h p −∇x L
=
−∇T h 0 ν h
or
    
.08702 .08936 1.4 p1 −.008340
 .08936 .13915 1.4   p2  =  −.015786  .
1.4 1.4 0 ν −.020000

The solution of the quadratic program is

 
.14196
p0 = , ν0 = −.004808,
−.15624
and the new estimates of the solutions are
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
70 METHODS

 
−.55804
x1 = x 0 + p 0 = ,
−.85624

λ1 = λ0 + ν0 = −.014808.
The complete iteration is given in Table 4.1. Notice the rapid convergence
rate, as expected for Newton’s method.

k xk λk || ∇x L || || h ||
0 -.70000 -.70000 -.010000 2 ×10−2 2 ×10−2
1 -.55804 -.85624 -.014808 2 ×10−3 5 ×10−2
2 -.60779 -.79780 -.016469 3 ×10−4 6 ×10−3
3 -.59988 -.80013 -.016839 8 ×10−6 7 ×10−5
4 -.60000 -.80000 -.016845 2 ×10−9 3 ×10−8
5 -.60000 -.80000 -.016845 3 ×10−16 4 ×10−15

Table 4.1: SQP method


Computing the updates pk and νk by solving the quadratic program cor-
responding to applying Newton’s method to the optimality conditions for the
original problem. As a result, this method will have a quadratic convergence
rate provided that ∇2 L(x∗ , λ∗ ) is nonsingular. This rapid convergence rate is
observed in the example. The Hessian of the Lagrangian will be nonsingular
if the regularity condition and the second-order sufficiency conditions for the
original optimization problem are satisfied, i.e., if ∇h(x∗ ) is of full rank and
if Z T ∇2xx L(x∗ , λ∗ )Z is positive definite, where Z is a basis matrix for the
null-space of ∇h(x∗ )T .
For the following two reasons, the SQP method outlined above is not
often used in this simple form:

(i) Convergence to a local solution is not guaranteed.

(ii) It is too expensive. The method outlined above requires the Hessians of
the objective function f and the constraint functions g. Once these have
been obtained, a quadratic program must be solved to determine the
updates pk and νk . It would be desirable to reduce these requirements
for derivatives, and to reduce the number of arithmetic operations re-
quired at each iteration of the algorithm. If larger problems are being
solved, it would also be desirable to reduce the storage requirements
by using computational techniques that do not require matrix storage.
4.4. SEQUENTIAL QUADRATIC PROGRAMMING 71

One way to reduce this expense is to use a quasi-Newton approximation to

the Hessian of the Lagrangian. If this is done, then, second derivatives need
not be computed. in addition, the quasi-Newton matrix can be maintained
in the form of a factorization, and in this way the arithmetic costs of the
method can also be reduced.
Choosing an update formula for the quasi-Newton approximation can be
a more complicated decision than in the unconstrained case.
Let us come back to (ii): Convergence results for the SQP programming
method are obtained by insisting that (xk+1 , λk+1 ) be a better estimate of
the solution than (xk , λk ). In the unconstrained case, progress is measured in
terms of the objective function. In our constrained case, progress is measured
in terms of an auxiliary merit function. Usually, a merit function is in the sum
of terms that include the objective function and the amount of infeasibility of
the constraints. If the new point reduces the objective function and reduces
infeasibility, then the value of the merit function will decrease. In many
instances, however, improvements in the objective value come at the expense
of feasibility, and vice versa, so the merit function must balance these two
aims. One example of a merit function is the quadratic penalty function

=:hT (x)
z }| {  m
X 
T 2
M(x) := f (x) + ρ (h(x)) h(x) = f (x) + ρ hi (x) ,
i=1

where ρ is some positive number (ρ > 0). The grater the value of ρ, the
grater the penalty for infeasibility. This merit function is a function of x
only; other examples of merit functions may be functions of both x and λ.
Example 4.9. Let us use the merit function

M(x) := f (x) + 10hT (x)h(x)

to measure the progress of the SQP method in the previous example. At the
first iteration

x0 = (−.7000, −.7000)T ,
f (x0 ) = .007447,
h(x0 ) = −.02,
M(x0 ) = .011447,

and
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
72 METHODS

x1 = (−.55804, −.85624)T ,
f (x1 ) = .006102,
h(x1 ) = .044563,
M(x1 ) = .025961,
so, in terms of this merit function, the point x1 is worse than the point x0 .
At the next iteration, though,

x2 = (−.60779, −.79780)T ,
f (x1 ) = .006641,
h(x1 ) = .005890,
M(x1 ) = .006988,
indicating that x2 is better than both x1 and x0 for this merit function. 

Ideally, the merit function would be chosen so that (x∗ , λ∗ ) would be a

local minimizer of the merit function if and only if it were a local solution of
the optimization problem. If this were true then a line search with respect
to the merit function could be performed:
     
xk+1 xk pk
= +α ,
λk+1 λk νk
where α is chosen so that

M(xk+1 , λk+1 ) < M(xk , λk ).

For this to be successful, the search direction from the quadratic program
would have to be descent direction for the merit function. Unfortunately, it
is rarely possible to guarantee that the local minimizers of the merit function
and the local solutions of the optimization problem coincide. For the merit
function M(x) := f (x) + ρ hT (x)h(x) (quadratic penalty) some of the local
minimizers of M approach local solutions of the constrained problem in the
limit as ρ → ∞. In addition, M may have local minima at points where
h(x) 6= 0, i.e., at infeasible points. Other merit functions have analogous
deficiencies that can limit the applicability of convergence theorems, or can
complicate the development of SQP methods.
In the unconstrained case we assume that search directions are descent
directions with respect to the objective function. Here, we assume that
either pk or the combined vector (pk , νk ) is a descent direction with respect
to the merit function. A common way to guarantee this is to insist that the
reduced Hessian for the quadratic program be positive definite. Come back
4.4. SEQUENTIAL QUADRATIC PROGRAMMING 73

to the constraints. If a quasi-Newton approximation to the Hessian is used

to define the quadratic program, then, positive definiteness of the reduced
Hessian is often guaranteed by the choice of quasi-Newton update formula.
If Newton’s method is used so that the Hessian in the quadratic program is
∇2xx L(xk , λk ), then, it is necessary to test if the reduced Hessian is positive
definite, and to modify if it is not. This testing is more complicated than in
the unconstrained case, particularly if the quadratic program is solved via
the linear system (⊕). In this case, the reduced Hessian may not available.
and more elaborate tests for positive definiteness must be used.
Near the solution of the constrained problem, we would normally like to
take a step of α = 1 in the line search, so that the quadratic convergence rate
of Newton’s method could be achieved. Hence, the merit function should be
chosen so that a step of α = 1 is guaranteed to be accepted in the limit as the
solution is approached. For certain merit functions this is not true. In such
cases it is possible to give examples where a step of α = 1 is unacceptable at
every iteration, no matter how close the current point is to the solution.
We will ignore many of these difficulties and only state that pk is a descent
direction for the quadratic penalty merit function. Even more, we will make
the simplifying assumption that the full Hessian of the Lagrangian (or some
approximation to it) is positive definite. If a quasi-Newton approximation
to the Hessian is used, this assumption is not reasonable, but for Newton’s
method it is restrictive. The lemma stated subsequently can be used to show
that

lim ∇M(xk ) = 0.
k→∞

For large values of ρ, local solutions of the constrained problem are ap-
proximative local minimizers of M(xk ) to within O(1/ρ), so this argument
provides a rough outline of a convergence theorem.
Lemma 4.1. Assume that (pk , νk ) is computed as the solution to the quadratic
problem

 1 T e

minimize p Hp+ pT ∇x L(xk , λk )
(QP)kapprox 2
T
subject to ∇ h(xk ) p + h(xk ) = 0,
e is some positive-definite approximation to ∇2xx L(xk , λk ). If pk 6= 0,
where H
then

pTk ∇M(xk ) < 0

CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
74 METHODS

M(x) := f (x) + ρhT (x)h(x);

i.e., pk is a descent direction with respect to this merit function.

Proof. See Nash, Sofer (1996). 

Finally, let us consider the following problem with inequality constraints

 minimize f (x)
(NLP) subject to h(x) = 0,

g(x) ≥ 0.

Here, we can develop an SQP method for this problem as well. Our earlier
quadratic program (QP)k was based on a quadratic approximation to the
Lagrangian function and a linear approximation to the constraints. If the
same approximation is used here, we obtain:

 1 T
 minimize 2
p (∇2xx L(xk , λk )) p + pT ∇x L(xk , λk )
(QP) k
subject to ∇T h(xk )p + h(xk ) = 0,

∇T g(xk )p + g(xk ) ≥ 0.

The solution (pk , νk )to this quadratic program provide the step to the next
estimate of the solution of the original constrained problem. This problem
can be solved using the active-set method discussed in Section 4.3.

Reduced-gradient methods try to maintain feasibility at every iteration. This
approach has several advantages:

(i) If each estimate of the solution is feasible, the algorithm can be stopped
before it converges and the approximate solution may still be useful.

(ii) Guaranteeing convergence is simpler because progress can be measured

directly using the value of the objective function, rather than with an
auxiliary merit function.
4.5. REDUCED-GRADIENT METHODS 75

The disadvantage of reduced-gradient methods is the computational ex-

pense of ensuring that nonlinear constraints remain satisfied at every itera-
tion.
We apply a version of the reduced-gradient method to our problem


minimize f (x)
(NLP)=
subject to h(x) = 0.

Our derivation begins in the same way as for SQP. The Lagrangian is

L(x, λ) := f (x) − λT h(x),

and, if Newton’s method is applied to the first-order optimality condition,
     
xk+1 xk pk
= + ,
λk+1 λk νk
where pk and νk are obtained as the solution to the linear system
 
2 pk
∇ L(xk , λk ) = −∇L(xk , λk ).
νk
This linear system has the form (cf. Section 4.4):

    
∇2xx L(xk , λk ) −∇h(xk ) pk −∇x L(xk , λk )
(⊕) = .
−∇T h(xk ) 0 νk h(xk )

At this point, the derivations for the two methods diverge. In the reduced-
gradient method, we use these formulas to derive a portion of the search
direction, the portion that lies in the null-space of the constraint gradients.
If Zk is a basis matrix for the null-space of ∇T h(xk ), and Yk is a basis matrix
for the range-space of ∇h(xk ), then

pk = Z k pZ + Yk pY ,
where pZ is the solution to the reduced problem


ZkT ∇2xx L(xk , λk ) Zk pZ = −ZkT ∇x L(xk , λk ).
This formula determines pZ in the reduced-gradient method. If all constraints
are linear, then the formula is equivalent to the formula for the reduced
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
76 METHODS

Newton method derived in Section 4.1. In that case, the matrix Zk will be
the same at every iteration, and

ZkT ∇x L(xk , λk ) = ZkT ∇f (xk ).

Hence, the reduced-gradient method is a generalization of the reduced New-
ton method for linearly constrained problems. Come back to the general
nonlinear case.
The remaining portion of the search direction is determined from the
condition that the new estimate of the solution must be feasible: h(xk+1 ) = 0.
Since xk+1 = xk + pk , this condition has the form

h(xk + Zk pZ + Yk pY ) = 0.
This is a system of m nonlinear equations in the m variables pY . (We assume
that ∇h(xk ) is a matrix of null rank.) If the constraints are linear, then
pY = 0. If they are nonlinear, some auxiliary algorithm must be applied to
this nonlinear system to determine pY . For example, Newton’s method could
be used.
Example 4.10. Let again be given (cf. Example 4.8):


minimize f (x) := e3x1 +4x2
(NLP)=
subject to h(x) := x21 + x22 − 1 = 0.

The solution being x∗ = (− 35 , − 54 )T with λ∗ = − 25 e−5 ≈ −.016845. Again, let

x0 = (−.7, −.7)T , even though this point is infeasible. There, ∇f, ∇2 f, h, ∇h
and ∇2 h were computed in Example 4.8. An estimate of the Lagrange mul-
tiplier is needed to determine the gradients and Hessians of the Lagrangian.
We compute a multiplier estimate λk by solving

λ

At x0 the multiplier estimate is λ0 = −.018616. Using this value, we obtain

 
.008340
∇x L = ∇f − λ∇h = ,
.015786
4.5. REDUCED-GRADIENT METHODS 77

 
.08702 .08936
∇2xx L 2
= ∇ f − λ∇ h = 2
.
.08936 .13915
We use variable-reduction to compute the null-space matrix Zk and the range-
space matrix Yk based on

∇T h = (2x1 , 2x2 ) = (B, N ).

Then,
   
−B −1 N − xx21
Zk = = ,
I 1
 
1 x1
Yk = p 2 .
x1 + x22 x2
At this iteration
 
−1
Z0 = .
1
The null-space portion of the search direction is obtained by solving


ZkT ∇2xx L Zk pk = −ZkT ∇x L,
or

(.081912)pk = −.0074466,
so that pZ = −.090909.
The remaining portion of the search direction, Yk pY , is determined by
solving

h(xk + Zk pZ + Yk pY ) = 0.
using Newton’s method. In this example,
 
y1
Yk p Y = γ
y2
for some unknown γ, where y1 and y2 are the components of Yk . If we define
x
bk := xk + Zk pZ , then the condition for Yk pY has the form
CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
78 METHODS

x1 + γy1 )2 + (b
ϕ(γ) := (b x2 + γy2 )2 − 1 = 0.
Applying Newton’s method to this equation gives the iteration

ϕ(γi )
γi+1 = γi −
ϕ0 (γi )
x1 + γi y1 )2 + (b
(b x 2 + γ i y 2 )2 − 1
= γi − .
2y1 (b
x1 + γi y1 ) + 2y2 (b
x2 + γ i y 2 )
In this example, we initialize the iteration with γ0 = 0. Then,

γ0 = 0, ϕ(γ0 ) = −3 × 10−3 ,
γ1 = .001753157308727014, ϕ(γ1 ) = 3 × 10−6 ,
γ2 = .001751607670563711, ϕ(γ2 ) = 2 × 10−12 ,
γ3 = .001751607669352957, ϕ(γ0 ) = −2 × 10−16 .
The overall search direction is

p0 = Z
0 pZ +Y0 pY  
−1 −.70711
= (−.090909) + (.001752)
 1  −.70711
.089671
=
−.092148
and
 
−.61033
x1 = x 0 + p 0 = .
−.79215
The complete iteration is given in Table 4.2. The initial guess of the solution
is not feasible, but all later estimate of the solution are feasible (to 16 digits).
The method converges rapidly, as expected for Newton’s method. 

k xk λk || ∇x L || || h ||
0 -.70000 -.70000 -.018616 5 ×10−3 2 ×10−2
1 -.61033 -.79215 -.016851 4 ×10−4 0
2 -.60001 -.79999 -.016845 3 ×10−7 0
3 -.60000 -.80000 -.016845 1 ×10−16 1 ×10−16

The reduced-gradient method corresponds to using Newton’s method in

the null-space of the constraints, so it can be expected to converge quadrat-
ically in nondegenerate cases. (This quadratic convergence rate is observed
4.5. REDUCED-GRADIENT METHODS 79

in the example above.) As before, Newton’s method is not guaranteed to

converge, and even if it does converge, it may converge to a maximum or
stationary point and not a minimum. Hence, some globalization strategy
must be employed to ensure convergence to a local solution of (NLP)= .
A line search can be used to guarantee convergence, as was done with the
SQP method. If all the solution estimates xk are feasible points, then the
value of the quadratic penalty merit function is

M(xk ) = f (xk ) + ρhT (xk )h(xk ) = f (xk ),

so a line search can be performed using the objective function f itself as a
merit function.
The line search for a reduced-gradient method is more complicated than
in the unconstrained case. for each trial value of α, the trial point must
satisfy

h(xk + αZk pZ + Yk pY ) = 0.
Hence, pY depends on α, and must be computed by solving a nonlinear
system of equations. For large value of α there may not be pY that satisfies
the constraints, further complicating the line search algorithm. For these
reasons, it is not entirely correct to say that the reduced-gradient method
produces a search direction, since in fact the method must search along an
arc defined by pZ and α.
Concerning descent, we state:

Lemma 4.2. Assume that the reduced-gradient method is applied to the prob-
lem (NLP)= . Let xk be the kth estimate of the solution, with h(xk ) = 0. Let
∇h(x) be of full rank for all x in a neigbourhood of Nk of xk . Also assume
that the null-space portion pZ of the search direction is computed from


ZkT ∇2xx L(xk , λk ) Zk pk = −ZkT ∇x L(xk , λk ),
where λk is an estimate of the Lagrange multipliers, and where


ZkT ∇2xx L(xk , λk ) Zk
is positive definite. Define pY () as the solution to

h(xk + Zk pZ + Yk pY ()) = 0.

CHAPTER 4. NONLINEAR PROGRAMMING: FEASIBLE-POINT
80 METHODS

f (xk + Zk pZ + Yk pY ()) < f (xk ),

for all sufficiently small positive values of .

Proof. See Nash, Sofer (1996). 

Variants of the reduced-gradient method have been developed that are
more flexible in that they allow some violation of the constraints. such meth-
ods might be considered a compromise between reduced-gradient and SQP
methods with their different advantages and disadvantages.
Chapter 5

Nonlinear Programming:
Penalty and Barrier Methods

For the feasible-point methods, considered in the previous chapter, there are

(i) As the number of constraints increases, the number of potential sub-

problems increases exponentially. While the hope is that the algorithm
will consider only a small proportion of these subproblems, there is no
known method to guarantee that this indeed will be the case.

(ii) The idea of keeping the constraints satisfied exactly, although easily
achieved in the case of linear constraints, is much more difficult to
accomplish in the case of nonlinear constraints, and in some cases may
not be desirable.

In this chapter, we discuss a group of methods, referred to as penalizing

methods, that remove some of these difficulties. Namely, these methods solve
a constrained optimization problem by solving a sequence of unconstrained
optimization problems. The hope is that in the limit, the latter solutions
will converge to the solution of the constrained problem. The unconstrained
problems involve an auxiliary function that incorporates the objective func-
tion or the Lagrangian function, together with penalty terms that measure
violations of the constraints. The auxiliary function also includes one or more
parameters that determine the relative importance of the constraints in the
auxiliary function. By changing these parameters appropriately, a sequence
of problems is generated where the effect of the constraints becomes increas-
ingly pronounced. In contrast to active-set methods, the auxiliary function
CHAPTER 5. NONLINEAR PROGRAMMING: PENALTY AND
82 BARRIER METHODS

takes into account all constraints, even when inequalities are present; thus the
combinatorial difficulties of guessing a correct active set are avoided. Further,
since they do not attempt to keep the constraints satisfied exactly, penaliza-
tion techniques can be more suitable for handling nonlinear constraints.

5.1 Classical Penalty and Barrier Methods

Here, we have two groups of methods:

(α) one group imposes a penalty for violating a constraint,

(β) the other group imposes a penalty for reaching the boundary of an
inequality constraint.

We call (α) penalty methods, and (β) barrier methods.

We start with a geometrical motivation. Let us consider


minimize f (x)
(P)con
subject to x ∈ M,

where M is the set of feasible points. Define


0, if x ∈ M,
σ(x) :=
+∞, if x 6∈ M,
being such to say, an infinite penalty for violating feasibility. Hence, the
problem (P)con with its constraint can be transformed into an equivalent
unconstrained problem

Conceptually, if we could solve (P)uncon , we would be done: a point x∗ solves

(P)con if and only if it solves (P)uncon . Unfortunately, this is not a practical
idea, since the objective function of the unconstrained minimization is not
defined outside of the feasible region. Even if we were to replace the “∞”
by a large number, the resulting unconstrained problem would be difficult to
solve because of its discontinuities.
5.1. CLASSICAL PENALTY AND BARRIER METHODS 83

Instead, barrier and penalty methods solve a sequence of unconstrained

subproblems that are more manageable, and that gradually approximate the
problem (P)uncon . This is achieved by replacing the ideal penalty σ by a
continuous function which gradually approaches σ.
In barrier methods, this function (called barrier term) approaches σ from
the interior of the feasible region. It creates a barrier that prevents the iter-
ates from reaching the boundary of the feasible region. in penalty methods,
this function (called a penalty term) approach σ from the exterior of the
feasible region. It serves as a penalty for being infeasible.
Barrier methods generate a sequence of strictly feasible iterates which con-
verge to a solution of the problem from the interior of the feasible region. For
this reason, they are called interior-point methods; in Chapter 3 we studied
them in the special case of linear programming! Since these methods require
the interior of the feasible region to be 6= ∅, they are not more appropriate
for problems without equality constraints (see, however, our discussion in the
linear case, including infeasibility, in Chapter 3).
In contrast, penalty methods permit the iterates to be infeasible. A gradu-
ally increasing penalty is imposed for violation of feasibility, however. Penalty
methods usually generate a sequence of points that that converge to a solu-
tion of the problem from the exterior of the feasible region. These methods
are usually more convenient on problems with equality constraints.
Despite their apparent differences, barrier and penalty methods have
much in common. Their convergence theories are similar, and the underly-
ing structure of their unconstrained problems is similar. Much of the theory
for barrier methods can be replicated for penalty methods and vice versa.
It is common to use the generic name “penalty methods” to describe both
methods, with interior penalty methods referring to (β) (barrier methods),
and exterior penalty methods referring to (α) (penalty methods).

Barrier Methods:
Consider the problem


minimize f (x)
(NLP)≥
subject to g(x) ≥ 0,

where g := (g1 , g2 , . . . , gs )T ; f and g being of class C 2 . Barrier methods are

strictly feasible methods, i.e., the iterates lie in the interior of the feasible re-
gion. We assume, therefore, that the feasible set M has a nonempty interior:
CHAPTER 5. NONLINEAR PROGRAMMING: PENALTY AND
84 BARRIER METHODS

M 0 6= ∅; i.e., ∃x0 ∈ Rn : gj (x0 ) > 0 ∀j ∈ {1, . . . , s}. We also assume that it

is possible to reach any boundary point by approaching it from the interior.
Barrier methods maintain feasibility by creating a barrier that keeps the
iterates away from the boundary of the feasible region. These methods use
a barrier term that approaches the infinite penalty function σ. Let ϕ(x) be
a function which is continuous on M 0 , and that becomes unbounded as the
boundary of the set is approached from its interior:

ϕ(x) → ∞ as gj (x) → 0+ for some j ∈ {1, 2, . . . , s} .

Two examples for such a function ϕ:

s
X
(a) ϕ(x) = − ln(gj (x)),
j=1

s
X 1
(b) ϕ(x) = − .
j=1
gj (x)

Figure 5.1: Effect of barrier term.

Now, let some τ > 0 be given. Then, τ ϕ(x) will approach σ(x) as τ → 0+ .
This is demonstrated in Figure 5.1 for a one-dimensional problem with bound
constraints. By adding a barrier term of the form τ ϕ(x) to the objective, we
obtain a barrier function

β(x, τ ) := f (x) + τ ϕ(x),

5.1. CLASSICAL PENALTY AND BARRIER METHODS 85

where τ is referred to as the barrier parameter. The best known barrier

function is the logarithmic barrier function

s
X
β(x, τ ) := f (x) + τ ln(gj (x)),
j=1

but the inverse barrier function

s
X 1
β(x, τ ) := f (x) + τ ,
j=1
gj (x)

lems:

(NLP)kuncon minimize β(x, τk )

x

for a sequence τk k∈N0 of barrier parameters τk > 0 that decrease monoton-
ically to zero: τk+1 ≤ τk ∀k ∈ N0 , τk → 0+ (k → ∞).
As the barrier term is ∞ at the boundary ∂M of the feasible set M , it
acts as a repelling force which derives the minimizers of the barrier function
away from the boundary into the interior M 0 of the feasible region. Thus,
any minimizer of the barrier function will be strictly feasible. As the barrier
parameter is decreased, however, the effect of the barrier term is diminished,
so that the iterates can gradually approach the boundary ∂M .
Why do we solve a sequence of problems? It might seem better to solve a
single unconstrained problem using a small µ –but this is usually not practi-
cal: When τ is small, the problems are difficult to solve (τ ϕ is close in shape
to the infinite penalty σ). If τ is decreased gradually, and if the solution of
one unconstrained problem is used as the starting point of the next problem,
these unconstrained problems tend to be much easier to solve.
Example 5.1. Consider the nonlinear program

 minimize f (x) := x1 − 2x2
(NLP)≥ subject to 1 + x1 − x22 ≥ 0,

x22 ≥ 0.

Then, the logarithmic barrier function gives the unconstrained problems

CHAPTER 5. NONLINEAR PROGRAMMING: PENALTY AND
86 BARRIER METHODS

x

for a sequence of decreasing µk . For a specific µk , the first-order necessary

optimality conditions are:

 τk

 1− = 0,
1 + x1 − x22
(NOC)k 2τk x2 τk
 −2 −
 − = 0.
2
1 + x 1 − x 2 x2

If the constraints are strictly satisfied, then, the denominators are positive.
We obtain an equation for x2 :

1
x22 − x2 − τ = 0,
2
giving

1 + 2τ
1+
x2 (τ ) := (τ := τk )
2
(why is the positive root taken?). Since x1 = x22 − 1 + τ , we get

1 + 2τ + 3τ
x1 (τ ) := .
2
The unconstrained objective is strictly convex, hence, this solution is the
unique local minimizer in the feasible region. As µ approaches 0 (i.e., when
k → ∞), we obtain
 √ 
  1+2(0)+3(0)−1  
x1 (τ ) 0
lim = √2 = =: x∗ .
τ →0+ x2 (τ ) 1+ 1+2(0) 1
2

It is easy to verify that this point x∗ is indeed the solution of (NLP)≥ .

Table 5.1 shows the values of x1 (τ ) and x2 (τ ) for a sequence of barrier
parameters τ = τk . The initial parameter is related to be τ0 := 1, and
consecutive parameters are decreased by a factor of 10. Observe that x(τ ) =
(x1 (τ ), x2 (τ ))T exhibits a linear rate of convergence to the optimal solution.

5.1. CLASSICAL PENALTY AND BARRIER METHODS 87

τ x1 (τ ) x2 (τ )
100 1.8660254 1.3660254
10−1 0.1977226 1.0477226
10−2 0.0199752 1.0049752
10−3 0.0019998 1.0004998
10−4 0.0002000 1.0000500
10−5 0.0000200 1.0000050
10−6 0.0000020 1.0000005

The previous examples illustrates a number of features that typically

occur in a barrier method:

(i) x(τk ) k∈N0 converges to a optimal solution x∗ . (It is possible to prove
convergence for barrier methods under mild conditions.)
(ii) τ 7→ x(τ ) (τ > 0) is a differentiable function, called the barrier tra-
jectory. If the logarithmic or inverse barrier method are used, it ex-
ists, provided x∗ is a regular point of the constraints that satisfies the
second-order sufficiency conditions. As well as the strict complemen-
tarity conditions. The existence of a trajectory can be used to develop
algorithms such as path-following algorithms for LP (similar to Chap-
ter 3). It can also be used to develop techniques that accelerate the
convergence of a barrier method.
(iii) Under a regularity condition again, for both the logarithmic and the
inverse barrier functions, we find estimates µ(τ ) of the Lagrange mul-
tiplier µ∗ at the optimal point x∗ . We explain this subsequently.

To (iii): We discuss it for the case of the logarithmic barrier function. Con-
sider a point x = x(τ ) which is a minimizer of the logarithmic barrier function
for a specific parameter τ . Setting the gradient of the barrier function (with
respect to x) to zero, we obtain

s
X 1
∇f (x) − τ ∇gj (x) = 0
j=1
gj (x)
if and only if
s
X τ
∇f (x) − ∇gj (x) = 0
gj (x)
j=1 | {z }
=:µj (τ )
CHAPTER 5. NONLINEAR PROGRAMMING: PENALTY AND
88 BARRIER METHODS

(at x = x(τ )). We therefore have a feasible point x(τ ) and a vector µ(τ )
satisfying the first-order necessary optimality conditions, except τ (instead
of 0) on the right-hand side of (cf. (NOC)τsf in Section 3.2.):

 Ps

 ∇f (x(τ )) − j=1 µj (τ )∇gj (x(τ )) = 0

 

 µj (τ )gj (x(τ )) = τ 
(NOC)τsf ≥ ∀j ∈ {1, 2, . . . , s}

 µj (τ ) 0 

 (>)

and (strict) feasibility.

As τ → 0+ , µ(τ ) can be viewed as an estimate of µ∗ , the Lagrange

multiplier at x∗ . (Indeed, if x∗ is a regular point, then: x(τ ) → x∗ (τ →
0) =⇒ µ(τ ) → µ∗ (τ → 0).)
Example 5.2. Consider the nonlinear problem

 minimize f (x) := x21 + x22
(NLP)≥ subject to x1 − 1 ≥ 0,

x2 + 1 ≥ 0,

the solution being x∗ = (1, 0)T . The inequality j = 1 is active at x∗ , and the
corresponding Lagrange multiplier is µ∗1 = 2. Suppose the problem is solved
via a logarithmic barrier method. Then, this method solves

x

for a decreasing sequence of barrier parameters τ = τk which converge to 0.

Now, (NOC)τsf gives

τ
2x1 − = 0,
x1 − 1
τ
2x2 − = 0,
x2 + 1
yielding
√ √
1+ 1 + 2τ −1 − 1 + 2τ
x1 (τ ) = , x2 (τ ) = .
2 2
5.1. CLASSICAL PENALTY AND BARRIER METHODS 89

The corresponding Lagrange multiplier estimates are:

τ √
µ1 (τ ) = = 1 + 2τ + 1,
x1 (τ ) − 1
τ √
µ2 (τ ) = = 1 + 2τ − 1,
x2 (τ ) + 1
When τ → 0, we obtain
       
x1 (τ ) 1 µ1 (τ ) 2
lim = and lim = .
τ →0 x2 (τ ) 0 τ →0 µ2 (τ ) 0
| {z } | {z } | {z } | {z }
=x(τ ) =:x∗ =µ(τ ) =:µ∗


Despite their attractive features, barrier methods also have potential dif-
ficulties; we mention the most important one: The unconstrained problems
become increasingly difficult to solve as the barrier parameter decreases. The
reason is that (with the exception of some special cases) the condition number
of the Hessian matrix of the barrier function at its minimum point becomes
increasingly large, tending to ∞ as τ → 0.
Example 5.3. We come back to our Example 5.2. There,
 τ 
2+ 0
 (x1 − 1)2 
∇2xx β(x, τ ) =  τ .
0 2+
(x2 + 1)2
Suppose now, that x(τ ) is a minimizer of β(·, τ ) for some τ > 0. Recall from
Example 5.2 that

τ τ
µ1 (τ ) = and µ2 (τ ) = .
x1 (τ ) − 1 x2 (τ ) + 1
When τ ≈ 0, then, µ1 (τ ) ≈ 0. Therefore,

 
µ21 (τ ) 4
!
 2+ τ 0  2+ 0
∇2xx β(x(τ ), τ ) =  ≈ .
µ22 (τ )  0
τ
2
0 2+
τ
The condition number of this matrix (i.e., the product of both its matrix
norm and the one of its inverse; exercise) is approximately equal to
CHAPTER 5. NONLINEAR PROGRAMMING: PENALTY AND
90 BARRIER METHODS

4
2+ τ 2 1
= 1+ = O( ),
2 τ τ
hence, the matrix is ill conditioned (highly sensitive against smallest per-
turbations). Although the calculations were performed for a point on the
barrier trajectory, the same will hold at all points in a neighbourhood of the
solution.
The contours of the barrier function in this example are shown in Fig-
ure 5.2 for τ = 1 and τ = 0.1. We see that for the smaller τ , the contours
(level sets, i.e., the “finger print” of β(·, τ )) are almost parallel to the line
x1 = 1. More precisely, the contours are almost parallel to the null-space of
the gradient of the active constraint at the solution. This is characteristic of
barrier functions.


Figure 5.2: Contours of the logarithmic barrier function.

The ill conditioning of the Hessian matrix of β(·, τ ) has several ramnifi-
cations. It rules out the use of an unconstrained method whose convergence
rate depends on the condition number of the Hessian matrix at the solution.
Therefore, Newton-type methods are usually the method of choice. The so-
lution to the Newton equations is also sensitive to the ill condition of the
Hessian matrix, however. The numerical errors can result in a poor search
direction.
In the early 1970s, the ill conditioning of the barrier functions led to
their abandonment. Interest in barrier methods was renewed in 1984, with
the announcement of Karmarkar’s method for LP and the discovery shortly
thereafter that this method is just a special case of a barrier method (cf. our
discussions above, and Chapter 3)!
5.1. CLASSICAL PENALTY AND BARRIER METHODS 91

Recently, special attention has been given to the development of special-

ized linear algebra techniques that compute a numerically stable, approxi-
mate solution to the Newton equations for a barrier function. We discussed
them in Section 3.2 (path-following methods).

Penalty Methods:
In contrast to barrier methods, penalty methods solve a sequence of un-
constrained optimization problem whose solution is usually infeasible to the
original constrained problem. A penalty for violation of the constraints is
incurred, however. As this penalty is increased, the iterates are forced to-
wards the feasible region. An advantage of penalty methods is that they do
not require the iterates to be strictly feasible. Thus, unlike barrier methods,
they are suitable for problems with equality constraints.
Consider first the equality-constrained problem


minimize f (x)
(NLP)=
subject to h(x) = 0,

where f and h := (h1 , h2 , . . . , hm )T are assumed to be of class C 2 . The

penalty for constraint violation will be a continuous function ψ with the
following property:

ψ(x) = 0, if x ∈ M (feasibility),
<∗>
ψ(x) > 0, otherwise.
The best-known such penalty is the quadratic-loss function

m
1X 2 1
(a) ψ(x) = hi (x) = hT (x)h(x).
2 i=1 2

Also possible is a penalty of the form

m
1X γ
(b) ψ(x) = | h (x) | , where γ ≥ 1.
γ i=1 i

The wight of the penalty is controlled by a penalty parameter ρ. As ρ in-

creases, ρψ approaches the “ideal penalty” σ: ρψ(x) → σ(x) (ρ → ∞). By
adding ρψ to f we obtain the penalty function
CHAPTER 5. NONLINEAR PROGRAMMING: PENALTY AND
92 BARRIER METHODS

π(x, ρ) := f (x) + ρψ(x).

The penalty method consists in solving a sequence of unconstrained problems
of the form

(NLP)kuncon minimize π(x, ρk )

x

for an increasing sequence ρk k∈N0
of ρk > 0 tending to infinity:

ρk ≤ ρk+1 ∀k ∈ N0 ,

ρk → ∞ (k → ∞).

In general, the minimizers of the penalty function violate the constraints

h(x) = 0. The growing penalty gradually forces these minimizers towards
the feasible region.
Penalty methods share many of the advantages of the barrier methods:

(ii) Under mild conditions, the sequence of penalty function minimizers

defines a continuous trajectory.

(iii) In the latter case, it is possible to get estimates of the Lagrange mul-
tipliers at the solution.

For example, consider the quadratic-loss penalty function

m
1 X 2
π(x, ρ) = f (x) + ρ h (x).
2 i=1 i
Its minimizer x(ρ) satisfies

m
X
∇x π(x(ρ), ρ) = ∇f (x(ρ)) + ρ ∇hi (x(ρ))hi (x(ρ))
i=1
= 0,
that is,
5.1. CLASSICAL PENALTY AND BARRIER METHODS 93

m 
X 
∇f (x(ρ)) − − ρhi (x(ρ)) ∇hi (x(ρ)) = 0.
i=1 | {z }
=:λi (ρ)

If x(ρ) converges to x∗ being a solution and regular with respect to the

constraints, then, λ(ρ) converges to the Lagrange multiplier λ∗ associated
with x∗ :

x(ρ) → x∗ (ρ → ∞) =⇒ λ(ρ) → λ∗ (ρ → ∞).

Penalty functions suffer from the same problems of the ill conditioning
as do barrier functions! As ρ increases, the condition number of the Hes-
sian matrix of π(x(ρ), ρ) increases, tending to ∞ as ρ → ∞. Therefore the
unconstrained minimization can become increasingly difficult to solve.
It is possible to apply penalty methods to problems with inequality con-
straints. Consider in class C 2 , e.g.,


minimize f (x)
(NLP)≥
subject to g(x) ≥ 0,

where g := (g1 , g2 , . . . , gs )T . Any continuous function which satisfies < ∗ >

can serve as a penalty. Thus, for example, the quadratic-loss penalty in this
case is

s
1X 2
ψ(x) = min {gj (x), 0} .
2 j=1

This function has continuous first derivatives

s
X 
∇ψ(x) = min {gj (x), 0} ∇gj (x),
j=1

but its second derivatives can be discontinuous at points where some con-
straint gj is satisfied exactly. The same observation holds for other simple
forms of the penalty function. Thus, we cannot safely use Newton’s method
to minimize the function. For this reason, straightforward penalty methods
have not been widely used for solving general inequality-constrained prob-
lems.
CHAPTER 5. NONLINEAR PROGRAMMING: PENALTY AND
94 BARRIER METHODS

Convergence of Penalization Methods:

Let us focus on the convergence of barrier methods when applied to the
inequality-constrained problem


minimize f (x)
(NLP)≥
subject to g(x) ≥ 0.

We denote the feasible set by M again, and write


Msf := x ∈ Rn gj (x) > 0 (j ∈ {1, 2, . . . , s}) ;
generically, under some regularity condition: Msf = M 0 .
In preparing our convergence theorem, we make the following assump-
tions:

(A1) f and g are continuous.


(A2) ∀α ∈ R : x ∈ Rn x ∈ M, f (x) ≤ α

(A3) Msf 6= ∅.

(A4) M = Msf (topological closure).

Assumptions (A1-2) imply that f has a minimum value on the set M (why?).
Assumption (A3) is necessary to define the barrier subproblems. Assumption
(A4) is necessary to avoid situations where the minimum point is isolated,
and does not have neighbouring interior points.
For example, consider f (x) := x to be minimized subject to the con-
straints x2 − 1 ≥ 0 and x + 1 ≥ 0. Then

M = [1, ∞) ∪ {−1},
where the isolated point x∗ = −1 is the minimizer, but because it is isolated
it is not possible to approach it from the interior of the feasible region, and
a barrier method could not converge to this solution. 
The barrier function will be of the form

β(x, τ ) = f (x) + τ ϕ(x),

where ϕ can be any function that is continuous on M 0 , and that satisfies
5.1. CLASSICAL PENALTY AND BARRIER METHODS 95

ϕ(x) → ∞ as gj (x) → 0+ .

We will state that under mild conditions, the sequence of barrier mini-
mizers has a convergent subsequence, and the limit of any such convergent
subsequence is a solution to the problem. Although in practice convergence
of the entire sequence of minimizers is observed, from a theoretical point of
view it is not always possible to guarantee convergence of the entire sequence,
but only convergence of some subsequence:
Example 5.4. Consider the one-variable problem


minimize − x2
(NLP)≥
subject to 1 − x2 ≥ 0.

The logarithmic barrier function is β(x, τ ) = −x2 − τ ln(1 √

− x2 ). It has a
single minimizer x∗ = 0 if τ = 1, and two minimizers x∗ = ∓ 1 − τ if τ < 1.
(The point x∗ = 0 is a local minimizer if τ < 1) Suppose that τ0 , τ1 , τ2 , . . . is
a sequence of decreasing barrier parameters < 1. Then, a possible sequence
of minimizers of β(·, τk ) is

xk = (−1)k 1 − τk .
This sequence oscillates between neighbourhoods of −1 and +1, and,
 hence,
nonconvergent! However, the subsequences x2ν ν∈N0 and x2ν+1 ν∈N0 both
converge to solutions of (NLP)≥ !


Theorem 5.1. Suppose that our given problem (NLP)≥ satisfies Assump-
tions (A1-4), and that a sequence of unconstrained minimization problems

x

is solved for 0 < τk+1 < τk ∀k ∈ N0 , lim τk = 0. Furthermore, suppose

k→∞
that the functions β(·, τk ) have a minimum in Msf ∀k ∈ N0 . Let this global
minimizers be called xk .
Then,

(a) f (xk+1 ) ≤ f (xk ) ∀k ∈ N0 (descent);

CHAPTER 5. NONLINEAR PROGRAMMING: PENALTY AND
96 BARRIER METHODS

(b) ϕ(xk+1 ) ≥ ϕ(xk ) ∀k ∈ N0 ;


(c) xk k∈N0 has a convergent subsequence;

(d) if xκν ν∈N0 is any convergent subsequence of unconstrained minimizers
of β(·, τκν ), then, its limit point is a global solution x∗ of (NLP)≥ .

Proof. See Nash, Sofer (1996). 

A key assumption in this theorem is that the barrier functions have a
minimizer. We now state a condition which guarantees this:

Lemma 5.1. Assume that M is compact. Then, for any fixed positive value
τ , there exists a point x(τ ) ∈ Msf which minimizes β(·, τ ).

Proof. See Nash, Sofer (1996), or exercise. 

Another case where the barrier minimizers always exist, is when the prob-
lem is convex, even when M is unbounded.
Convergence results for penalty methods can be developed in a similar
manner.
Bibliography

 G. Golub and C.V. Loan, Matrix Computations. The Johns Hopkins Uni-
versity Press, Baltimore, 1989.

 N.I.M. Gould and S. Lyffer, An introduction to algorithms for nonlinear

optimization, RAL-TR-2002-031, Computational Science and Engineer-
ing Department, Atlas Centre, Rutherford Appleton Laboratory, Oxford-
shire, 2003.

 H. Hamacher and K. Klamroth, Lineare und Netzwerkoptimierung. Lin-

ear and Network Optimization. Ein bilinguales Lehrbuch, Viewegs Ver-
lagsgesellschaft, 2000.

forthcoming.

Hill, 1996.

 J. Nocedal and S. J. Wright. Numerical Optimization, Springer Series in

Operations Research. Springer, 1999.

 P. Spellucci, Numerische Verfahren der nichtlinearen Optimierung, ISNM

Lehrbuch, Birkhäuser Verlag, 1993.
{"ft_lang_label":"__label__en","ft_lang_prob":0.8375545,"math_prob":0.99581105,"size":63499,"snap":"2019-43-2019-47","text_gpt3_token_len":21434,"char_repetition_ratio":0.13955429,"word_repetition_ratio":0.065832175,"special_character_ratio":0.33983213,"punctuation_ratio":0.14885929,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9986623,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2019-10-22T11:51:31Z\",\"WARC-Record-ID\":\"<urn:uuid:a839b855-b05e-453d-8b20-ccac411859a0>\",\"Content-Length\":\"561895\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:242e0379-56b2-41bb-99a0-5c15eb48bea2>\",\"WARC-Concurrent-To\":\"<urn:uuid:70de17b1-dc9a-4990-a438-b6b92f037908>\",\"WARC-IP-Address\":\"151.101.202.152\",\"WARC-Target-URI\":\"https://www.scribd.com/document/385637782/lecture-notes-series-optimisation-constrained-pdf\",\"WARC-Payload-Digest\":\"sha1:7FZCNLKHWLNYQN4IRG6BJDJN5MRLUU5W\",\"WARC-Block-Digest\":\"sha1:YD2GPYK5QVNGY4WXCJL3IROVTDA7XE75\",\"WARC-Identified-Payload-Type\":\"application/xhtml+xml\",\"warc_filename\":\"/cc_download/warc_2019/CC-MAIN-2019-43/CC-MAIN-2019-43_segments_1570987817685.87_warc_CC-MAIN-20191022104415-20191022131915-00397.warc.gz\"}"}