Stability of Augmented System Factorizations in Interior-Point Methods.
Some implementations of interior-point algorithms obtain their search directions by solving symmetric indefinite systems of linear equations. The conditioning of the coefficient matrices in these so-called augmented systems deteriorates on later iterations, as some of the diagonal elements grow without bound. Despite this apparent difficulty, the steps produced by standard factorization procedures are often accurate enough to allow the interior-point method to converge to high accuracy. When the underlying linear program is nondegenerate, we show that convergence to arbitrarily high accuracy occurs, at a rate that closely approximates the theory. We also explain and demonstrate what happens when the linear program is degenerate, where convergence to acceptable accuracy (but not arbitrarily high accuracy) is usually obtained.
Introduction
. We focus on the core linear algebra operation in primal-dual
interior-point methods for linear programming: solution of a system of linear equations
whose coefficient matrix is large, sparse, and symmetric. In existing codes, the linear
system is formulated in two different ways. One formulation, usually called the augmented
system formulation, has a symmetric indefinite coefficient matrix. The other
involves a more compact (but generally denser) symmetric positive-definite matrix.
A diagonal matrix D is involved in both formulations, where D has the disconcerting
property that some of its elements grow to 1 as the iterates approach the solution set.
This blowup in D can produce ill conditioning in the coefficient matrix of the linear
system. In this paper, we examine the augmented system and look at how various
factorization algorithms for this system behave as this ill conditioning develops.
We restrict our study to three standard factorization algorithms - the Bunch-
Parlett, Bunch-Kaufman, and sparse Bunch-Parlett algorithms. The last of these
has been used in at least one practical interior-point code for linear programming
(see Fourer and Mehrotra [5]). We assume that no attempt is made to improve the
conditioning of the underlying linear systems by guessing whether each component
of the solution is at a bound. Preprocessing of this kind detracts from the intuitive
appeal of interior-point algorithms, namely, that they avoid explicit guessing about
the contents of the basis.
In numerical experiments with feasible linear programs, we find that two distinct
scenarios arise.
1. Even when the iterates are very close to the solution set, the computed search
directions are good enough to produce rapid convergence of the algorithm at
nearly the rates predicted by the theory. This performance is a little sur-
prising. Since the matrix is poorly conditioned, we might have expected the
computed directions to be too inaccurate to allow the algorithm to make much
progress. This scenario usually occurs when the underlying linear program
has a unique primal-dual solution.
2. Near the solution, calculation of the search direction fails because of break-down
of the matrix factorization, or else the computed search direction is so
inaccurate that the interior-point method can move only a tiny distance along
Mathematics and Computer Science Division, Argonne National Laboratory, 9700 South Cass
Avenue, Argonne, IL 60439. This work was supported by the Mathematical, Information, and Computational
Sciences Division subprogram of the Office of Computational and Technology Research,
U.S. Department of Energy, under Contract W-31-109-Eng-38.
it before violating a bound. This scenario usually occurs when the underlying
linear program is degenerate.
Our analysis in this paper explains these observations through a close examination of
the behavior of factorization algorithms on the highly structured matrices that arise
in our application. The effects of roundoff error are tracked by using fairly standard
techniques from backward error analysis.
The most successful interior-point methods for practical linear programming problems
are primal-dual methods. The best-known potential-reduction algorithm in this
class was devised by Kojima, Mizuno, and Yoshise [9]; the review paper of Todd [17]
contains a wealth of historical information on potential-reduction methods. Early
developments in path-following methods are surveyed by Gonzaga [7], while Mizuno,
Todd, and Ye [15] describe an important variant of these methods that does not require
the iterates to stay within a cramped neighborhood of the central path. Zhang
[25] extended the path-following approach further, allowing the iterates to be infeasible
while retaining global convergence and polynomial complexity; see also Wright
[21]. Some of these developments took place in the context of linear complementarity,
a class of problems that includes linear programming as a special case.
On the computational side, the OB1 code described by Lustig, Marsten, and
Shanno [10] generated search directions of the type described in this paper. They
compute the maximum step ff   that could be taken along this direction without violating
the positivity bounds, then set the actual step length to .995 ff   . Mehrotra's [14]
predictor-corrector search direction differs from the one analyzed in this paper, but
under our assumptions below, the difference vanishes as the solution is approached.
Newer codes, such as those described by Mehrotra [14], Fourer and Mehrotra [5],
Lustig, Marsten, and Shanno [12], Vanderbei [18], and Xu, Hung, and Ye [23] all implement
Mehrotra's predictor-corrector strategy. These newer codes continue to use
step lengths based on ff   ; hence, we pay particular attention to the effect of roundoff
error on this quantity.
Previous analysis of the ill-conditioned linear systems that arise in interior-point
and barrier methods has been carried out by Poncele'on [16] and Wright [22]. Pon-
cele'on [16] showed that these systems are not too sensitive to structured perturbations
from a certain class provided that the underlying optimization problem is well con-
ditioned. Wright [22] analyzed Gaussian elimination in the context of interior-point
algorithms for linear complementarity problems.
Simultaneously with the original version of this paper, and independently, Fors-
gren, Gill, and Shinnerl [4] performed an analysis of the augmented system in barrier
algorithms. Their analysis tends to be more detailed than ours, and a few of the
results overlap. However, they assume that the factorization algorithms select the
large diagonal elements as pivots before any others, a pattern that does not generally
occur in practice.
Vavasis [19] gives an illuminating discussion of the augmented system in contexts
other than optimization. He presents a solution method that is provably stable in
a certain sense, but which is not guaranteed to produce "useful" steps in the sense
of this paper. Duff [3] also discusses augmented systems in a general context and
describes a sparse factorization procedure.
2. Interior-Point Methods. We consider the linear program in standard form:
. The dual of (1) is
(2)
solution if x
is feasible for (1), (-   ; s   ) is feasible for (2), and s   and x   are complementary; that
is,
We denote the set of primal-dual solutions by S.
Each iterate (-; x; s) of a primal-dual interior-point method satisfies the strict
inequality directions are found by applying a modification of
Newton's method to the following system of nonlinear equations:
Specifically, the search
direction (\Delta-; \Deltax; \Deltas) satisfies the linear equations4
\Deltax
\Delta-
\Deltas5 =4
where oe 2 [0; 1] is known as the centering parameter and the important quantity - is
defined by
The step length ff along the search direction is determined by various factors; mini-
mally, the updated x and s components are required to stay strictly positive:
At least half the components of (x; s) - the critical components - become very close
to their lower bound of zero during the later stages of the algorithm. Despite this
property, the step length ff can be quite close to one without violating the property
(6), when the search direction (\Delta-; \Deltax; \Deltas) is an exact solution of (5). If perturbations
caused by roundoff are present in the critical components of (\Delta-; \Deltax; \Deltas),
the requirement (6) can severely curtail the allowable step length and slow the con-
vergence. Hence, it is important that the critical components of (\Delta-; \Deltax; \Deltas) be
computed to high relative accuracy. This point provides the focus for much of our
error analysis.
Throughout the paper we use u to denote unit roundoff, which we define implicitly
by the statement that when x and y are any two floating-point numbers, op denotes
+; \Gamma; \Theta; =, and fl(z) denotes the floating-point approximation of any real number z,
we have
Since our concern is with the internal workings of a single interior-point iterate, we
omit iteration counters from all quantities. For this reason, we use the order notation
O(\Delta) in a slightly unconventional way. When - and j are two nonnegative numbers,
we there is a positive constant C (not too large) such that - Cj.
We say that a matrix or vector is O(j) if its norm is O(j). We say that
For the purposes of this paper, we are mainly interested in how the factorizations
behave relative to - and u. The dimensions m and n are ignored in our use of the
notation O(\Delta).
If G is a matrix, G \Deltaj denotes its j-th column, while G i\Delta denotes the i-th row. The
matrix whose elements are jG ij j is denoted by jGj.
We use k \Delta k to denote any one of the equivalent matrix norms k
k \Delta k1 . When G is rectangular, the 2-norm condition number is defined as follows.
Definition 1. Let G be a rectangular matrix with full rank, and suppose that
svmax (G) and svmin (G) denote the largest and smallest singular values of G,
respectively. The 2-norm condition number of G is
If G is square and nonsingular, this definition coincides with the usual definition
3. Definitions and Assumptions. We assume throughout that the problems
(1), (2) are feasible; that is, there exists at least one triple (-; x; s) satisfying the
constraints implies existence of solutions
to (1), (2). The following theorem gives another consequence of feasibility.
Theorem 3.1. Suppose that (1) and (2) are feasible and that (-; x; s) is any
point with there exists a solution (\Delta-; \Deltax; \Deltas) to (5).
Proof. The proof follows from Section 6 of Wright [21]. See, in particular, Lemma
6.2, Theorem 6.3, and the remarks in the last two paragraphs of [21].
Note that A need not have full rank for Theorem 3.1 to hold.
The set of basic indices ng can be defined as
while the nonbasic set N is
It is well known that B and N form a partition of f1; ng and that there is
at least one solution (-   ; x   ; s   ) that is strictly complementary, that is, x
(Goldman and Tucker [6]). The cardinality of B is denoted by jBj. By partitioning
the columns of A according to B and N , we define
so that B is m \Theta jBj and N is m \Theta jN j. We say that the linear program is nondegenerate
and the primal-dual solution is unique. We assume also that B is reasonably
well conditioned in nondegenerate problems.
We do not confine our analysis to one specific primal-dual algorithm. Rather, we
rely on a set of assumptions that is satisfied by a variety of algorithms. The first of
these assumptions concerns the iterates, the search directions, and the relationship
between - and the current infeasibility.
Assumption 1. The sequence of iterates (-; x; s) generated by the interior-point
algorithm satisfies the following properties when - becomes sufficiently small:
(11a)
In addition, the infeasibility is O(-); that is,
Assumption 1 is not very strong. G-uler and Ye [8] study algorithms in which all
iterates are strictly feasible; that is
In fact they require that x and s be slightly separated from the boundary of the
positive orthant, in the sense that
for some constant fl 2 (0; 1). They show that all limit points of such algorithms
are strictly complementary solutions of (1), (2) and that most path-following and
potential-reduction algorithms do in fact satisfy (14). It is easy to infer from their
results that (11) holds for all subsequences that approach these limit points. Moreover,
(12) is trivially satisfied for all feasible algorithms.
The infeasible-interior-point algorithm described by Wright [20] satisfies Assumption
1. So does the algorithm in [21], provided that the sequence or iterates (x; s) is
bounded. Implemented algorithms such as those of Vanderbei [18], Lustig, Marsten,
and Shanno [10, 11], and Xu, Hung, and Ye [23] usually step a fixed multiple of the
distance to the boundary rather than enforce a potential reduction condition or a condition
like (14). Nevertheless, the iteration sequence usually satisfies the properties
of Assumption 1 for most practical problems.
Finally, we state without proof a technical lemma for use in later sections.
Lemma 3.2. Let H be a square matrix partitioned as
where H 11 and H 22 are also square. Suppose that H 11 and H 22 \Gamma H
are
nonsingular. Then H is nonsingular, and
4. Exact and Approximate Search Directions. By defining r
in (5), we obtain4
\Delta-
\Deltax
\Deltas5 =4
By eliminating \Deltas from this system, we obtain the augmented system formulation:
A
\Delta-
\Deltax
(16a)
In Wright [22], we performed an error analysis on a system like (16a), but in the
context of a specific path-following algorithm for the monotone linear complementarity
problem. Some of our results from [22] are relevant to the present case of (16), as we
discuss later.
Potential difficulties with the formulation (16) arise from two sources - possible
rank-deficiency in certain submatrices of A, and the fact that some diagonal elements
of X approach zero while others approach +1. Despite the effects of
ill conditioning and finite precision, we find that the approximate search directions
obtained from (16) by using standard factorization procedures are often remarkably
good. They allow the interior-point algorithm to take near-unit steps and to make
substantial improvements in the duality measure -. In the following theorem, we
specify a set of conditions for which this happy situation holds. In later sections, we
identify situations under which these conditions hold.
In the remainder of the paper, we use ff   to denote the largest number in [0; 1]
such that
(17a)
decreasing for ff 2 [0; ff   ]:
Theorem 4.1. Suppose that Assumption 1 holds. Let (\Delta-; \Deltax; \Deltas) be the exact
solution of (5) (equivalently, (16)), and let ( c
\Delta-; c
\Deltax; c
\Deltas) be an approximation to this
step. Suppose that the centering parameter oe in (5) lies in the range [0; 1=2] and that
the following conditions hold:
(18a)
\Deltas
(18c)
as in (17), and suppose -
ff   is obtained by replacing (\Deltax; \Deltas) with ( c
\Deltax; c
\Deltas)
in (17). Then for all - sufficiently small, we have
and
ff   c
\Deltax)
Proof. From (11a) and (18a), we have
so these components do not restrict the value of ff   . Since u is much smaller than 1,
we use (18b) as well to deduce that
\Deltas
\Deltas
Similarly, we can show that xB
For the decrease condition (17b) we show that the duality gap actually decreases
over the entire interval [0; 1] for both exact and approximate search directions, so that
this condition does not play a role in determining ff   or -
ff   . For the exact direction,
we have from (5), (18a), and oe 2 [0; 1=2] that
d
for all ff 2 [0; 1]. Hence, for - sufficiently small, the duality gap is decreasing over
[0; 1]. For the approximate direction ( c
\Deltax; c
\Deltas), this bound can be modified slightly
to account for the inexactness. We omit the details, which are straightforward but
messy, and state the conclusion as
d
dff
\Deltax)
Again, we find that the duality gap is decreasing over the whole interval ff 2 [0; 1].
Hence, the only condition that can bound ff   and -
ff   away from 1 is (17a), and
then only for the N-components of x and the B-components of s. In fact, ff   satisfiesff
From (5), we have x i \Deltas
we have
For we have from (11a) and (18a) that j\Deltax i =x
An identical argument can be used for the other term in (22), so we haveff   -
proving (19).
For the maximum step length -
ff   along the approximate direction ( c
\Deltax; c
\Deltas), we
have from (18c) and (11b) that
c
\Deltas i
c
Hence, from (22), we have-
c
\Deltas i
c
For all sufficiently small -, the estimates (20) follow immediately from this last expression

Finally, for the estimate of potential decrease (21), we have from (5) that
\Deltax)
\Deltas)
where we have used Assumption 1 and (18) to estimate the remainder terms. Finally,
we obtain (21) by substituting
ff
5. The Augmented System. In the remainder of the paper, we focus on the
procedure based on (16) for finding the search directions. In this section, we define a
generalized form of the matrix in (16a) which we call a canonical matrix. We show that
if the backward error analysis of the solution procedure satisfies a certain condition
Condition 1 below - then the approximate step ( c
\Delta-; c
\Deltax; c
\Deltas) obtained from (16)
in a finite-precision environment is "useful" in the sense of Theorem 4.1.
In later sections, we define conditions under which these standard algorithms
for solving symmetric indefinite systems satisfy Condition 1 and hence yield useful
search directions. Our sharpened, specialized error analysis yields much stronger
results than naive application of the standard results. We also gain insight into how
the algorithms work even when the nondegeneracy assumptions of Sections 6, 7, and
8 fail to hold, and why they continue to generate useful search directions even for
degenerate problems until - is quite small.
Given a symmetric matrix T of order - n, the factorization procedures yield
where P is a permutation matrix, L is unit lower triangular, and D is a block-diagonal
matrix with 1 \Theta 1 and 2 \Theta 2 diagonal blocks. We denote the counterparts of these
matrices that are actually computed in the finite-precision environment by -
L and -
D,
respectively.
Given the system T z = d and the data P , -
L, and -
D from the factorization,
we find the computed solution -
z by performing two vector permutations with P ,
triangular substitutions with -
L and -
L T , and a blockwise inversion of -
D. Each of
these operations (except the permutations) may introduce additional roundoff error,
which must be accounted for in the error analysis.
For each of the methods, we focus on a single step of the factorization procedure
applied to a matrix T with properties like those of our given system (16a), which we
now define.
Definition 2. A matrix T is a canonical matrix if it is a symmetric permutation
where
are small;
-   is diagonal with all diagonal elements of
=\Omega\Gamma16 and
We call T a degenerate canonical matrix if it has the form
where the zero blocks are nonvacuous. In keeping with our particular application
(16a), we use m and n to denote the number of rows and columns in the composite
matrix [B j N ], respectively, and -
to denote the total dimension of T .
Corresponding to our canonical matrix, we define a canonical error matrix. We
prove that for each of the factorizations, the error matrix has this form.
Definition 3. Let T be a canonical matrix. The corresponding canonical error
matrix \Delta is a matrix of the same dimension as T such that
where ffi u and the elements of \Delta u are O(u).
An important role in the pivot selection process is played by the quantities - i ,
which denote the magnitude of the largest off-diagonal element in column i, that is,
A sufficient condition for useful steps. The following condition states the
common goal of our backward error analysis of the three factorization procedures.
When this condition is satisfied along with nondegeneracy of the linear program, the
result of Theorem 4.1 holds.
Condition 1. Given the system T z = d, where T is a canonical matrix, the
symmetric factorization and solution process yields a computed solution -
z that satisfies
d;
where \Delta is a canonical error matrix associated with T and -
O(u). We allow
for a perturbed right-hand side -
d because of the nature of our particular system (16a).
The residuals r b and r c are computed as the difference of O(1) quantities, so O(u)
perturbations will appear when they are evaluated in the obvious way. Addition of
the terms s N and -X \Gamma1
may give rise to errors of similar magnitude.
Theorem 5.1. Suppose that Assumption 1 holds and that the problem is nonde-
generate, that is, moderate. Suppose that the procedure for solving
Condition 1, and denote the approximate solution to (16a) by ( c
\Delta-; c
\Deltax).
Then for all sufficiently small -, we have
and
Proof. We prove (31) by appealing to (5). By partitioning A into B and N
according to (10), and partitioning the diagonal matrices S and X accordingly, we
see that the matrix in (5) is a permutation of6 6 6 6 4
Because of (11), the diagonal elements in XB and SN
are\Omega\Gammae/3 while the matrices SB
and XN are O(-). In addition, B is square and well conditioned, so the matrix (33)
is an O(-) perturbation of a uniformly nonsingular matrix. From (5), we then have
so the result (31) follows from (11) and (12).
To derive the relative error estimate (32), consider the system (16a). By permuting
the matrix in accord with the B [ N partition, we can rewrite (16a) as follows:4
\Delta-
\Deltax N5 =4
From (11), we have for sufficiently small - that the diagonals in X \Gamma1
while the diagonals of X \Gamma1
coefficient matrix is canonical.
By defining
\Delta-
we can restate the system as
MB MN
z N
dB
dN
From our assumption on B, we have
Because of Condition 1, the computed solution -
z of (35) satisfies
'-
MB MN
'-
z N
dB
dN
O(u) and the canonical error matrix \Delta satisfies
By combining this estimate with (35) and (36), we obtain
'-
MB MN
'-
z N
O(-) from (31), we have
z N
O(u)
'-
O(-)
O(-)
O(-u)
O(u)
so when we add the effect of -
d, we find that the right-hand side of (37) is O(u).
For the coefficient matrix in (37) we use Lemma 3.2 with
Lemma 3.2 yields the following estimates:
By combining these observations with (38), we obtain
O(u)
O(-u)
giving (32).
Next, we examine the accuracy of c
\Deltas, which is calculated by substituting c
\Delta- and
c
\Deltax into (16b).
Theorem 5.2. Suppose that the assumptions of Theorem 5.1 are satisfied and
that c
\Deltas is evaluated in floating-point arithmetic from the formula (16b), with c
\Deltax
replacing \Deltax. We then have
(39a)
\Deltas
\Deltas
Proof. Standard roundoff error analysis applied to (16b) shows that
c
\Deltaxj
By differencing (16b) and (40), we obtain
\Deltaxj
we have from (11) that
By combining these estimates with (32) and (41), we obtain the desired result (39a).
For we have from (11) again that
while from (31) and (32) we have
By substituting in (41), we obtain (39b).
The last two results show that the requirements of Theorem 4.1 are satisfied, so
that the algorithm can make significant progress along these search directions. We
summarize the combination of Theorems 4.1, 5.1, and 5.2 as a corollary.
Corollary 5.3. Suppose that Assumption 1 holds and that the problem is non-
degenerate, that is moderate. Suppose that the procedure for solving
Condition 1. If the approximate step is computed with oe 2 [0; 1=2], then
for all sufficiently small -, the formulae (19), (20), and (21) are satisfied.
6. The Bunch-Kaufman Factorization. We show in this section that a procedure
for solving (16a) based on the Bunch-Kaufman factorization satisfies Condition
1, so that the conclusion of Corollary 5.3 applies. Since much of the analysis of this
section can be reused in the analysis of the Bunch-Parlett and sparse Bunch-Parlett
algorithms, we give the details here and refer to them in later sections.
It is sufficient to describe just the first stage of the procedure. Later stages apply
the same technique recursively to the remaining submatrix.
The pivot selection procedure for Bunch-Kaufman [2] is as follows.
Choose find r such that -
else
1 \Theta 1 pivot; choose P 1 so that (P 1 TP T
else
2 \Theta 2 pivot; choose P 1 so that (P 1 TP T
end.
If we denote the 1 \Theta 1 or 2 \Theta 2 pivot block by E and write
the first step of the factorization yields
I
I
C. The algorithm continues by applying this procedure to
T . Note that the - i are generally changed by each stage of the factorization. The
contains the subdiagonals in the first one or two columns of the L
factor.
Bunch and Kaufman [2] show that for the particular choice
17)=8, we
have
so there is a modest bound on element growth during each stage of the factorization.
When applied to canonical matrices, the Bunch-Kaufman procedure selects pivots
of specific types and produces a reduced submatrix that is also canonical. We state
these results in the following two theorems, whose proofs are tedious and are relegated
to the

Appendix

.
Theorem 6.1. Let one step of the Bunch-Kaufman factorization be applied to a
canonical matrix that is not degenerate. Then
(a) The pivot block E will be either
(i) a 1 \Theta 1 block, chosen from among the diagonal elements of  ; or
(ii) a 2 \Theta 2 block, in which the off-diagonal element E 12 is one of the elements
of B;
(b) The matrix remaining after the elimination is canonical, and the absolute
change in the elements of   is at most O(1);
(c) Using the notation from (42), we have that
pivot.
Theorem 6.2. Let one step of the Bunch-Kaufman factorization be applied to a
degenerate canonical matrix. Then
(a) The pivot block E will be either
(i) a 1 \Theta 1 block, chosen from any of the diagonals (large or small); or
(ii) a 2 \Theta 2 block, in which all the elements are O(-
(b) The matrix remaining after the elimination is canonical (not necessarily de-
generate), and the absolute change to the remaining matrix is O(-
Because of Assumption 1, our initial matrix in (16a) is canonical. Barring pathological
growth in the remaining submatrices, one of Theorems 6.1 and 6.2 applies at
every stage of the Bunch-Kaufman factorization.
If B is square in the original matrix (corresponding to a nondegenerate linear
program), then the remaining matrices encountered at every stage of the factorization
are not degenerate. After a 1 \Theta 1 pivot, the dimensions of B are unchanged, while a
2 \Theta 2 pivot shrinks B by exactly one row and column, so it remains square. When a
pivot causes B to disappear altogether, the reduced matrix has the form
u). It follows that in the case of square B, Theorem 6.1 is sufficient to analyze the
entire factorization. The following result gives the backward error analysis for the
factorization in this case.
Corollary 6.3. Let the Bunch-Kaufman factorization be applied to a canonical
matrix T in which B is square. Then, for all sufficiently small -, we obtain computed
factors -
L and -
D such that
\Delta is a canonical error matrix associated with T .
Proof. We prove the result by an induction argument on the dimension -
of the matrix T . The induction is made slightly more complex than usual by the form
of the canonical matrix, notably, the presence of the square matrix B of dimension
For
Therefore
(45) holds with -
For -
we have two cases there are two
elements of
on the diagonal, while the off-diagonals are O(- +u).
Hence, a 1 \Theta 1 pivot is chosen. If there is no pivoting, the first step of elimination
yields
L has unit diagonals, we obtain by expanding the factors that
\Delta;
where
so -
\Delta is a canonical error matrix associated with T . The same logic applies if pivoting
occurs.
In the remaining case the pivot is 2 \Theta 2, we have -
holds trivially with -
We now examine a canonical matrix of dimension - n ? 2 in which B is square,
and examine the first stage of the factorization. Because the matrix is canonical and
nondegenerate, Theorem 6.1 applies. For some permutation matrix P 1 , we have from
(42) and (43) that the first stage yields partial factors -
I 0
where
Note that \Delta D is a canonical error matrix corresponding to -
. By the proof of Theorem
6.1, the (2; 2) submatrix of -
D 1 is canonical, so we use the inductive hypothesis to
deduce that the -
L, -
D factors of this submatrix satisfy
for some permutation matrix P 2 and some canonical error matrix -
corresponding
to ( -
We compose the overall factors of T as follows:
I 0
I 0
where
and so
By substituting (47) and (46) into (48), we obtain
\Theta -

where
are canonical error matrices
corresponding to -
T , we have
Hence, -
\Delta is a canonical error matrix corresponding to T .
We complete the proof by noting that Theorem 6.1 can be applied to the remaining
because it is also canonical and nondegenerate.
Given the system T z = d and the data P , -
L, and -
D from the factorization,
the computed solution - z is found by performing two vector permutations with P ,
triangular substitutions with -
L and -
L T , and a blockwise inversion of -
D. The 2 \Theta 2
diagonal blocks in -
D can be handled by the Gaussian elimination procedure outlined
in the following technical lemma, which is proved in Appendix A.3. It is easy to show
that the elements of the pivot block E satisfy the condition (49).
Lemma 6.4. Consider the 2 \Theta 2 linear system in which E is symmetric
with
for some ffi 2 (0; 1). Then if we compute the solution by applying Gaussian elimination
to the permuted system
then the computed solution - y satisfies
where
The additional error that is introduced during recovery of the solution with the
computed factors -
L, -
D, and -
L T is quantified in the next result.
Lemma 6.5. Suppose the assumptions and notation of Corollary 6.3 hold. Then
the computed solution -
z to the system -
\DeltaP
\Delta is a canonical error matrix associated with T .
Proof. From standard results for triangular substitution, the computed solution
of -
A similar result holds for triangular substitution with the transpose -
For solution of -
z a , we note that -
D is block-diagonal with 1 \Theta 1 and 2 \Theta 2
blocks. For the 2 \Theta 2 pivot blocks that arise in the Bunch-Kaufman procedure, the
assumptions of Lemma 6.4 hold, so the computed solution -
y of a 2 \Theta 2 subsystem
When E is a 1 \Theta 1 block, the estimate (53) holds trivially. Hence, the computed
solution -
z b of -
z a
By combining the error expressions for the three component systems, we find that
our computed solution - z satisfies
Multiplying the matrix products, we find that (52) is satisfied with
From our earlier discussions on the composition of -
L and -
D, it is easy to see that the
absolute matrix product j -
Lj T contains all O(1) elements, except for the large
diagonals, which occur in the same positions as in PTP T . Hence P -
\DeltaP T is a canonical
error matrix corresponding to PTP T , and our proof is complete.
We can now summarize the effects of roundoff error on the entire solution process
for (16) in the following theorem.
Theorem 6.6. Suppose T is a canonical matrix in which B is square. Then,
for all sufficiently small -, the Bunch-Kaufman factorization followed by the solution
process outlined above satisfies Condition 1.
Proof. As we noted immediately following Condition 1, the actual right-hand side
may differ by terms of O(u) from its "theoretical" value d. From (52), the computed
solution -
z to T z = d satisfies
\DeltaP
d;
Substituting from (45), we obtain
\DeltaP
d;
so Condition 1 follows when we set
\Delta.
We have shown that in the case of a nondegenerate linear program, the procedure
based on applying Bunch-Kaufman to (16a) leads to approximate steps ( c
\Delta-; c
\Deltax; c
\Deltas)
that satisfy the conditions of Theorem 4.1. The estimate (20) implies that during
the final iterations of a primal-dual algorithm, near-unit steps can be taken along
these directions without leaving the nonnegative orthant. Moreover, if the centering
parameter oe is small or zero, a large reduction in the duality gap - can be expected.
In the extreme case linear convergence with a rate
constant of O(u) can be attained if the actual step length is close to -
ff   . Most practical
algorithms choose the step length to be a fixed multiple - typically :95 or :9995 -
of -
ff   , and indeed these methods often converge rapidly during their final stages. For
algorithms that use a more theoretically justifiable definition of step length the story
is not, unfortunately, this simple. In [22, Section 4], for instance, extra restrictions
are applied to ff to ensure that (12) and (14) continue to hold at the next iterate.
These restrictions may result in ff being much smaller than one. This case is analyzed
in [22, Section 4], so we do not repeat it here.
Finally, we note that the lower triangle L produced by the Bunch-Kaufman factorization
may contain elements that are much larger than those of the original matrix
T . This phenomenon has been closely scrutinized in a recent report by Ashcraft,
Grimes, and Lewis [1], who observe that it leads to convergence difficulties in a non-linear
programming code. In the context of our canonical matrix of Theorem 6.1,
this blowup problem does not occur. As we show in part (c) of the theorem, the
contribution CE \Gamma1 made by one step of Bunch-Kaufman is either O(-) or O(1). The
blowup problem may occur, however, when we have a degenerate canonical matrix as
in Theorem 6.2. We only have to deal with matrices like this when the linear program
itself is degenerate, and in this case there are other, more serious difficulties to face,
as we discuss in Section 9.
7. The Bunch-Parlett Factorization. The Bunch-Parlett searches the entire
remaining matrix for each pivot, not just one or two columns. The pivot selection
procedure is as follows.
Choose rs
choose P 1 so that (P 1 TP T
else
choose P 1 so that (P 1 TP T
rs
end.
The elimination step is identical to Bunch-Kaufman, and the process of using the
LDL T factorization to solve the system T z = d is the same as in the preceding
section. As in Bunch-Kaufman, the value
17)=8 leads to the modest bound
of 2:57 on element growth at each stage.
When applied to canonical matrices, the Bunch-Parlett factorization proceeds in
three stages:
1. All the diagonal elements of   are selected as 1 \Theta 1 pivots;
2. 2 \Theta 2 pivots of the type described in Theorem 6.1(a) are chosen;
3. When no more 2 \Theta 2 pivots like this are available and the remaining matrix
contains only elements of size O(- + u), a combination of small 1 \Theta 1 and
pivots is used to complete the factorization process.
We prove this assertion in the following lemma.
Theorem 7.1. Suppose that the Bunch-Parlett procedure is applied to a canonical
matrix. Then the factorization proceeds according to the three-stage outline above. If
the canonical matrix has B square and is nonvacuous, the factorization is completed
by stages 1 and 2; stage 3 is vacuous.
Proof. Assuming that   is not vacuous, we have at the pivot selection step that
The pivot element will therefore be one of the
large diagonals corresponding to  . The remaining matrix is updated by subtracting
O(-). Hence, the remaining matrix
retains canonical form.
We can apply this argument inductively until all the diagonals in   are exhausted.
At the end of stage 1, the remaining matrix has the form
Stage 2 now begins. If B is not vacuous, we have -
In fact, by the assumption B
=\Omega\Gamma245 we have - off
=\Omega\Gamma338 and the element T rs that
achieves the maximum comes from B. The 2 \Theta 2 block with off-diagonal element T rs
is selected as the pivot. After the elimination step, the size of B is reduced by one
row and column. The proof of Theorem 6.1(b) can be applied again here to show that
the remaining matrix is also canonical, so 2 \Theta 2 pivots of this type will continue to be
selected until B vanishes.
The number of steps in stage 2 is min(rows(B); columns(B)). At the end of this
stage, the remaining matrix is square with dimension
all its elements have size O(- + u). In stage 3, both 1 \Theta 1 and 2 \Theta 2 pivots may be
used to factor this matrix. If B is square, the factorization is complete after stage 2.
The other major results of Section 6 continue to hold when the Bunch-Parlett
algorithm is used instead of Bunch-Kaufman; only trivial adjustments to the analysis
in Section 6 and Appendix A.1 are necessary. We summarize the conclusions in the
following theorem.
Theorem 7.2. Suppose T is a canonical matrix in which B is square. Then,
for all sufficiently small -, the Bunch-Parlett factorization followed by the solution
process outlined in Section 6 satisfies Condition 1.
8. The Sparse Bunch-Parlett Factorization. Several authors (notably Fourer
and Mehrotra [5]) have proposed a sparse variant of the Bunch-Parlett factorization
that compromises between maintaining sparsity and limiting element growth in the
remaining matrix. We outline the pivot selection procedure as described by [5], with
a slight modification noted below.
For each index we define the degree n i to be the number of off-diagonal
nonzeros in row i. We also define an estimate of the joint nonzero content
of rows i and j by
is termed oxo if both of T ii and T jj are zero, tile if one of T ii and T jj is zero, and full
if both of T ii and T jj are nonzero. We define a cost associated with using (55) as the
pivot block in each of these three cases by
tile: (n
The cost is an estimate of the fill-in associated with using (55) as the pivot block.
For prospective pivots, we define stability criteria in terms of the usual constant
and the off-diagonal norms - i defined in (29). Any 1 \Theta 1 pivot must satisfy
ii
while a 2 \Theta 2 pivot (55) must have
The pivot selection procedure is as follows.
for
for i with
consider T ii with degree
if any of these elements
accept as a 1 \Theta 1 pivot and exit;
else label it as unstable;
for unstable pivots T ii from the previous loop
consider 2 \Theta 2 pivots involving T ii , with costs at most
for oxo, tile, and full pivots, respectively;
if any of these blocks satisfy (57)
accept as a 2 \Theta 2 pivot and exit;
end.
The pivot selection pattern for the sparse Bunch-Parlett algorithm is essentially
the same as for the Bunch-Kaufman algorithm, as described in Theorems 6.1 and
6.2. We prove this result in the appendix, since the analysis differs a little from the
Bunch-Kaufman case.
Theorem 8.1. The results of Theorems 6.1 and 6.2 hold when the sparse Bunch-
Parlett factorization is used in place of the Bunch-Kaufman procedure. To obtain this
result, we modified the acceptance condition (56) for 1 \Theta 1 pivots. In the description
of [5], the right-hand side is 1=ffi rather than 2=ffi. With the original choice, the
sparse Bunch-Parlett algorithm applied to a degenerate canonical matrix could allow
another type of pivot: a 2 \Theta 2 pivot in which one diagonal is from   and the other
has size O(- + u). A pivot of this type is poorly conditioned and will generally lead
to instability during the blockwise inversion of -
D.
The other major results of Section 6 also continue to hold when the sparse Bunch-
Parlett algorithm is used in place of Bunch-Kaufman. We summarize the conclusions
in the following theorem.
Theorem 8.2. Suppose T is a canonical matrix in which B is square. Then, for
all sufficiently small -, the sparse Bunch-Parlett factorization followed by the solution
process outlined in Section 6 satisfies Condition 1.
9. The Degenerate Case. When the linear program (1), (2) is degenerate -
the three factorization procedures can no longer run to completion with
just the two kinds of pivots described in Theorem 6.1. The nonsquare shape of B in
the matrix (34) means that pivots of size O(- + u) - either 1 \Theta 1 or 2 \Theta 2 - are
used at some point in the factorization process. The factorizations fail only if these
pivots are exactly zero, which happens often on small problems but not otherwise.
The more common outcome is that the interior-point algorithm makes only slow or
erratic progress after - has achieved a certain (small) value. In this section we sketch
the reasons for this outcome.
In all the factorizations above, the large diagonal elements in X \Gamma1
N SN are used
as 1 \Theta 1 pivots. Even though these pivots are not necessarily used before any others
(except in the Bunch-Parlett algorithm), the factorizations behave as if they were
solving the system (16) in the equivalent, partitioned form
\Delta-
\Theta
The coefficient matrix in (58) is an O(-) perturbation of the matrix
Since B is well conditioned by Definition 2, the matrix in (60) has 2 min(jBj; m)
nonzero singular values of
magnitude\Omega\Gammagni In the nondegenerate case, (60) is well
conditioned. Otherwise, it has jm \Gamma jBjj zero singular values. When jBj ! m, the null
space of (60) is spanned by
Z-
Z is an m \Theta (m \Gamma jBj) matrix of full rank such that B T -
the null space of (60) is spanned by the matrix
--
Z
Z spans the null space of B. For small -, these null spaces are not altered
much by the perturbation of size O(-) that is present in the matrix (58), because
the nonzero singular values of (60) are well separated from zero. Perturbations in
the solution of (58) due to roundoff will occur mainly in the space of small singular
values. Hence, when jBj ! m, the perturbations occur mostly in the range space of
the matrix (61), that is, in the components of \Delta-. Similarly, when jBj ? m, the
perturbations occur in the range space of the matrix (61), that is, in the components
of \Deltax B .
The main source of difficulty is inaccuracy in the computed residual vectors r b
and r c which, as mentioned above, contain errors of O(u). In the case jBj ? m,
these perturbations are magnified by the inverses of the small singular values, usually
leading to errors of size about O(u=-) in the components of \Deltax B . The large relative
errors in \Deltax B induce large relative errors in \Deltas B through the formula (16b). The
step length to the boundary ff   may therefore be sharply curtailed because of the
nonnegativity requirements (17a). In the case jBj ! m, the large relative errors in \Delta-
induce errors in \Deltax N through the formula (59), while in turn induce large relative
errors in \Deltas N through (16b). The step length may again be curtailed as a result.
Errors from sources other than the vector r are less significant.
If we have a strictly feasible starting point (see (13)), then we can simply set
throughout the algorithm. In this case, we can fix r at zero in the computations and
avoid the problem above. It is usually not easy to find such a starting point, however,
so some thought should be given to other ways of dealing with the problem.
One option is to simply terminate the algorithm when it stalls, declaring success
if both - and r are small. This option works well for most purposes, since stalling
usually occurs only after - is reduced to O(u), by which time the problem has usually
converged to acceptable accuracy. Fourer and Mehrotra [5] report that the convergence
criteria are usually satisfied before the ill effects of roundoff are seen. Our
testing in Section 10 allows a similar conclusion.
A second option is to switch to a termination procedure when the interior-point
algorithm stalls. A finite termination procedure (see, for example, Ye [24]) or crossover
to the simplex method (Meggido [13]) could be activated.
A third option is simply to fix r at zero in the computations once it has reached
the O(u) level, because at this stage our current point is feasible to within the limits of
floating-point arithmetic. By doing so, we are effectively introducing a perturbation
into the problem to freeze the infeasibility at its current level. This perturbation has
an interesting effect: It moves the solution to a particular vertex of the previously
optimal face, changing the B [ N partition appropriately. If we continue to run the
interior-point algorithm to higher accuracy, it eventually converges to this vertex,
but only after going through many more iterates (and taking some sharp turns in
the process). The result of this process is similar to what we would achieve with a
crossover to simplex, but the computational cost would generally be much higher.
10. Computational Experiments. We report here on some computational experiments
that demonstrate the effects described above. Our testbed algorithm is the
path-following algorithm described in Wright [21]. In exact
arithmetic, this algorithm achieves superlinear convergence because it eventually always
takes affine-scaling steps in (5)) with step length ff approaching 1. This
algorithm performs well on practical problems, but is not as fast as codes that use the
Mehrotra predictor-corrector heuristic, for which no solid convergence theory exists,
except in the nondegenerate case. The asymptotic behavior in finite precision is quite
similar for the two algorithms.
To show that the finite precision effects are not confined to "nice" problems, we
generate problems with fairly wide variations in the components of A, x
N .
The matrix A is dense and random, with elements defined by
where every instance of - is selected from a uniform distribution on the interval [0; 1].
(We choose all the elements in the first row of A to be positive to ensure that the
feasible region is bounded.) We control the size of the index sets B and N (to control
the amount of degeneracy) and set
We choose a particular solution (- setting
s
where each - is as before. The vectors b and c are determined by the choices of A
and (-
The LAPACK Bunch-Kaufman factorization routines dsytrf and dsytrs are used
to solve (16a). These routines (and the rest of our code) use double-precision arith-
metic, giving on the SPARC-5 on which these results were obtained.
We report on problems with problems smaller than this, exactly
zero pivots often occur in degenerate cases, leading to breakdown.) Termination
occurs an artificially stringent criterion, chosen to give us a clear
look at asymptotic effects.
The first result is for a nondegenerate problem, for which
1 shows the sizes of - and krk on each iterate. For the reasons that we outlined
immediately following Condition 1, krk stabilizes at a magnitude of O(u). The duality
gap - does not converge subquadratically (as it would in exact arithmetic) but rather
exhibits extremely fast linear convergence, with a rate constant of about 10 \Gamma10 . This
is exactly the effect predicted by formula (21) for the affine-scaling steps that are
taken on the last four iterations.
To see that the pivots have the properties predicted by Theorem 6.1, we examine
the matrix D from the Bunch-Kaufman factorization. Table 2 shows D at iteration
As expected, there are six 1 \Theta 1 pivots of
pivots in which the diagonals are tiny and the off-diagonals
are\Omega\Gammae/4 The
same structure is present in D at every iteration after iteration 15.
Our second example is for a dual degenerate problem with As
can be seen from Table 3, the algorithm achieves fairly high accuracy after about 20
iterations, but no further improvement can be made after that point. The behavior
is consistent with the discussion of Section 9. It suggests that the results of Section 6
are "tight," in that we cannot prove that "useful" search directions are obtained for
arbitrarily small -.
Examination of the D factor for the second example (Table 4) shows that the
pivot pattern is in line with the predictions of Theorems 6.1 and 6.2. Together, these
results imply that there are exactly min(m; jBj) of the stable 2 \Theta 2 pivots with an
3 4.3 1.6
5 3.1 -12.0
19 -22.1 -13.8
20 -33.3 -14.2 termination

Table
The D factor at iteration 17 of the nondegenerate test problem (
Row/Column Pivot Block

Table
Dual degenerate problem:
19 -6.0 -13.8
22 -14.8 -13.8
The D factor at iteration 17 of the degenerate test problem with
magnitude less than
Row/Column Pivot Block
9
off-diagonal from B, and jN of the large 1 \Theta 1 pivots. Together, these
stable pivots account for
stages of the factorization, so unstable pivots are used on the remaining submatrix
whose dimension is jm \Gamma jBjj. In Table 4, we see that the last two 1 \Theta 1 pivots are
unstable, as expected. As we described in the first part of Section 9, the errors in
c
\Deltax B and c
\Deltas B are preventing further progress. On iteration 100, the computed affine
step has k c
its exact counterpart would have k\Deltax B k1 = O(-).
By comparing components of c
\Deltas B with s B , we find that the step to the boundary
is sharply curtailed by the restriction s
\Deltas (23)). The remaining
components of the step do not contain deleterious errors; we have
\Deltas
Finally, we consider a primal degenerate problem with m. The iteration
schedule in Table 5 shows similar behavior to the dual degenerate problem. The D
factor from iteration 100 is shown in Table 6. All pivots are stable except for the last
two 1 \Theta 1 blocks, which again matches the prediction (63). As discussed in Section 9,
the deleterious errors occur in the subvector c
\Delta-, so errors are induced in c
\Deltas N and
c
\Deltax N through formulas 59 and (16b). On iteration 100, we have k c
\Deltas for the affine scaling step. The components c
\Deltax B and c
\Deltas B are not
affected; their 1-norms are :17(\Gamma18) and :51(\Gamma12), respectively.
Primal degenerate problem:
100 -17.6 -14.0

Table
The D factor at iteration 17 of the degenerate test problem with
Row/Column Pivot Block
14
A. Proofs of Theorems from Sections 6 and 8.
A.1. Proof of Theorem 6.1. We prove (a) by systematically excluding the
other possible choices for pivots:
(iii) The pivot is 1 \Theta 1 and is a diagonal element from either the (1; 1) or (2; 2)
blocks of the canonical matrix. Inspection of the Bunch-Kaufman algorithm
shows that T 11 is chosen as pivot if either
r
Now, since - r is the maximum off-diagonal in some column of (26), we have
O(1), while since T 11 comes from either the (1; 1) or (2; 2) block of (26),
we have jT fixed, we have from (64) that
is the magnitude of the largest off-diagonal in some row/column of
(26), we have that - 1 is the 1-norm of some row or column of B. But (65)
is incompatible with Hence jT 11 j from the (1; 1)
or (2; 2) blocks cannot be used as a pivot.
A similar argument holds when T rr is chosen as pivot, where T rr is one of the
small diagonals.
(iv) The pivot is 2 \Theta 2 and involves at least one element from  . Since all the
off-diagonals in (26) are O(1), the quantities - i , are all O(1).
pivot with diagonal elements T 11 and T rr must have
which implies that T 11 and T rr are both O(1). Since all the diagonals of
cannot be candidates for T 11 and T rr .
(v) The pivot is 2 \Theta 2, and the pivot block is drawn either entirely from the (1; 1)
block of (26) or entirely from the (2; 2) block. In this case, T 1r - the element
for which jT 1r 1r has the largest magnitude in
its column of (26), and since its column includes either a row or column of
B, we have that one of the rows or columns of B is O(- u). As in (iii), we
have a contradiction, since this estimate is incompatible with B
=\Omega\Gamma43 and
This completes the proof of part (a).
We turn to (b), examining the effects of one step of elimination performed with
pivot selection corresponding to the two cases (i) and (ii). For (i), suppose the (i; i)
element of   is chosen as the pivot. After symmetric permutation of the canonical
matrix, to place the pivot in the (1; 1) position, we obtain
-~
-~
where ~
P is some permutation matrix, N \Deltai denotes the i-th column of N , ~
N is obtained
from N by removing N \Deltai , and ~
is obtained from   by removing its i-th row and
column. Since j( +O(-+u)) \Gamma1
O(-), the submatrix that remains after elimination
is
~
~
ii
\Theta
~
~
It is easy to see that (66) is canonical, so our result is proved for case (i).
For case (ii), the proof is a little messier. Suppose the diagonals of the 2 \Theta 2
pivot are the (i; i) element of E 1 and the (j; j) element of E 2 . After symmetric
rearrangement to put this pivot in the upper left corner, (26) becomes
I
\Deltaj
I
where
P is some permutation matrix;
is the i-th row of N ;
N is N with N i\Delta removed;
is the i-th row of B, with its j-th element removed;
;i is the j-th column of B, with its i-th element removed;
its i-th and j-th column removed.
By the choice of B ij , it is either the largest element in its row or the largest element in
its column of B. From our assumptions on B, we deduce that jB ij j
Denoting
the pivot block by E, we have
Therefore the elimination step yields the remaining matrix
\Deltaj
where
It is obvious that (68) satisfies Definition 2, except possibly for the conditioning of
the remaining matrix -
B. This matrix is obtained by pivoting the (i; j) element of B
to the (1; 1) position and then doing one step of Gaussian elimination. In fact, we are
doing partial pivoting since, as noted above, B ij is the largest element in either its
row or its column. Hence, the conditioning of the reduced submatrix -
B is unlikely to
differ much from -(B), so it is reasonable to assert that -
We have shown that our stated result holds for both cases (i) and (ii), so our
proof of part (b) is complete.
For part (c), note that whether the pivot block is 1 \Theta 1 or 2 \Theta 2. For
pivots, we have pivots, we have from
=\Omega\Gamma17 that
A.2. Proof of Theorem 6.2. Again, we prove (a) by excluding the other possible
choice for a pivot:
(iii) The pivot is 2 \Theta 2 and contains at least one element from  . In a degenerate
canonical matrix, we have - pivot with
diagonal elements T 11 and T rr must have
which implies that both diagonals are O(-+u), so neither element can come
from  .
In the case of either a 1 \Theta 1 or 2 \Theta 2 pivot made up of elements of size O(-+u), we
can use the standard argument about element growth in Bunch-Kaufman (that is, the
argument that leads to (44)) to deduce the result (b). In the remaining case, where
the pivot is a single diagonal element from  , we have in the notation of (42) that
Hence, the update to the remaining submatrix is
bounded by
which certainly has size O(- + u).
A.3. Proof of Lemma 6.4. Proof. In floating-point arithmetic, the LU factorization
of (50) yields the following approximate LU factors:
where
It is well known that for triangular substitution applied to any triangular system
h, the computed solution - z satisfies (U h, where jEU
applying this observation to each of the matrices in (69), we find that the computed
solution -
y of (50) satisfies
where
By multiplying out the coefficient matrix in (70), we obtain
where
(The last equality follows from (49).) Hence, (71) can be written as
satisfies the bound (51).
A.4. Proof of Theorem 8.1. Proof. We start by proving the analog of Theorem
6.1(a). As in the earlier proof, we systematically exclude the three other possible
choices of pivots.
(iii) The pivot is 1 \Theta 1 and is a diagonal element from either the (1; 1) or (2; 2)
blocks of (26). Then this pivot (T ii , say) will be O(- u). According to the
stability criterion (56) we then have which implies that one of
the rows or columns of B is O(-+u). However, this estimate is incompatible
with
O(1), so this kind of pivot cannot occur.
(iv) The pivot is 2 \Theta 2 and involves at least one diagonal element from  . First,
we show that we cannot have both diagonals from  . If this were the case,
then at least one of these diagonals (T ii , say) would have been considered as
a 1 \Theta 1 pivot at an earlier point in the algorithm. But if it was considered, it
would have been accepted, since
ii
for sufficiently small -. Hence, at most one of the diagonals is from  .
Without loss of generality, suppose in (57) that T ii is from   while the remaining
diagonal T jj is O(- + u). In fact, we have
and so fi fi fi fi fi
Hence, from (57), we have
From the second row of this inequality, we have
is the 1-norm of one of the rows or columns of B, so this estimate
contradicts our assumptions on B. Hence, this type of pivot cannot occur.
(v) The pivot is 2 \Theta 2, and the pivot block E is drawn either entirely from the
of (26) or entirely from the (2; 2) block. In this case, all elements
of E are O(- + u). From (57), we have as above that
Taking the second row of this relation, we obtain
where, by definition, - i and - j are both nonnegative. Consider two cases.
we have from (72) that
For the reasons outlined earlier, the assumptions on B are inconsistent with
this bound on - i , so this case cannot hold. For the other case jT ij
we have
which is also disallowed by our assumptions. Hence, pivots of this type cannot
occur.
The proof of the remaining parts (b) and (c) of Theorem 6.1 is identical in this
case.
Turning now to the case of a degenerate canonical matrix and the analog of
Theorem 6.2, we start by showing that no 2 \Theta 2 pivots may contain diagonal elements
from  .
Note that for a degenerate matrix, the off-diagonals, and hence the quantities - i ,
all have size O(- + u). If the pivot is a 2 \Theta 2 block in which both diagonals are from
, then one of them (T ii , say) must have been considered previously as a 1 \Theta 1 pivot.
But if it was considered, it would have been accepted, since
ii
Hence, this type of pivot cannot occur.
If just one of the diagonals is from  , this diagonal element (T jj , say) must not
have been considered earlier as a 1 \Theta 1 pivot, since then it would have been accepted
for the reason described above. Hence, the other pivot T ii , which has size O(- + u),
must have been considered as a 1 \Theta 1 pivot and rejected. Because of (56), T ii must
On the other hand, since the 2 \Theta 2 pivot is accepted, we must have
Consider first the case of T 2
j. Then from the first block row in (74), this
inequality implies that
which contradicts our assumption that T jj has
remaining case has
From (74) and (73), we have
which is a contradiction. Hence this kind of pivot - in which exactly one of the
diagonals comes from   - cannot occur either, and we are done.
For the analog of part (b) of Theorem 6.2, we have from (56) and (57) and the
definition of C and E in (42) that
Hence, the update matrix CE \Gamma1 C T is bounded as follows:
giving the result.

Acknowledgments

. I thank the editor, Linda Kaufman, for the care she took
with this paper, and two anonymous referees for their meticulous reports which improved
both presentation and content. I also thank John Lewis for interesting discussions
about the Bunch-Kaufman factorization during his visit to Argonne in April,
1995.



--R

Accurate symmetric indefinite linear equation solvers.
Some stable methods for calculating inertia and solving symmetric linear systems
The solution of augmented systems
Stability of symmetric ill-conditioned systems arising in interior methods for constrained optimization
Solving symmetric indefinite systems in an interior-point method for linear programming
Theory of linear programming


An O( p nL) iteration potential reduction algorithm for linear complementarity problems
Computational experience with a primal-dual interior point method for linear programming


On finding primal- and dual-optimal bases
On the implementation of a primal-dual interior point method
On adaptive step primal-dual interior-point algorithms for linear programming

Potential reduction methods in mathematical programming
Technical Report SOR 92-5
Stable numerical algorithms for equilibrium systems



A simplified homogeneous and self-dual linear programming algorithm and its implementation
On the finite convergence of interior-point algorithms for linear programming
On the convergence of a class of infeasible-interior-point methods for the horizontal linear complementarity problem
--TR

--CTR
Stephen J. Wright, Superlinear Convergence of a Stabilized SQP Method to a Degenerate Solution, Computational Optimization and Applications, v.11 n.3, p.253-275, Dec. 1998
Csaba Mszros, Detecting "dense" columns in interior point methods for linear programs, Computational Optimization and Applications, v.36 n.2-3, p.309-320, 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
symmetric indefinite matrices;interior-point methods