On the Stability of Null-Space Methods for KKT Systems.
This paper considers the numerical stability of null-space methods for Karush--Kuhn--Tucker (KKT) systems, particularly in the context of quadratic programming. The methods we consider are based on the direct elimination of variables, which is attractive for solving large sparse systems. Ill-conditioning in a certain submatrix A in the system is shown to adversely affect the method insofar as it is commonly implemented. In particular, it can cause growth in the residual error of the solution, which would not normally occur if Gaussian elimination or related methods were used. The mechanism of this error growth is studied and is not due to growth in the null-space basis matrix Z, as might have been expected, but to the indeterminacy of this matrix. When LU factors of A are available it is shown that an alternative form of the method is available which avoids this residual error growth. These conclusions are supported by error analysis and Matlab experiments on some extremely ill-conditioned test problems. These indicate that the alternative method is very robust in regard to residual error growth and is unlikely to be significantly inferior to the methods based on an orthogonal basis matrix. The paper concludes with some discussion of what needs to be done when LU factors are not available.
Introduction
A Karush-Kuhn-Tucker (KKT) system is a linear system
y
version of this paper was presented at the Dundee Biennial Conference in Numerical
Analysis, June, 1995 and the Manchester IMA Conference on Linear Algebra, July 1995.
R. Fletcher and T. Johnson
involving a symmetric matrix of the form
Such systems are characteristic of the optimization problem
subject to A T
in which there are linear equality constraints, and the objective is a quadratic func-
tion. The KKT system (1.1) represents the first order necessary conditions for a locally
minimum solution of this problem, and y is a vector of Lagrange multipliers (see [3]
for example). Problems like (1.3) arise in many fields of study, such as in Newton's
method for nonlinear programming, and in the solution of partial differential equations
involving incompressible fluid flows, incompressible solids, and the analysis of plates and
shells. Also problems with inequality constraints are often solved by solving a sequence of
equality constrained problems, most particularly in the active set method for quadratic
programming.
In (1.2) and (1.3), G is the symmetric n \Theta n Hessian matrix of the objective function,
A is the n \Theta m Jacobian matrix of the linear constraints, and m  n. We assume that
A has full rank, for otherwise K would be singular. In some applications, A does not
immediately have full rank, but can readily be reduced to a full rank matrix by a suitable
transformation.
There are various ways of solving KKT systems, most of which can be regarded as
symmetry-preserving variants of Gaussian elimination with pivoting (see for example
Forsgren and Murray [4]). This approach is suitable for a one-off solution of a large
sparse KKT system, by incorporating a suitable data structure which permits fill-in in
the resulting factors. Our interest in KKT systems arises in a Quadratic Programming
(QP) context, where we are using the so-called null-space method to solve the sequence
of equality constrained problems that arise. This method is described in Section 2. An
important feature of QP is that the successive matrices K differ only in that one column
is either added to or removed from A. The null-space method allows this feature to be
used advantageously to update factors of the reduced Hessian matrix that arises when
solving the KKT system. However in this paper we do not consider the updaing issue, but
concentrate on the solution of a single problem like (1.3), but in a null-space context. In
fact the null-space method is related to one of the above mentioned variants of Gaussian
Elimination, and this point is discussed towards the end of Section 3.
In this paper we study the numerical stability of the null-space method when the
matrix K is ill-conditioned. This arises either when the matrix A is close to being rank
deficient or when the reduced Hessian matrix is ill-conditioned. It is well known however
that Gaussian elimination with pivoting usually enables ill-conditioned systems to be
solved with small backward error (that is the computed solution is the exact solution of
Stability of Null-Space Methods 3
a nearby problem). As Wilkinson [6] points out, the size of the backward error depends
only on the growth in certain reduced matrices, and the amount of growth is usually
negligible for an ill-conditioned matrix. Although it is possible for exponential growth
to occur (we give an example for a KKT system), this is most unlikely in practice. A
consequence of this is that if the computed solution is substituted into the system of
equations, a very accurate residual is obtained. Thus variants of Gaussian elimination
with pivoting usually provide a very stable method for solving ill-conditioned systems.
However this argument does not carry over to the null-space method and we indicate
at the end of Section 2 that there are serious concerns about numerical stability when A is
nearly rank deficient. We describe some Matlab experiments in Section 6 which support
these concerns. In particular the residual of the KKT system is seen to be proportional
to the condition number of A. We present some error analysis in Section 4 which shows
how this arises.
When LU factors of A are available, we show in Section 3 that there is an alternative
way of implementing a null-space method, which avoids the numerical instability. This is
also supported by Matlab experiments. The reasons for this are described in Section 5,
and we present some error analysis which illustrates the difference in the two approaches.
In practice, when solving large sparse QP problems, LU factors are not usually available
and it is more usual to use some sort of product form method. We conclude with some
remarks about what can be done in this situation to avoid numerical instability.
Null-Space Methods
A null-space method (see e.g. [3]) is an important technique for solving quadratic programming
problems with equality constraints. In this section we show how the method
can be derived as a generalised form of constraint elimination. The key issue in this
procedure is the formation of a basis for the null space of A. We determine the basis
in such a way that we are able to solve large sparse problems efficiently. When A is
ill-conditioned we argue that there is serious concern for the numerical stability of the
method.
The null space of A may be defined by
and has dimension when A has full rank. Any matrix
whose columns are a basis for N (A) will be referred to as a null-space matrix for A. Such
a matrix satisfies A T and has linearly independent columns. A general specification
for computing a null-space matrix is to choose an n \Theta (n \Gamma m) matrix V such that the
matrix
4 R. Fletcher and T. Johnson
is non-singular. Its inverse is then partitioned in the following way
n\Gammam
It follows from the properties of the inverse that A T . By construc-
tion, the columns of Z are linearly independent, and it follows that these columns form
a basis for N (A).
The value of this construction is that it enables us to parametrize the solution set of
the (usually) underdetermined system A T in (1.3) by
Here Y b is one particular solution of A T any other solution x differs from
Y b by a vector, Zv say, in N (A). Thus (2.2) provides a general way of eliminating the
constraints, by expressing the problem in terms of the reduced variables v. Hence if (2.2)
is substituted into the objective function of (1.3), we obtain the reduced problem
We refer to Z T GZ as the reduced Hessian matrix and Z T (GY b\Gammac) as the reduced gradient
vector (at the point sufficient condition for (2.3) to have a unique minimizer
is that Z T GZ is positive definite. In this case there exist Choleski factors Z T
and (2.3) can be solved by finding a stationary point, that is by solving the linear system
Then substitution of v into (2.2) determines the solution x of (1.3). The vector Gx \Gamma c is
the gradient of the objective function at the solution, so a vector y of Lagrange multipliers
satisfying can then be obtained from
by virtue of the property that Y T I. The vectors x and y also provide the solution
to (1.1) as can readily be verified.
In practice, when A is a large sparse matrix, the matrices Y and Z are usually
substantially dense and it is impracticable to store them explicitly. Instead, products
with Y and Z or their transposes are obtained by solving linear systems involving A.
For example the vector could be computed by solving the linear
system
A
by virtue of (2.1). Likewise solving the system
A
Stability of Null-Space Methods 5
provides the products u 1
t. These computations require an invertible
representation of the matrix A to be available.
Solving systems involving A is usually a major cost with the null-space method. To
keep this cost as low as possible, it is preferable to choose the matrix V to be sparse.
Other choices (for example based on the QR factors of A, see [3]) usually involve significantly
more fill-in and computational expense. In particular it is attractive to choose the
columns of V to be unit vectors, using some form of pivoting to keep A as well conditioned
as possible. In this case, assuming that the row permutation has been incorporated into
A, it is possible to write
"I
where A 1
is an m \Theta m nonsingular submatrix. Then (2.1) becomes
and provides an explicit expression for Y and Z. In particular we see that
We refer to this choice of V as direct elimination as it corresponds to directly using the
first m variables to eliminate the constraints (see [3]). We shall adopt this choice of V
throughout the rest of the paper.
The reduced Hessian matrix Z T GZ is also needed for use in (2.3), and can be calculated
in a similar way. The method is to compute the vectors Z T GZe k for
of the unit matrix I n\Gammam . The computation
is carried out from right to left by first computing the vector z by solving the
system
A T z
Then the product Gz k is computed, followed by the solution of
The partition u 2
is then column k of Z T GZ as required. The lower triangle of Z T GZ
is then used to calculate the Choleski factor L. A similar approach is essentially used
in an active set method for QP, in which the Choleski factor of Z T GZ is built up over
a sequence of iterations. (If indefinite QP problems are solved, it may be required to
solve KKT systems in which Z T GZ is indefinite. We note that such systems can also be
solved in a numerically stable way which preserves symmetry, see Higham [5] in regard
to the method of Bunch and Kaufmann [1]).
6 R. Fletcher and T. Johnson
An advantage of the null-space approach is that we only need to have available a
subroutine for the matrix product Gv. Thus we can take full advantage of sparsity or
structure in G, without for example having to allow for fill-in as Gaussian elimination
would require. The approach is most convenient when Z T GZ is sufficiently small to
allow it to be stored as a dense matrix. In fact there is a close relationship between
the null-space method and a variant of Gaussian elimination, as we shall see in the next
section, and the matrix Z T GZ is the same submatrix in both methods. Thus it would
be equally easy (or difficult) to represent Z T GZ in a sparse matrix format with either
method.
To summarize the content of this section we can enumerate the steps implied by (2.2)
through (2.5)
1. Calculate Z T GZ as in (2.10) and (2.11).
2. Calculate by a solve with A T as in (2.6) with
3. Calculate requiring a product with G.
4. Calculate u 2
by a solve with A as in (2.7).
5. Solve Z T
to determine v as in (2.4).
6. Calculate by a solve with A T as in (2.6).
7. Calculate requiring a product with G.
8. Calculate by a solve with A, which also provides z g.
When direct elimination based on (2.9) is used, we shall refer to this as Method 1. Step 1
requires m) solves with either A or A T and products with G to set up the
reduced Hessian matrix. The remaining steps require 4 solves and 2 products, plus a
solve with Z T GZ. In some circumstances these counts can be reduced. If
steps 2 and 3 are not required. If the multiplier part y of the solution is not of interest
then steps 7 and 8 are not needed.
We now turn to the concerns about the numerical stability of the null-space method
when A (and hence A 1
and is ill-conditioned. In this case A is close to a rank deficient
say, which has a null space of higher dimension. When we solve systems
like (2.10) and (2.11), the matrix Z that we are implicitly using is badly determined.
Therefore, because of round-off error, we effectively get a significantly different Z matrix
each time we carry out a solve. Thus the computed reduced Hessian matrix Z T GZ does
not correspond to any one particular Z matrix. As we shall see in the rest of the paper,
this can lead to solutions with significant residual error.
Stability of Null-Space Methods 7
3 Using LU factors of A
In this section we consider the possibility that we can readily compute LU factors of A
given by
is unit lower triangular and U is upper triangular. We can assume that a row
permutation has been made which enables us to bound the elements of L 1
and L 2
by
1. As we shall see, these factors permit us to circumvent the difficulties caused by
ill-conditioning to a large extent. (Unfortunately, LU factors are not always available,
and some indication is given in Section 7 as to what might be done in this situation.)
We also describe how the steps in the null-space method are changed. Finally we explore
some connections with Gaussian elimination and other methods, which provide some
insight into the likelihood of growth in Z.
A key observation is that if LU factors of A are available, then it is possible to express
Z in the alternative form
in which the UU \Gamma1 factors arising from (2.9) and (3.1) are cancelled out. A minor
disadvantage, compared to (2.9), is that L 2
is needed, which is likely to be less sparse
than A 2
and also requires additional storage. However if A is ill-conditioned, this is
manifested in U (but not usually L) being ill-conditioned, so that (3.2) enables Z to be
defined in a way which is well-conditioned. In calculating the reduced Hessian matrix it
is convenient to define
I
and replace equations (2.10) and (2.11) by
and
I
The steps of the resulting null-space method are as follows (using subscript 1 to denote
the first m rows of a vector or matrix, and subscript 2 to denote the last
1. Calculate Z T GZ as in (3.4) and (3.5).
2. Calculate s 1
8 R. Fletcher and T. Johnson
3. Calculate requiring a product with G.
4. Calculate u 2
5. Solve Z T
for v.
6. Calculate
7. Calculate
8. Calculate requiring a product with G.
9. Calculate
10. Calculate
In the above, inverse operations involving L 1
and U are done by forward or backward
substitution. The method is referred to as Method 2 in what follows. (For comparability
with Method 1, we have also included the calculation of the reduced gradient z, although
this would not normally be required.) Note that all solves with the n \Theta n matrix A
are replaced by solves with smaller m \Theta m matrices. Also steps 1, 4, 6 and 10 use the
alternative definition (3.2) of Z and so avoid a potentially ill-conditioned calculation with
A (or A 1
We consider the numerical stability of both Method 1 and Method 2 in more
detail in the next section.
In the rest of this section, we explore some connections between this method and some
variants of Gaussian elimination, and we examine the factored forms that are provided by
these methods. It is readily observed (but not well known) that there are block factors of
K corresponding to any null-space method in this general format. These are the factors
I
A T
(using blanks to denote a zero matrix). This result is readily verified by using the
equation I derived from (2.1). This expression makes it clear that
inverse representations of the matrices A and Z T GZ will be required to solve (1.1).
However these factors are not directly useful as a method of solution as they also involve
the matrices Y T GY and Y T GZ whose computation we wish to avoid in a null-space
method. Equation (3.6) also shows that K \Gamma1 will become large when either A or Z T GZ
is ill-conditioned, and we would expect the spectral condition number to behave like
A M where
When using direct elimination (2.8) we may partition K in the form
G 11
G 21
G 22
A
Stability of Null-Space Methods 9
When A has LU factors (3.1) then it is readily verified that another way of factorizing
K is given by6 4
G 11
G 21
G 22
A
I
Z U
U T7 56 4
where Z is defined by (3.2) and G 1
]. Note that the matrix U occurs on
the reverse diagonal of the middle factor, but that no operations with U \Gamma1 are required
in the calculation of the factors. Thus any ill-conditioning associated with U does not
manifest itself until the factors are used in solving the KKT system (1.1). If there is no
growth in Z then the backward error in (3.7) will be small, indicating the potential for
a small residual solution of the KKT system. We show in Section 5 how this can come
about. Another related observation is that if A is rank deficient, then the factors (3.6)
do not exist (since the calculation of Y involves A \Gamma1and hence U
be calculated without difficulty.
The factorization (3.7) of K is closely related to some symmetry preserving variants
of Gaussian elimination. Let us start by eliminating A 2
and the sub-diagonal elements
of A 1
by row operations. (As before we can assume that row pivoting has been used.)
The outcome of these row operations is that6 4
G 21
G 22
A
I
G 22
]. Note that these row operations are exactly those used by Gaussian
elimination to form (3.1). To restore symmetry in the factors, we repeat the above
procedure in transposed form, that is we make column operations on A Tand A T, which
gives rise to (3.7).
We can also interleave these row and column operations without affecting the final
result. If we pair up the first row and column operation, then the second row and
column operation, and so on, then we get the method of 'ba' pivots described by Forsgren
and Murray [4]. Thus these methods essentially share the same matrix factors. The
difference is that in the null-space method, Z T GZ is calculated by matrix solves with
A, as described in Section 2, whereas in these other methods it is obtained by row and
column operations on the matrix K.
This association with Gaussian elimination enables us to bound the growth in the
R. Fletcher and T. Johnson
factors of K. The bound is attained for the critical case typified by the matrix
for which Row operations with pivots in the (1,7), (2,8), (3,9) and
positions leads to the
Then column operations with pivots in the (7,1), (8,2), (9,3) and (10,4) positions gives
rise to 2
which corresponds to the middle factor in (3.7). In this case and the
corresponding matrix Z is given by
Stability of Null-Space Methods 11
In general it is readily shown that when m ! n, growth of 2 2m in the maximum modulus
element of K can occur. For the null-space method based on (3.2), this example also
illustrates the maximumpossible growth of 2  1. In practice however
such growth is most unlikely and it is usual not to get any significant growth in Z.
4 Numerical Stability of Method 1
In this and the next section we consider the effect of ill-conditioning in the matrix K on
the solutions obtained by null-space methods based on direct elimination. In particular
we are interested to see whether or not we can establish results comparable to those for
Gaussian elimination. We shall show that the forward error in x is not as severe as would
be predicted by the condition number of K. We also look at the residual errors in the
solution and show that Method 2 is very satisfactory in this respect, whereas Method 1
is not.
In order to prevent the details of the analysis from obscuring the insight that we
are trying to provide, we shall adopt the following simple convention. We imagine that
we are solving a sequence of problems in which either A or M (the spectral condition
numbers of A and increasing without bound. We then use the notation
O(h) to indicate a quantity that is bounded in norm by ckhk on this sequence, where
there exists an implied constant c that is independent of A or M , but may contain a
modest dependence on n. Also we shall assume that the system is well scaled so that
1. This enables us for example to deduce that multiplication of an
error bound O(") by A \Gamma1 causes the bound to be increased to O(A "). We also choose
to assume that the KKT system models a situation in which the exact solution vectors
x and y exist and are not unreasonably large in norm, that is
A similar assumption is needed in order to show that Gaussian elimination provides
accurate residuals, so we cannot expect to dispense with this assumption. Sometimes it
may be possible to argue that we are solving a physical problem which is known to have
a well behaved solution.
Another assumption that we make is that the choice of the matrix V in (2.8) (and
hence the partitioning of A) is made using some form of pivoting. Now the exact solution
for Z is given by
from (3.3), using the factors of A defined in (3.1). It follows that
where L is the spectral condition number of L. Assuming that partial pivoting is used,
so that jl ij j  1, and that negligible growth occurs in L \Gamma1, it then follows that negligible
growth occurs in Z and we can assert that
12 R. Fletcher and T. Johnson
Another consequence of this assumption is that we are able to neglect terms of O(L ")
relative to terms of O(A ") when assessing the propagation of errors for Method 2.
We shall now sketch some properties (Wilkinson [6]) of floating point arithmetic of
relative precision ". If a nonsingular system of n linear equations solved by
Gaussian elimination, the computed solution b
x is the exact solution of a perturbed system
referred to as the backward error. E can be bounded by an
expression of the form aeOE(n)" in which ae measures the growth in A during the elimination
and OE(n) is a modest quadratic in n. For ill-conditioned systems, and assuming that
partial pivoting is used, growth is rare and can be ignored. Also this bound usually
overstates the dependence on n which is unlikely to be a dominant factor. Hence for the
backward error
We can measure the accuracy of the solution either by the forward error b
x
or by computing the residual
x. Using
where A is some condition number of A. Since assuming that A "  1,
it follows that b
Likewise we can deduce that
These bounds are likely to be realistic and tell us that for Gaussian elimination, ill-conditioning
affects the forward error in x but not the residual r, as long as b
x is of
reasonable magnitude.
Wilkinson also gives expressions for the backward error in a scalar product and hence
in the product +Ax. The computed product b s is the exact product of a system in
which the relative perturbation in each element of b and A is no more than n" where n
is the dimension of x. We can express this as
if we make the assumption that b and A are O(1).
The first stage in a null-space calculation is the determination of Z T GZ, which we
denote by M . In Method 1, this is computed as in (2.10) and (2.11). In (2.10) a column
z k of the matrix Z is computed which, by applying (4.4), satisfies
Stability of Null-Space Methods 13
where A is the spectral condition number of A. The product with G introduces negligible
error, and the solution of (2.11) together with (4.5) shows that
Ab
Multiplying by L \Gamma1 and extracting the b
partition gives
using (4.7) and then (4.2). Hence we have established that
c
The argument has been given in some detail as it is important to see why the error in
M is O(A ") and not O( 2
A "). We also observe that hence that
c
We now turn to the solution of the KKT system using Method 1. We shall assume
that systems involving A and M are solved in such a way that (4.5) applies. Using (4.6),
and assuming that the computed quantities b s; b t; are O(1), the residual errors in
the sequence of calculations are then
Ab
c
A
y
z
These results, together with (4.8), may be combined to get the forward errors in the solution
vectors b
x and b
y. Multiplying through equations (4.9) and (4.13) by A \GammaT magnifies
the errors by a factor A (since we are assuming that A = O(1)), giving
We can get a rather better bound from (4.11) and (4.15) by first multiplying through by
using to give
14 R. Fletcher and T. Johnson
from the second partition of the solution. However the first partition of (4.15) gives
Combining (4.8) and (4.12) gives
We can now chain through the forward errors, noting that a product with Z or Z T does
not magnify the error in a previously computed quantity (by virtue of (4.2). However
the product M \Gamma1 b
in (4.21) magnifies the error in b
by a factor M and the product
in (4.20) magnifies the error in b g by a factor A . The outcome is that
and
A M "): (4.23)
As we would expect, the forward errors are affected by the condition numbers of A and
M . However although the condition number of K is expected to be of the order  2
A M ,
we see that this factor only magnifies the error in the y part of the solution, with the x
part being less badly affected.
When K is ill-conditioned we must necessarily expect that the forward errors are
adversely affected. A more important question is to ask whether the solution satisfies the
equations of the problem accurately. There are three measures of interest, the residuals
of the KKT system (1.1), and the reduced gradient
Gx is the negative gradient vector at the solution. We note that
the vector z is computed as a by-product of step 8 of Method 1.
If we compute r we obtain b
as in (4.6), and it follows from (4.13)
and the definition of A that A T b
When computing q we obtain
b z
from (4.14) and (4.15). Thus the accuracy of b
q depends on that of b z. From (4.19) and
it follows that
Stability of Null-Space Methods 15
from (4.17). (Notice that it is important not to use (4.22) here which would give an
unnecessary factor of M .) Then (4.8), (4.12), (4.11) and (4.10) can be used, giving
Thus we are able to predict under our assumptions that the reduced gradient b z and the
residual b q are adversely affected by ill-conditioning in A, but not by ill-conditioning in
M . However the residual b r is unaffected by ill-conditioning either in A or M .
Simulations are described in Section 6 which indicate that these error bounds reliably
predict the actual effects of ill-conditioning. Method 1 is seen to be unsatisfactory in
that an accurate residual q cannot be obtained when A is ill-conditioned. We shall show
in the next section that Method 2 does not share this disadvantage.
The main results of this section and the next are summarised and discussed in Section
7.
5 Numerical Stability of Method 2
In this section we assess the behaviour of Method 2 in the presence of ill-conditioning
in K. Although we cannot expect any improvement for the forward errors, we are able
to show that Method 2 is able to give accurate residuals that are not affected by ill-
conditioning. The relationship between Method 2 and Gaussian elimination described
towards the end of Section 3 gives some hope of proving this result. However this is not
immediate because Method 2 does not make direct use of the factors (3.7) in the same
way that Gaussian elimination does.
A fundamental difficulty with the analysis of Method 2 is that we can deduce from
(4.7) that
and this result cannot be improved if LU factors are available. To see this, we know that
the computed factors of any square matrix A satisfy
when there is no growth in b
U . If are the exact factors, it follows that
U
say, where Q is the strict lower triangular part of L
U \Gamma1 and R is the upper triangular
part. Because L
L is unit lower triangular and U b
U \Gamma1 is upper triangular we can deduce
that
involves an inverse operation
with b
U we can expect that b
L and L differ by O(A "). This result has been confirmed
by computing the LU factors of a Hilbert matrix in single and double precision Fortran.
On applying the result to our matrix A, it follows that (5.1) holds.
R. Fletcher and T. Johnson
Fortunately all is not lost because we are still able to compute a null-space matrix
which accurately satisfies the equation Z T
Z denote the null-space matrix
obtained from b
L in exact arithmetic. It follows that b
and hence from (5.2) that
We also have b
long as A "  1. Our analysis will express the errors that
arise in Method 2 in terms of b
Z rather than Z and this enables us to avoid the A factor
in the residuals.
The first step in Method 2 is to compute GZ as in (3.4) and (3.5). In this
section we denote c
Z as the value computed from b
Z in exact arithmetic and
use c c
M to denote the computed value of c
M . It readily follows, using results like (4.2),
that
c c
We now consider the solution of the KKT system using Method 2. As in equations
through (4.15) we assume that the computed quantities b s; b t; are O(1). Then
the residual errors in the sequence of calculations are
A Tb s 1
c c
It is readily seen from these equations that the forward errors will propagate in a similar
way to Method 1.
Turning to the residual errors, the computed value of the residual r is
from (5.10), (5.9), (5.5) and (5.3). When computing q we obtain b
for Method 1, and it follows from (5.12) that b q 1
O("). From (5.3) we can deduce that
so it follows that
Stability of Null-Space Methods 17
Thus the accuracy of the residual b
q depends on that of b z, as for Method 1. For b
z we can
use (5.13), (5.11), (5.10) and (5.9) to get
Now we can invoke (5.4) and (5.8) giving
from (5.7) and (5.6). Thus we have established under our assumptions that all three
measures of accuracy for the KKT system are O(") for Method 2 and are not affected by
ill-conditioning in either A or M . These results are again supported by the simulations
in the next section.

Figure

1. Condition numbers of K, A and L
6 Numerical Experiments
In order to check the predictions of Sections 4 and 5, some experiments have been carried
out on artifically generated KKT systems. These experiments have been carried out in
Matlab for which the machine precision is They suggest that the upper bounds
given by the error analysis accurately reflect the actual behaviour of an ill-conditioned
system. Another phenomenon that occurs when the ill-conditioning is very extreme is
also explained.
The KKT systems have been constructed in the following way. To make A ill-conditioned
we have chosen it as the first m columns of the n \Theta n Hilbert matrix, where
provides a sequence of problems for which the condition
number of A increases exponentially. Factors are calculated by the
Matlab routine lu which uses Gaussian Elimination with partial pivoting, and A is replaced
by PA. In the first instance the matrix G is generated by random numbers in
the range [\Gamma1; 1]. However to make positive definite, a multiple of the unit
matrix is added to the G 22
partition of G, chosen so that the smallest eigenvalue of M
R. Fletcher and T. Johnson
is changed to 10 1\Gammak for some positive integer k. The assumptions of the analysis require
that the KKT system has a solution that is O(1). To achieve this, exact solutions x and
y are generated by random numbers in [\Gamma1; 1], and the right hand sides c and b are
calculated from (1.1). For each value of m, 10 runs are made with a different random
number seed and the statistics are averaged over these 10 runs.
First of all we examine the effect of increasing the condition number of A whilst
keeping M well-conditioned. To do this we increase m from 2 up to 10, whilst fixing
1. The resulting condition numbers of K, A and L are plotted in Figure 1. It can
be seen that the slope of the unbroken line (K ) is about twice that of the dashed line
1, this is consistent with the estimate K   2
A M that we deduced
in Section 3. The condition number of L (dotted line) shows negligible increase, showing
that there is no growth in L \Gamma1, thus enabling us to assert that O(1). The levelling
out of the K graph for is due to round-off error corrupting the least
eigenvalue of K.

Figure

2. Error growth vs. A for Method 1 Figure 3. Error growth vs. A for Method 2
The effect of the conditioning of A on the different types of error is illustrated in

Figures

2 and 3. The forward error is shown by the two unbroken lines, the upper line
being the error in y and the lower line being the error in x. The upper line has a slope
of about 2 on the log-log scale, and the lower line has a slope of about 1, and both have
an intercept with the y-axis of about 10 \Gamma16 . This is precisely in accordance with (4.23)
and (4.22). It can also be seen that both methods exhibit the same forward error. The
computed value of the residual error shown by the dashed line and both
methods show the O(") behaviour as predicted by (4.24) and (5.14), with the increasing
condition number having no effect.
The difference between Methods 1 and 2 is shown by the computed values of the
residual (dotted line) and the reduced gradient
line). As we would expect from (4.27), these graphs are superimposed, and they clearly
show the influence of A on the error growth for Method 1, as predicted by (4.28).
Negligible error growth is observed for Method 2 as predicted by (5.16), except for an
Stability of Null-Space Methods 19
increase in q for A greater than about 10 9 . This feature is explained later in the section.

Figure

4. Error growth vs. M for Method 1 Figure 5. Error growth vs. M for Method 2
We now turn to see the influence of ill-conditioning in M on the errors. To do
this we fix carry out a sequence of calculations with
causes M to increase exponentially. Each calculation is the
average of ten runs as above. The results are illustrated in Figures 3 and 4, using the
same key. The forward errors are again seen to be the same for both methods and they
both have a slope of about 1 on the log-log scale, corresponding to the M factor in (4.22)
and (4.23). The upper line for the forward error in y lies about 10 5 units above that for
the forward error in x, as the extra factor of A in (4.23) would predict. The residual r is
seen to be unaffected by the conditioning of M as above. The residual q and the reduced
gradient z are also unaffected by M , but the graphs for Method 1 lie above those for
Method 2, due to the A factor in (4.28). All these effects are in accordance with what
the error analysis predicts.
To examine the anomalous behaviour of q in Figure 3 in more detail, we turn to
a sequence of more ill-conditioned test problems obtained by using the last m columns
of the Hilbert matrix to define A. The results for Method 2 are illustrated in Figure
6 and the anomalous behaviour (dotted line) is now very evident. The reason for this
becomes apparent when it is noticed that it sets in when the forward error in y, and
hence the value of b
y, becomes greater than unity. This possibility has been excluded in
our error analysis by the assumption that b
O(1). The anomalous behaviour sets in
when  2
that is A ' (M ") \Gamma1=2 , or in this case A ' 10 8 , much as Figures 3
and 6 illustrate. For greater values of A there is a term O(b y") in the expression for b
indicating that the error is of the form  2
. The fact that this part of the graph of
q is parallel to the graph of the forward error in y supports this conclusion.
The above calculations have also been carried out using a Vandermonde matrix in
place of the Hilbert matrix and very similar results have been obtained.
R. Fletcher and T. Johnson

Figure

6. Error growth for Method 2 for a more ill-conditioned matrix
7 Summary and Discussion
In this paper we have examined the effect of ill-conditioning on the solution of a KKT
system by null-space methods based on direct elimination. Such methods are important
because they are well suited to take advantage of sparsity in large systems. However they
have often been criticised for a lack of numerical stability, particularly when compared
to methods based on QR factors. We have studied two methods: Method 1 in which
an invertible representation of A in (2.8) is used to solve systems, and Method 2 in
which LU factors (3.1) of A are available. We have presented error analysis backed up by
numerical simulations which, under certain assumptions on growth, provide the following
conclusions
ffl Both methods have the same forward error bounds, with b
A M ").
ffl Both methods give accurate residuals if A is well conditioned, even if M is ill-conditioned

gives an accurate residual
for Method 1.
ffl Both methods give an accurate residual if A is ill-conditioned.
These conclusions do indicate that Method 1 is adversely affected by ill-conditioning in
A, even though the technique for solving systems involving A is able to provide accurate
residuals. The reasons for this are particularly interesting. For example one might expect
that when A is ill-conditioned, then A \Gamma1 would be large and we might therefore expect
from (2.1) that Z would be large. In fact we have seen that as long as V is chosen
suitably, then growth in Z is very unlikely (the argument is similar to that for Gaussian
elimination). Of course if V is badly chosen then Z can be large and this will cause
Stability of Null-Space Methods 21
significant error. One might also expect that because the forward error in computing Z
is necessarily of order O(A "), it would follow that no null-space method could provide
accurate residuals.
The way forward, which is exploited in the analysis for Method 2, is that Method 2
determines a matrix b
Z for which b
O("). Thus, although the null-space is inevitably
badly determined when A is ill-conditioned, Method 2 fixes on one particular basis matrix
Z that is well behaved. This basis is an exact basis for an O(") perturbation to A.
Method 2 is able to solve this perturbed problem accurately. On the other hand Method 1
essentially obtains a different approximation to Z for every solve with A. Thus the
computed reduced Hessian matrix Z T GZ does not correspond accurately to any one
particular b
Z matrix.
In passing, it is interesting to remark that computing the factors
and defining
, also provides a stable approach, not so much because it avoids the
growth in Z (we have seen that this is rarely a problem), but because it also provides a
fixed null-space reference basis, which is an exact basis for an O(") perturbation to A.
In the context of quadratic programming, a common solution method for large sparse
systems is to use some sort of product form method (Gauss-Jordan, Bartels-Golub-Reid,
Forrest-Tomlin etc. It is not clear that such methods provide O(") solutions to the
systems involving A that are solved in Method 1 (although B-G-R may be stable in
this respect). However the main difficulty comes when the product form becomes too
unweildy and is re-inverted. If A is ill-conditioned, the refactorization of A is likely to
determine a basis matrix Z that differs by O(A ") from that defined by the old product
form. Thus the old reduced Hessian matrix Z T GZ would not correspond accurately to
that defined by the new Z matrix after re-inversion. The only recourse would be to
re-evaluate Z T GZ on re-inversion, which might be very expensive. Thus we do not see a
product form method on its own as being suitable. Our paper has shown that if a fixed
reference basis is generated then accurate residuals are possible. It is hoped to show how
this might be done in a subsequent paper by combining a product form method with
another method such as LU factorization.



--R

Parlett B.
Erisman A.
Practical Methods of Optimization
Newton methods for large-scale linear equality- constrained minimization
Stability of the Diagonal Pivoting Method with Partial Pivot- ing
The Algebraic Eigenproblem
--TR
null-space method;KKT system;ill-conditioning