Constraint Preconditioning for Indefinite Linear Systems.
The problem of finding good preconditioners for the numerical solution of indefinite linear systems is considered. Special emphasis is put on preconditioners that have a 2  2 block structure and that incorporate the (1,2) and (2,1) blocks of the original matrix.  Results concerning the spectrum and form of the eigenvectors of the preconditioned matrix and its minimum polynomial are given. The consequences of these results are considered for a variety of Krylov subspace methods. Numerical experiments validate these conclusions.
Introduction
In this paper we are concerned with investigating a new class of preconditioners
for indefinite systems of linear equations of a sort which arise in constrained
optimization as well as in least-squares, saddle-point and Stokes problems. We
attempt to solve the indefinite linear system
A
where A 2 IR n\Thetan is symmetric and B 2 IR m\Thetan . Throughout the paper we shall
assume that m - n and that A is non-singular, in which case B must be of full
rank.
Example 1. Consider the problem of minimizing a function of n variables
subject to m linear equality constraints on the variables, i.e.
minimize
subject to
Any finite solution to (1.2) is a stationary point of the Lagrangian function
where the - i are referred to as Lagrangian multipliers. By differentiating L with
respect to x and - the solution to (1.2) is readily seen to satisfy linear
equations of the form (1.1) with x d. For this
application these are known as the Karush-Kuhn-Tucker (KKT) conditions. 2
Example 2 (The Stokes Problem). The Stokes equations in compact form are
defined by
div
Discretising equations (1.3) together with the boundary conditions
defines a linear system of equations of the form (1.1), where b
Carsten Keller, Nick Gould and Andy Wathen
Among the most important iterative methods currently available, Krylov sub-space
methods apply techniques that involve orthogonal projections onto sub-spaces
of the form
The most common schemes that use this idea are the method of conjugate
gradients (CG) for symmetric positive definite matrices, the method of minimum
residuals (MINRES) for symmetric and possibly indefinite matrices and
the generalised minimum residual method (GMRES) for unsymmetric matrices,
although many other methods are available-see for example Greenbaum [12].
One common feature of the above methods is that the solution of the linear
system (1.1) is found within n +m iterations in exact arithmetic-see Joubert
and Manteuffel [14, p. 152]. For very large (and possibly sparse) linear systems
this upper limit on the number of iterations is often not practical. The idea
of preconditioning attempts to improve on the spectral properties, i.e. the
clustering of the eigenvalues, such that the total number of iterations required
to solve the system to within some tolerance is decreased substantially.
In this paper we are specifically concerned with non-singular preconditioners
of the form
n\Thetan approximates, but is not the same as A. The inclusion of the
exact representation of the (1; 2) and (2; 1) matrix blocks in the preconditioner,
which are often associated with constraints (see Example 1), leads one to hope
for a more favourable distribution of the eigenvalues of the (left-)preconditioned
linear system
Since these blocks are unchanged from the original system, we shall call G a
constraint preconditioner. A preconditioner of the form G has recently been
used by Luk-san and Vl-cek [16] in the context of constrained non-linear programming
problems-see also Coleman [4], Polyak [18] and Gould et al. [11].
Here we derive arguments that confirm and extend some of the results in [16]
and highlight the favourable features of a preconditioner of the form G. Note
that Golub and Wathen [10] recently considered a symmetric preconditioner of
the form (1.4) for problems of the form (1.1) where A is non-symmetric.
In Section 2 we determine the eigensolution distribution of the preconditioned
system and give lower and upper bounds for the eigenvalues of G \Gamma1 A
Constraint Preconditioning 3
in the case when the submatrix G is positive definite. Section 3 describes the
convergence behaviour of a Krylov subspace method such as GMRES, Section 4
investigates possible implementation strategies, while in Section 5 we give numerical
results to support the theory developed in this paper.
Preconditioning A
For symmetric (and in general normal) matrix systems, the convergence of an
applicable iterative method is determined by the distribution of the eigenvalues
of the coefficient matrix. In particular it is desirable that the number of
distinct eigenvalues, or at least the number of clusters, is small, as in this case
convergence will be rapid. To be more precise, if there are only a few distinct
eigenvalues then optimal methods like CG, MINRES or GMRES will terminate
(in exact arithmetic) after a small and precisely defined number of steps.
We prove a result of this type below. For non-normal systems convergence as
opposed to termination is not so readily described-see Greenbaum [12, p. 5].
2.1 Eigenvalue Distribution
The eigenvalues of the preconditioned coefficient matrix G \Gamma1 A may be derived
by considering the general eigenvalue problem
x
y
x
y
Y Z
be an orthogonal factorisation of B T , where
n\Thetam and Z 2 IR n\Theta(n\Gammam) is a basis for
the nullspace of B. Premultiplying (2.6) by the non-singular and square matrix6 6 4
and postmultiplying by its transpose gives6 6 4
x z
x y
x z
x y
with and where we made use of the equalities
Performing a simultaneous sequence of row and column inter-
4 Carsten Keller, Nick Gould and Andy Wathen
changes on both matrices in (2.7) reveals two lower block-triangular matrices
~
and thus the preconditioned coefficient matrix G \Gamma1 A is similar to
I
\Theta (Z T GZ)
Here the precise forms of \Theta, \Upsilon and \Gamma are irrelevant for the argument that
they are in general non-zero. We just proved the following theorem.
Theorem 2.1. Let A 2 IR (n+m)\Theta(n+m) be a symmetric and indefinite matrix
of the form
where A 2 IR n\Thetan is symmetric and B 2 IR m\Thetan is of full rank. Assume Z is
an n \Theta (n \Gamma m) basis for the nullspace of B. Preconditioning A by a matrix
of the form
n\Thetan is symmetric, G 6= A and B 2 IR m\Thetan is as above, implies
that the matrix G \Gamma1 A has
(1) an eigenvalue at 1 with multiplicity 2m; and
eigenvalues which are defined by the generalised eigenvalue
problem Z T AZx
Note that the indefinite constrained preconditioner applied to the indefinite linear
system (1.1) yields the preconditioned matrix P which has real eigenvalues.
Remark 1. In the above argument we assumed that B has full row rank and
consequently applied an orthogonal factorisation of B T which resulted in a
Constraint Preconditioning 5
upper triangular matrix R 2 IR m\Thetam . If B does not have full row rank, i.e.
rows and columns can be
deleted from both matrices in (2.7), thus giving a reduced system of dimension
This removal of the redundant information does not
impose any restriction on the proposed preconditioner, since all mathematical
arguments equivalently apply to the reduced system of equations.
2.2 Eigenvector Distribution
We mentioned above that the termination for a Krylov subspace method is
related to the location of the eigenvalues and the number of corresponding
linearly independent eigenvectors. In order to establish the association between
eigenvectors and eigenvalues we expand the general eigenvalue problem (2.7),
yielding
From (2.11) it may be deduced that either In the former case
equations (2.9) and (2.10) simplify to
which can consequently be written as
Y Z
and
z
. Since Q is orthogonal, the general
eigenvalue problem (2.12) is equivalent to considering
with w 6= 0 if and only if oe = 1. There are m linearly independent eigenvectors
corresponding to linearly
independent eigenvectors (corresponding to eigenvalues
Now suppose - 6= 1, in which case x Equations (2.9) and (2.10) yield
6 Carsten Keller, Nick Gould and Andy Wathen
The general eigenvalue problem (2.14) defines
m) of these are not equal to 1 and for which two cases have to be
distinguished. If x z 6= 0, y must satisfy
from which follows that the corresponding eigenvectors are defined by
If x we deduce from (2.15) that
and hence that As
z x T
0 in this case, no
extra eigenvectors arise.
Summarising the above, it is evident that P has
We now show that, under realistic assumptions, these eigenvectors are in fact
linearly independent.
Theorem 2.2. Let A 2 IR (n+m)\Theta(n+m) be a symmetric and indefinite matrix
of the form
where A 2 IR n\Thetan is symmetric and B 2 IR m\Thetan is of full rank. Assume the
preconditioner G is defined by a matrix of the form
n\Thetan is symmetric, G 6= A and B 2 IR m\Thetan is as above. Let Z
denote an n \Theta (n \Gamma m) basis for the nullspace of B and suppose that Z T GZ
is positive definite. The preconditioned matrix G \Gamma1 A has n+m eigenvalues
as defined by Theorem 2.1 and linearly independent eigenvectors.
There are
eigenvectors of the form
that correspond to the
case
Constraint Preconditioning 7
eigenvectors of the form
z x T
arising from
z
linearly independent,
m) eigenvectors of the form
that
correspond to the case - 6= 1.
Proof. To prove that the m eigenvectors of P are linearly independent
we need to show
y (1)
a (1)a (1)
x z
(2)
(2)
x y
(2)
(2)
y (2)
a (2)a (2)
x z
y (3)
a (3)a (3)
implies that the vectors a are zero vectors. Multiplying
by A and G \Gamma1 , and recalling that in the previous equation the first
matrix arises from the case - m), the second matrix from
the case - the last matrix arises
from - k
y (1)
a (1)a (1)
x z
(2)
(2)
x y
(2)
(2)
y (2)
a (2)a (2)
x z
y (3)
a (3)a (3)
Subtracting Equation (2.16) from (2.17) we obtain6 6 4
x z
y (3)
a (3)a (3)
8 Carsten Keller, Nick Gould and Andy Wathen
which simplifies to6 6 4
x z
y (3)
a (3)a (3)
since - k
The assumption that Z T GZ is positive definite implies that x z
in (2.18) are linearly independent and thus a (3)
Similarly, a (2)
follows from the linear independence of
x z
(2)
x y
(2)
thus (2.16) simplifies to6 6 4
y (1)
a (1)a (1)
But y (1)
are linearly independent and thus a (1)
Remark 2. Note that the result of Theorem 2.2 remains true if Z T (flA
is positive definite for some scalars fl and oe-see Parlett [17, p. 343] for details.
To show that the eigenvector bounds of Theorem 2.2 can in fact be attained,
consider the following two examples.
Example 3 (Minimum bound). Consider the matrices
so that 2. The preconditioned matrix P has an eigenvalue at
1 with multiplicity 3, but only one eigenvector arising from case (1) in Theorem
2.2. This eigenvector may be taken to be
. 2
Example 4 (Maximum bound). Let A 2 IR 3\Theta3 be defined as in Example 3,
but assume A. The preconditioned matrix P has an eigenvalue at 1 with
multiplicity 3 and clearly a complete set of eigenvectors. These may be taken
to be
, and
Constraint Preconditioning 9
2.3 Eigenvalue Bounds
It is apparent from the calculations in the previous section that the eigenvalue
at 1 with multiplicity 2m is independent of the choice of G in the preconditioner.
On the contrary, the eigenvalues that are defined by (2.14) are highly
sensitive to the choice of G. If G is a close approximation of A, we can expect a
more favourable distribution of eigenvalues and consequently may expect faster
convergence of an appropriate iterative method. In order to determine a good
factorisation of A it will be helpful to find intervals in which the
are located. If G is a positive definite matrix one possible approach is
provided by Cauchy's interlace theorem.
Theorem 2.3 (Cauchy's Interlace Theorem, see [17, Theorem 10.1.1]).
Suppose n\Thetan is symmetric, and that
Label the eigenpairs of T and H as
Then
Proof. See Parlett [17, p. 203]. 2
The applicability of Theorem 2.3 is verified by recalling the definitions of Q and
Z given in the previous section, and by considering the generalised eigenvalue
problems
and
Carsten Keller, Nick Gould and Andy Wathen
Since G is positive definite so is Q T GQ, and we may therefore write
RR T . Rewriting
and (2.20) gives
and
where
Now, since the matrix M \Gamma1 Q T AQM \GammaT is similar to G
defines the same eigenvalues ff i A. We may therefore
apply Theorem 2.3 directly. The result is that the
In particular, the - i are
bounded by the extreme eigenvalues of G \Gamma1 A so that the - i will necessarily be
clustered if G is a good approximation of A. Furthermore, a good preconditioner
G for A implies that Z T GZ is at least as good a preconditioner for Z T AZ. To
show that the preconditioner Z T GZ can in fact be much better, consider the
following example, taken from the CUTE collection [3].
Example 5. Consider the convex quadratic programming problem BLOWEYC
which may be formulated as
subject to
Z
Selecting a size parameter of 500 discretisation intervals defines a set of linear
equations of the form (1.1), where Letting G be the
diagonal of A, we may deduce by the above theory that the extreme eigenvalues
of G \Gamma1 A give a lower and upper bound for the defined by
the general eigenvalue problem (2.14). In Figure 2.1 (a) the 1002 eigenvalues
of G \Gamma1 A are drawn as vertical lines, whereas Figure 2.1 (b) displays the 500
eigenvalues of (Z T GZ)
The spectrum of Figure 2.1 (a) is equivalent to a graph of the entire spectrum
of P, but with an eigenvalue at 1 and multiplicity 502 removed. Rounded to two
decimal places the numerical values of the two extreme eigenvalues of G \Gamma1 A are
Constraint Preconditioning 11
-0.4
size of eigenvalue
(a) Eigenvalues of G \Gamma1 A
-0.4
size of eigenvalue
(b) Eigenvalues of (Z T GZ)

Figure

2.1: Continuous vertical lines represent the eigenvalues of (a) G \Gamma1 A and
0:02 and 1:98, whereas the extreme eigenvalues of (Z T GZ) are given
by 0:71 and 1. Note that for this example a large number of eigenvalues of
are clustered in the approximate intervals [0:02; 0:38] and [1:65; 1:97]. The
eigenvalue distribution in Figure (2.1) (b) reveals that there is one eigenvalue
near 0:71 and a group of eigenvalues near 1. It follows that any appropriate
iterative method that solves (1.5) can be expected to converge in a very small
number of steps; this is verified by the numerical results presented in Section 5.It is readily seen from Example 5 that in this case the bounds provided by
Theorem 2.3 are not descriptive in that there is significantly more clustering of
the eigenvalues than implied by the theorem.
Convergence
In the context of this paper, the convergence of an iterative method under pre-conditioning
is not only influenced by the spectral properties of the coefficient
matrix, but also by the relationship between the dimensions n and m. In par-
ticular, it follows from Theorem 2.1 that in the special case when the
preconditioned linear system (1.5) has only one eigenvalue at 1 with multiplicity
2n. For gives an eigenvalue at 1 with multiplicity 2m and
(generally distinct) eigenvalues whose value may or may not be equal to
1. Before we examine how these results determine upper bounds on the number
12 Carsten Keller, Nick Gould and Andy Wathen
of iterations of an appropriate Krylov subspace method, we recall the definition
of the minimum polynomial of a matrix.
Definition 1. Let A 2 IR (n+m)\Theta(n+m) . The monic polynomial f of minimum
degree such that is called the minimum polynomial of A.
The importance of this definition becomes apparent when considering subsequent
results and by recalling that similar matrices have the same minimum
polynomial.
The Krylov subspace theory states that the iteration with any method with
an optimality property such as GMRES will terminate when the degree of the
minimum polynomial is attained-see Axelsson [1, p. 463] (To be precise, the
number may be less in special cases where b is a combination of a few eigenvectors
that affect the 'grade' of A with respect to b). In particular, the degree of
the minimum polynomial is equal to the dimension of the corresponding Krylov
subspace (for general b) and so the following theorems are relevant.
Theorem 3.1. Let A 2 IR (n+m)\Theta(n+m) be a symmetric and indefinite matrix
of the form
where A 2 IR n\Thetan is symmetric and B 2 IR m\Thetan is of full rank. Let
If A is preconditioned by a matrix of the form
where G 2 IR n\Thetan , G 6= A and B 2 IR m\Thetan is as above, then the Krylov
subspace K(P; b) is of dimension at most 2 for any b.
Proof. Writing the preconditioned system (2.8) in its explicit form we
observe that P is in fact given by
I 0
\Upsilon I
where \Upsilon is non-zero if and only if A 6= G. To show that the dimension
of the corresponding Krylov subspace is at most 2 we need to determine
the minimum polynomial of the system. It is evident from (3.23) that the
Constraint Preconditioning 13
eigenvalues of P are all 1 and so the
minimum polynomial is of order 2. 2
Remark 3. It is of course possible in the case to solve the (square)
constrained equation Bx and then to obtain x
gives motivation for why the result of Theorem 3.1 is independent of G.
Remark 4. The important consequence of Theorem 3.1 is that termination
of an iteration method such as GMRES will occur in at most 2 steps for any
choice of b, even though the preconditioned matrix is not diagonalisable (unless
Theorem 3.2. Let A 2 IR (n+m)\Theta(n+m) be a symmetric and indefinite matrix
of the form
where A 2 IR n\Thetan is symmetric and B 2 IR m\Thetan is of full rank. Assume
n and that A is non-singular. Furthermore, assume A is preconditioned by
a matrix of the form
n\Thetan is symmetric, G 6= A and B 2 IR m\Thetan is as above. If Z T GZ
is positive definite, where Z is an n \Theta (n \Gamma m) basis for the nullspace of B,
then the dimension of the Krylov subspace K(P; b) is at most n \Gamma m+ 2.
Proof. From the eigenvalue derivation in Section 2.1 it is evident that the
characteristic polynomial of the preconditioned linear system (1.5) is
To prove the upper bound on the dimension of the Krylov subspace we need
to show that the order of the minimum polynomial is less than or equal to
2.
14 Carsten Keller, Nick Gould and Andy Wathen
Expanding the polynomial
obtain a matrix of the form6 6 4
Here \Psi n\Gammam is defined by the recursive formula
\Theta
with base cases \Psi
Note that the (2; 1), (2; 2) and (3; 2) entries of matrix (3.24) are in fact zero,
since the - i m) are the eigenvalues of S, which is similar to
a symmetric matrix and is thus diagonalisable. Thus (3.24) may be written
as 2
and what remains is to distinguish two different cases for the value of \Phi n\Gammam ,
that is \Phi In the former case the order the minimum
polynomial of P is less than or equal to thus the dimension
of the Krylov subspace K(P; b) is of the same order. In the latter case the
dimension of K(P; b) is less than or equal to n \Gamma m+ 2 since multiplication
of (3.25) by another factor (P \Gamma I) gives the zero matrix.
The upper bound on the dimension of the Krylov subspace, as stated in Theorem
3.2, can be reduced in the special case when (Z repeated
eigenvalues. This result is stated in Theorem 3.3. The following (ran-
domly generated) example shows that the bound in Theorem 3.2 is attainable.
Example 6. Let A 2 IR 6\Theta6 and B T 2 IR 6\Theta2 be given by
2:69 1:62 1:16 1:60 0:81 \Gamma1:97
1:62 6:23 \Gamma1:90 1:89 0:90 0:05
0:81 0:90 \Gamma0:16 0:01 1:94 0:38
Constraint Preconditioning 15
and assume that diag(A). For the above matrices the (3; 1) entry of (3.25)
is
\Gamma0:22 \Gamma0:02
It follows that the minimum polynomial is of order 6 and thus the bound given
in Theorem 3.2 is sharp. 2
Theorem 3.3. Let A 2 IR (n+m)\Theta(n+m) be a symmetric and indefinite matrix
of the form
where A 2 IR n\Thetan is symmetric and B 2 IR m\Thetan is of full rank. Assume
non-singular and that A is preconditioned by a matrix of the
n\Thetan is symmetric, G 6= A and B 2 IR m\Thetan is as above.
Furthermore, let Z be an n \Theta (n \Gamma m) basis for the nullspace of B and
assume (Z T GZ) m) distinct eigenvalues
of respective multiplicity - i , where
the dimension of the Krylov subspace K(P; b) is at most k 2.
Proof. The proof is similar to the one for Theorem 3.2. In the case
when (Z distinct eigenvalues of multiplicity - i we
may, without loss of generality, write the characteristic polynomial of P as
Y
Y
Expanding (y) we obtain the matrix6 6 4
Carsten Keller, Nick Gould and Andy Wathen
Y
Here \Psi k is given by the recursive formula
\Theta
with base cases \Psi
Note that the (2; 1), (2; 2) and (3; 2) blocks of matrix (3.26) are in fact zero.
It follows that, for \Phi k 6= 0, a further multiplication of (3.26) by
the zero matrix and thus the dimension of Krylov subspace K(P; b) is less
then or equal to k 2. 2
To verify that the bound in Theorem 3.3 is attainable consider the following
example.
Example 7. Let A 2 IR 4\Theta4 , G 2 IR 4\Theta4 and B T 2 IR 4\Theta1 be given by
. Then two of the eigenvalues that are
defined by the generalised eigenvalue problem (2.14) are distinct and given by
[2; 4]. It follows that the (3; 1) entry of (3.26) is non-zero with
and so the minimum polynomial is of order 4. 2
Implementation
There are various strategies that can be used to implement the proposed pre-
conditioner, two of which are used in the numerical results in Section 5. The
first strategy applies the standard (preconditioned) GMRES algorithm [20],
where the preconditioner step is implemented by means of a symmetric indefinite
factorisation of (1.4). Such a factorisation of the preconditioner may be
much less demanding than the factorisation of the initial coefficient matrix if
G is a considerably simpler matrix than A. The second approach, discussed in
the next section, is based on an algorithm that solves a reduced linear system
Constraint Preconditioning 17
4.1 Conjugate Gradients on a Reduced Linear System
In [11] Gould et al. propose a Conjugate Gradient like algorithm to solve
equality constrained quadratic programming problems such as the one described
in Example 1. The algorithm is based on the idea of computing an implicit basis
Z which spans the nullspace of B. The nullspace basis is then used to remove
the constraints from the system of equations, thus allowing the application of
the Conjugate Gradients method to the (positive definite) reduced system.
Assume that W GZ is a symmetric and positive definite preconditioner
matrix of dimension (n \Gamma m) \Theta (n \Gamma m) and Z is an n \Theta (n \Gamma m) matrix.
The algorithm can then be stated as follows.
Algorithm 4.1: Preconditioned CG for a Reduced System.
(1) Choose an initial point x satisfying
(2) Compute
\Gammag
(3) Repeat the following steps until
x
The computation of the preconditioned residual in (4.29) is often the most
expensive computational factor in the algorithm. Gould et al. suggest avoiding
the explicit use of the nullspace Z, but instead to compute g + by applying a
Carsten Keller, Nick Gould and Andy Wathen
symmetric indefinite factorisation of
r +#
In practice (4.30) can often be factored efficiently by using the MA27 package
of the Harwell Subroutine Library when G is a simple matrix block, whereas
the direct application of MA27 to the original system (1.1) is limited by space
requirements as well as time for large enough systems [6]. In this context the
factorisation consists of three separate routines, the first two of which analyse
and factorise the matrix in (4.30). They need to be executed only once in Step
(1) of Algorithm 4.1. Repeated calls to the third routine within MA27 apply
forward- and backward-substitutions to find the initial point x in Step (1), solve
for g in (4.27) and also to find g + in (4.29).
Remark 5. The computation of the projected residual g + is often accompanied
by significant roundoff errors if this vector is much smaller than the residual
Iterative refinement is used in (4.28) to redefine r + so that its norm is closer
to that of g + . The result is a dramatic reduction of the roundoff errors in the
projection operation-see Gould et al. [11].
5 Numerical Results
We now present the results of numerical experiments that reinforce the analysis
given in previous sections. The test problems we use are partly randomised
sparse matrices (Table 5.1) and partly matrices that arise in linear and non-linear
optimization (Table 5.2)-see Bongartz et al. [3]. As indicated through-
out, all matrices are of the form
where A 2 IR n\Thetan is symmetric, B 2 IR m\Thetan has full rank and m - n.
Four different approaches to finding solutions to (1.1) are compared-three
iterative algorithms based on Krylov subspaces, and the direct solver MA27
which applies a sparse variant of Gaussian elimination-see Duff and Reid [6].
To investigate possible favourable aspects of preconditioning it makes sense to
compare unpreconditioned with preconditioned solution strategies. The indefinite
nature of matrix (5.31) suggests the use of MINRES in the unpreconditioned
case. As outlined in Section 4 we employ two slightly different strategies
in order to implement the preconditioner G. The first method applies a standard
(full) GMRES(A) code (PGMRES in Tables 5.1 and 5.2 below), which is
Constraint Preconditioning 19
mathematically equivalent to MINRES(A) for symmetric matrices A, whereas
the second approach implements Algorithm 4.1 (RCG in Tables 5.1 and 5.2 be-
low). The choice in the preconditioner is made for both PGMRES
and RCG.
Random I Random II Random III Random IV
non-zero entries in A 2316 9740 39948 39948
non-zero entries in B 427 1871 3600 686
MINRES # of iterations 174 387 639 515
time in seconds 0:4 3:1 17:5 13:1
PGMRES # of iterations 46 87 228 242
time in seconds 0:2 3:9 96:0 108:9
RCG # of iterations 36 67 197 216
time in seconds 0:1 1:0 5:9 5:9
time in seconds 0:1 0:9 5:4 2:9

Table

5.1: Random test problems
All tests were performed on a SUN Ultra SPARCII-300MhZ (ULTRA-30)
workstation with 245 MB physical RAM and running SunOS Release 5:5:1.
Programs were written in standard Fortran 77 using the SUN WorkShop f77
compiler (version 4:2) with the -0 optimization flag set. In order to deal with
large sparse matrices we implemented an index storage format that only stores
non-zero matrix elements-see Press et al. [19]. The termination criterion for
all iterative methods was taken to be a residual vector of order less than 10 \Gamma6
in the 2-norm.
As part of its analysis procedure, MA27 accepts the pattern of some coefficient
matrix and chooses pivots for the factorisation and solution phases of
subsequent routines. The amount of pivoting is controlled by the special parameter
Modifying u within its positive range influences
the accuracy of the resulting solution, whereas a negative value prevents any
pivoting-see Duff and Reid [6]. In this context, the early construction of some
of the test examples with the default value accompanied by difficulties
in the form of memory limitations. We met the trade-off between less use
of memory and solutions of high enough accuracy by choosing the parameter
value Tables 5.1 and 5.2 .
The time measurements for the eight test examples indicate that the itera-
Carsten Keller, Nick Gould and Andy Wathen
non-zero entries in A 3004 672 1020 13525
non-zero entries in B 2503 295 148 50284
MINRES # of iterations 363 no convergence 51 180
time in seconds 1:8 no convergence 0:1 14:2
PMGREM # of iterations
time in seconds 0:6 0:1 0:2 13:2
RCG # of iterations
time in seconds 0:6 0:1 0:2 16:2
time in seconds 0:5 0:1 0:1 15:6

Table

5.2: CUTE test problems
tion counts for each of the three proposed iterative methods are comparable as
far as operation counts, i.e. work, is concerned. The numerical results suggest
that the inclusion of the (1; 2) and (2; 1) block of A into the preconditioner,
together with results in a considerable reduction of iterations,
where the appropriate bounds of Theorems 3.1, 3.2 and 3.3 are attained in all
cases. Specifically, Theorem 3.1 applies in context of problem CVXQP1.
Test problems RANDOM III and RANDOM IV in Table 5.1 emphasise the
storage problems that are associated with the use of long recurrences in the
PGMRES algorithm. The time required to find a solution to both RANDOM
III and RANDOM IV via the PGMRES algorithm is not comparable to any
of the other methods, which is due to the increased storage requirements and
the data trafficking involved. A solution to the memory problems is to restart
PGMRES after a prescribed number of iterations, but the iteration counts for
such restarts would not be comparable with those of full PGMRES.
The relevance of the time measurements for MA27 are commented on in the
next section.
6 Conclusion
In this paper we investigated a new class of preconditioner for indefinite linear
systems that incorporate the (2; 1) and (2; 2) blocks of the original matrix.
These blocks are often associated with constraints. In our numerical results
we used a simple diagonal matrix G to approximate the (1; 1) block of A,
Constraint Preconditioning 21
even though other approximations, such as an incomplete factorisation of A,
are possible. We first showed that the inclusion of the constraints into the
preconditioner clusters at least 2m eigenvalues at 1, regardless of the structure
of G. However, unless G represents A exactly, P does not have a complete set
of linearly independent eigenvectors and thus the standard convergence theory
for Krylov subspace methods is not readily applicable.
To find an upper bound on the number of iterations, required to solve linear
system of the form (1.1) by means of appropriate subspace methods, we used
a minimum polynomial argument. Theorem 3.1 considers the special condition
m, in which case termination is guaranteed in two iterations. For
Theorem 3.2 gives a general (sharp) upper bound on the dimension of the Krylov
subspace, whereas Theorem 3.3 defines a considerably stronger result if some
of the defined by (Z T GZ) are repeated.
In the special case when G is a positive definite matrix block we were able to
apply Cauchy's interlacing theorem in order to give an upper and lower bound
for the eigenvalues that are defined by the (2; 2) block of matrix (2.8).
To confirm the analytical results in this paper we used three different sub-space
methods, MINRES of Paige and Saunders for the unpreconditioned matrix
system and RCG of Gould et al. and also PGMRES of Saad and Schultz for
the preconditioned case. Overall, the results show that the number of iterations
is decreased substantially if preconditioning is applied. The Krylov subspaces
that are built during the execution of the two preconditoned implementations
are in theory of equal dimension for any of the eight test examples, and thus
PGMRES and RCG can be expected to terminate in the same number of steps.
However, convergence to any prescribed tolerance may occur for a different
number of steps since PGMRES and RCG minimize different quantities. This
can be seen in some of the examples. Nevertheless, we note that convergence
for both methods is attained much earlier than suggested by the bounds in
Theorems 3.1, 3.2 and 3.3.
The time measurements for MA27 in the last section suggest that the preconditioned
conjugate gradients algorithm, discussed in Section 4.1, is a suitable
alternative to the direct solver. Whereas both MINRES and especially PGM-
RES are considerably slower than MA27, the timings for RCG are in virtually
all cases comparable. For problems of large enough dimension or bandwidth
the resources required by MA27 must become prohibitive in which case RCG
becomes even more competitive.
22 Carsten Keller, Nick Gould and Andy Wathen

Acknowledgments

The authors would like to to thank Gene. H. Golub for his insightful comments
during the process of the work.



--R

Cambridge University Press

CUTE: Constrained and Unconstrained Testing Environment
Linearly constrained optimization and projected preconditioned conjugate gradients
Direct Methods for Sparse Matrices
The multifrontal solution of indefinite sparse symmetric linear equations
Perturbation of eigenvalues of preconditioned Navier-Stokes operators
Polynomial Based Iteration Methods for Symmetric Linear Systems
Matrix Computations
An iteration for indefinite system and its application to the Navier-Stokes equations
On the solution of equality constrained quadratic programming problems arising in optimiza- tion
Iterative Methods for Solving Linear Systems

Iterative methods for non-symmetric linear systems
Iterative Methods for Linear and Nonlinear Equations

The Symmetric Eigenvalue Problem

Numerical Recipes in Fortran: The Art of Scientific Computing
GMRES: a generalised minimal residual algorithm for solving nonsymmetric linear systems
--TR

--CTR
Joo-Siong Chai , Kim-Chuan Toh, Preconditioning and iterative solution of symmetric indefinite linear systems arising from interior point methods for linear programming, Computational Optimization and Applications, v.36 n.2-3, p.221-247, April     2007
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
Z. Dostl, An optimal algorithm for a class of equality constrained quadratic programming problems with bounded spectrum, Computational Optimization and Applications, v.38 n.1, p.47-59, September 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. Cafieri , M. D'Apuzzo , V. Simone , D. Serafino, Stopping criteria for inner iterations in inexact potential reduction methods: a computational study, Computational Optimization and Applications, v.36 n.2-3, p.165-193, April     2007
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
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
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
indefinite matrices;preconditioning;krylov subspace methods