On the Solution of Equality Constrained Quadratic Programming Problems Arising in Optimization.
We consider the application of the conjugate gradient method to the solution of large equality constrained quadratic programs arising in nonlinear optimization. Our approach is based implicitly on a reduced linear system and generates iterates in the null space of the constraints. Instead of computing a basis for this null space, we choose to work directly with the matrix of constraint gradients, computing projections into the null space by either a normal equations or an augmented system approach. Unfortunately, in practice such projections can result in significant rounding errors. We propose iterative refinement techniques, as well as an adaptive reformulation of the quadratic problem, that can greatly reduce these errors without incurring high computational overheads. Numerical results illustrating the efficacy of the proposed approaches are presented.
Introduction
A variety of algorithms for nonlinearly constrained optimization [7, 8, 12, 29, 31] use
the conjugate gradient (CG) method [25] to solve subproblems of the form
minimize
x
subject to
In nonlinear optimization, the n-vector c usually represents the gradient rf of the objective
function or the gradient of the Lagrangian, the n \Theta n symmetric matrix H stands for either
the Hessian of the Lagrangian or an approximation to it, and the solution x represents a
search direction. The equality constraints (1.2) are obtained by linearizing the constraints
of the optimization problem at the current iterate. We will assume here that A is an m \Theta n
that A has full row rank so that the constraints (1.2) constitute
linearly independent equations. We also assume for convenience that H is positive
definite in the null space of the constraints, as this guarantees that (1.1)-(1.2) has a unique
solution. This positive definiteness assumption is not needed in trust region methods, but
our discussion will also be valid in that context because trust region methods normally
terminate the CG iteration as soon as negative curvature is encountered (see [36, 38], and,
by contrast, [23]).
The use of an iterative method such as CG is attractive in large scale optimization
because, when the number of variables is large, it can be cost effective to solve (1.1)-
(1.2) approximately, and only increase the accuracy of the solution as the iterates of the
optimization algorithm approach the minimizer. In addition, the properties of the CG
method merge very well with the requirements of globally convergent optimization methods
(see e.g. [36]). In this paper we study how to apply the preconditioned CG method to (1.1)-
(1.2) so as to keep the computational cost at a reasonable level while ensuring that rounding
errors do not degrade the performance of the optimization algorithm.
The quadratic program (1.1)-(1.2) can be solved by computing a basis Z for the null
space of A, using this basis to eliminate the constraints, and then applying the CG method
to the reduced problem. We will argue, however, that due to the form of the preconditioners
used in practice, the explicit use of Z will cause the iteration to be very expensive, and that
significant savings can be achieved by means of approaches that bypass the computation of
Z altogether. The price to pay for these alternatives is that they can give rise to excessive
roundoff errors that can slow the optimization iteration and may even prevent it from
converging.
As we shall see, these errors cause the constraints (1.2) not to be satisfied to the desired
accuracy. We describe iterative refinement techniques that can improve the accuracy of the
solution in highly ill-conditioned problems. We also propose a mechanism for redefining
the vector c adaptively that does not change the solution of the quadratic problem but that
has more favorable numerical properties.
Notation. Throughout the paper k \Delta k stands for the ' 2 matrix or vector norm, while the
G-norm of the vector x is defined to be
x T Gx, where G is a given positive-definite
matrix. We will denote the floating-point unit roundoff (or machine precision) by ffl m . We
let -(A) denote the condition number of A, i.e.
are the nonzero singular values of A.
2. The CG method and linear constraints
A common approach for solving linearly constrained problems is to eliminate the constraints
and solve a reduced problem (c.f. [17, 20]). More specifically, suppose that Z is
an n \Theta (n \Gamma m) matrix spanning the null space of A. Then the columns of A T
together with the columns of Z span R n , and any solution x   of the linear equations (1.2)
can be written as
for some vectors x A
. The constraints (1.2) yield
which determines the vector x A
. Substituting (2.1) into (1.1), and omitting constant terms
is a constant now) we see that x Z
solves the reduced problem
minimize
x Z2
where
As we have assumed that the reduced Hessian H ZZ is positive definite, (2.3) is equivalent
to the linear system
We can now apply the conjugate gradient method to compute an approximate solution of
the problem (2.3), or equivalently the system (2.4), and substitute this into (2.1) to obtain
an approximate solution of the quadratic program (1.1)-(1.2).
This strategy of computing the normal component A T x A exactly and the tangential
component Zx Z inexactly is compatible with the requirements of many nonlinear optimization
algorithms which need to ensure that, once linear constraints are satisfied, they remain
so throughout the remainder of the optimization calculation (cf. [20]).
Let us now consider the practical application of the CG method to the reduced system
(2.4). It is well known that preconditioning can improve the rate of convergence of the
CG iteration (c.f. [1]), and we therefore assume that a preconditioner W ZZ is given. W ZZ
is a symmetric, positive definite matrix of dimension which might be chosen to
reduce the span of, and to cluster, the eigenvalues of W \Gamma1
or could be the result
of an automatic scaling of the variables [7, 29]. Regardless of how W ZZ is defined, the
preconditioned conjugate gradient method applied to (2.4) is as follows (see, e.g. [20]).
Algorithm I. Preconditioned CG for Reduced Systems.
Choose an initial point x Z , compute r
ZZ
r Z and p \Gammag Z .
Repeat the following steps, until a termination test is satisfied:
r Z
Z
ZZ
r Z
\Gammag Z
Z / g Z
and r Z / r Z
This iteration may be terminated, for example, when r Z
ZZ
r Z is sufficiently small.
Once an approximate solution is obtained, it must be multiplied by Z and substituted in
(2.1) to give the approximate solution of the quadratic program (1.1)-(1.2). Alternatively,
we may rewrite Algorithm I so that the multiplication by Z and the addition of the term
A T x A
is performed explicitly in the CG iteration. To do so, we introduce, in the following
algorithm, the n-vectors x;
Algorithm II Preconditioned CG (in Expanded Form) for Reduced Systems.
Choose an initial point x satisfying (1.2), compute
\Gammag. Repeat the following steps, until a convergence test is satisfied:
x
This will be the main algorithm studied in this paper. Several types of stopping tests
can be used, but since their choice depends on the requirements of the optimization method,
we shall not discuss them here. In the numerical tests reported in this paper we will use
the quantity r T
ZZ
r Z to terminate the CG iteration.
Note that the vector g, which we call the preconditioned residual, has been explicitly
defined to be in the range of Z. As a result, in exact arithmetic, all the search directions
generated by Algorithm II will also lie in the range of Z, and thus the iterates x will all
satisfy (1.2). Rounding errors when computing (2.17) may cause p to have a component
outside the range of Z, but this component will normally be too small to cause difficulties.
3. Implementation of the Projected CG Method
Algorithm II constitutes an effective method for computing the solution to (1.1)-(1.2)
and has been successfully used in various algorithms for large scale optimization (cf. [16, 28,
39]). The main drawback is the need for a null-space basis matrix Z, whose computation
and manipulation can be costly, and which can sometimes give rise to unnecessary ill-conditioning
[9, 10, 18, 24, 33, 37]. These difficulties will become apparent when we describe
practical procedures for computing Z and when we consider the types of preconditioners
W ZZ used in practice. Let us begin with the first issue.
3.1. Computing a basis for the null space
There are many possible choices for the null-space matrix Z. Possibly the best strategy
is to choose Z so as to have orthonormal columns, for this provides a well conditioned
representation of the null space of A. However computing such a null-space matrix can be
very expensive when the number of variables is large; it essentially requires the computation
of a sparse LQ factorization of A and the implicit or explicit generation of Q, which has
always been believed to be rather expensive when compared with the alternatives described
in [24]. Recent research [30, 35] has suggested that it is in fact possible to generate Q
as a product of sparse Householder matrices, and that the cost of this may, after all,
be reasonable. We have not experimented with this approach, however, because, to our
knowledge, general purpose software implementing it is not yet available.
Another possibility is to try to compute a basis of the null-space which involves as
few nonzeros as possible. Although this problem is computationally hard [9], sub-optimal
heuristics are possible but still rather expensive [10, 18, 33, 37].
A more economical alternative is based on simple elimination of variables [17, 20]. To
define Z we first group the components of x into m basic or dependent variables (which for
simplicity are assumed to be the first m variables) and
and partition A as
where the m \Theta m basis matrix B is assumed to be nonsingular. Then we define
\GammaB
I
which clearly satisfies and has linearly independent columns. In practice Z is not
formed explicitly; instead we compute and store sparse LU factors [13] of B, and compute
products of the form Zv and Z T v by means of solves using these LU factors. Ideally we
would like to choose a basis B that is as sparse as possible and whose condition number is
not significantly worse than that of A, but these requirements can be difficult to achieve. In
simply ensuring that B is well conditioned can be difficult when the task of choosing
a basis is delegated to a sparse LU factorization algorithm such as MA48 [15]. Some recent
codes (see, e.g., [19]) have been designed to compute a well-conditioned basis, but it is not
known to us to what extent they reach their objective.
3.2. Preconditioning
These potential drawbacks of the null-space basis (3.1) are not sufficiently serious to
prevent its effective use in Algorithm II. However, when considering practical choices for
the preconditioning matrix W ZZ , one exposes the weaknesses of this approach. Ideally, one
would like to choose W ZZ so that W \Gamma1
thus
ZZ
is the perfect preconditioner. However, it is unlikely that Z T HZ or its inverse are sparse
matrices, and even if Z T HZ is of small dimension, forming it can be quite costly. Therefore
operating with this ideal preconditioner is normally out of the question.
In this paper we consider preconditioners of the form
ZZ
where G is a symmetric matrix such that Z T GZ is positive definite. Some suggestions on
how to choose G have been made in [32]. Two particularly simple choices are
The first choice is appropriate when H is dominated by its diagonal. This is the case, for
example, in barrier methods for constrained optimization that handle bound constraints
l - x - u by adding terms of the form \Gamma-
to the objective
function, for some positive barrier parameter -. The choice I arises in several trust
region methods for constrained optimization [7, 12, 29], where the preconditioner (which
derives from a change of variables) is thus given by
ZZ
Regardless of the choice of G, the preconditioner (3.3) requires operations with the
inverse of the matrix Z T GZ. In some applications [16, 39] Z, defined by (3.1), has a simple
enough structure that forming and factorizing the (n \Gamma m) \Theta (n \Gamma m) matrix Z T GZ is not
expensive when G has a simple form. But if the LU factors of B are not very sparse and
the number of constraints m is large, forming Z T GZ may be rather costly, even if
as it requires the solution of 2m triangular systems with these LU factors. In this case it is
preferable not to form Z T GZ, but rather compute products of the form (Z
solving (Z T using the CG method. This inner CG iteration has been employed
in [29] with I, and can be effective on some problems-particularly if the number
of degrees of freedom, very small. But it can fail when Z is badly conditioned
and tends to be expensive. Moreover, since the matrix Z T GZ is not known explicitly, it is
difficult to construct effective preconditioners for accelerating this inner CG iteration.
In summary when the preconditioner has the form (3.3), and when Z is defined by means
of (3.1), the computation (2.15) of the preconditioned residual g is often so expensive as
to dominate the cost of the optimization algorithm. The goal of this paper is to consider
alternative implementations of Algorithm II whose computational cost is more moderate
and predictable. Our approach is to avoid the use of the null-space basis Z altogether.
3.3. Computing Projections
To see how to bypass the computation of Z, let us begin by considering the simple case
when so that the preconditioner W ZZ is given by (3.4). If P Z denotes the orthogonal
projection operator onto the null space of A,
then the preconditioned residual (2.15) can be written as
This projection can be performed in two alternative ways.
The first is to replace P Z by the equivalent formula
and thus to replace (3.6) with
We can express this as
is the solution of
Noting that (3.10) are the normal equations, it follows that v + is the solution of the least
squares problem
minimize
and that the desired projection g + is the corresponding residual. This approach can be
implemented using a Cholesky factorization of AA T .
The second possibility is to express the projection (3.6) as the solution of the augmented
system /
I A T
!/
r +!
This system can be solved by means of a symmetric indefinite factorization that uses 1 \Theta 1
and 2 \Theta 2 pivots [21].
Let us suppose now that the preconditioner has the more general form (3.3). The
preconditioned residual (2.15) now requires the computation
This may be expressed as
if G is non-singular, and can be found as the solution of
G A T
!/
r +!
whenever Z T GZ is non-singular (see, e.g., [20, Section 5.4.1]). While (3.14) is far from
appealing when G \Gamma1 does not have a simple form, (3.15) is a useful generalization of (3.12).
Clearly the system (3.12) may be obtained from (3.15) by setting I, and the perfect
preconditioner results if other choices for G are also possible; all that is required
is that Z T GZ be positive definite. The idea of using the projection (3.7) in the CG method
dates back to at least [34]; the alternative (3.15), and its special case (3.12), are proposed
in [8], although [8] unnecessarily requires that G be positive definite. A more recent study
on preconditioning the projected CG method is [11].
Hereafter we shall write (2.15) as
where P is any of the projection operators we have mentioned above.
Note that (3.8), (3.12) and (3.15) do not make use of the null space matrix Z and only
require factorization of matrices involving A. Unfortunately they can give rise to significant
round-off errors, particularly as the CG iterates approach the solution. The difficulties are
caused by the fact that as the iterations proceed, the projected vector
increasingly small while r does not. Indeed, the optimality conditions of the quadratic
program (1.1)-(1.2) state that the solution x   satisfies
for some Lagrange multiplier vector -. The vector Hx + c, which is denoted by r in
Algorithm II, will generally stay bounded away from zero, but as indicated by (3.16), it
will become increasingly closer to the range of A T . In other words r will tend to become
orthogonal to Z, and hence, from (3.13), the preconditioned residual g will converge to zero
so long as the smallest eigenvalue of Z T GZ is bounded away from zero.
That this discrepancy in the magnitudes of will cause numerical difficulties
is apparent from (3.9), which shows that significant cancellation of digits will usually
take place. The generation of harmful roundoff errors is also apparent from (3.12)/(3.15)
because will be small while the remaining components v + remain large. Since the magnitude
of the errors generated in the solution of (3.12)/(3.15) is governed by the size of the
large component v + , the vector g + will contain large relative errors. These arguments will
be made more precise in the next section.
Example 1.
We applied Algorithm II to solve problem CVXEQP3 from the CUTE collection [4],
with In this and all subsequent experiments, we use the simple
preconditioner (3.4) corresponding to the choice used both the normal equations
(3.8) and augmented system (3.12) approaches to compute the projection. The results are
given in Figure 1, which plots the residual
r T g as a function of the iteration number. In
both cases the CG iteration was terminated when r T g became negative, which indicates
that severe errors have occurred since r T
must be positive-continuing the
iteration past this point resulted in oscillations in the norm of the gradient without any
significant improvement. At iteration 50 of both runs, r is of order 10 5 whereas its projection
g is of

Figure

also plots the cosine of the angle between the preconditioned residual g and
the rows of A. More precisely, we define
A T
where A i is the i-th row of A. Note that this cosine, which should be zero in exact arithmetic,
increases indicating that the CG iterates leave the constraint manifold
Severe errors such as these are not uncommon in optimization calculations; see x7 and
[27]. This is of grave concern as it may cause the underlying optimization algorithms to
behave erratically or fail.
In this paper we propose several remedies. One of them is based on an adaptive redefinition
of r that attempts to minimize the differences in magnitudes between
. We also describe several forms of iterative refinement for the projection operation. All
these techniques are motivated by the roundoff error analysis given next.
4. Analysis of the Errors
We now present error bounds that support the arguments made in the previous section,
particularly the claim that the most problematic situation occurs in the latter stages of the
PCG Augmented System
Iteration
resid
cos
PCG Normal Equations
Iteration
resid
cos

Figure

1: Conjugate gradient method with two options for the projection
iteration when g + is converging to zero, but r + is not. For simplicity, we shall assume
henceforth that A has been scaled so that shall only consider the
simplest possible preconditioner, as opposed to exact, quantity will
be denoted by a subscript c.
Let us first consider the normal equations approach. Here is given by (3.9)
where (3.10) is solved by means of the Cholesky factorization of AA T . In finite precision,
instead of the exact solution v + of the normal equations we obtain v
the error \Deltav
with . Recall that ffl m denotes unit roundoff and -(A) the condition number of
A.
We can now study the total error in the projection vector g + . To simplify the analysis,
we will ignore the errors that arise in the computation of the matrix-vector product A T v
and in the subtraction given in (3.9), because these errors will be dominated by
the error in v + whose magnitude is estimated by (4.1). Under these assumptions, we have
from (3.9) that the computed projection
and the exact projection
1 The bound (4.1) assumes that there are no errors in the formation of AA T and Ar + , or in the backsolves
using the Cholesky factors; this is a reasonable assumption in our context. We should also note that (4.1)
can be sharpened by replacing the term possible
diagonal scalings D.
and thus the error in the projection lies entirely in the range of A T . We then have from
(4.1) that the relative error in the projection satisfies
This error can be significant when -(A) is large or when
is large.
Let us consider the ratio (4.4) in the case when kr much larger than its projection
We have from (3.9) that kr and by the assumption that
Suppose that the inequality above is achieved. Then (4.4) gives
which is simpler to interpret than (4.4). We can thus conclude that the error in the projection
(4.3) will be large when either -(A) or the ratio kr large.
When the condition number -(A) is moderate, the contribution of the ratio (4.4) to
the relative error (4.3) is normally not large enough to cause failure of the optimization
calculation. But as the condition number -(A) grows, the loss of significant digits becomes
severe, especially since -(A) appears squared in (4.3). In Example 1,
and we have mentioned that the ratio (4.4) is of order O(10 6 ) at iteration 50. The bound
(4.3) indicates that there could be no correct digits in g + , at this stage of the CG iteration.
This is in agreement with our test, for at this point the CG iteration could make no further
progress.
Let us now consider the augmented system approach (3.15). Again we will focus on the
choice I, for which the preconditioned residual is computed by solving
I A T
!/
r +!
using a direct method. There are a number of such methods, the strategies of Bunch and
Kaufman [5] and Duff and Reid [14] being the best known examples for dense and sparse
matrices, respectively. Both form the LDL T factorization of the augmented matrix (i.e.
the matrix appearing on the left hand side of (4.5)), where L is unit lower triangular and
D is block diagonal with 1 \Theta 1 or 2 \Theta 2 blocks.
This approach is usually (but not always) more stable than the normal equations ap-
proach. To improve the stability of the method, Bj-orck [2] suggests introducing a parameter
ff and solving the equivalent system
!/
r +!
An error analysis [3] shows that
where j depends on n and m and in the growth factor during the factorization, and oe 1 -
are the nonzero singular values of A. It is important to notice that now
-(A)-and not - 2 (A)-enters in the bound. If ff - oe m (A), this method will give a solution
that is never much worse than that obtained by a tight perturbation analysis, and therefore
can be considered stable for practical purposes. But approximating oe m (A) can be difficult,
and it is common to simply use
In the case which concerns us most, when kg converges to zero while kv
the term inside the last square brackets in (4.7) is approximately kv + k, and we obtain
where we have assumed that ff = 1. It is interesting to compare this bound with (4.3). We
see that the ratio (4.4) again plays a crucial role in the analysis, and that the augmented
system approach is likely to give a more accurate solution than the method of normal
equations in this case. This cannot be stated categorically, however, since the size of the
factor j is difficult to predict.
The residual update strategy described in x6 aims at minimizing the contribution of
the ratio (4.4), and as we will see, has a highly beneficial effect in Algorithm II. Before
presenting it, we discuss various iterative refinement techniques designed to improve the
accuracy of the projection operation.
5. Iterative Refinement
Iterative refinement is known as an effective procedure for improving the accuracy of a
solution obtained by a method that is not backwards stable. We will now consider how to
use it in the context of our normal equations and augmented system approaches.
5.1. Normal Equations Approach
Let us suppose that we choose I and that we compute the projection P A r
the normal equations approach (3.9)-(3.10). An appealing idea for trying to improve the
accuracy of this computation is to apply the projection repeatedly. Therefore rather than
computing in (2.15), we let where the projection is
applied as many times as necessary to keep the errors small. The motivation for this multiple
projections technique stems from the fact that the computed projection
have only a small component, consisting entirely of rounding errors, outside of the null space
of A, as described by (4.2). Therefore applying the projection P A to the first projection
c will give an improved estimate because the ratio (4.4) will now be much smaller. By
repeating this process we may hope to obtain further improvement of accuracy.
The multiple projection technique may simply be described as setting g +
performing the following steps:
solve L(L
set
where L is the Cholesky factor of AA T . We note that this method is only appropriate when
although a simple variant is possible when G is diagonal.
Example 2.
We solved the problem given in Example 1 using multiple projections. At every CG
iteration we measure the cosine (3.17) of the angle between g and the columns of A. If this
cosine is greater than 10 \Gamma12 , then multiple projections are applied until the cosine is less
than this value. The results are given in Figure 2, and show that the residual
r T g was
reduced much more than in the plane CG iteration (Figure 1). Indeed the ratio between
the final and initial values of
r T g is 10 \Gamma16 , which is very satisfactory.
It is straightforward to analyze the multiple projections strategy (5.1)-(5.2) provided
that, as before, we make the simplifying assumption that the only rounding errors we make
are in forming L and solving (5.1). We obtain the following result which can be proved by
induction. For
where as in (4.1)
A simple consequence of (5.3)-(5.4) and the assumption that A has norm one is that
and thus that the error converges R-linearly to zero with rate
Of course, this rate can not be sustained indefinitely as the other errors we have ignored
in (5.1)-(5.2) become important. Nonetheless, one would expect (5.5) to reflect the true
behaviour until k(g
small multiple of the unit roundoff ffl m . It should
Iteration
residual

Figure

2: CG method using multiple projections in the normal equations approach.
be stressed, however, that this approach is still limited by the fact that the condition number
of A appears squared in (5.5); improvement can be guaranteed only if
We should also note that multiple projections are almost identical in their form and
numerical properties to fixed precision iterative refinement to the least squares problem [3,
p.125]. Fixed precision iterative refinement is appropriate because the approach we have
chosen to compute projections is not stable. To see this, compare (4.3) with a perturbation
analysis of the least squares problem [3, Theorem 1.4.6]), which gives
Here the dependence on the condition number is linear-not quadratic. Moreover, since
is multiplied by kg is small the effect of the condition number of A is
much smaller in (5.7) than in (4.3).
We should mention two other iterative refinement techniques that one might consider,
but that are either not effective or not practical in our context.
The first is to use fixed-precision iterative refinement [3, Section 2.9] to attempt to
improve the solution v + of the normal equations (3.10). This, however, will generally
be unsuccessful because fixed-precision iterative refinement only improves a measure of
backward stability [21, p.126], and the Cholesky factorization is already a backward stable
method. We have performed numerical tests and found no improvement from this strategy.
However, as is well known, iterative refinement will often succeed if extended-precision
is used to evaluate the residuals. We could therefore consider using extended precision
iterative refinement to improve the solution v + of the normal equations (3.10). So long as
and the residuals of (3.10) are smaller than one in norm, we can expect that
the error in the solution of (3.10) will decrease by a factor ffl m-(A) 2 until it reaches O(ffl m ).
But since optimization algorithms normally use double precision arithmetic for all their
computations, extending the precision may not be simple or efficient, and this strategy is
not suitable for general purpose software.
For the same reason we will not consider the use of extended precision in (5.1)-(5.2) or
in the iterative refinement of the least squares problem.
5.2. Augmented System Approach
We can apply fixed precision iterative refinement to the solution obtained from the
augmented system (3.15). This gives the following iteration.
Compute
solve
G A T
!/
\Deltag
ae g
ae v
and update
Note that this method is applicable for general preconditioners G. When
an appropriate value of ff is in hand, we should incorporate it in this iteration, as described
in (4.6). The general analysis of Higham [26, Theorem 3.2] indicates that, if the condition
number of A is not too large, we can expect high accuracy in v + and good accuracy in g +
in most cases.
Example 3.
We solved the problem given in Example 1 using this iterative refinement technique. As
in the case of multiple projections discussed in Example 2, we measure the angle between
g and the columns of A at every CG iteration. Iterative refinement is applied as long as
the cosine of this angle is greater than 10 \Gamma12 . The results are given in Figure 3.
We observe that the residual
r T g is decreased almost as much as with the multiple
projections approach, and attains an acceptably small value. We should point out, however,
that the residual increases after it reaches the value 10 \Gamma10 , and if the CG iteration is
continued for a few hundred more iterations, the residual exhibits large oscillations. We
will return to this in x6.1.
In our experience 1 iterative refinement step is normally enough to provide good accu-
racy, but we have encountered cases in which 2 or 3 steps are beneficial.
6. Residual Update Strategy
We have seen that significant roundoff errors occur in the computation of the projected
residual vector is much smaller than the residual r + . We now describe a procedure
Iteration
residual

Figure

3: CG method using iterative refinement in the augmented system approach.
for redefining r + so that its norm is closer to that of g + . This will dramatically reduce the
roundoff errors in the projection operation.
We begin by noting that Algorithm II is theoretically unaffected if, immediately after
computing r + in (2.14), we redefine it as
for some y This equivalence is due to the condition and the fact that r
is only used in (2.15) and (2.16). It follows that we can redefine r + by means of (6.1) in
either the normal equations approach (3.8)/(3.13) or in the augmented system approach
(3.12)/(3.15) and the results would, in theory, be unaffected.
Having this freedom to redefine r + , we seek the value of y that minimizes
where G is any symmetric matrix for which Z T GZ is positive definite, and G \Gamma1 is the
generalized inverse of G. The vector y that solves (6.2) is obtained as
This gives rise to the following modification of the CG iteration.
Algorithm III Preconditioned CG with Residual Update.
Choose an initial point x satisfying (1.2), compute find the
vector y that minimizes kr\GammaA T y, compute
and set \Gammag. Repeat the following steps, until a convergence test is satisfied:
x
This procedure works well in practice, and can be improved by adding iterative refinement
of the projection operation. In this case, at most 1 or 2 iterative refinement steps
should be used. Notice that there is a simple interpretation of Steps (6.6) and (6.7). We
first obtain y by solving (6.2), and as we have indicated the required value is
(3.15). But (3.15) may be rewritten as
G A T
!/
and thus when we obtain g + in Step (6.7), it is as if we had instead found it by solving (6.11).
The advantage of using (6.11) compared to (3.15) is that the solution in the latter may be
dominated by the large components v + , while in the former g + are the large componentsof
course, in floating point arithmetic, the zero component in the solution of (6.11) will
instead be tiny rounded values provided (6.11) is solved in a stable fashion. Viewed in this
way, we see that Steps (6.6) and (6.7) are actually a limited form of iterative refinement in
which the computed v + , but not the computed g + which is discarded, is used to refine the
solution. This "iterative semi-refinement" has been used in other contexts [6, 22].
There is another interesting interpretation of the reset r / r \Gamma A T y performed at the
start of Algorithm III. In the parlance of optimization, c is the gradient of
the objective function (1.1) and r \Gamma A T y is the gradient of the Lagrangian for the problem
(1.1)-(1.2). The vector y computed from (6.2) is called the least squares Lagrange multiplier
estimate. (It is common, but not always the case, for optimization algorithms to set
in (6.2) to compute these multipliers.) Thus in Algorithm III we propose that the initial
residual be set to the current value of the gradient of the Lagrangian, as opposed to the
gradient of the objective function.
One could ask whether it is sufficient to do this resetting of r at the beginning of
Algorithm III, and omit step (6.6) in subsequent iterations. Our computational experience
shows that, even though this initial resetting of r reduces its magnitude sufficiently to avoid
errors in the first few CG iteration, subsequent values of r can grow, and rounding errors
may reappear. The strategy proposed in Algorithm III is safe in that it ensures that r
is small at every iteration, but one can think of various alternatives. One of them is to
monitor the norm of r and only apply the residual update when it seems to be growing.
6.1. The Case
There is a particularly efficient implementation of the residual update strategy when
I. Note that (6.2) is precisely the objective of the least squares problem (3.11) that
occurs when computing via the normal equations approach, and therefore the desired
value of y is nothing other than the vector v + in (3.10) or (3.12). Furthermore, the first
block of equations in (3.12) shows that r . Therefore, in this case (6.6) can
be replaced by r and (6.7) is In other words we have applied the
projection operation twice, and this is a special case of the multiple projections approach
described in the previous section.
Based on these observations we propose the following variation of Algorithm III that
requires only one projection per iteration. We have noted that (6.6) can be written as
. Rather than performing this projection, we will define r where g is the
projected residual computed at the previous iteration. The resulting iteration is given by
Algorithm III with the following two changes:
Omit
Replace (6.10) by g / g + and r /
This strategy has performed well in our numerical experiments and avoids the extra
storage and computation required by Algorithm III. We now show that it is mathematically
equivalent to Algorithm III - which in turn is mathematically equivalent to Algorithm II.
The arguments that follow make use of the fact that, when we have that
The first iteration is clearly the same as that of Algorithm III, except that the value
we store in r in the last step is not r us consider the effect that
this has on the next iteration. The numerator in the definition (6.3) of ff now becomes
T g which equals r T P g. Thus the formula of ff is theoretically unchanged, but the
has the advantage that it can never be negative, as is the
case with (6.3) when rounding errors dominate the projection operation. Next, the step
which is different from the value calculated in Algorithm III.
Step (6.6) is omitted in the new variant of Algorithm III. The projected residual calculated
in (6.7) is now P (P r +ffHp) which is mathematically equivalent to the value PP (r +ffHp)
calculated in Algorithm III (recall that (6.6) can be written as that the new
strategy applies the double projection only to r. Finally let us consider the numerator in
(6.8). In the new variant, it is given by
whereas in Algorithm III it is given by
By expanding these expressions we see that the formula for fi is mathematically equivalent
in both cases, but that in the new variant the projection is applied selectively.
Example 4.
We solved the problem given in Example 1 using this residual update strategy with
I. The results are given in Figure 4 and show that the normal equations and augmented
system approaches are equally effective in this case. We do not plot the cosine (3.17) of
the angle between the preconditioned residual and the columns of A because it was very
small in both approaches, and did not tend to grow as the iteration progressed. For the
normal equations approach this cosine was of order 10 \Gamma14 throughout the CG iteration; for
the augmented system approach it was of order 10 \Gamma15 . Note that we have obtained higher
accuracy than with the iterative refinement strategies described in the previous section;
compare with Figures 2 and 3.
Augmented System
Iteration
residual
Normal Equations
Iteration
residual

Figure

4: Conjugate gradient method with the residual update strategy.
To obtain a highly reliable algorithm for the case when I we can combine the
residual update strategy just described with iterative refinement of the projection operation.
This gives rise to the following iteration which will be used in the numerical tests reported
in x7.
Algorithm IV Residual Update and Iterative Refinement for
Choose an initial point x satisfying (1.2), compute
where the projection is computed by the normal equations (3.8) or augmented
system (3.12) approaches, and set \Gammag. Choose a tolerance ' max . Repeat
the following steps, until a convergence test is satisfied:
x
Apply iterative refinement to P r
until (3.17) is less than '
We conclude this discussion by elaborating on the point made before Example 4 concerning
the computation of the steplength parameter ff. We have noted that the formula
is preferable to (6.12) since it cannot give rise to cancellation. Similarly the
stopping test should be based on g T g rather than on g T r. The residual update implemented
in Algorithm IV does this change automatically, but we believe that these expressions are to
be recommended in other implementations of the CG iteration, provided the preconditioner
is based on
To test this, we repeated the computation reported in Example I using the augmented
system approach; see Figure 1. The only change is that Algorithm II now used the new
for ff and for the stopping test. The CG iteration was now able to continue
past iteration 70 and was able to reach the value
We also repeated the
calculation made in Example 3. Now the residual reached the level
and the
large oscillations in the residual mentioned in Example 3 no longer took place. Thus in
both cases these alternative expressions for ff and for the stopping test were beneficial.
6.2. General G
We can also improve upon the efficiency of Algorithm III for general G, using slightly
outdated information. The idea is simply to use the obtained when computing g + in
(6.7) as a suitable y rather than waiting until after the following step (6.5) to obtain a
slightly more up-to-date version. The resulting iteration is given by Algorithm III, with
the following two changes:
Omit
Replace (6.10) by g / g + and r / r obtained as a
bi-product from (6.7).
Notice, however, that for general G, the extra matrix-vector product A T v + will be required,
since we no longer have the relationship that we exploited when
Although we have not experimented on this idea here, it has proved to be beneficial in
other, similar circumstances [22].
7. Numerical Results
We now test the efficacy of the techniques proposed in this paper on a collection of
quadratic programs of the form (1.1)-(1.2). The problems were generated during the last
iteration of the interior point method for nonlinear programming described in [7], when this
method was applied to a set of test problems from the CUTE [4] collection. We apply the
CG method with preconditioner (3.4) (i.e. with to solve these quadratic programs.
We use the augmented system and normal equations approaches to compute projections,
and for each we compare the standard CG iteration (stand) with the iterative refinement
(ir) techniques described in x5 and the residual update strategy combined with iterative
refinement (update) as given in Algorithm IV. The results are given in Table 1. The first
column gives the problem name, and the second, the dimension of the quadratic program.
To test the reliability of the techniques proposed in this paper we used a very demanding
stopping test: the CG iteration was terminated when
In these experiments we included several other stopping tests in the CG iteration, that
are typically used by trust region methods for optimization. We terminate if the number
of iterations exceeds 2(n \Gamma m) where denotes the dimension of the reduced system
(2.4); a superscript 1 in Table 1 indicates that this limit was reached. The CG iteration
was also stopped if the length of the solution vector is greater than a "trust region radius"
that is set by the optimization method (see [7]). We us a superscript 2 to indicate that this
safeguard was activated, and note that in these problems only excessive rounding errors
can trigger it. Finally we terminate if p T Hp ! 0, indicated by 3 or if r T g ! 0, indicated by
4 . Note that the standard CG iteration was not able to meet the stopping test for any of
the problems in Table 1, but that iterative refinement and update residual were successful
in most cases.

Table

2 reports the CPU time for the problems in Table 1. Note that the times for the
standard CG approach (stand) should be interpreted with caution, since in some of these
problems it terminated prematurely. We include the times for this standard CG iteration
only to show that the iterative refinement and residual update strategies do not greatly
increase the cost of the CG iteration.
Next we report on 3 problems for which the stopping test
could not be
met by any of the variants. For these three problems, Table 3 provides the least residual
norm attained for each strategy.
As a final, but indirect test of the techniques proposed in this paper, we report the
results obtained with the interior point nonlinear optimization code described in [7] on 29
nonlinear programming problems from the CUTE collection. This code applies the CG
method to solve a quadratic program at each iteration. We used the augmented system
Augmented System Normal Equations
Problem dim stand ir update stand ir update
CORKSCRW 147
COSHFUN
OPTCTRL6

Table

1: Number of CG iterations for the different approaches. A 1 indicates that the
iteration limit was reached, 2 indicates termination from trust region bound, 3 indicates
negative curvature was detected and 4 indicates that r T
Augmented System Normal Equations
Problem dim stand ir update stand ir update
COSHFUN
OPTCTRL6

Table

2: CPU time in seconds. 1 indicates that the iteration limit was reached, 2 indicates
termination from trust region bound, 3 indicates negative curvature was detected and 4
indicated that r T
Augmented System Normal Equations
Problem dim stand ir update stand ir update
OBSTCLAE 900 2.3D-07 1.5D-07 5.5D-08 2.3D-07 9.9D-08 4.2D-08

Table

3: The least residual norm:
r T g attained by each option.
and normal equations approaches to compute projections, and for each of these strategies
we tried the standard CG iteration (stand) and the residual update strategy (update) with
iterative refinement described in Algorithm IV. The results are given in Table 4, where
"fevals" denotes the total number of evaluations of the objective function of the nonlinear
problem, and "projections" represents the total number of times that a projection operation
was performed during the optimization. A * indicates that the optimization algorithm
was unable to locate the solution.
Note that the total number of function evaluations is roughly the same for all strategies,
but there are a few cases where the differences in the CG iteration cause the algorithm to
follow a different path to the solution. This is to be expected when solving nonlinear
problems. Note that for the augmented system approach, the residual update strategy
changes the number of projections significantly only in a few problems, but when it does
the improvements are very substantial. On the other hand, we observe that for the normal
equations approach (which is more sensitive to the condition number -(A)) the residual
update strategy gives a substantial reduction in the number of projections in about half
of the problems. It is interesting that with the residual update, the performance of the
augmented system and normal equations approaches is very similar.
8. Conclusions
We have studied the properties of the projected CG method for solving quadratic programming
problems of the form (1.1)-(1.2). Due to the form of the preconditioners used
by some nonlinear programming algorithms we opted for not computing a basis Z for the
null space of the constraints, but instead projecting the CG iterates using a normal equations
or augmented system approach. We have given examples showing that in either case
significant roundoff errors can occur, and have presented an explanation for this.
We proposed several remedies. One is to use iterative refinement of the augmented
system or normal equations approaches. An alternative is to update the residual at every
iteration of the CG iteration, as described in x6. The latter can be implemented particularly
efficiently when the preconditioner is given by I in (3.3).
Our numerical experience indicates that updating the residual almost always suffices
to keep the errors to a tolerable level. Iterative refinement techniques are not as effective
by themselves as the update of the residual, but can be used in conjunction with it, and
the numerical results reported in this paper indicate that this combined strategy is both
economical and accurate.
9.

Acknowledgements

The authors would like to thank Andy Conn and Philippe Toint for their helpful input
during the early stages of this research.
Augmented System Normal Equations
f evals projections f evals projections
Problem n m stand update stand update stand update stand update
CORKSCRW 456 350 64 61 458 422
COSHFUN
GAUSSELM 14 11 25 26 92 93 28 41 85 97
HAGER4 2001 1000
OBSTCLAE 1024 0 26 26 6233 6068 26 26 6236 6080
OPTCNTRL
OPTCTRL6 122

Table

4: Number of function evaluations and projections required by the optimization
method for the different implementations of the CG iteration.



--R

Iterative solution methods.


CUTE: Constrained and unconstrained testing environment.
Some stable methods for calculating inertia and solving symmetric linear equations.
Linear least squares solutions by Housholder transfor- mations
Primal and primal-dual methods for nonlinear programming
Linearly constrained optimization and projected preconditioned conjugate gradients.
The null space problem I: Complexity.
The null space problem II: Algorithms.
A preconditioned conjugate gradient approach to linear equality constrained minimization.
A global convergence theory for general trust-region based algorithms for equality constrained optimization
Direct methods for sparse matrices.
The multifrontal solution of indefinite sparse symmetric linear equations.
The design of MA48

Practical Methods of Optimization.
Computing a sparse basis for the null-space
SNOPT: an SQP algorithm for large-scale constrained optimization
Practical Optimization.
Matrix Computations.
Iterative methods for ill-conditioned linear systems from optimiza- tion
Solving the trust-region subproblem using the Lanczos method
Sparse orthogonal schemes for structural optimization using the force method.
Methods of conjugate gradients for solving linear systems.
Iterative refinement and LAPACK.

Implicit nullspace iterative methods for constrained least squares problems.
On the implementation of an algorithm for large-scale equality constrained optimization
Multifrontal computation with the orthogonal factors of sparse matrices.
Indefinitely preconditioned inexact newton method for large sparse equality constrained nonlinear programming problems.
Preconditioning reduced matrices.
Substructuring methods for computing the null space of equilibrium matrices.
The conjugate gradient method in extremal problems.
QR Factorization of Large Sparse Overdetermined and Square Matrices with the Multifrontal Method in a Multiprocessing Environment.
The conjugate gradient method and trust regions in large scale optimiza- tion
Nested dissection for sparse nullspace bases.
Towards an efficient sparsity exploiting Newton method for minimization.
On large-scale nonlinear network optimization
--TR

--CTR
Luca Bergamaschi , Jacek Gondzio , Manolo Venturin , Giovanni Zilli, Inexact constraint preconditioners for linear systems arising in interior point methods, Computational Optimization and Applications, v.36 n.2-3, p.137-147, April     2007
H. S. Dollar , N. I. Gould , W. H. Schilders , A. J. Wathen, Using constraint preconditioners with regularized saddle-point problems, Computational Optimization and Applications, v.36 n.2-3, p.249-270, April     2007
Luca Bergamaschi , Jacek Gondzio , Giovanni Zilli, Preconditioning Indefinite Systems in Interior Point Methods for Optimization, Computational Optimization and Applications, v.28 n.2, p.149-171, July 2004
S. Bocanegra , F. F. Campos , A. R. Oliveira, Using a hybrid preconditioner for solving large-scale linear systems arising from interior point methods, Computational Optimization and Applications, v.36 n.2-3, p.149-164, April     2007
Nicholas I. M. Gould , Dominique Orban , Philippe L. Toint, GALAHAD, a library of thread-safe Fortran 90 packages for large-scale nonlinear optimization, ACM Transactions on Mathematical Software (TOMS), v.29 n.4, p.353-372, December
S. Cafieri , M. D'Apuzzo , V. Simone , D. Serafino, On the iterative solution of KKT systems in potential reduction software for large-scale quadratic problems, Computational Optimization and Applications, v.38 n.1, p.27-45, September 2007
Nicholas I. M. Gould , Philippe L. Toint, An iterative working-set method for large-scale nonconvex quadratic programming, Applied Numerical Mathematics, v.43 n.1-2, p.109-128, October 2002
Meizhong Dai , David P. Schmidt, Adaptive tetrahedral meshing in free-surface flow, Journal of Computational Physics, v.208 n.1, p.228-252, 1 September 2005
Silvia Bonettini , Emanuele Galligani , Valeria Ruggiero, Inner solvers for interior point methods for large scale nonlinear programming, Computational Optimization and Applications, v.37 n.1, p.1-34, May       2007
conjugate gradient method;nonlinear optimization;quadratic programming;iterative refinement;preconditioning