Stability of Linear Equations Solvers in Interior-Point Methods.
Primal-dual interior-point methods for  linear complementarity and linear programming problems solve a linear system of equations to obtain a  modified Newton step at each iteration. These linear systems become increasingly ill-conditioned in the later stages of the algorithm, but the computed steps are often sufficiently accurate to be useful. We use error analysis techniques tailored to the special structure of these linear systems to explain this observation and examine how theoretically superlinear convergence of a path-following  algorithm is affected by the roundoff errors.
Introduction
. The monotone linear complementarity problem (LCP) is the
problem of finding a vector pair (x; y) 2 R l n \Theta R l n such that
(1)
where M (a real, n \Theta n positive semidefinite matrix) and q (a real vector with n
are given. Note that M need not be symmetric. It is well known that (1)
includes the linear programming problem as a special case. Specifically, for the linear
programming formulation
min z
c T z subject to Az - b; z - 0;
(2)
where A 2 R l m\Thetap , we can introduce the dual variable -
l m for the constraint
and obtain the following necessary and sufficient conditions for optimality of
the primal-dual pair (z; -):
z
c
\Gammab
z
(3a)
z 0:
For appropriate definitions of M and q, (3) has the form (1). Little is lost from either
the practical or theoretical point of view by applying interior-point algorithms for (1)
to the special cases of linear and convex quadratic programming, provided that the
special structure of each problem is exploited in the solution of the linear systems at
each iteration.
Interior-point methods for (1) generate a sequence of iterates are
strictly positive. Many such methods require a linear system of the form
- u
where
This work was based on research supported by the Office of Scientific Computing, U.S. Department
of Energy, under Contract W-31-109-Eng-38.
y Mathematics and Computer Science Division, Argonne National Laboratory, 9700 South Cass
Avenue, Argonne, IL 60439.
to be solved for a search direction (u; v) at each iteration. Affine-scaling methods solve
(4) with to find a search direction, then step a fraction of the distance along this
direction to the boundary of the nonnegative orthant defined by (x; y) - 0. Affine-
scaling steps (u; v) are simply Newton steps for the system of nonlinear equations
Path-following methods (see, for example, Monteiro and Adler [10], Zhang [19], Wright
[14]) generate steps by using generally positive values of oe in (4). (As we see later, the
algorithm of [14] allows some iterates in an attempt to attain the rapid local
convergence associated with Newton's method.) Potential-reduction methods (see,
for example, Kojima, Mizuno, and Yoshise [6], Kojima, Kurita, and Mizuno [5]) also
determine search directions by solving systems like (4), but they refer to a logarithmic
potential function to decide how far to move along the computed direction. Predictor-corrector
methods (see, for example, Ye and Anstreicher [18], Ji, Potra, and Huang
[4], Potra [12]) take steps with either
The system (4) is highly structured; since the diagonals of X and Y are strictly
positive, we can rearrange the system to obtain
(6a)
e:
In the case of linear programming (2), equation (6a) contains even more structure.
Its form is
A
u-
where Y z and Y - are positive diagonal matrices and
The matrix in (7) can be made symmetric indefinite by multiplying the first block
row by \Gamma1. This system can be reduced even further by eliminating either u z or u - .
For instance, if u z is eliminated, we obtain
z
z
Some interior-point codes for linear programming use the formulation (8), with
modifications for handling dense columns in A and for dealing with non-standard linear
programming formulations (see Lustig, Marsten, and Shanno [7, 8] and Xu, Hung,
and Ye [17]). Other codes, notably those of Fourer and Mehrotra [1] and Vanderbei
[13], handle the formulation (7). Analysis of algorithms for these formulations are
discussed in another preprint [16]. In this paper, we focus on the system arising from
general monotone LCP (6), and analyze the behavior of Gaussian elimination with
pivoting applied to this system.
Since for any index least one of x i and y i is zero at the solution, we
would expect some of the diagonal elements of X \Gamma1 Y to approach zero and some to
approach +1 as the solution set is approached. Hence the coefficient matrix in (6a)
tends to become increasingly ill-conditioned during the later stages of the algorithm.
From the standard error analysis of linear systems, we might therefore expect that
LINEAR EQUATIONS IN INTERIOR-POINT METHODS 3
rounding errors in the step (u; v) make it useless in advancing the algorithm towards
convergence. In this paper, we show that while the theoretical superlinear properties
suggested by the exact analysis are not generally observed, implementations of the
algorithms can still exhibit rapid convergence if the parameters are set to appropriate
values. For a particular path-following infeasible-interior-point algorithm with strong
theoretical convergence properties, these conclusions are presented in Section 4 and
confirmed by computational experiments in Section 5. Section 3 lays the groundwork
by deriving bounds on the rounding errors in the computed values of (u; v), for a wide
class of algorithms that includes the algorithm of Sections 4 and 5. Section 2 presents
the assumptions and a fundamental result from error analysis.
Linear systems that arise in logarithmic barrier methods for constrained optimization
methods were analyzed by Poncele'on [11]. The Newton equations for each
logarithmic subproblem are similar to (6a) in that the large elements occur only on the
diagonal. Despite the apparent ill-conditioning of these systems, Poncele'on showed
that their sensitivity to structured perturbations from a certain class is governed by
the conditioning of the underlying problem and does not depend on the current value
of the barrier parameter. Poncele'on's analysis is somewhat different from that of
Section 3 - she looks at the relative error in the components of the solution, rather
than starting with the absolute error - but her conclusions are consistent with those
obtained in Section 3.
In subsequent sections, subscripts denote components of a vector, while iteration
indices (usually appear as subscripts on scalars and as superscripts on vectors and
matrices. The sets B and N form a partition of the index set f1; ng defined in
Assumption 1 below. If x 2 R l
and so on. For the matrix M 2 R l n\Thetan , we have
and similarly for MBB , MNB , and MNN . Given any matrix we define
and denote the j-th column of H by H \Delta;j . The notation k \Delta k denotes
the 1-, 2-, or 1-norm of a vector or matrix, while -(\Delta) denotes the corresponding
condition number.
For any two nonnegative numbers - and -, we there is a
moderate constant - such that -. When W is a matrix or vector, we write
O(-) to denote O(-). We use to indicate that both
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 representation of any real number z, we have
(See Golub and Van Loan [2, Section 2.4.2].) We assume throughout that u is small
enough that O(u) - 1, where O(\Delta) is the order notation defined above.
2. Assumptions and Basic Results. In the remainder of the paper, we focus
on path-following interior-point methods. "Infeasible" variants of these methods are
among the most widely used practical algorithms for linear programming and LCPs
[7, 8] and have also been the subject of extensive theoretical investigations, which
have shown that they can have strong convergence properties under weak assumptions
[19, 15, 14]. The path-following infeasible-interior-point framework stated below also
includes the class of predictor-corrector methods, for appropriate choices of the initial
point and parameters.
Path-following algorithms restrict their iterates to neighborhoods of the
denotes the nonnegative orthant in R l n . All iterates generated by the algorithms
lie in N (fl min ), where fl min 2 (0; 1=2) is a constant. Other quantities needed
to define the general algorithm include
We define the algorithmic framework as follows:
Algorithm PFI
given
choose oe k 2 [0; 1] and find
\GammaI
choose
and
The decrease in kr k k at each iteration is linear - in fact, r
so the last condition in Algorithm PFI is equivalent to
Hence, kr k so the infeasibility is always bounded by a multiple of
the complementarity gap -.
When the initial point is feasible (r predictor-corrector algorithms such as
that of Ji, Potra, and Huang [4] are special cases of Algorithm PFI. This framework
also includes the infeasible-interior-point algorithms of Zhang [19] and Wright [15, 14].
These algorithms choose fl k+1 and oe k so that a step ff k of nontrivial length can always
be taken without violating the required conditions.
In practical implementations of interior-point methods, the framework of Algorithm
PFI is usually modified slightly. In linear programming codes, different step
lengths are usually chosen for the primal and dual components of x, as experience has
shown that this strategy tends to reduce the number of iterations slightly. Moreover,
explicit membership of the neighborhood (10) is usually not enforced.
strategy, for which there is no supporting theory, is to find the largest value of
LINEAR EQUATIONS IN INTERIOR-POINT METHODS 5
ff in [0; 1] that keeps in the nonnegative orthant and then choose
ff k to be a fixed fraction of this length.) The predictor-corrector strategy of Mehrotra
[9], used also in the codes of Lustig, Marsten, and Shanno [8], Vanderbei [13], and
Xu, Hung, and Ye [17], adds extra terms to the lower part of the right-hand side on
"corrector" iterations. Nevertheless, the coefficient matrices used in these practical
algorithms are the same as in (11), and our conclusions about the accuracy of the
computed steps continue to hold, with minor modifications to the analysis of Section
3.
For most of our analysis, we make the following assumptions about the data for
problem (1) and its solution set:
Assumption 1.
a) Problem (1) has a unique solution strict
complementarity holds). We can define an associated partition B; N of the
index set ng such that x
b) The quantities
are all moderate in size.
Assumption 1 implies that the coefficient matrix in (11) approaches a nonsingular
limit, since there are 2n \Theta 2n permutation matrices P and \Pi such that
Y   X
\GammaI MBN MBB 0
and each of the submatrices on the diagonal of (13) is nonsingular.
When the problem (1) is derived from a linear program as in (3), existence of a
solution implies existence of a strictly complementary solution. However, for both this
special case and the general case of M symmetric positive semidefinite, uniqueness
of the solution and well-conditioning of MBB are often not satisfied in practice, so
Assumption 1 is quite strong. As we see, however, this Assumption plays an important
role in showing that the errors in the computed solutions are not disastrous for the
interior-point algorithm, just as well-conditioning of the square coefficient matrix A
in a linear system is needed to ensure that the relative errors in the computed
version of z are not too large. Our computational experience (Section 5) tends to
indicate that Assumption 1 is necessary as well as sufficient for rapid local convergence
of the algorithm.
We make one further assumption on the iterates generated by the basic algorithm.
Assumption 2. The iterates generated by Algorithm PFI satisfy
lim
Of course, it is not necessary to make this assumption for any reasonable instance of
Algorithm PFI, since convergence to a solution should be one of the properties implied
by the particular schemes for choosing oe k , ff k , and fl k . We make this assumption here
merely to divorce the error estimates of the next section from any particular variant
of Algorithm PFI.
By the implicit function theorem, the nonsingularity of the matrix (13), equation
(11), and kr k
6 STEPHEN J. WRIGHT
We also have the following simple result.
Lemma 2.1. There are positive constants C 1 and C 2 such that for all k sufficiently
large, we have
(15a)
Proof. Because of Assumption 2, we can define an index K and positive constants
C 2 such that
Therefore, since
Also,
Therefore (15a) holds with
. The proof of (15b) is
similar.
Finally, we state a result from the roundoff error analysis of Gaussian elimination,
for reference in the next section.
Theorem 2.2. Suppose that the m \Theta m linear system solved by using
Gaussian elimination, possibly with row and/or column pivoting. Let us denote the
row permutation matrix by P , the column permutation matrices by \Pi, the computed
unit lower triangular factor by -
L, and the computed upper triangular factor by -
U .
Then the computed solution - z solves the perturbed system
U
Proof. Follows immediately from Theorem 6.4 of Higham [3].
During Gaussian elimination, the size of the largest element in each column of
the remaining submatrix may grow as multiples of the pivot rows are added to later
rows in the matrix. We quantify this growth by the growth factor ae, defined as the
smallest positive number such that
We then have the following simple corollary of Theorem 2.2.
Corollary 2.3. Let the system be as in Theorem 2.2, and suppose the
pivots -
U jj are chosen so that j -
LINEAR EQUATIONS IN INTERIOR-POINT METHODS 7
3. Error Bounds for the Steps. In this section, we derive estimates for the
difference between the step actually computed by solving (6), which we denote by
and the corresponding exact values, denoted by (u; v). We treat the cases in
which (-u; -
v) is determined by Gaussian elimination with row partial pivoting and with
complete pivoting.
We start with a purely technical result.
Lemma 3.1. Let G be a square matrix partitioned as
G
where G 11 and G 22 are also square. Suppose that G 11 and G 22 \Gamma G
are
nonsingular. Then G is nonsingular and G \Gamma1 has the form
Our first main result concerning components of -
u is the following.
Theorem 3.2. Let -
u be computed by applying Gaussian elimination with row
partial pivoting to (6a), and suppose that the growth factor ae is not too large. Then
for all - sufficiently small, we have
Proof. We assume that - is much smaller than all the quantities in Assumption
1(b), so that - 1. We retain only the lowest-order terms in u and - in the analysis
since, by our assumptions, higher order terms are small enough to be absorbed into
lower-order terms with minor perturbations of the coefficients.
From Theorem 2.2 and (18), we have, by permuting the rows and columns of (6a),
that
uN
where
(20a)
(20c)
(20d)
Now from Assumption 1(b) and Lemma 2.1, we have
and kMBB k, kMBN k, kMNB k, and kMNN k are all O(1). Combining these observations
with (20), we obtain
8 STEPHEN J. WRIGHT
Therefore (19) can be rewritten as
EBN
uB
where
and
e:
If we denote the coefficient matrix in (21) by G, with G
EBB and so on,
we have from Assumption 1(b) and Lemma 3.1 that
Since
and since kXN Y \Gamma1
O(-) from (15a), we have
\Theta I

\Theta I

Hence
Substitution in Lemma 3.1, together with Assumption 1(b) and some manipulation,
yields
We have exact arithmetic, but roundoff errors in the calculation of
restrict us to assuming that O(-+u). Using this fact together
with Assumption 1(b), formula (14), and Lemma 2.1, we have
Therefore, from (21), we have
as required.
Our next result bounds the difference between u and -
u.
LINEAR EQUATIONS IN INTERIOR-POINT METHODS 9
Theorem 3.3. Suppose the assumptions of Theorem 3.2 hold. Then for all -
sufficiently small, we have
(23a)
and, for all i 2 N ,
(24a)
Proof. The expression (23a) follows immediately from (14) and Theorem 3.2,
since
For (23b), note first from (6) and (21) that
ENN
Hence, from (14), (22), (23a), Assumption 1(b), and Theorem 3.2, we have
From (15a), we have
so (23b) is proved.
From (4), we have
and therefore
Because x i y O(1). Also from Assumption 2 and (14), we
have for i 2 N that v i =y Hence, (24a) is obtained by using these estimates
in (25). For (24b), we have from (15a), (23b), and (24a) that for all i 2 N ,
Note that in Theorems 3.2 and 3.3, we have ignored possible errors that are
introduced into the computation during the formation of the right-hand side b from
the vectors r, x, and y, and the scalars oe and -. Since the formation process introduces
a relative perturbation of O(u) into each component of b, we lose nothing by ignoring
the perturbations.
We now turn to recovery of the step - v. From the exact formula (6b), we have
In the actual computation of -
v, we have only the computed value - u available to us.
Moreover, errors are introduced when each of the five or six floating-point operations
on the right-hand side of (26) are performed. The exact nature of these errors will
depend on the order in which the operations in (26) are performed. Two possibilities
are suggested by the parentheses in the expressions
\Gammay
We can, however, perform an analysis that takes all the possibilities into account, as
we show in the following theorem.
Theorem 3.4. Suppose the assumptions of Theorem 3.3 hold and that - v is
computed from the formula (6b) (equivalently, (26)) with - u replacing u. Then we
have
(27a)
and, for all i 2 B,
(28a)
Proof. In all the formulae of this proof, we use the notation
to represent scalar quantities of order u (We certainly have
the proof.) Recall that relative errors of O(u) are incurred whenever a real number is
approximated by a floating-point number and when an arithmetic operation involving
two floating-point numbers is performed (cf. (9)). Therefore, regardless of the order
in which the operations required to recover -
are performed, we have
\Gammay i
By rearranging and combining terms in (29), we obtain
By substituting from (26), we obtain
Consider first the case of i 2 B. From (30) together with Assumption 2, Theorem
3.2, expressions (15) and (23a), and O(u), we have
proving (27a). For i 2 N , we have from Assumption 2, Theorem 3.2, and expressions
and (23b) that
LINEAR EQUATIONS IN INTERIOR-POINT METHODS 11
Therefore we obtain from (30) that
as required.
The inequalities (28) are proved in the same way as (24).
Gaussian elimination with complete pivoting is possibly more relevant to practical
algorithms, since sparse elimination algorithms rearrange both rows and columns and
hence can be regarded as approximations to the complete pivoting strategy. The main
error results for complete pivoting are the same as those for partial pivoting. To justify
this claim, we note first that the nonbasic indices will eventually be used as pivots
before any of the basic indices are used, because of the large sizes of y
Moreover, the error matrices EBN and ENN are O(u) rather than O(-
the elements cannot appear in a pivot row except on the diagonal, so
they cannot "contaminate" other elements in the nonbasic columns of M +X \Gamma1 Y . In
other words, -
actually solves the system
uN
where
(32a)
(32c)
(32d)
By defining G as the coefficient matrix in (31) and partitioning as before, we obtain
after some manipulation that
Therefore, using kb O(1), we have
which is the same error result as the one obtained for partial pivoting in Theorem
3.2. The other results in the section also hold for complete pivoting, with minor
modifications to the proofs.
4. Effect of Roundoff Error on Local Convergence. We now consider the
algorithm in [14], which can be described as follows. Given parameters fl k+1 2
1), the step ff k is chosen as
subject to
(34a)
for all ff 2 [0; ff k ]. The choices of oe k , fl k+1 , and fi k at each iteration are made according
to the following scheme.
fl), and
Find k from (11), (33), and (34) with
then accept this step; t k+1 / to next k;
Find k from (11), (33), and (34) with
accept this step; t k+1 / t k ; go to next k;
This algorithm takes two types of steps - "safe" steps, for which oe k -
oe, and "fast"
steps, for which oe Theoretically, the safe steps ensure good global convergence
properties and complexity, while the fast steps ensure asymptotic superlinear conver-
gence. The counter t k keeps track of the number of fast steps taken prior to iteration
k.
The choice of step length ff k ensures that kr k+1 To see this, note
from condition (34b) that
Y
Since
Y
Y
we have
fikr
so condition (12) holds with
LINEAR EQUATIONS IN INTERIOR-POINT METHODS 13
We focus on this algorithm because of its strong theoretical properties, namely,
global convergence from any positive starting point
when properly initialized, and superlinear local convergence. Also, the method performs
well in computational tests and is quite similar (at least in its "non-superlinear"
phase where oe k -
oe) to the algorithm implemented by Lustig, Marsten, and Shanno
[8]. We assume throughout that finite termination does not occur, that is, the algorithm
generates an infinite sequence of strictly positive iterates
In this section, we examine how the behavior of this algorithm is affected when
the computed steps (-u are used in place of the exact steps We start by
showing that near-unit steplengths can eventually be taken by this algorithm without
violating the positivity condition Consequently, there exists the possibility
of rapid convergence of the sequence of complementarity gaps - k to zero, even
in the presence of roundoff error. We refine the results to show that for the safe steps
oe), we actually have ff sufficiently small.
In all the analysis below, our convention is to use the iteration index k in the
statement of each result, but omit it in the proofs.
Lemma 4.1. For all k sufficiently large, we have
for all ff 2 [0; - ff k ], where
Proof. We consider first the indices i 2 N . From (27b), we have
while from Assumption 2 we have for large k that y k
for all ff 2 [0; 1] and all k sufficiently large. On the other hand, we have from (24b)
that
Therefore, if x
then we must have
Hence
ff satisfying (35).
The case of i 2 B is proved in a similar way by using (28b).
We now show that near-unit steps produce fast linear convergence of the complementarity
gap to zero.
Theorem 4.2. If k is sufficiently large, then
Proof. For any we have from the second part of equation (4) that
14 STEPHEN J. WRIGHT
Now, by Assumption 2 and relations (15a) and (27), we have
A similar result holds for jy i (-u )j. For the last term in (36), we have from (14),
(27), and Theorem 3.2 that
We also have -=(x i y Using these estimates in (36), we obtain
By summing over i, we obtain
which yields the desired result.
We now examine the safe steps, for which oe k -
oe, and show that ff
the criteria (33) and (34) for large enough k, even when the computed search direction
used in place of the exact direction That is, a unit step is taken.
Theorem 4.3. Suppose that -
oe is substantially larger than u, in a sense to be
defined below. Then for all sufficiently large k, if a safe step (with oe k - oe,
and is computed, the step length parameter satisfying (33) and (34) will be
Proof. From (38), we have
Therefore (34b) will hold for all ff 2 [0; 1] (with (u; v) replaced by (-u; - v) and fi
provided that the term in square brackets in (39) is nonnegative. But nonnegativity
is guaranteed for u - oe and - sufficiently small, so ff inequality.
Consider now (34a), with From (37) and (38), we have
for some positive constants C
that (34a) is satisfied if
or, equivalently,
oe
LINEAR EQUATIONS IN INTERIOR-POINT METHODS 15
This last inequality, and therefore (34a), holds provided that - and u are small enough
with respect to -
oe, as we have assumed.
Finally, we show that
is decreasing on the interval ff 2 [0; 1]. Since from (38), we have
-(ff)
the derivative -
nonpositive provided that
and, by our assumptions in the first part of this proof, -
oe dominates
the O(- we have that (41) holds.
We have shown that all ff 2 [0; 1] satisfy the conditions (34) with (u; v) replaced
by (-u; - v). Moreover, the function in (33) is decreasing over this interval. We conclude
that ff is the step chosen by the line search procedure, giving the result.
We turn now to fast step, for which oe
The (exact) analysis in Wright [14] shows that fast steps are eventually always
taken by the algorithm. Note that if the fast step is accepted, we have
The following theorem gives an estimate for the length of a fast step.
Theorem 4.4. If a fast step is attempted at iteration k, where k is sufficiently
large and u - 1, then
Proof. As in (37), we have for oe
From (38) we have
Putting (43) and (44) together, we deduce that (34a) holds
provided that
which in turn is true if
This last inequality is implied by the following bound on ff:
where we have used (42) to derive the final inequality.
Using (44) again, we have
so the inequality (34b) is satisfied when
that is, when
Providing we can show that -
is decreasing on [0; 1], we have from (45) and (46)
that the result holds for defined by
However,
so - k (ff) is certainly decreasing on [0; 1], and the result is proved.
This result accurately indicates the behavior of fast steps on later iterations of the
algorithm. The quantity j k is typically either extremely small or else quite significant
(that is,
\Omega\Gamma1699 depending on the sign of certain products such as -
and so on. When j k is tiny and the value of ff k is very close to 1,
and the fast step is accepted with a large reduction in -. When j k is larger, or when
O(u), the fast step may not lead to a very large decrease in - and may even be
rejected in favor of a safe step.
We summarize the results of this section in the following theorem.
Theorem 4.5. Suppose that u is much smaller than -
oe, in the sense of Theorem
4.3. Then for all sufficiently large k we have that either
(i) a fast step is taken, and
with
or
(ii) a safe step is taken, with
Proof. The condition for fast-step acceptance yields (47). The estimate (48)
follows from Theorems 4.2 and 4.4 and the identity oe The safe-step estimate
(49) follows from Theorem 4.2 when we use ff
LINEAR EQUATIONS IN INTERIOR-POINT METHODS 17
5. Computational Results. The algorithm of Section 4 was implemented in
double-precision Fortran, using the LAPACK routines dgetrf and dgetrs to solve the
linear system (6a). Our test problems are of two types.
(i) The matrix M has the form dense with all
elements drawn from a uniform distribution on [\Gamma1; 1], while D is diagonal
with diagonal elements of the form 10 - , where - is drawn from a uniform
distribution on [0; 1]. Since we choose n r - n, the rank of M is n r . (Rank-
deficiency of M is not an artificial feature; in certain applications of (1),
including (3), M is structurally rank-deficient.) The solutions x   and y   are
chosen so that the even-numbered components of x   and the odd-numbered
components of y   are zero, while the remaining components are uniform on
(ii) The matrix M has the form (3a), where the matrix A is dense with elements
of the form are drawn from uniform distributions on
respectively. If p and m denote the dimensions of of z and
-, respectively, we choose the even-numbered components
of both z   and -   to be nonzero. (It is a consequence of nondegeneracy of the
solution of (2) that the same number of components of z   and -   be nonzero.
This requirement is also necessary for nonsingularity of MBB .) The nonzero
components of z   , -   and the complementary vector pair in (3) are all drawn
from [0; 1].
We use the following values for the constants:
We choose the centering parameter oe k at each safe iteration according to the formula
Though we made no special effort to tune these constants to their optimal values, our
experience indicates that the choices (50), (51) are efficient for these and other types
of problems.
In

Tables

1-4 we tabulate the behavior of the algorithm of Section 4. Many of
the uninteresting safe iterates are omitted. An asterisk in the last column indicates
that a fast step was taken from this iterate. We terminate the algorithm when - falls
below
These tables indicate rapid convergence of the algorithm during its final stages.
Typically, the algorithm takes only fast steps after it has decreased - below a certain
threshold. (The experience of the author and others indicates that this threshold is
quite small for linear and quadratic problems, so that superlinear convergence does not
set in until quite late in the process. Preliminary experience with nonlinear problems
indicates that fast steps are typically taken at an earlier stage, that is, the threshold
is not so small.)
The behavior observed in Tables 1-4 certainly confirms the efficacy of Gaussian
elimination with partial pivoting in the context of this interior-point method. The
linear algebra continues to produce good steps even when - is extremely small. The
convergence of - to zero appears to be superlinear in each case (even quadratic, in
the case of Table 4). These tables do not, however, show the asymptotic behavior
suggested by Theorem 4.5. To see it, we must continue to run the algorithm past
the point of convergence. Table 5 shows what happens when we continue to iterate
on the problem of Table 3 until - is reduced below 10 \Gamma100 . (The late asymptotic

Table
Convergence of the algorithm: Problem type (i),
Fast step?
4.3
3.2 3.3
3 2.7 -11.7
19 -9.5 -13.2

Table
Convergence of the algorithm: Problem type (i),
Fast step?
4.1 5.5
4 3.3 -9.8
22 -6.7 -11.6
LINEAR EQUATIONS IN INTERIOR-POINT METHODS 19

Table
Convergence of the algorithm: Problem type (i),
Fast step?
4.3 5.8
4.1 5.5
4 3.2 3.9
26 -4.4 -11.7
28 -6.7 -11.7
29 -8.1 -11.7

Table
Convergence of the algorithm: Problem type (ii),
Fast step?
3 3.1 4.2
4 2.1 3.1
6 1.4 2.2
9

Table
Later iterates on the problem of Table 3
Fast step?
34 -42.1 -11.9
36 -50.9 -11.8
42 -102.4 -11.9 terminate
convergence was qualitatively similar on all the problems we tried, so we report just
this one instance.) Note that fast steps are taken on each iteration with decrease
iteration - the 38th - on which a
safe step is taken with a decrease ratio of almost exactly oe . The existence of
these two kinds of steps and their effects on - k are in close accord with the predictions
of Theorem 4.5.
Note that in all the tables the residual norm kr k k decreases to O(u) but no further.
As discussed in the proof of Theorem 3.2, this behavior is due to roundoff error in the
calculation of r k via the formula r
We experimented with a version of the code in which a modified complete pivoting
strategy was used for solving (6a). The columns of the coefficient matrix were ordered
by decreasing value of k \Delta k1 before Gaussian elimination with partial pivoting was
applied. Asymptotically, this strategy has the effect of ordering the nonbasic columns
first, so the analysis at the end of Section 3 still applies. As predicted in that analysis,
this version of the algorithm behaves only slightly differently from the partial pivoting
version described above.
The assumption that MBB is nonsingular (indeed, well conditioned) plays an
important role in the analysis of Sections 3 and 4. Theoretically, the algorithm of
Section 4 is known to have fast local convergence even when MBB is singular and
the solution is not unique. We tested to see whether fast convergence was attainable
in practice by forming a problem from the class (i) with n contains
50 indices, the submatrix MBB is certainly rank deficient. The result of this run is
summarized in Table 6. It is clear that the behavior indicated in Theorem 4.5 does
not occur. After taking two fast steps and converging to - iteration 14,
the algorithm stalls and makes very little progress from that point on. This and other
similar examples suggest that the assumption of MBB nonsingular probably cannot
be relaxed.

Acknowledgement

. I thank the editor and referees for their insightful comments
on an earlier draft.
LINEAR EQUATIONS IN INTERIOR-POINT METHODS 21

Table
Convergence in the case of MBB rank-deficient: Problem type (i),
Fast step?
4.2 5.9
3 2.7 4.2
14 -7.3 -7.0
100 -8.5 -8.3



--R

Solving symmetric indefinite systems in an interior-point method for linear programming
Matrix Computations

A predictor-corrector method for linear complementarity problems with polynomial complexity and superlinear convergence

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 the implementation of a primal-dual interior point method
Interior path-following primal-dual algorithms

An o(nl) infeasible-interior-point algorithm for LCP with quadratic convergence
Technical Report SOR 92-5



A simplified homogeneous and self-dual linear programming algorithm and its implementation
On quadratic and O( p nL) convergence of a predictor-corrector algorithm for LCP
On the convergence of a class of infeasible-interior-point methods for the horizontal linear complementarity problem
--TR
primal-dual interior-point methods;error analysis;stability