Optimization I; Chapter 3 56
Chapter 3 Quadratic Programming
3.1 Constrained quadratic programming problems
A special case of the NLP arises when the objective functional f is quadratic
and the constraints h, g are linear in x ∈ lRn. Such an NLP is called a Quadratic
Programming (QP) problem. Its general form is
minimize f(x) :=
1
2
xT Bx − xT b (3.1a)
over x ∈ lRn
subject to A1x = c (3.1b)
A2x ≤ d , (3.1c)
where B ∈ lRn×n is symmetric, A1 ∈ lRm×n , A2 ∈ lRp×n, and b ∈ lRn , c ∈
lRm , d ∈ lRp.
As we shall see in this chapter, the QP (3.1a)-(3.1c) can be solved iteratively
by active set strategies or interior point methods where each iteration requires
the solution of an equality constrained QP problem.
3.2 Equality constrained quadratic programming
If only equality constraints are imposed, the QP (3.1a)-(3.1c) reduces to
minimize f(x) :=
1
2
xT Bx − xT b (3.2a)
over x ∈ lRn
subject to Ax = c , (3.2b)
where A ∈ lRm×n , m ≤ n. For the time being we assume that A has full row
rank m.
The KKT conditions for the solution x∗ ∈ lRn of the QP (3.2a),(3.2b) give rise
to the following linear system
(
B AT
A 0
)
︸ ︷︷ ︸
=: K
(
x∗
λ∗
)
=
(
b
c
)
, (3.3)
where λ∗ ∈ lRm is the associated Lagrange multiplier.
We denote by Z ∈ lRn×(n−m) the matrix whose columns span KerA, i.e., AZ = 0.
Optimization I; Chapter 3 57
Definition 3.1 KKT matrix and reduced Hessian
The matrix K in (3.3) is called the KKT matrix and the matrix ZT BZ is
referred to as the reduced Hessian.
Lemma 3.2 Existence and uniqueness
Assume that A ∈ lRm×n has full row rank m ≤ n and that the reduced Hessian
ZT BZ is positive definite. Then, the KKT matrix K is nonsingular. Hence,
the KKT system (3.3) has a unique solution (x∗, λ∗).
Proof: The proof is left as an exercise.
•
Under the conditions of the previous lemma, it follows that the second order
sufficient optimality conditions are satisfied so that x∗ is a strict local minimizer
of the QP (3.2a),(3.2b). A direct argument shows that x∗ is in fact a global
minimizer.
Theorem 3.3 Global minimizer
Let the assumptions of Lemma 3.2 be satisfied and let (x∗, λ∗) be the unique
solution of the KKT system (3.3). Then, x∗ is the unique global solution of the
QP (3.2a),(3.2b).
Proof: Let x ∈ F be a feasible point, i.e., Ax = c, and p := x∗ − x. Then,
Ap = 0. Substituting x = x∗ − p into the objective functional, we get
f(x) =
1
2
(x∗ − p)T B(x∗ − p) − (x∗ − p)T b =
=
1
2
pT Bp − pT Bx∗ + pT b + f(x∗) .
Now, (3.3) implies Bx∗ = b− AT λ∗. Observing Ap = 0, we have
pT Bx∗ = pT (b − AT λ∗) = pT b − (Ap)T λ∗︸ ︷︷ ︸
= 0
,
whence
f(x) =
1
2
pT Bp + f(x∗) .
In view of p ∈ Ker A, we can write p = Zu , u ∈ lRn−m, and hence,
f(x) =
1
2
uT ZT BZu + f(x∗) .
Since ZT BZ is positive definite, we deduce f(x) > f(x∗). Consequently, x∗ is
the unique global minimizer of the QP (3.2a),(3.2b).
•
Optimization I; Chapter 3 58
3.3 Direct solution of the KKT system
As far as the direct solution of the KKT system (3.3) is concerned, we distinguish
between symmetric factorization and the range-space and null-space approach.
3.3.1 Symmetric indefinite factorization
A possible way to solve the KKT system (3.3) is to provide a symmetric fac-
torization of the KKT matrix according to
P T KP = LDLT , (3.4)
where P is an appropriately chosen permutation matrix, L is lower triangular
with diag(L) = I, and D is block diagonal.
Based on (3.4), the KKT system (3.3) is solved as follows:
solve Ly = P T
(
b
c
)
, (3.5a)
solve Dŷ = y , (3.5b)
solve LT ỹ = ŷ , (3.5c)
set
(
x∗
λ∗
)
= P ỹ . (3.5d)
3.3.2 Range-space approach
The range-space approach applies, if B ∈ lRn×n is symmetric positive definite.
Block Gauss elimination of the primal variable x∗ leads to the Schur complement
system
AB−1AT λ∗ = AB−1b − c (3.6)
with the Schur complement S ∈ lRm×m given by S := AB−1AT . The range-
space approach is particularly effective, if
• B is well conditioned and easily invertible (e.g., B is diagonal or block-
diagonal),
• B−1 is known explicitly (e.g., by means of a quasi-Newton updating for-
mula),
• the number m of equality constraints is small.
Optimization I; Chapter 3 59
3.3.3 Null-space approach
The null-space approach does not require regularity of B and thus has a wider
range of applicability than the range-space approach.
We assume that A ∈ lRm×n has full row rank m and that ZT BZ is positive
definite, where Z ∈ lRn×(n−m) is the matrix whose columns span Ker A which
can be computed by QR factorization (cf. Chapter 2.4).
We partition the vector x∗ according to
x∗ = Y wY + ZwZ , (3.7)
where Y ∈ lRn×m is such that [Y Z] ∈ lRn×n is nonsingular and wY ∈ lRm , wZ ∈
lRn−m.
Substituting (3.7) into the second equation of (3.3), we obtain
Ax∗ = AY wY + AZ︸︷︷︸
= 0
wZ = c , (3.8)
i.e., Y wY is a particular solution of Ax = c.
Since A ∈ lRm×n has rank m and [Y Z] ∈ lRn×n is nonsingular, the product
matrix A[Y Z] = [AY 0] ∈ lRm×m is nonsingular. Hence, wY is well determined
by (3.8).
On the other hand, substituting (3.7) into the first equation of (3.3), we get
BY wY + BZwZ + A
T λ∗ = b .
Multiplying by ZT and observing ZT AT = (AZ)T = 0 yields
ZT BZwZ = Z
T b − ZT BY wY . (3.9)
The reduced KKT system (3.9) can be solved by a Cholesky factorization of the
reduced Hessian ZT BZ ∈ lR(n−m)×(n−m). Once wY and wZ have been computed
as the solutions of (3.8) and (3.9), x∗ is obtained according to (3.7).
Finally, the Lagrange multiplier turns out to be the solution of the linear system
arising from the multiplication of the first equation in (3.7) by Y T :
(AY )T λ∗ = Y T b − Y T Bx∗ . (3.10)
3.4 Iterative solution of the KKT system
If the direct solution of the KKT system (3.3) is computationally too costly,
the alternative is to use an iterative method. An iterative solver can be ap-
plied either to the entire KKT system or, as in the range-space and null-space
approach, use the special structure of the KKT matrix.
Optimization I; Chapter 3 60
3.4.1 Krylov methods
The KKT matrix K ∈ lR(n+m)×(n+m) is indefinite. In fact, if A has full row rank
m, K has n positive and m negative eigenvalues. Therefore, for the iterative
solution of (3.3) Krylov subspace methods like GMRES (Generalized Minimum
RESidual) and QMR (Quasi Minimum Residual) are appropriate candidates.
3.4.2 Transforming range-space iterations
We assume B ∈ lRn×n to be symmetric positive definite and suppose that B̃ is
some symmetric positive definite and easily invertible approximation of B such
that B̃−1B ∼ I.
We choose KL ∈ lR(n+m)×(n+m) as the lower triangular block matrix
KL =
(
I 0
−AB̃−1 I
)
, (3.11)
which gives rise to the regular splitting
KLK =
(
B̃ AT
0 S̃
)
︸ ︷︷ ︸
=: M1
−
(
B̃(I − B̃−1B) 0
A(I − B̃−1B) 0
)
︸ ︷︷ ︸
=: M2 ∼ 0
, (3.12)
where S̃ ∈ lRm×m is given by
S̃ := − AB̃−1AT . (3.13)
We set
ψ := (x, λ)T , α := (b, c)T .
Given an iterate ψ(0) ∈ lRn+m, we compute ψ(k) , k ∈ lN, by means of the
transforming range-space iterations
ψ(k+1) = ψ(k) + M−11 KL(α − Kψ(k)) = (3.14)
= (I − M−11 KLK)ψ(k) + M−11 KLα , k ≥ 0 .
The transforming range-space iteration (3.14) will be implemented as follows:
d(k) = (d
(k)
1 , d
(k)
2 )
T := α − Kψ(k) , (3.15a)
KLd
(k) = (d
(k)
1 , − AB̃−1d(k)1 + d(k)2 )T , (3.15b)
M1ϕ
(k) = KLd
(k) , (3.15c)
ψ(k+1) = ψ(k) + ϕ(k) . (3.15d)
Optimization I; Chapter 3 61
3.4.3 Transforming null-space iterations
We assume that x ∈ lRn and λ ∈ Rm admit the decomposition
x = (x1, x2)
T , x1 ∈ lRm1 , x2 ∈ lRn−m1 , (3.16a)
λ = (λ1, λ2)
T , λ1 ∈ lRm1 , λ2 ∈ lRm−m1 , (3.16b)
and that A ∈ lRm×n and B ∈ lRn×n can be partitioned by means of
A =
(
A11 A12
A21 A22
)
, B =
(
B11 B12
B21 B22
)
, (3.17)
where A11, B11 ∈ lRm1×m1 with nonsingular A11.
Partitioning the right-hand side in (3.3) accordingly, the KKT system takes the
form
B11 B12 | AT11 AT21
B21 B22 | AT12 AT22
−− −− | −− −−
A11 A12 | 0 0
A21 A22 | 0 0
x∗1
x∗2
−
λ∗1
λ∗2
=
b1
b2
−
c1
c2
. (3.18)
We rearrange (3.18) by exchanging the second and third rows and columns
B11 A
T
11 | B12 AT21
A11 0 | A12 0
−− −− | −− −−
B21 A
T
12 | B22 AT22
A21 0 | A22 0
x∗1
λ∗1
−
x∗2
λ∗2
=
b1
c1
−
b2
c2
. (3.19)
Observing B12 = B
T
21, in block form (3.19) can be written as follows
(
A BT
B D
)
︸ ︷︷ ︸
=: K
(
ψ∗1
ψ∗2
)
︸ ︷︷ ︸
=: ψ∗
=
(
α1
α2
)
︸ ︷︷ ︸
=: α
, (3.20)
where ψ∗i := (x
∗
i , λ
∗
i )
T , αi := (bi, ci)
T , 1 ≤ i ≤ 2.
We note that block Gauss elimination in (3.20) leads to the Schur complement
system
(
A BT
0 S
)(
ψ∗1
ψ∗2
)
=
(
α1
α2 − Ba−1α1
)
. (3.21)
The Schur complement S is given by
S = D − BA−1BT , (3.22)
Optimization I; Chapter 3 62
where
A−1 =
(
0 A−111
A−T11 −A−T11 B11A−111
)
. (3.23)
We assume that Ã11 ∈ lRm1×m1 is an easily invertible approximation of A11 and
define
à =
(
B11 Ã
T
11
Ã11 0
)
. (3.24)
We remark that the inverse of à is given as in (3.23) with A−111 , A−T11 replaced
by Ã−111 and Ã
−T
11 , respectively.
We introduce the right transform
KR =
(
I −Ã−1BT
0 I
)
, (3.25)
which gives rise to the regular splitting
KKR =
(
à 0
B S̃
)
︸ ︷︷ ︸
=: M1
−
(
(I −AÃ−1)Ã (−I +AÃ−1)BT
0 0
)
︸ ︷︷ ︸
=: M2 ∼ 0
, (3.26)
where
S̃ := D − BÃ−1BT . (3.27)
Given a start iterate ψ(0) = (ψ
(0)
1 , ψ
(0)
2 )
T , we solve the KKT system (3.20) by
the transforming null-space iterations
ψ(k+1) = ψ(k) + KRM
−1
1 (α − Kψ(k)) = (3.28)
= (I − KRM−11 K)ψ(k) + KRM−11 α , k ≥ 0 .
3.5 Active set strategies for convex QP problems
3.5.1 Primal active set strategies
We consider the constrained QP problem
minimize f(x) :=
1
2
xT Bx − xT b (3.29a)
over x ∈ lRn
subject to Cx = c (3.29b)
Ax ≤ d , (3.29c)
Optimization I; Chapter 3 63
where B ∈ lRn×n is symmetric positive definite, C ∈ lRm×n , A ∈ lRp×n, and
b ∈ lRn , c ∈ lRm , d ∈ lRp.
We write the matrices A and C in the form
A =
a1
·
ap
, ai ∈ lRn , C =
c1
·
cm
, ci ∈ lRn . (3.30)
The inequality constraints (3.29c) can be equivalently stated as
aTi x ≤ di , 1 ≤ i ≤ m . (3.31)
The primal active set strategy is an iterative procedure:
Given a feasible iterate x(ν) , ν ≥ 0, we determine an active set
Iac(x(ν)) ⊂ {1, …, p} (3.32)
and consider the corresponding constraints as equality constraints, whereas the
remaining inequality constraints are disregarded. Setting
p = x(ν) − x , b(ν) = Bx(ν) − b , (3.33)
we find
f(x) = f(x(ν) − p) = 1
2
pT Bp − (b(ν))T p + g ,
where g := 1
2
(x(ν))T Bx(ν) − bT x(ν).
The equality constrained QP problem to be solved at the (ν+1)-st iteration
step is then:
minimize
1
2
pT Bp − (b(ν))T p (3.34a)
over p ∈ lRn
subject to Cp = 0 (3.34b)
aTi p = 0 , i ∈ Iac(x(ν)) , (3.34c)
We denote the solution of (3.34a)-(3.34c) by p(ν). The new iterate x(ν+1) is then
obtained according to
x(ν+1) = x(ν) − ανp(ν) , αν ∈ [0, 1] , (3.35)
where αν is chosen such that x
(ν+1) stays feasible.
In particular, for i ∈ Iac(x(ν)) we have
aTi x
(ν+1) = aTi x
(ν) − αν aTi p(ν) = aTi x(ν) ≤ di .
Optimization I; Chapter 3 64
On the other hand, if aTi p
(ν) ≥ 0 for some i /∈ Iac(x(ν)), it follows that
aTi x
(ν+1) = aTi x
(ν) − αν aTi p(ν) ≤ aTi x(ν) ≤ di .
Finally, if aTi p
(ν) < 0 for i /∈ Iac(x(ν)), we have
aTi x
(ν+1) = aTi x
(ν) − αν aTi p(ν) ≤ di ⇐⇒ αν ≤
aTi x
(ν) − di
aTi p
(ν)
.
Consequently, in order to guarantee feasibility, we choose
αν := min (1, min
i/∈Iax(x(ν))
aTi x
(ν) − di
aTi p
(ν)
) . (3.36)
aTi p
(ν)<0
We define
Ib`(p(ν)) := {i /∈ Iac(x(ν)) | aTi p(ν) < 0 , min
i/∈Iax(x(ν))
aTi x
(ν)−di
aTi p
(ν)
< 1} . (3.37)
Clearly,
aTi x
(ν+1) = aTi
(
x(ν) − ανp(ν)
)
= di , i ∈ Ib`(x(ν)) .
Therefore, Ib`(p(ν)) is referred to as the set of blocking constraints.
We specify Iac(x(ν+1)) by adding to Iac(x(ν)) the most restrictive blocking con-
straint:
Iac(x(ν+1)) := (3.38)
Iac(x(ν)) ∪ {j ∈ Ib`(x(ν) |
aTj x
(ν) − dj
aTj p
(ν)
= min
i/∈Iac(x(ν))
aTi x
(ν) − di
aTi p
(ν)
}.
aTi p
(ν)<0
Further information with respect to a proper specification of the set of active
constraints is provided by systematically checking the KKT conditions:
Assume that p(ν) = 0 is the solution of the QP problem (3.34a)-(3.34c). Since
p(ν) satisfies the KKT conditions associated with that QP problem, there exist
Lagrange multipliers λ(ν) ∈ lRm and µ(ν)i , i ∈ Iac(x(ν)), such that
−b(ν) = −(Bx(ν) − b) = −
m∑
i=1
λ
(ν)
i ci −
∑
i∈Iac(x(ν))
µ
(ν)
i ai . (3.39)
If we set
µ
(ν)
i := 0 , i ∈ {1, ..., p} \ Iac(x(ν)) ,
Optimization I; Chapter 3 65
it is clear that x(ν), λ(ν), and µ(ν) satisfy the first KKT condition with respect
to the original QP problem (3.29a)-(3.29c).
Since x(ν) is feasible, the second and third KKT conditions also hold true.
We check the fourth KKT condition in terms of the sign of the multiplier µ(ν):
If
µ
(ν)
i ≥ 0 , i ∈ Iac(x(ν)) ,
the fourth KKT condition holds true. Consequently, x(ν) is a strict local mini-
mum, since B is symmetric positive definite.
On the other hand, if there exists j ∈ Iac(x(ν)) such that
µ
(ν)
j < 0 ,
we remove that constraint from the active set. We show that this strategy
produces a direction p in the subsequent iteration step that is feasible with
respect to the dropped constraint:
Theorem 3.4 Feasible and descent direction for active set strategy
We assume that x̂ ∈ lRn satisfies the KKT conditions for the QP problem
(3.34a)-(3.34c), i.e., in particular (3.39) holds true. Further, assume that the
constraint gradients ci, 1 ≤ i ≤ m, ai, i ∈ Iac(x̂) are linearly independent. Fi-
nally, suppose there is j ∈ Iac(x̂) such that µ̂j < 0.
Let p is the solution of the QP problem
minimize
1
2
pT Bp − b̂T p (3.40a)
over p ∈ lRn
subject to Cp = 0 (3.40b)
aTi p = 0 , i ∈ Iac(x̂) \ {j} . (3.40c)
where b̂ := Bx̂− b.
Then, p is a feasible direction for the constraint j, i.e.,
aTj p ≤ 0 .
Moreover, if p satisfies the second order sufficient optimality conditions for
(3.40a)-(3.40c), then aTj p < 0, i.e., p is a descent direction.
Proof. Since p is a solution of the QP problem (3.40a)-(3.40c), there exist
Lagrange multipliers λ̃i, 1 ≤ i ≤ m, and µ̃i, i ∈ Iac(x̂), i 6= j, such that
Bp − b̂ = −
m∑
i=1
λ̃ici −
∑
i∈Iac(x̂),i 6=j
µ̃iai . (3.41)
Optimization I; Chapter 3 66
Let Z be the null space basis of the matrix
[
[ci]
T
1≤i≤m [ai]
T
i∈Iac(x̂),i 6=j
]
.
The second order necessary optimality conditions imply that
ZT BZ
is positive semidefinite. Since p has the form p = ZpZ for some vector pZ , we
deduce
pT Bp ≥ 0 .
Since we have assumed that (3.39) holds true (with b(ν), x(ν), λ
(ν)
i , µ
(ν)
i replaced
by b̂, x̂, λ̂i, µ̂i), subtraction of (3.39) from (3.41) yields
Bp = −
m∑
i=1
(λ̃i − λ̂i)ci −
∑
i∈Iac(x̂),i 6=j
(µ̃i − µ̂i)ai + µ̂jaj . (3.42)
Forming the inner product with p and observing cTi p = 0, 1 ≤ i ≤ m, and
aTi p = 0, i ∈ Iac(x̂), i 6= j, we find
pT Bp = µ̂ja
T
j p . (3.43)
Since µ̂j < 0, we must have a
T
j p ≤ 0.
If the second order sufficient optimality conditions are satisfied , it follows from
(3.43) that
aTj p = 0 ⇐⇒ pT Bp = pTZZT BZpZ = 0 ⇐⇒ pZ = 0 ,
which implies p = 0. Due to the linear independence of the constraint gradients,
then (3.42) gives us µ̂j = 0, which is a contradiction. Hence, we must have
pT Bp > 0, whence aTj p < 0.
•
Corollary 3.5 Strictly decreasing objective functional
Suppose that p(ν) 6= 0 is a solution of the quadratic subprogram (3.34a)-(3.34c)
and satisfies the second order sufficient optimality conditions. Then, the ob-
jective functional of the original QP problem is strictly decreasing along the
direction p(ν).
Proof. Let us denote by Z the null space basis matrix associated with (3.34a)-
(3.34c). Then ZT BZ is positive definite, and we find that p(ν) is the unique
global minimizer of (3.34a)-(3.34c).
On the other hand, p = 0 is also a feasible point for (3.34a)-(3.34c). Conse-
quently, the value of the objective functional at p = 0 must be larger, i.e.,
1
2
(p(ν))T Bp(ν) − (b(ν))T p(ν) < 0 .
Optimization I; Chapter 3 67
Since (p(ν))T Bp(ν) ≥ 0 and αν ∈ [0, 1], we obtain
1
2
αν(p
(ν))T Bp(ν) − (b(ν))T p(ν) < 0 .
It follows that
1
2
(x(ν) − ανp(ν))T )B(xν) − ανp(ν)) − bT (x(ν) + ανp(ν)) =
=
1
2
(x(ν))T Bx(ν) − bT x(ν) + αν
[1
2
αν(p
(ν))T Bp(ν) − (b(ν))T p(ν)
]
<
<
1
2
(x(ν))T Bx(ν) − bT x(ν) .
•
As far as the specification of the active set is concerned, after having com-
puted the Lagrange multipliers λ(ν), µ(ν), if p(ν) = 0, usually the most negative
multiplier µ
(ν)
i is removed from the active set. This leads us to the following
algorithm:
Primal active set strategy
Step 1: Compute a feasible starting point x(0) and determine the set Iac(x(0))
of active inequality constraints.
Step 2: For ν ≥ 0, proceed as follows:
Compute p(ν) as the solution of (3.34a)-(3.34c).
Case 1: If p(ν) = 0, compute the multipliers λ
(ν)
i , 1 ≤ i ≤ m,µ
(ν)
i , i ∈ Iac(x(ν))
that satisfy (3.39).
If µ
(ν)
i ≥ 0, i ∈ Iac(x(ν)), then stop the algorithm.
The solution is x∗ = x(ν).
Otherwise, determine j ∈ Iac(x(ν)) such that
µ
(ν)
j = min
i∈Iac(x(ν))
µ
(ν)
i .
Set x(ν+1) = x(ν) and Iac(x(ν+1)) := Iac(x(ν)) \ {j}.
Case 2: If p(ν) 6= 0, compute αν and set
x(ν+1) := x(ν) − ανp(ν).
In case of blocking constraints, compute Iac(x(ν+1)) according to (3.38).
There are several techniques to compute an initial feasible point x(0) ∈ F . A
common one requires the knowledge of some approximation x̃ of a feasible point
which should not be ”too infeasible”. It amounts to the solution of the following
Optimization I; Chapter 3 68
linear programming problem:
minimize eT z , (3.44a)
over (x, z) ∈ lRn × lRm+p
subject to cTi x + γizi = ci , 1 ≤ i ≤ m , (3.44b)
aTi x− γm+izm+i ≤ di , 1 ≤ i ≤ p , (3.44c)
z ≥ 0 , (3.44d)
where
e = (1, ..., 1)T , γi =
{
−sign(cTi x̃− ci) , 1 ≤ i ≤ m
1 , m + 1 ≤ i ≤ m + p . (3.45)
A feasible starting point for this linear problem is given by
x(0) = x̃ , z
(0)
i =
{
|cTi x̃− ci| , 1 ≤ i ≤ m
max(aTi−mx̃− di−m, 0) , m + 1 ≤ i ≤ m + p
.(3.46)
Obviously, the optimal value of the linear programming subproblem is zero, and
any solution provides a feasible point for the original one.
Another technique introduces a measure of infeasibility in the objective func-
tional in terms of a penalty parameter β > 0:
minimize
1
2
xT Bx− xT b + βt , (3.47a)
over (x, t) ∈ lRn × lR
subject to cTi x− ci ≤ t , 1 ≤ i ≤ m , (3.47b)
− (cTi x− ci) ≤ t , 1 ≤ i ≤ m , (3.47c)
aTi x− di ≤ t , 1 ≤ i ≤ p , (3.47d)
t ≥ 0 . (3.47e)
For sufficiently large penalty parameter β > 0, the solution of (3.47a)-(3.47e) is
(x, 0) with x solving the original quadratic programming problem.
3.5.2 Primal-dual active set strategies
We consider a primal-dual active set strategy which does not require feasibility
of the iterates. It is based on a Moreau-Yosida type approximation of the
indicator function of the convex set of inequality constraints
K := {v ∈ lRp | vi ≤ 0 , 1 ≤ i ≤ p} . (3.48)
The indicator function IK : K → lR of K is given by
IK(v) :=
{
0 , v ∈ K
+∞ , v /∈ K . (3.49)
Optimization I; Chapter 3 69
The complementarity conditions
aTi x
∗ − di ≤ 0 , µ∗i ≥ 0 ,
µ∗i (a
T
i x
∗ − di) = 0 , 1 ≤ i ≤ p
can be equivalently stated as
µ∗ ∈ ∂IK(v∗) , v∗i := aTi x∗ − di , 1 ≤ i ≤ p , (3.50)
where ∂IK denotes the subdifferential of the indicator function IK as given by
µ ∈ ∂IK(v) ⇐⇒ IK(v) + µT (w − v) ≤ IK(w) , w ∈ lRp . (3.51)
Using the generalized Moreau-Yosida approximation of the indicator function
IK , (3.50) can be replaced by the computationally more feasible condition
µ∗ = σ
[
v∗ + σ−1µ∗ − PK(v∗ + σ−1µ∗)
]
, (3.52)
where σ is an appropriately chosen positive constant and PK denotes the pro-
jection onto K as given by
PK(w) :=
{
wi , wi < 0
0 , wi ≥ 0 , 1 ≤ i ≤ p .
Note that (3.52) can be equivalently written as
µ∗i = σ max
(
0, aTi x
∗ − di +
µ∗i
σ
)
. (3.53)
Now, given startiterates x(0), λ(0), µ(0), the primal-dual active set strategy pro-
ceeds as follows: For ν ≥ 1 we determine the set Iac(x(ν)) of active constraints
according to
Iac(x(ν)) := {1 ≤ i ≤ p | aTi x(ν) − di +
µ
(ν)
i
σ
> 0} . (3.54)
We define
Iin(x(ν)) := {1, …, p} \ Iac(x(ν)) (3.55)
and set
p(ν) := card Iac(x(ν)) , (3.56)
µ
(ν+1)
i := 0 , i ∈ Iin(x(ν)) . (3.57)
We compute (x(ν+1), λ(ν+1), µ̃(ν+1)) ∈ lRn × lRm × lRp(ν) as the solution of the
KKT system associated with the equality constrained quadratic programming
Optimization I; Chapter 3 70
problem
minimize Q(x) :=
1
2
xT Bx − xT b (3.58a)
over x ∈ lRn
subject to cTi x = ci , 1 ≤ i ≤ m , (3.58b)
aTi x = di , i ∈ Iac(x(ν)) . (3.58c)
Since feasibility is not required, any startiterate (x(0), λ(0), µ(0)) can be chosen.
The following results show that we can expect convergence of the primal-dual
active set strategy, provided the matrix B is symmetric positive definite.
Theorem 3.3 Reduction of the objective functional
Assume B ∈ lRn×n to be symmetric, positive definite and refer to ‖ · ‖E :=
((·)T B(·))1/2 as the associated energy norm. Let x(ν), ν ≥ 0, be the iterates
generated by the primal-dual active set strategy. Then, for Q(x) := 1
2
xT Bx−xT b
and ν ≥ 1 there holds
Q(x(ν)) − Q(x(ν−1)) = (3.59)
= −1
2
‖x(ν) − x(ν−1)‖2E −
∑
i∈Iac(x(ν)
µ
(ν)
i (di − aTi x(ν−1)) ≤ 0 .
i/∈Iac(x(ν−1))
Proof. Observing the KKT conditions for (3.58a)-(3.58c), we obtain
Q(x(ν)) − Q(x(ν−1)) =
= − 1
2
‖x(ν) − x(ν−1)‖2E + (x(ν) − x(ν−1))T (Bx(ν) − b) =
= − 1
2
‖x(ν) − x(ν−1)‖2E −
∑
i∈Iac(x(ν))
µ
(ν)
i a
T
i (x
(ν) − x(ν−1)) .
For i ∈ Iac(x(ν)) we have aTi x(ν) = di, whereas aTi x(ν−1) = di for i ∈ Iac(x(ν−1)).
we thus get
Q(x(ν)) − Q(x(ν−1)) =
= − 1
2
‖x(ν) − x(ν−1)‖2E −
∑
i∈Iac(x(ν))
µ
(ν)
i (di − aTi x(ν−1)) .
i/∈Iac(x(ν−1))
But µ
(ν)
i ≥ 0, i ∈ Iac(x(ν)) and aTi x(ν−1) ≤ di, i /∈ Iac(x(ν−1)) which gives the
assertion.
•
Corollary 3.4 Convergence of a subsequence to a local minimum
Optimization I; Chapter 3 71
Let x(ν), ν ≥ 0, be the iterates generated by the primal-dual active set strategy.
Then, there exist a subsequence lN′ ⊂ lN and x∗ ∈ lRn such that x(ν) → x∗, ν →
∞, ν ∈ lN′. Moreover, x∗ is a local minimizer of (3.29a)-(3.29c).
Proof. The sequence (x(ν))lN is bounded, since otherwise we would have
Q(x(ν)) → +∞ as ν → ∞ in contrast to the result of the previous theo-
rem. Consequently, there exist a subsequence lN′ ⊂ lN and x∗ ∈ lRn such
that x(ν) → x∗, ν → ∞, ν ∈ lN′. Passing to the limit in the KKT system for
(3.58a)-(3.58c) shows that x∗ is a local minimizer.
•
The following result gives an a priori estimate in the energy norm.
Theorem 3.5 A priori error estimate in the energy norm
Let x(ν), ν ≥ 0, be the iterates generated by the primal-dual active set strategy
and lN′ ⊂ lN such that x(ν) → x∗, ν →∞, ν ∈ lN′. Then, there holds
‖x(ν) − x∗‖2E ≤ 2
∑
i∈Iac(x∗)
µ∗i (di − aTi x(ν)) . (3.60)
i/∈Iac(x(ν))
Proof. The proof is left as an exercise.
3.6 Interior-point methods
Interior-point methods are iterative schemes where the iterates approximate a
local minimum from inside the feasible set. For ease of exposition, we restrict
ourselves to inequality constrained quadratic programming problems of the form
minimize Q(x) :=
1
2
xT Bx − xT b (3.61a)
over x ∈ lRn
subject to Ax ≤ d, (3.61b)
where B ∈ lRn×n is symmetric, positive semidefinite, A = [ai]1≤i≤p ∈ lRp×n , b ∈
lRn, and d ∈ lRp.
We already know from Chapter 2 that the KKT conditions for (3.61a)-(3.61b)
can be stated as follows:
If x∗ ∈ lRn is a solution of (3.61a)-(3.61b), there exists a multiplier µ∗ ∈ lRp
such that
Bx∗ + AT µ∗ − b = 0 , (3.62a)
Ax∗ − d ≤ 0 , (3.62b)
µ∗i (Ax
∗ − d)i = 0 , 1 ≤ i ≤ p , (3.62c)
µ∗i ≥ 0 , 1 ≤ i ≤ p . (3.62d)
Optimization I; Chapter 3 72
By introducing the slack variable z := d − Ax, the above conditions can be
equivalently formulated as follows
Bx∗ + AT µ∗ − b = 0 , (3.63a)
Ax∗ + z∗ − d = 0 , (3.63b)
µ∗i zi = 0 , 1 ≤ i ≤ p , (3.63c)
z∗i , µ
∗
i ≥ 0 , 1 ≤ i ≤ p . (3.63d)
We can rewrite (3.63a)-(3.63d) as a constrained system of nonlinear equations.
We define the nonlinear map
F (x, µ, z) =
Bx + AT µ− b
Ax + z − d
ZDµe
, z, µ ≥ 0 , (3.64)
where
Z := diag(z1, …, zp) , Dµ := diag(µ1, …, µp) , e := (1, …, 1)
T .
Given a feasible iterate (x, µ, z), we introduce a duality measure according to
κ :=
1
p
p∑
i=1
ziµi =
zT µ
p
. (3.65)
Definition 3.6 Central path
The set of points (xτ , µτ , zτ ) , τ > 0, satisfying
F (xτ , µτ , zτ ) =
0
0
τe
, zτ , µτ > 0 (3.66)
is called the central path.
The idea is to apply Newton’s method to (3.64) to compute (xσκ, µσκ, zσκ) on
the central path, where σ ∈ [0, 1] is a parameter chosen by the algorithm. The
Newton increments (∆x, ∆z, ∆µ) are the solution of the linear system
B AT 0
A 0 I
0 Z Dµ
∆x
∆µ
∆z
=
−rb
−rd
−ZDµe + σκe
, (3.67)
where
rb := Bx + A
T µ − b , rd := Ax + z − d .
The new iterate (x, µ, z) is then determined by means of
(x, µ, z) = (x, µ, z) + α (∆x, ∆µ, ∆z) (3.68)
Optimization I; Chapter 3 73
with α chosen such that (x, µ, z) stays feasible.
Since Dµ is a nonsingular diagonal matrix, the increment ∆z in the slack variable
can be easily eliminated resulting in
(
B AT
A −D−1µ Z
)(
∆x
∆µ
)
=
(
−rb
−rd −D−1µ g
)
, (3.69)
where g := −ZDµe + σκe.
Optimization I; Chapter 3 74
3.7 Logarithmic barrier functions
We consider the inequality constrained quadratic programming problem (3.62a)-
(3.62b). Algorithms based on barrier functions are iterative methods where the
iterates are forced to stay within the interior
F int := { x ∈ lRn | aTi x − ci < 0 , 1 ≤ i ≤ p } . (3.70) Barrier functions have the following properties: • They are infinite outside F int. • They are smooth within F int. • They approach ∞ as x approaches the boundary of F int. Definition 3.7 Logarithmic barrier function For the quadratic programming problem (3.62a)-(3.62b) the objective function- als B(β)(x) := Q(x) − β p∑ i=1 log(di − aTi x) , β > 0 (3.71)
are called logarithmic barrier functions. The parameter β is referred to as the
barrier parameter.
Theorem 3.8 Properties of the logarithmic barrier function
Assume that the set S of solutions of (3.62a)-(3.62b) is nonempty and bounded
and that the interior F int of the feasible set is nonempty. Let {βk}lN be a
decreasing sequence of barrier parameters with βk → 0 as k →∞. Then there
holds:
(i) For any β > 0 the logarithmic barrier function B(β)(x) is convex in F Int
and attains a minimizer x(β) on F int. Any local minimizer x(β) is also a global
minimizer of B(β)(x).
(ii) If {x(βk)}lN is a sequence of minimizers, then there exists lN′ ⊂ lN such
that x(βk) → x∗ ∈ S, k ∈ lN′.
(iii) If Q∗ is the optimal value of the objective functional Q in (3.62a)-(3.62b),
then for any sequence {x(βk)}lN of minimizers there holds
Q(x(βk)) → Q∗ , B(βk)(x(βk)) → Q∗ as k →∞ .
Proof. We refer to M.H. Wright; Interior methods for constrained optimiza-
tion. Acta Numerica, 341-407, 1992.
•
Optimization I; Chapter 3 75
We will have a closer look at the relationship between a minimizer of B(β)(x)
and a point (x, µ) satisfying the KKT conditions for (3.62a)-(3.62b).
If x(β) is a minimizer of B(β)(x), we obviously have
∇xB(β)(x(β)) = ∇xQ(x(β)) +
p∑
i=1
β
di − aTi x(β)
ai = 0 . (3.72)
Definition 3.9 Perturbed complementarity
The vector z(β) ∈ lRp with components
z
(β)
i :=
β
di − aTi x(β)
, 1 ≤ i ≤ p (3.73)
is called perturbed or approximate complementarity.
The reason for the above definition will become obvious shortly.
Indeed, in terms of the perturbed complementarity, (3.72) can be equivalently
stated as
∇xQ(x(β)) +
p∑
i=1
z
(β)
i ai = 0 . (3.74)
We have to compare (3.74) with the first of the KKT conditions for (3.62a)-
(3.62b) which is given by
∇xL(x, µ) = ∇xQ(x) +
p∑
i=1
µiai = 0 . (3.75)
Obviously, (3.75) looks very much the same as (3.74).
The other KKT conditions are as follows:
aTi x − di ≤ 0 , 1 ≤ i ≤ p , (3.76a)
µi ≥ 0 , 1 ≤ i ≤ p , (3.76b)
µi (a
T
i x − di) = 0 , 1 ≤ i ≤ p . (3.76c)
Apparently, (3.76a) and (3.76b) are satisfied by x = x(β) and µ = z(β).
However, (3.76c) does not hold true, since it follows readily from (3.73) that
z
(β)
i (di − aTi x(β)) = β > 0 , 1 ≤ i ≤ p . (3.77)
On the other hand, as β → 0 a minimizer x(β) and the associated z(β) come
closer and closer to satisfying (3.76c). This is reason why z(β) is called perturbed
(approximate) complementarity.
Optimization I; Chapter 3 76
Theorem 3.10 Further properties of the logarithmic barrier function
Assume that F int 6= ∅ and that x∗ is a local solution of (3.62a)-(3.62b) with
multiplier µ∗ such that the KKT conditions are satisfied.
Suppose further that the LICQ, strict complementarity, and the second order
sufficient optimality conditions fold true at (x∗, µ∗). Then there holds:
(i) For sufficiently small β > 0 there exists a strict local minimizer x(β) of
B(β)(x) such that the function x(β) is continuously differentiable in some neigh-
borhood of x∗ and x(β) → x∗ as β → 0.
(ii) For the perturbed complementarity z(β) there holds:
z(β) → µ∗ as β → 0 . (3.78)
(iii) For sufficiently small β > 0, the Hessian ∇xxB(β)(x) is positive definite.
Proof. We refer to M.H. Wright; Interior methods for constrained optimiza-
tion. Acta Numerica, 341-407, 1992.
•