Stability of the Diagonal Pivoting Method with Partial Pivoting.
LAPACK and LINPACK both solve symmetric indefinite linear systems using the diagonal pivoting method with the partial pivoting strategy of Bunch and Kaufman [Math. Comp., 31 (1977), pp. 163--179]. No proof of the stability of this method has appeared in the literature. It is tempting to argue that the diagonal pivoting method is stable for a given pivoting strategy if the growth factor is small. We show that this argument is false in general and give a sufficient condition for stability. This condition is not satisfied by the partial pivoting strategy because the multipliers are unbounded. Nevertheless, using a more specific approach we are able to prove the stability of partial pivoting, thereby filling a gap in the body of theory supporting LAPACK and LINPACK.
Introduction
. LAPACK is renowned for the numerical reliability of the algorithms
it employs. The LAPACK Users' Guide [1] states that "almost all the
algorithms in LAPACK (as well as LINPACK and EISPACK) are [normwise back-
ward] stable" [1, p. 74], and the algorithms not covered by this statement are known
to be stable in appropriately weakened senses. The analyses to back up these claims
of stability are spread throughout the research literature of the last 35 years. While
writing the book Accuracy and Stability of Numerical Algorithms [14] we realised that
there is no proof in the literature of the stability of the method used in LAPACK and
LINPACK for solving symmetric indefinite linear systems. Furthermore, the stability
is not a direct consequence of existing results. The purpose of this paper is to prove
the stability of the method and thereby to fill a gap in the body of theory supporting
LAPACK and LINPACK.
In the remainder of the introduction we briefly describe the method to be anal-
ysed: the diagonal pivoting method with the partial pivoting strategy of Bunch and
Kaufman [5].
Let A 2 IR n\Thetan be symmetric. If A is nonzero, we can find a permutation \Pi and
an integer so that
with E nonsingular. Then we can compute the factorization
I n\Gammas
This process can be repeated recursively on the (n \Gamma s) \Theta (n \Gamma s) Schur complement
Department of Mathematics, University of Manchester, Manchester, M13 9PL, England
(na.nhigham@na-net.ornl.gov). This work was supported by Engineering and Physical Sciences
Research Council grants GR/H/52139 and GR/H/94528.
N. J. HIGHAM
The result is a factorization
where L is unit lower triangular and D is block diagonal with each diagonal block
having dimension 1 or 2. This factorization is essentially a symmetric block form
of Gaussian elimination, with pivoting, and it costs n 3 =3 flops 1 (the same cost as
Cholesky factorization of a positive definite matrix) plus the cost of determining the
permutations \Pi . This method for computing a block LDL T factorization is called the
diagonal pivoting method. Given the factorization (1.2) of a nonsingular A, a linear
system readily solved by substitution, and by solving 2 \Theta 2 linear systems
corresponding to any 2 \Theta 2 diagonal blocks of D.
The strategy for choosing \Pi is crucial for achieving stability. Bunch and Parlett
proposed a complete pivoting strategy, which requires the whole active submatrix
to be searched on each stage of the factorization and therefore requires up to n 3 =6
comparisons. Bunch [3] proved that the diagonal pivoting method with complete pivoting
satisfies a backward error bound almost as good as that for Gaussian elimination
with complete pivoting. Bunch and Kaufman [5] devised a partial pivoting strategy
that searches at most two columns at each stage and so requires only O(n 2 ) com-
parisons. The LAPACK driver routines xSYSV (simple) and xSYSVX (expert) and the
routines xSIFA/xSISL all use the diagonal pivoting method with partial
pivoting to solve a linear system with a symmetric (indefinite) coefficient matrix.
To describe the partial pivoting strategy it suffices to define the pivot choice for
the first stage of the factorization. Recall that s denotes the size of the pivot block.
Algorithm 1 (Bunch-Kaufman Partial Pivoting Strategy). This algorithm determines
the pivot for the first stage of the diagonal pivoting method with partial
pivoting applied to a symmetric matrix A 2 IR n\Thetan .
there is nothing to do on this stage of the elimination.
r := minfi - 2: ja i1
if ja 11 j - ff-
else
if ja 11 joe - ff- 2
else if ja rr j - ffoe
choose \Pi to swap rows and columns 1 and r.
else
choose \Pi to swap rows and columns 2 and r,
so that j(\PiA\Pi T )
1 A flop is a floating point addition, subtraction, multiplication or division.
To understand the partial pivoting strategy it helps to consider the
a
.
and to note that the pivot is one of a 11 , a rr and
\Theta a11
arr
(or, rather, since
this matrix with - replaced by a r1 ).
The value of the constant
17)=8 is determined by regarding ff as a free
parameter and equating a bound for the element growth over two stages to a
bound for the element growth over one see [5] or [14] for the details.
A growth factor can be defined for the diagonal pivoting method in just the same
way as for Gaussian elimination:
where the a (k)
ij are the elements of the Schur complements arising in the course of
the factorization. From the derivation of the constant ff it is easy to show that
partial pivoting, which is larger than the bound
elimination with partial pivoting (GEPP). But, it seems that as
for GEPP, large element growth is rare in practice [5], [9].
2. Stability of the Diagonal Pivoting Method. Since the growth factor for
the diagonal pivoting method with partial pivoting is bounded, and is usually small
in practice, does it not follow that the method is stable in the same sense as for
GEPP? This is a tempting argument, and one that is neither used nor warned against
in the existing literature. However, it is easy to show that the argument is false, by
exhibiting an example where the diagonal pivoting method has a small growth factor
but is unstable. An example (not produced by partial pivoting) is, with
with a 2 \Theta 2 pivot followed by a 1 \Theta 1 pivot,
\Gammaffl \Gammaffl \Gamma15
The growth factor ae n is 1, yet kLk1=kAk1 is unbounded as ffl ! 0,
which suggests that the factorization, however it is computed, may not provide a stable
way to solve linear systems in finite precision arithmetic. The instability is
confirmed by a Matlab experiment, in which the unit roundoff
We solved a linear system in two different ways.
First, we computed the factorization in (2.1) using the diagonal pivoting method, as
specified in (1.1) (with taking a 2 \Theta 2 pivot on the first step and using GEPP
4 N. J. HIGHAM

Table
Backward error for computed solution of indefinite system of order 3.
Diagonal Explicit
pivoting factors
to solve linear systems involving this pivot. For comparison, we evaluated the explicit
formulae for the LDL T factors in (2.1), and used the explicit inverse of D(1: 2; 1: 2)
when solving the linear system involving D. Table 2.1 shows the normwise relative
backward error of the computed solution b x,
(see [16] or [14, Th. 7.1] for a proof of the latter equality), which would be of order
u for a stable solution method. As ffl decreases the computations become unstable.
We note that stability is obtained if, in (1.1), we take the natural 1 \Theta 1 pivot a 11
instead of the ill conditioned 2 \Theta 2 pivot A(1: 2; 1: 2); interestingly, though, the 2 \Theta 2
pivot shares with those chosen by the Bunch-Kaufman partial pivoting strategy the
property that it is indefinite. Partial pivoting is stable on this example.
We conclude that a small growth factor is not, by itself, enough to guarantee
stability of the diagonal pivoting method. A sufficient condition for stability can
be obtained by regarding the block LDL T factorization computed by the diagonal
pivoting method as a special case of a block LU factorization. Error analysis for
block LU factorization is given by Demmel, Higham and Schreiber [8], and a suitable
modification of this analysis gives the following result: if linear systems involving 2 \Theta 2
pivots are solved in a normwise backward stable fashion then the condition
for a modest constant c n , is sufficient to ensure that the diagonal pivoting method
produces a factorization with a small relative residual and provides computed solutions
to linear systems that have a small backward error. Unfortunately, condition (2.2)
does not hold for the partial pivoting strategy of Bunch and Kaufman, as is shown by
the following example. For ffl ? 0, the diagonal pivoting method with partial pivoting
produces the factorization, with
As ffl ! 0, kLk1kDk1kL T k1 =kAk1 !1, and indeed the multipliers are unbounded.
Even 1 \Theta 1 pivots can lead to arbitrarily large elements in L, as the following example
with shows (again, partial pivoting selects
A
It is worth emphasizing that large elements in a factor of a matrix do not necessarily
imply that the factorization is unstable. For example, in the (point) LDL T
factorization of a symmetric positive definite matrix A with
the ratio kLk1=kAk1 can be arbitrarily large, yet the factorization is guaranteed to
be stable. One such example is, with ffl ? 0,
Our conclusion is that existing results for LU factorization and block LU factorization
do not directly imply the stability of the diagonal pivoting method with
partial pivoting. Any proof of stability must make use of the particular properties of
the partial pivoting strategy.
The only claims of stability that we have found in the literature are in the paper
by Bunch, Kaufman and Parlett [6] and in the LINPACK Users' Guide [9, p. 5.19];
in both cases, residual bounds of the form
are stated
without proof, where p is a polynomial; we prove a result of this form and, in Theorem
4.2, a backward error result for the computed solution of We note that
much of Bunch's analysis of the diagonal pivoting method in [3] is specific to complete
pivoting, so his analysis does not readily yield results for partial pivoting.
In the rest of the paper we present a new analysis to show that partial pivoting
is indeed a stable pivoting strategy for the diagonal pivoting method.
3. Background Results from Error Analysis. We collect in this section
some standard error analysis results that will be needed later. For our model of
floating point arithmetic we take
where u is the unit roundoff. All the results we quote remain true under a weaker
model that accommodates machines without a guard digit [14, x2.4], provided some
of the constants are increased slightly.
We introduce the constant
which carries with it the implicit assumption that nu ! 1. Useful properties are (a)
Proofs of the following results can be found in [14]. First, for matrix multiplication

n\Thetan is a nonsingular triangular matrix and the system
solved by substitution then
6 N. J. HIGHAM
Third, if a linear system n\Thetan , is solved without breakdown by
Gaussian elimination without pivoting, then the computed solution satisfies
U
L and b
U are the computed LU factors.
We will use the norm defined by
(for which kABk M - is the best bound of this form that holds for all
4. Error Analysis.
4.1. 2 \Theta 2 Linear Systems . Crucial to the error analysis that follows is a
backward error result for the solution of linear systems involving 2 \Theta 2 pivots. Note
that, in the notation of Algorithm 1, the pivot is
a 11 a r1
a r1 a rr
For this subsection and the later analysis, it is convenient to tabulate the conditions
that must hold for a 2 \Theta 2 pivot to be selected:
(4.1a)
(4.1b)
(4.1c)
(4.1d)
where the fourth inequality is a consequence of the previous two (note that (4.1c)
implies oe 6= 0).
Suppose, first, that linear systems are solved by GEPP. By (4.1a),
computes the LU
factorization
- a r1 a rr
a 11 a r1
a 11
a r132
4 a r1 a rr
0 a r1 \Gamma a 11 a rr
a
From (3.3), we have the backward error result
U j:
Now
a r1
a r1
using (4.1d). It follows that
using the numerical value of ff specified in Algorithm 1. Strictly, we should append
)" to this bound, to account for replacing j b
U j by a bound for jLjjU j; we
omit the second order term for the moment and reinstate it later. Note that the result
(4.2) holds trivially for a 1 \Theta 1 pivot E.
The main alternative to using GEPP to solve the systems is to use the
explicit inverse of E, as is done in the implementations of the diagonal pivoting method
with partial pivoting in LAPACK and LINPACK (see the auxiliary routine xLASYF
in LAPACK and xSIFA in LINPACK). In both LAPACK and LINPACK,
solved by evaluating
a r1
' a 11
a r1
\Delta a rr
a r1
a rr
a r1
a r17 5 b;
which corresponds to using an explicit formula for the inverse of a 2 \Theta 2 matrix (or,
equivalently, Cramer's rule), with scaling to avoid overflow. The term
a r1
\Delta a rr
a r1
appears to be a potential source of instability, since for arbitrary a 11 , a r1 and a rr the
relative error in the computed b - is unbounded. However, by exploiting the condition
(4.1d) for a 2 \Theta 2 pivot, which we rewrite as
a 2
we can obtain a very satisfactory error bound for b
-. Using the model (3.1) we have
' a 11
a r1
\Delta a rr
a r1
where 1: 4, which implies [14, Lemma 3.1]
a r1
\Delta a rr
a r1
Hence
a 2
It is then straightforward to show that, denoting the matrix in (4.3) by Z,
so that
8 N. J. HIGHAM
using (A.3). The Oettli-Prager theorem [15], [14, Th. 7.3] then implies that
Again, strictly a second order term should be added to the bound, this time to account
for the fact that jxj rather than jbxj appears on the right-hand side of (4.4).
The conclusion is that whether the linear system involving the 2 \Theta 2 pivot
is solved by GEPP or by using the explicit inverse, we have
for an integer constant c. It is worth stressing that such a result does not hold for an
arbitrary 2 \Theta 2 (symmetric) matrix E-we have fully exploited the pivoting conditions
in the derivation.
4.2. Componentwise Backward Error Analysis. Now we carry out a componentwise
backward error analysis of the diagonal pivoting method. We make only
one assumption about the pivoting strategy: that (4.5) holds for the 2 \Theta 2 pivots. For
convenience, we assume, without loss of generality, that no interchanges are needed,
which amounts to redefining A := PAP T in (1.2).
To begin, we consider the first stage of the factorization, using the notation of
(1.1). The submatrix L
l j is the jth column of L T
is the jth column of C T , then, from (4.5),
Hence, overall,
We assume that the Schur complement is computed as
The remaining stages of the diagonal pivoting method factorize the Schur complement
as
S , and we assume, inductively, that the computed factors
s; u) is a constant depending on n \Gamma s and u. We therefore have computed
factors b
L and b
D of A that satisfy
I 0
2 If the Schur complement is computed as
21 then the same bound (4.9) ensues.
Now, from (4.6) we have the inequalities
and
Using (4.7) and (4.8) we have
Overall, then, we have
Hence
where d(n; u) is clearly of the form p(n)u +O(u 2 ), where p is a linear polynomial.
Now we analyse the substitution stages when the LDL T factorization is used to
solve a linear system From (3.2) and (4.5), the computed solutions to the
three systems Ly
Thus
where j\DeltaAj is bounded in (4.9) and
On bringing back into account the row and column interchanges, we obtain the following
result.
Theorem 4.1. Let A 2 IR n\Thetan be symmetric and let b
x be a computed solution to
the linear system produced by the diagonal pivoting method with any pivoting
strategy. If for all linear systems involving 2 \Theta 2 pivots (4.5) holds, then
N. J. HIGHAM
where p is a linear polynomial and PAP T - b
is the factorization computed by
the diagonal pivoting method.
The bound in (4.10) is analogous to the bound in (3.3) that holds for Gaussian
elimination. We have already seen that the assumption (4.5) in Theorem 4.1 holds
for the partial pivoting strategy of Bunch and Kaufman, provided linear systems
are solved by GEPP or by using the explicit inverse. It is easy to show that
this assumption also holds for the complete pivoting strategy of Bunch and Parlett
under the same conditions (interestingly, for the 2 \Theta 2 pivots E that arise with
the Bunch-Parlett strategy, GEPP applied to a is identical to Gaussian
elimination with complete pivoting).
4.3. Normwise Analysis for Partial Pivoting. To show that the diagonal
pivoting method is stable for a particular pivoting strategy, we need to show that
the matrix j b
suitably bounded. We now specialise to partial pivoting.
For partial pivoting, b
L can be arbitrarily large, so stability is not an immediate
consequence of Theorem 4.1. We therefore need to look closely at the elements of
the matrix j b
j. For simplicity, we bound the matrix jLjjDjjL T j containing the
exact factors, which makes only a second order change to the overall bounds, since
Initially, we examine the contribution from the blocks of L and D produced by
the first stage of the factorization. For this more delicate part of the analysis we take
full account of the interchanges in our notation. Note that
I
I jL T
We first bound
For a 1 \Theta 1 pivot, F is a vector with elements
each of which is trivially
bounded by max
Now consider a 2 \Theta 2 pivot. Algorithm 1 dictates that \Pi in (1.1) swaps rows and
columns 2 and r so that, as noted earlier,
a 11 a r1
a r1 a rr
Using (A.1) and (4.1a), we have
Next, we need to bound
G := jL 21 jjEjjL T
First, consider a 1 \Theta 1 pivot. In cases (1) and (2) of Algorithm 1 we have
ff
case (1),
oe
ff
case (2).
In case (3),
ff
For a 1 \Theta 1 pivot,
For a 2 \Theta 2 pivot (case (4) of Algorithm 1), using (A.2) we have
oe
(using (4.1b))
The remaining blocks of jLjjDjjL T j are composed of blocks of L and D that make
factors of Schur complements of A. But every Schur complement satisfies
where ae n is the growth factor. Hence, applying the bounds above recursively to the
(2; 2) block in (4.11), we deduce the (pessimistic) bound
We mention in passing that in early drafts of this paper we had a weaker version
of (4.5) in which jEj in the bound was replaced by jEj
2 . We were still able
to obtain a satisfactory bound for k jLjjDjjL T indicating that partial pivoting is
somewhat more tolerant of how the 2 \Theta 2 systems are solved than might be thought
from the analysis above.
N. J. HIGHAM
Using the bound (4.14) in Theorem 4.1 we obtain the following normwise backward
stability result for partial pivoting.
Theorem 4.2. Let A 2 IR n\Thetan be symmetric and let b
x be a computed solution to
the linear system produced by the diagonal pivoting method with the partial
pivoting strategy of Bunch and Kaufman, where linear systems involving 2 \Theta 2 pivots
are solved by GEPP or by use of the explicit inverse. Then
where p is a quadratic.
Theorem 4.2 has the same form as Wilkinson's result for GEPP applied to a
nonsymmetric system (see, e.g., [14, x9.2]), though of course the numerical value of
ae n is usually different for the two methods.
5. Discussion. The backward error matrix \DeltaA in (4.9) is necessarily symmetric,
but that in (4.15) is not, in general. However, we can take \DeltaA in (4.15) to be
symmetric, at the cost of increasing the bound by a factor n, because of the following
result of Bunch, Demmel and Van Loan [4]: if
such that
.
Sorensen and Van Loan [10, x5.3.2] modify the Bunch-Kaufman partial pivoting
strategy by redefining, in Algorithm 1,
This small change has the pleasing effect of ensuring that for a positive definite matrix
no interchanges are done (and that, as for the Bunch-Kaufman strategy, only 1 \Theta 1
pivots are used in this case). At the same time it leaves the growth factor bound
unchanged, and all our analysis remains valid for this variant.
For sparse symmetric matrices, Duff, Reid and co-authors compute the block
using a pivoting strategy very different from that of Bunch and
Kaufman [11], [12], [13]. We describe the strategy in [13] as it applies to the first stage
of the factorization: a 11 is defined to be an acceptable 1 \Theta 1 pivot, from the point of
view of numerical stability, if
ja i1 j;
is a tolerance; the matrix
a 11 a r1
a r1 a rr
is an acceptable 2 \Theta 2 pivot if
From among the acceptable pivots one is chosen that best preserves sparsity, according
to some particular sparsity criterion. The conditions (5.1) and (5.2) ensure that kLk1
is bounded by a multiple of ' \Gamma1 , which then implies bounds on the growth factor,
and hence on kDk1 . The stability of this pivoting strategy is therefore immediate,
since (2.2) is satisfied. An interesting contrast is that the Bunch-Kaufman strategy
involves a fixed amount of searching for a pivot, and the reasons for its stability are
subtle, whereas the Duff et al. strategy more directly forces stability by bounding
the multipliers, but gives up the fixed amount of searching of the Bunch-Kaufman
strategy.
We emphasize that the aim of this work was to obtain a rigorous backward error
bound for the diagonal pivoting method with partial pivoting. The actual performance
of the method is affected by the size of the growth factor. More work is needed
to investigate the behaviour of the growth factor, about which less is known than
the growth factor for Gaussian elimination with partial pivoting. Although the unboundedness
of kLk1 does not preclude backward stability, it does have implications
for the practical behaviour of the method; see Ashcraft, Grimes and Lewis [2] for a
thorough study for both dense and sparse matrices. Finally, we mention that the
implementation of the diagonal pivoting method with partial pivoting in LAPACK
2.0 can be unstable when kLk1 is large, as pointed out and explained in [2]. The
potential instability stems from replacing a symmetric rank-2 update by two rank-1
updates, via the use of an eigendecomposition. This problem will be corrected in a
future release of LAPACK.

Acknowledgements

. It is a pleasure to thank Philip Gill and Michael Saunders
for valuable comments, particularly at early stages of this work. I also thank Jim
Bunch, Des Higham and John Lewis for suggesting improvements to draft manuscripts.


Appendix

.
In this appendix we bound three matrix expressions involving a 2 \Theta 2 pivot from
partial pivoting,
a 11 a r1
a r1 a rr
First, we note that
using (4.1d). Hence
using (4.1d) again. Next,
14 N. J. HIGHAM
Finally,



--R


Accurate symmetric indefinite linear equation solvers.
Analysis of the diagonal pivoting method
The strong stability of algorithms for solving symmetric linear systems
Some stable methods for calculating inertia and solving symmetric linear systems
Decomposition of a symmetric matrix
Direct methods for solving symmetric indefinite systems of linear equations
Stability of block LU factorization
Society for Industrial and Applied Mathematics
Solving Linear Systems on Vector and Shared Memory Computers
The factorization of sparse symmetric indefinite matrices

Direct solution of sets of linear equations whose matrix is sparse
Accuracy and Stability of Numerical Algorithms
Compatibility of approximate solution of linear equations with given error bounds for coefficients and right-hand sides
On the compatibility of a given solution with the data of a linear system
--TR
numerical stability;rounding error analysis;diagonal pivoting method;LAPACK;partial pivoting;growth factor;LINPACK;symmetric indefinite matrix;LDLT factorization