STABILITY OF LINEAR EQUATIONS SOLVERS IN INTERIOR-POINT METHODS^{∗}

STEPHEN J. WRIGHT^{†}

SIAM J. MAT. ANAL. APPLS \begin{matrix}c \\ ◯ \end{matrix} 1994 Society for Industrial and Applied Mathematics Vol. 0, No. 0, pp. 000–000, Month 1994 012

Abstract. 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.

Key words. primal-dual interior-point methods, error analysis, stability

AMS(MOS) subject classifications. 65G05, 65F05, 90C33

1. Introduction. The monotone linear complementarity problem (LCP) is the problem of finding a vector pair (x, y) ∈ R l ^{n} × R l ^{n} such that

y = Mx + q,  (x, y) ≥ 0,  x^{T} y = 0, (1)

where M (a real, n × n positive semidefinite matrix) and q (a real vector with n elements) 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

(2)  min_{z} c^{T} z subject to Az ≥ b, z ≥ 0,

where A ∈ R l ^{m×p}, we can introduce the dual variable λ ∈ R l ^{m} for the constraint Az ≥ b and obtain the following necessary and sufficient conditions for optimality of the primal-dual pair (z, λ):

[ \begin{matrix}0 & −A^{T} \\ A & 0 \end{matrix} ] [ \begin{matrix}z \\ λ \end{matrix} ] + [ \begin{matrix}c \\ −b \end{matrix} ] ≥ [ \begin{matrix}0 \\ 0 \end{matrix} ] , [ \begin{matrix}z \\ λ \end{matrix} ] ≥ [ \begin{matrix}0 \\ 0 \end{matrix} ] , (3a)

z^{T} (c − A^{T} λ) + λ^{T} (Az − b) = 0. (3b)

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 (x^{k}, y^{k}) that are strictly positive. Many such methods require a linear system of the form

[ \begin{matrix}M & −I \\ Y & X \end{matrix} ] [ \begin{matrix}u \\ v \end{matrix} ] = [ \begin{matrix}r \\ −XY e + σμe \end{matrix} ] , (4)

where

X = diag(x_{1}, x_{2}, …, x_{n}), Y = diag(y_{1}, y_{2}, …, y_{n}), e = (1, 1, …, 1)^{T} ,

μ = x^{T} y/n,  r = y − Mx − q,  σ ∈ [0, 1],

∗ 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.

† 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 σ = 0 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. Affinescaling steps (u, v) are simply Newton steps for the system of nonlinear equations

F (x, y) = [ \begin{matrix}Mx + q − y \\ XY e \end{matrix} ] = [ \begin{matrix}0 \\ 0 \end{matrix} ] . (5)

Path-following methods (see, for example, Monteiro and Adler [10], Zhang [19], Wright [14]) generate steps by using generally positive values of σ in (4). (As we see later, the algorithm of [14] allows σ = 0 on 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. Predictorcorrector methods (see, for example, Ye and Anstreicher [18], Ji, Potra, and Huang [4], Potra [12]) take steps with either σ = 0 or σ = 1.

The system (4) is highly structured; since the diagonals of X and Y are strictly positive, we can rearrange the system to obtain

(6a)  (M + X^{−1}Y)u = r − y + σμX^{−1}e,
(6b)  v = −X^{−1}Y u − y + σμX^{−1}e.

In the case of linear programming (2), equation (6a) contains even more structure. Its form is

[ \begin{matrix}Z^{−1}Y_{z} & −A^{T} \\ A & Λ^{−1}Y_{λ} \end{matrix} ] [ \begin{matrix}u_{z} \\ u_{λ} \end{matrix} ] = [ \begin{matrix}A^{T} λ − c + σμZ^{−1}e \\ −Az + b + σμΛ^{−1}e \end{matrix} ] , (7)

where Y_{z} and Y_{λ} are positive diagonal matrices and

Z = diag(z_{1}, …, z_{p}), Λ = diag(λ_{1}, …, λ_{m}).

The matrix in (7) can be made symmetric indefinite by multiplying the first block row by −1. This system can be reduced even further by eliminating either u_{z} or u_{λ}. For instance, if u_{z} is eliminated, we obtain

(AY_{z}^{−1}ZA^{T} + Λ^{−1}Y_{λ})u_{λ} = (b − Az + σμΛ^{−1}e) − AY_{z}^{−1}Z(A^{T} λ − c + σμZ^{−1}e). (8)

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 i = 1, …, n, at least one of x_{i} and y_{i} is zero at the solution, we would expect some of the diagonal elements of X^{−1}Y to approach zero and some to approach +∞ 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 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ón [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ón 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ón’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 k) appear as subscripts on scalars and as superscripts on vectors and matrices. The sets B and N form a partition of the index set {1, 2, …, n} defined in Assumption 1 below. If x ∈ R l ^{n}, then

x_{B} = [x_{i}]_{i∈B},  X = diag(x) = diag(x_{1}, x_{2}, …, x_{n}),  X_{B} = diag(x_{B}),

and so on. For the matrix M ∈ R l ^{n×n}, we have

M_{BN} = [M_{ij}]_{i∈B,j∈N},

and similarly for M_{BB}, M_{NB}, and M_{NN}. Given any matrix H = [h_{ij}], we define |H| = [|h_{ij}|] and denote the j-th column of H by H_{·,j}. The notation ∥ · ∥ denotes the 1-, 2-, or ∞-norm of a vector or matrix, while κ(·) denotes the corresponding condition number.

For any two nonnegative numbers ν and χ, we write ν = O(χ) if there is a moderate constant τ such that ν ≤ τχ. When W is a matrix or vector, we write W = O(χ) to denote ∥W ∥ = O(χ). We use ν = Ω(χ) to indicate that both ν = O(χ) and χ = O(ν).

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 +, −, ×, /, and fl(z) denotes the floating point representation of any real number z, we have

fl(x op y) = (x op y)(1 + δ), |δ| ≤ u. (9)

(See Golub and Van Loan [2, Section 2.4.2].) We assume throughout that u is small enough that O(u) ≪ 1, where O(·) 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 (x^{k}, y^{k}) to neighborhoods of the form

N (γ) = {(x, y) > 0 | x_{i}y_{i} ≥ γ(x^{T} y/n)} ⊂ R l _{+}^{n} × R l _{+}^{n}, (10)

where R l _{+}^{n} denotes the nonnegative orthant in R l ^{n}. All iterates generated by the algorithms lie in N (γ_{min}), where γ_{min} ∈ (0, 1/2) is a constant. Other quantities needed to define the general algorithm include

μ_{k} = (x^{k})^{T} y^{k}/n,  r^{k} = y^{k} − Mx^{k} − q,  γ_{max} ∈ (γ_{min}, 1/2],
X^{k} = diag(x^{k}), Y^{k} = diag(y^{k}).

We define the algorithmic framework as follows:

Algorithm PFI given γ_{0} ∈ [γ_{min}, γ_{max}] and (x^{0}, y^{0}) ∈ N (γ_{0})

for k = 0, 1, 2, … choose σ_{k} ∈ [0, 1] and find (u^{k}, v^{k}) that satisfies

[ \begin{matrix}M & −I \\ Y^{k} & X^{k} \end{matrix} ] [ \begin{matrix}u^{k} \\ v^{k} \end{matrix} ] = [ \begin{matrix}r^{k} \\ −X^{k}Y^{k}e + σ_{k}μ_{k}e \end{matrix} ]; (11)

choose γ_{k+1} ∈ [γ_{min}, γ_{max}] and α_{k} > 0 such that
(x^{k+1}, y^{k+1}) = (x^{k}, y^{k}) + α_{k}(u^{k}, v^{k}) ∈ N (γ_{k+1})
and ∏_{j=0}^{k}(1 − α_{j}) ≤ Kμ_{k+1}/μ_{0} when r^{0} ≠ 0, for some constant K > 0;

end (for)

The decrease in ∥r^{k}∥ at each iteration is linear — in fact, r^{k+1} = ∏_{j=0}^{k}(1 − α_{j})r^{0} — so the last condition in Algorithm PFI is equivalent to

∥r^{k+1}∥/∥r^{0}∥ ≤ Kμ_{k+1}/μ_{0}. (12)

Hence, ∥r_{k}∥ = O(μ_{k}) for all k, so the infeasibility is always bounded by a multiple of the complementarity gap μ.

When the initial point is feasible (r^{0} = 0), 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 γ_{k+1} and σ_{k} so that a step α_{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. (A more common strategy, for which there is no supporting theory, is to find the largest value of α in [0, 1] that keeps (x^{k}, y^{k}) + α(u^{k}, v^{k}) in the nonnegative orthant and then choose α_{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 (x^{∗}, y^{∗}) such that x^{∗} + y^{∗} > 0 (i.e. strict complementarity holds). We can define an associated partition B, N of the index set {1, …, n} such that x_{i}^{∗} > 0 for all i ∈ B and y_{i}^{∗} > 0 for all i ∈ N.

b) The quantities

∥M∥,  ∥M_{BB}^{−1}∥,  ∥X_{B}^{∗}∥,  ∥(X_{B}^{∗})^{−1}∥,  ∥Y_{N}^{∗}∥,  ∥(Y_{N}^{∗})^{−1}∥

are all moderate in size. Assumption 1 implies that the coefficient matrix in (11) approaches a nonsingular

limit, since there are 2n × 2n permutation matrices P and Π such that

P [ \begin{matrix}M & −I \\ Y^{∗} & X^{∗} \end{matrix} ] Π = [ \begin{matrix}X_{B}^{∗} & 0 & 0 & 0 \\ 0 & Y_{N}^{∗} & 0 & 0 \\ −I & M_{BN} & M_{BB} & 0 \\ 0 & M_{NN} & M_{NB} & −I \end{matrix} ] , (13)

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 M_{BB} 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 Az = b 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_{k→∞} (x^{k}, y^{k}) = (x^{∗}, y^{∗}).

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 σ_{k}, α_{k}, and γ_{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 ∥r^{k}∥ = O(μ_{k}), we have

∥(u^{k}, v^{k})∥ = O(μ_{k}), ∥(x^{k}, y^{k}) − (x^{∗}, y^{∗})∥ = O(μ_{k}). (14)

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

C_{1}μ_{k} ≤ x_{i}^{k} ≤ C_{2}μ_{k}, ∀i ∈ N, (15a) C_{1}μ_{k} ≤ y_{i}^{k} ≤ C_{2}μ_{k}, ∀i ∈ B. (15b)

Proof. Because of Assumption 2, we can define an index K and positive constants C̄_{1} and C̄_{2} such that

x_{i}^{k} ∈ [C̄_{2}, C̄_{1}], ∀k ≥ K, ∀i ∈ B,
y_{i}^{k} ∈ [C̄_{2}, C̄_{1}], ∀k ≥ K, ∀i ∈ N.

Therefore, since (x^{k}, y^{k}) > 0, we have for k ≥ K, i ∈ N, that

x_{i}^{k}y_{i}^{k} < (x^{k})^{T} y^{k} = nμ_{k} ⇒ x_{i}^{k} < \frac{nμ_{k}}{y_{i}^{k}} ≤ \frac{nμ_{k}}{C̄_{2}} .

Also,

x_{i}^{k}y_{i}^{k} ≥ γ_{k}μ_{k} ≥ γ_{min}μ_{k} ⇒ x_{i}^{k} ≥ \frac{γ_{min}μ_{k}}{y_{i}^{k}} ≥ \frac{γ_{min}μ_{k}}{C̄_{1}} .

Therefore (15a) holds with C_{1} = γ_{min}/C̄_{1} and C_{2} = n/C̄_{2}. 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 × m linear system Az = b is 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 Π, the computed unit lower triangular factor by L̂, and the computed upper triangular factor by Û. Then the computed solution ẑ solves the perturbed system (A + H)ẑ = b, where

|P HΠ| ≤ ϵ_{m}(2 + ϵ_{m})|L̂||Û|, (16)

and ϵ_{m} = mu/(1 − mu) = O(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 ρ, defined as the smallest positive number such that

(17)  max_{i=1,...,m} |Û_{ij}| ≤ ρ max_{i=1,...,m} |(P AΠ)_{ij}|,  for all j = 1, 2, …, m.

We then have the following simple corollary of Theorem 2.2.

Corollary 2.3. Let the system Az = b be as in Theorem 2.2, and suppose the pivots Û_{jj} are chosen so that |L̂_{ij}| ≤ 1 for all i = 2, …, m, j = 1, …, i − 1. Then

∥H_{·,j}∥ ≤ mϵ_{m}(2 + ϵ_{m})ρ∥A_{·,j}∥u. (18)

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 (û, v̂), and the corresponding exact values, denoted by (u, v). We treat the cases in which (û, 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 = [ \begin{matrix}G_{11} & G_{12} \\ G_{21} & G_{22} \end{matrix} ] ,

where G_{11} and G_{22} are also square. Suppose that G_{11} and G_{22} − G_{21}G_{11}^{−1}G_{12} are nonsingular. Then G is nonsingular and G^{−1} has the form

[ \begin{matrix}G_{11}^{−1} + G_{11}^{−1}G_{12}(G_{22} − G_{21}G_{11}^{−1}G_{12})^{−1}G_{21}G_{11}^{−1} & −G_{11}^{−1}G_{12}(G_{22} − G_{21}G_{11}^{−1}G_{12})^{−1} \\ −(G_{22} − G_{21}G_{11}^{−1}G_{12})^{−1}G_{21}G_{11}^{−1} & (G_{22} − G_{21}G_{11}^{−1}G_{12})^{−1} \end{matrix} ] .

Our first main result concerning components of û is the following.

Theorem 3.2. Let û be computed by applying Gaussian elimination with row partial pivoting to (6a), and suppose that the growth factor ρ is not too large. Then for all μ sufficiently small, we have

∥û_{B}∥ = O(u + μ), ∥û_{N}∥ = O(μ).

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

[ \begin{matrix}M_{BB} + X_{B}^{−1}Y_{B} + E_{BB} & M_{BN} + E_{BN} \\ M_{NB} + E_{NB} & M_{NN} + X_{N}^{−1}Y_{N} + E_{NN} \end{matrix} ] [ û_{B} ; û_{N} ]
(19)  = [ \begin{matrix}r_{B} − y_{B} − σμX_{B}^{−1}e \\ r_{N} − y_{N} − σμX_{N}^{−1}e \end{matrix} ] ,

where

(20a)  ∥E_{BB}∥ ≤ ρ O(∥M_{BB} + X_{B}^{−1}Y_{B}∥ + ∥M_{NB}∥)u,
(20b)  ∥E_{NB}∥ ≤ ρ O(∥M_{BB} + X_{B}^{−1}Y_{B}∥ + ∥M_{NB}∥)u,
(20c)  ∥E_{BN}∥ ≤ ρ O(∥M_{NN} + X_{N}^{−1}Y_{N}∥ + ∥M_{BN}∥)u,
(20d)  ∥E_{NN}∥ ≤ ρ O(∥M_{NN} + X_{N}^{−1}Y_{N}∥ + ∥M_{BN}∥)u.

Now from Assumption 1(b) and Lemma 2.1, we have

∥X_{B}^{−1}Y_{B}∥ = O(μ), ∥X_{N}^{−1}Y_{N}∥ = O(μ^{−1}),

and ∥M_{BB}∥, ∥M_{BN}∥, ∥M_{NB}∥, and ∥M_{NN}∥ are all O(1). Combining these observations with (20), we obtain

∥E_{BB}∥ = O(u),  ∥E_{NB}∥ = O(u),  ∥E_{BN}∥ = O(μ^{−1}u),  ∥E_{NN}∥ = O(μ^{−1}u).

Therefore (19) can be rewritten as

[ \begin{matrix}M_{BB} + Ē_{BB} & M_{BN} + Ē_{BN} \\ M_{NB} + Ē_{NB} & M_{NN} + X_{N}^{−1}Y_{N} + Ē_{NN} \end{matrix} ] [ û_{B} ; û_{N} ] = [ \begin{matrix}b_{B} \\ b_{N} \end{matrix} ] , (21)

where

∥Ē_{BB}∥ = O(μ + u),  ∥Ē_{NB}∥ = O(u),
(22)  ∥Ē_{BN}∥ = O(μ^{−1}u),  ∥Ē_{NN}∥ = O(μ^{−1}u),

and

b_{B} = r_{B} − y_{B} − σμX_{B}^{−1}e, b_{N} = r_{N} − y_{N} − σμX_{N}^{−1}e.

If we denote the coefficient matrix in (21) by G, with G_{11} = M_{BB} + Ē_{BB} and so on, we have from Assumption 1(b) and Lemma 3.1 that

G_{11}^{−1} = [I + M_{BB}^{−1}Ē_{BB}]^{−1}M_{BB}^{−1} = M_{BB}^{−1} + O(μ + u).

Since

∥G_{21}∥ = O(1),  ∥G_{12}∥ = O(1 + μ^{−1}u),  ∥G_{11}^{−1}∥ = O(1),

and since ∥X_{N}Y_{N}^{−1}∥ = O(μ) from (15a), we have

G_{22} − G_{21}G_{11}^{−1}G_{12} = X_{N}^{−1}Y_{N} [I + X_{N}Y_{N}^{−1}(M_{NN} + Ē_{NN}) − X_{N}Y_{N}^{−1}G_{21}G_{11}^{−1}G_{12}]
= X_{N}^{−1}Y_{N} [I + O(μ)O(1 + μ^{−1}u) + O(μ)O(1)O(1)O(1 + μ^{−1}u)]
= X_{N}^{−1}Y_{N}[I + O(μ + u)].

Hence

(G_{22} − G_{21}G_{11}^{−1}G_{12})^{−1} = (I + O(μ + u)) Y_{N}^{−1}X_{N} = O(μ).

Substitution in Lemma 3.1, together with Assumption 1(b) and some manipulation, yields

(G^{−1})_{11} = M_{BB}^{−1} + O(μ + u),  (G^{−1})_{12} = O(μ + u),
(G^{−1})_{21} = O(μ),  (G^{−1})_{22} = O(μ).

We have ∥r∥ = O(μ) in exact arithmetic, but roundoff errors in the calculation of r = y − Mx − q restrict us to assuming that ∥r∥ = O(μ + u). Using this fact together with Assumption 1(b), formula (14), and Lemma 2.1, we have

∥b_{B}∥ = O(μ + u), ∥b_{N}∥ = O(1).

Therefore, from (21), we have

û_{B} = (G^{−1})_{11}b_{B} + (G^{−1})_{12}b_{N} = O(μ + u),
û_{N} = (G^{−1})_{21}b_{B} + (G^{−1})_{22}b_{N} = O(μ),

as required.

Our next result bounds the difference between u and û.

Theorem 3.3. Suppose the assumptions of Theorem 3.2 hold. Then for all μ sufficiently small, we have

(23a)  ∥û_{B} − u_{B}∥ = O(μ + u),
(23b)  ∥û_{N} − u_{N}∥ = O(μ(μ + u)),

and, for all i ∈ N,

(24a)  u_{i}/x_{i} = −1 + O(σ + μ),
(24b)  û_{i}/x_{i} = −1 + O(σ + μ + u).

Proof. The expression (23a) follows immediately from (14) and Theorem 3.2, since

∥û_{B} − u_{B}∥ ≤ ∥û_{B}∥ + ∥u_{B}∥ = O(μ + u).

For (23b), note first from (6) and (21) that

M_{NB}u_{B}+(M_{NN} +X_{N}^{−1}Y_{N})u_{N} = b_{N} = (M_{NB}+Ē_{NB})û_{B}+(M_{NN} +X_{N}^{−1}Y_{N} +Ē_{NN})û_{N}.

Hence, from (14), (22), (23a), Assumption 1(b), and Theorem 3.2, we have

X_{N}^{−1}Y_{N}(u_{N} − û_{N}) = M_{NB}(û_{B} − u_{B}) + M_{NN}(û_{N} − u_{N}) + Ē_{NB}û_{B} + Ē_{NN} û_{N}
= O(μ + u) + O(μ) + O(μu) + O(u) = O(μ + u).

From (15a), we have X_{N}Y_{N}^{−1} = O(μ) and so (23b) is proved.

From (4), we have

y_{i}u_{i} + x_{i}v_{i} = −x_{i}y_{i} + σμ

and therefore

\frac{u_{i}}{x_{i}} = −1 − \frac{v_{i}}{y_{i}} + \frac{σμ}{x_{i}y_{i}} . (25)

Because x_{i}y_{i} ≥ γ_{min}μ, we have μ/(x_{i}y_{i}) = O(1). Also from Assumption 2 and (14), we have for i ∈ N that v_{i}/y_{i} = O(μ). Hence, (24a) is obtained by using these estimates in (25). For (24b), we have from (15a), (23b), and (24a) that for all i ∈ N,

\frac{û_{i}}{x_{i}} = \frac{u_{i}}{x_{i}} + \frac{1}{x_{i}} O(μ(μ + u)) = −1 + O(μ + σ) + O(μ + u).

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 σ 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

v_{i} = −y_{i}(1 + u_{i}/x_{i}) + σμ/x_{i}. (26)

In the actual computation of v̂, we have only the computed value û 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

−y_{i}[1 + (u_{i}/x_{i})] + (σμ)/x_{i}, [−y_{i} − (y_{i}u_{i})/x_{i}] + σ(μ/x_{i}).

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 û replacing u. Then we have

(27a)  ∥v̂_{B} − v_{B}∥ = O(μ(μ + u)),
(27b)  ∥v̂_{N} − v_{N}∥ = O(μ + u),

and, for all i ∈ B,

(28a)  v_{i}/y_{i} = −1 + O(σ + μ),
(28b)  v̂_{i}/y_{i} = −1 + O(σ + μ + u).

Proof. In all the formulae of this proof, we use the notation δ_{j}, j = 1, 2, …, to represent scalar quantities of order u (We certainly have δ_{j} ≤ 10u throughout 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 v̂_{i} are performed, we have

v̂_{i} = [−y_{i} (1 + δ_{1} + \frac{û_{i}}{x_{i}} (1 + δ_{2})) (1 + δ_{3}) + \frac{σμ}{x_{i}}(1 + δ_{4})] (1 + δ_{5}). (29)

By rearranging and combining terms in (29), we obtain

v̂_{i} = −y_{i}(1 + δ_{6}) − \frac{y_{i}û_{i}}{x_{i}} (1 + δ_{7}) + \frac{σμ}{x_{i}}(1 + δ_{8})
= −y_{i} − \frac{y_{i}u_{i}}{x_{i}} + \frac{σμ}{x_{i}} + \frac{y_{i}}{x_{i}} (u_{i} − û_{i}) − δ_{6}y_{i} − δ_{7} \frac{y_{i}û_{i}}{x_{i}} − δ_{8} \frac{σμ}{x_{i}}.

By substituting from (26), we obtain

|v_{i} − v̂_{i}| = |\frac{y_{i}}{x_{i}} | |u_{i} − û_{i}| + O(u) [|y_{i}| + |\frac{y_{i}û_{i}}{x_{i}} | + |\frac{σμ}{x_{i}} |] . (30)

Consider first the case of i ∈ B. From (30) together with Assumption 2, Theorem 3.2, expressions (15) and (23a), and δ_{j} = O(u), we have

|v_{i} − v̂_{i}| = O(μ(μ + u)) + O(u) [O(μ) + O(μ(μ + u)) + O(σμ)] = O(μ(μ + u)),

proving (27a). For i ∈ N, we have from Assumption 2, Theorem 3.2, and expressions (15) and (23b) that

|\frac{y_{i}}{x_{i}} | = O(μ^{−1}),  |u_{i} − û_{i}| = O(μ(μ + u)),  |\frac{y_{i}û_{i}}{x_{i}} | = O(1),  |\frac{σμ}{x_{i}} | = O(1).

Therefore we obtain from (30) that

|v_{i} − v̂_{i}| = O(μ + u), ∀i ∈ N,

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_{i}/x_{i}, i ∈ N. Moreover, the error matrices E_{BN} and E_{NN} are O(u) rather than O(μ^{−1}u), because the elements y_{i}/x_{i}, i ∈ N cannot appear in a pivot row except on the diagonal, so they cannot “contaminate” other elements in the nonbasic columns of M + X^{−1}Y. In other words, û actually solves the system

[ \begin{matrix}M_{BB} + X_{B}^{−1}Y_{B} + E_{BB} & M_{BN} + E_{BN} \\ M_{NB} + E_{NB} & M_{NN} + X_{N}^{−1}Y_{N} + E_{NN} \end{matrix} ] [ û_{B} ; û_{N} ]
(31)  = [ \begin{matrix}r_{B} − y_{B} − σμX_{B}^{−1}e \\ r_{N} − y_{N} − σμX_{N}^{−1}e \end{matrix} ] ,

where

(32a)  ∥E_{BB}∥ ≤ ρ O(∥M_{BB} + X_{B}^{−1}Y_{B}∥ + ∥M_{NB}∥)u = O(u),
(32b)  ∥E_{NB}∥ ≤ ρ O(∥M_{BB} + X_{B}^{−1}Y_{B}∥ + ∥M_{NB}∥)u = O(u),
(32c)  ∥E_{BN}∥ ≤ ρ O(∥M_{NN}∥ + ∥M_{BN}∥)u = O(u),
(32d)  ∥E_{NN}∥ ≤ ρ O(∥M_{NN}∥ + ∥M_{BN}∥)u = O(u).

By defining G as the coefficient matrix in (31) and partitioning as before, we obtain after some manipulation that

(G^{−1})_{11} = M_{BB}^{−1} + O(μ + u),  (G^{−1})_{12} = O(μ),
(G^{−1})_{21} = O(μ),  (G^{−1})_{22} = O(μ).

Therefore, using ∥b_{B}∥ = O(μ + u) and ∥b_{N}∥ = O(1), we have

û_{B} = (G^{−1})_{11}b_{B} + (G^{−1})_{12}b_{N} = O(μ + u),
û_{N} = (G^{−1})_{21}b_{B} + (G^{−1})_{22}b_{N} = O(μ),

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 γ_{k+1} ∈ (γ_{min}, γ_{k}] and β_{k} ∈ [0, 1), the step α_{k} is chosen as

α_{k} = arg min_{α∈(0,1]} μ_{k}(α) =^{△} (x^{k} + αu^{k})^{T} (y^{k} + αv^{k})/n, (33)

subject to

(x^{k}, y^{k}) + α(u^{k}, v^{k}) ∈ N (γ_{k+1}), (34a) μ_{k}(α) ≥ (1 − α)(1 − β_{k})μ_{k}, if r^{k} ≠ 0, (34b)

for all α ∈ [0, α_{k}]. The choices of σ_{k}, γ_{k+1}, and β_{k} at each iteration are made according to the following scheme.

given γ̄ ∈ (0, 1), γ_{min}, γ_{max} with 0 < γ_{min} < γ_{max} < 1/2, σ̄ ∈ (0, 1/2), ρ ∈ (0, γ̄), and (x^{0}, y^{0}) with x_{i}^{0}y_{i}^{0} ≥ γ_{max}μ_{0} > 0;

t_{0} ← 1, γ_{0} ← γ_{max};

for k = 0, 1, 2, …

if  μ_{k} = 0  then stop;

Find (u^{k}, v^{k}) and α_{k} from (11), (33), and (34) with

σ_{k} = 0, β_{k} = γ̄^{t_{k}}, γ_{k+1} = γ_{min} + γ̄^{t_{k}}(γ_{max} − γ_{min});

if μ_{k}(α_{k}) ≤ ρμ_{k}

then accept this step; t_{k+1} ← t_{k} + 1; go to next k; end if

Find (u^{k}, v^{k}) and α_{k} from (11), (33), and (34) with

σ_{k} ∈ [σ̄, 1/2], β_{k} = 0, γ_{k+1} = γ_{k};

accept this step; t_{k+1} ← t_{k}; go to next k;

end for. This algorithm takes two types of steps — “safe” steps, for which σ_{k} ≥ σ̄, and “fast” steps, for which σ_{k} = 0. Theoretically, the safe steps ensure good global convergence properties and complexity, while the fast steps ensure asymptotic superlinear convergence. The counter t_{k} keeps track of the number of fast steps taken prior to iteration k.

The choice of step length α_{k} ensures that ∥r^{k+1}∥ = O(μ_{k+1}). To see this, note from condition (34b) that

μ_{k+1} = (1 − α_{k})(1 − β_{k})μ_{k} = [∏_{j=0}^{k}(1 − α_{j})(1 − β_{j})] μ_{0}.

Since

∏_{j=0}^{k}(1 − β_{j}) ≥ ∏_{j=0}^{∞}(1 − γ̄^{j}) = β̂ > 0, r^{k+1} = ∏_{j=0}^{k}(1 − α_{j})r^{0},

we have

μ_{k+1}/μ_{0} ≥ β̂∥r^{k+1}∥/∥r^{0}∥,

so condition (12) holds with K = β̂^{−1}.

We focus on this algorithm because of its strong theoretical properties, namely, global convergence from any positive starting point (x^{0}, y^{0}), polynomial complexity 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 σ_{k} ≥ σ̄) 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 (x^{k}, y^{k}).

In this section, we examine how the behavior of this algorithm is affected when the computed steps (û^{k}, v̂^{k}) are used in place of the exact steps (u^{k}, v^{k}). We start by showing that near-unit steplengths can eventually be taken by this algorithm without violating the positivity condition (x^{k}, y^{k}) > 0. 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 (σ_{k} ≥ σ̄), we actually have α_{k} = 1 when μ_{k} is 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

(x^{k} + αû^{k}, y^{k} + αv̂^{k}) ≥ 0,

for all α ∈ [0, ᾱ_{k}], where

|1 − ᾱ_{k}| = O(σ_{k} + μ_{k} + u). (35)

Proof. We consider first the indices i ∈ N. From (27b), we have

|v̂_{i}| = |v_{i}| + O(μ + u) = O(μ + u),

while from Assumption 2 we have for large k that y_{i}^{k} ≈ y_{i}^{∗} > 0. Hence y_{i}^{k} + αv̂_{i}^{k} > 0 for all α ∈ [0, 1] and all k sufficiently large. On the other hand, we have from (24b) that

x_{i} + αû_{i} = x_{i} + αx_{i}(−1 + O(σ + μ + u)).

Therefore, if x_{i} + α_{i}û_{i} = 0 for some index i, then we must have

1 − α + αO(σ + μ + u) = 0 ⇒ |1 − α| = O(σ + μ + u).

Hence x_{i} + αû_{i} ≥ 0 for all α ∈ [0, ᾱ], for ᾱ satisfying (35).

The case of i ∈ 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

(x^{k} + α_{k}û^{k})^{T} (y^{k} + α_{k}v̂^{k}) = [1 − α_{k}(1 − σ_{k}) + O(μ_{k} + u)](x^{k})^{T} y^{k}.

Proof. For any i = 1, …, n, we have from the second part of equation (4) that

(x_{i} + αû_{i})(y_{i} + αv̂_{i})
= x_{i}y_{i} + αy_{i}u_{i} + αx_{i}v_{i} + αy_{i}(û_{i} − u_{i}) + αx_{i}(v̂_{i} − v_{i}) + α^{2}û_{i}v̂_{i}
(36)  = (1 − α)x_{i}y_{i} + ασμ + αy_{i}(û_{i} − u_{i}) + αx_{i}(v̂_{i} − v_{i}) + α^{2}û_{i}v̂_{i}.

Now, by Assumption 2 and relations (15a) and (27), we have

i ∈ N ⇒ |x_{i}(v̂_{i} − v_{i})| = O(μ)O(μ + u),
i ∈ B ⇒ |x_{i}(v̂_{i} − v_{i})| = O(1)O(μ(μ + u)).

A similar result holds for |y_{i}(û_{i} − u_{i})|. For the last term in (36), we have from (14), (27), and Theorem 3.2 that

i ∈ B ⇒ |û_{i}v̂_{i}| ≤ O(μ + u)(|v_{i}| + |v̂_{i} − v_{i}|) = O(μ(μ + u))
i ∈ N ⇒ |û_{i}v̂_{i}| = O(μ)(|v_{i}| + |v̂_{i} − v_{i}|) = O(μ(μ + u)).

We also have μ/(x_{i}y_{i}) ≤ 1/γ_{min} = O(1). Using these estimates in (36), we obtain

(x_{i} + αû_{i})(y_{i} + αv̂_{i}) = x_{i}y_{i}(1 − α) + ασμ + αO(μ(μ + u))
(37)  = x_{i}y_{i}(1 − α) + αμ [σ + O(μ + u)] .

By summing over i, we obtain

(x + αû)^{T} (y + αv̂) = x^{T} y [(1 − α + ασ) + αO(μ + u)] , (38)

which yields the desired result.

We now examine the safe steps, for which σ_{k} ≥ σ̄, and show that α_{k} = 1 satisfies the criteria (33) and (34) for large enough k, even when the computed search direction (û^{k}, v̂^{k}) is used in place of the exact direction (u^{k}, v^{k}). That is, a unit step is taken.

Theorem 4.3. Suppose that σ̄ is substantially larger than u, in a sense to be defined below. Then for all sufficiently large k, if a safe step (with σ_{k} ≥ σ̄, γ_{k+1} = γ_{k}, and β_{k} = 0) is computed, the step length parameter satisfying (33) and (34) will be α_{k} = 1.

Proof. From (38), we have

\frac{(x + αû)^{T} (y + αv̂)}{x^{T} y} = (1 − α) + α [σ + O(μ + u)] . (39)

Therefore (34b) will hold for all α ∈ [0, 1] (with (u, v) replaced by (û, v̂) and β_{k} = 0), provided that the term in square brackets in (39) is nonnegative. But nonnegativity is guaranteed for u ≪ σ̄ and μ sufficiently small, so α_{k} = 1 satisfies this inequality.

Consider now (34a), with γ_{k+1} = γ_{k}. From (37) and (38), we have

\frac{(x_{i} + αû_{i})(y_{i} + αv̂_{i})}{(x + αû)^{T} (y + αv̂)/n} = \frac{(1 − α)x_{i}y_{i} + αμ[σ + O(μ + u)]}{(1 − α)μ + αμ[σ + O(μ + u)]}
(40)  ≥ \frac{x_{i}y_{i}/μ + [σ − C_{10}(μ + u)]α/(1 − α)}{1 + [σ + C_{11}(μ + u)]α/(1 − α)},

for some positive constants C_{10} and C_{11}. Since x_{i}y_{i} ≥ γ_{k}μ for all i = 1, …, n, we find that (34a) is satisfied if

\frac{σ − C_{10}(μ + u)}{σ + C_{11}(μ + u)} ≥ γ_{max},

or, equivalently,

\frac{μ + u}{σ} ≤ \frac{1 − γ_{max}}{C_{10} + C_{11}} .

This last inequality, and therefore (34a), holds provided that μ and u are small enough with respect to σ̄, as we have assumed.

Finally, we show that

μ̂(α) =^{△} (x + αû)^{T} (y + αv̂)/n

is decreasing on the interval α ∈ [0, 1]. Since from (38), we have

\frac{μ̂(α)}{μ} = 1 − α + ασ + αO(μ + u),

the derivative μ̂′(α) is nonpositive provided that

− 1 + σ + O(μ + u) ≤ 0. (41)

Since σ̄ ≤ σ ≤ 1/2 and, by our assumptions in the first part of this proof, σ̄ dominates the O(μ + u) term, we have that (41) holds.

We have shown that all α ∈ [0, 1] satisfy the conditions (34) with (u, v) replaced by (û, v̂). Moreover, the function in (33) is decreasing over this interval. We conclude that α_{k} = 1 is the step chosen by the line search procedure, giving the result.

We turn now to fast step, for which σ_{k} = 0, β_{k} = γ̄^{t_{k}}, and γ_{k+1} = γ_{min}+γ̄^{t_{k}}(γ_{max}− γ_{min}). 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

γ_{k} − γ_{k+1} = (γ̄^{t_{k}−1} − γ̄^{t_{k}})(γ_{max} − γ_{min}) = γ̄^{t_{k}}(γ̄^{−1} − 1)(γ_{max} − γ_{min}) = Ω(β_{k}). (42)

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

α_{k} ≥ [1 + η_{k} \frac{μ_{k} + u}{γ̄^{t_{k}}}]^{−1} ,

where η_{k} satisfies 0 ≤ η_{k} ≤ O(1).

Proof. As in (37), we have for σ_{k} = 0 that

(x_{i} + αû_{i})(y_{i} + αv̂_{i}) = x_{i}y_{i}(1 − α) + μαO(μ + u) ≥ γ_{k}μ(1 − α) − μα(μ + u)η_{k}^{(1)}, (43)

where 0 ≤ η_{k}^{(1)} ≤ O(1). From (38) we have

(x + αû)^{T} (y + αv̂)/n = μ [1 − α + αO(μ + u)] ≤ μ [1 − α + α(μ + u)η_{k}^{(2)}] , (44)

where 0 ≤ η_{k}^{(2)} ≤ O(1). Putting (43) and (44) together, we deduce that (34a) holds provided that

γ_{k}(1 − α) − α(μ + u)η_{k}^{(1)} ≥ γ_{k+1}[1 − α + α(μ + u)η_{k}^{(2)}],

which in turn is true if

(γ_{k} − γ_{k+1})(1 − α) ≥ α(μ + u)(η_{k}^{(1)} + η_{k}^{(2)}).

This last inequality is implied by the following bound on α:

(45) α ≤ [1 + \frac{(μ + u)(η_{k}^{(1)} + η_{k}^{(2)})}{γ_{k} − γ_{k+1}}]^{−1} = [1 + \frac{(μ + u)(η_{k}^{(1)} + η_{k}^{(2)})}{γ̄^{t_{k}}(γ̄^{−1} − 1)(γ_{max} − γ_{min})}]^{−1} ,

where we have used (42) to derive the final inequality.

Using (44) again, we have

(x + αû)^{T} (y + αv̂)/n = μ [1 − α + αO(μ + u)] ≥ μ [1 − α − α(μ + u)η_{k}^{(3)}] ,

where 0 ≤ η_{k}^{(3)} ≤ O(1), so the inequality (34b) is satisfied when

1 − α − α(μ + u)η_{k}^{(3)} ≥ (1 − α)(1 − β_{k}),

that is, when

α ≤ [1 + \frac{(μ + u)η_{k}^{(3)}}{β_{k}} ]^{−1} = [1 + \frac{(μ + u)η_{k}^{(3)}}{γ̄^{t_{k}}} ]^{−1} . (46)

Providing we can show that μ̂_{k}(α) is decreasing on [0, 1], we have from (45) and (46) that the result holds for η_{k} defined by

η_{k} = max (η_{k}^{(3)}, \frac{η_{k}^{(1)} + η_{k}^{(2)}}{(γ̄^{−1} − 1)(γ_{max} − γ_{min})}) .

However,

μ̂′_{k}(α) ≤ −μ [1 − O(μ + u)] < 0,

so μ̂_{k}(α) 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 η_{k} is typically either extremely small or else quite significant (that is, η_{k} = Ω(1)), depending on the sign of certain products such as û^{T} v̂, û_{i}v̂_{i}, and so on. When η_{k} is tiny and μ_{k} + u ≪ γ̄^{t_{k}}, the value of α_{k} is very close to 1, and the fast step is accepted with a large reduction in μ. When η_{k} is larger, or when γ̄^{t_{k}} = 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 σ̄, in the sense of Theorem 4.3. Then for all sufficiently large k we have that either

(i) a fast step is taken, and

μ_{k+1} ≤ ρμ_{k}, (47)

with

μ_{k+1} = O (\frac{μ_{k} + u}{γ̄^{t_{k}}}) μ_{k}, (48)

or (ii) a safe step is taken, with

μ_{k+1} = [σ_{k} + O(μ_{k} + u)]μ_{k}. (49)

Proof. The condition for fast-step acceptance yields (47). The estimate (48) follows from Theorems 4.2 and 4.4 and the identity σ_{k} = 0. The safe-step estimate (49) follows from Theorem 4.2 when we use α_{k} = 1.

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 M = ADA^{T}, where A is n × n_{r} dense with all elements drawn from a uniform distribution on [−1, 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}. (Rankdeficiency 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 [0, 1].

(ii) The matrix M has the form (3a), where the matrix A is dense with elements of the form τ_{1}10^{τ_{2}}, where τ_{1} and τ_{2} are drawn from uniform distributions on [−\frac{1}{2}, \frac{1}{2}] and [0, 1], respectively. If p and m denote the dimensions of of z and λ, respectively, we choose the even-numbered components 2, 4, …, min(p, m) 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 M_{BB}.) 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:

γ_{min} = 10^{−5},  γ_{max} = 10^{−2},  ρ = γ̄/2,  γ̄ = 10^{−1}. (50)

We choose the centering parameter σ_{k} at each safe iteration according to the formula

(51)  σ_{k} = mid(σ̄, μ_{k}/\sqrt{n}, .2),  where σ̄ = .01.

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 10^{−20}.

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^{−100}. (The late asymptotic

Table 1

Convergence of the algorithm: Problem type (i), n = n_{r} = 20

k log_{10} μ_{k} log_{10} ∥r^{k}∥_{1} Fast step?
0  4.0  4.7
1  3.7  4.3
2  3.2  3.3
3  2.7  -11.7
⋮  ⋮  ⋮
15  -3.7  -13.4
16  -4.8  -13.3
17  -5.7  -13.2  ∗
18  -7.0  -13.1  ∗
19  -9.5  -13.2  ∗
20  -14.5  -13.5  ∗
21  -23.8  -13.4 terminate

Table 2

Convergence of the algorithm: Problem type (i), n = n_{r} = 100

k log_{10} μ_{k} log_{10} ∥r^{k}∥_{1} Fast step?
0  4.9  6.5
1  4.6  6.2
2  4.1  5.5
3  3.8  4.9
4  3.3  -9.8
⋮  ⋮  ⋮
21  -5.8  -11.6
22  -6.7  -11.6
23  -7.8  -11.6
24  -9.0  -11.6  ∗
25  -10.9  -11.5  ∗
26  -14.7  -11.6  ∗
27  -22.1  -11.8 terminate

Table 3

Convergence of the algorithm: Problem type (i), n = 100, n_{r} = 60.

k log_{10} μ_{k} log_{10} ∥r^{k}∥_{1} Fast step?
0  4.5  6.1
1  4.3  5.8
2  4.1  5.5
3  3.7  5.0
4  3.2  3.9
5  2.8  -9.5
⋮  ⋮  ⋮
24  -2.8  -11.8
25  -3.5  -11.7
26  -4.4  -11.7
27  -5.6  -11.7
28  -6.7  -11.7  ∗
29  -8.1  -11.7  ∗
30  -10.1  -11.8  ∗
31  -13.8  -11.8  ∗
32  -20.8  -11.8 terminate

Table 4

Convergence of the algorithm: Problem type (ii), n = 200, matrix A is 40 × 160

k log_{10} μ_{k} log_{10} ∥r^{k}∥_{1} Fast step?
0  4.0  5.2
1  3.8  4.9
2  3.5  4.7
3  3.1  4.2
4  2.1  3.1
5  1.7  2.6
6  1.4  2.2
7  1.0  1.7
8  0.5  1.2
9  0.2  0.8
10  -0.7  -0.3
11  -2.1  -1.8  ∗
12  -4.4  -4.1  ∗
13  -8.3  -7.9  ∗
14  -15.6  -10.2  ∗
15  -28.2  -10.2 terminate

Table 5 Later iterates on the problem of Table 3

k log_{10} μ_{k} log_{10} ∥r^{k}∥_{1} Fast step?
⋮  ⋮  ⋮
32  -20.8  -11.8  ∗
33  -31.4  -11.9  ∗
34  -42.1  -11.9  ∗
35  -46.8  -11.9  ∗
36  -50.9  -11.8  ∗
37  -61.7  -11.9  ∗
38  -71.9  -11.9
39  -73.9  -11.9  ∗
40  -78.4  -11.9  ∗
41  -90.5  -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 factors between 10^{−4} and 10^{−12}, except for one iteration — the 38th — on which a safe step is taken with a decrease ratio of almost exactly σ_{k} = 10^{−2}. 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 ∥r_{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^{k} = y^{k} − Mx^{k} − q.

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 ∥ · ∥_{∞} 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 M_{BB} 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 M_{BB} 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_{r} = 25. Since B contains 50 indices, the submatrix M_{BB} 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 μ_{k} ≈ 10^{−7} by iteration 14, the algorithm stalls and makes very little progress from that point on. This and other similar examples suggest that the assumption of M_{BB} nonsingular probably cannot be relaxed.

Acknowledgement. I thank the editor and referees for their insightful comments on an earlier draft.

Table 6

Convergence in the case of M_{BB} rank-deficient: Problem type (i), n = 100, n_{r} = 25

k log_{10} μ_{k} log_{10} ∥r^{k}∥_{1} Fast step?
0  4.2  5.9
1  4.0  5.5
2  3.6  5.1
3  2.7  4.2
4  2.1  3.4
⋮  ⋮  ⋮
11  -0.2  0.2
12  -1.2  -0.9  ∗
13  -3.7  -3.5  ∗
14  -7.3  -7.0
15  -7.4  -7.1
16  -7.4  -7.2
17  -7.6  -7.4
⋮  ⋮  ⋮
98  -8.4  -8.2
99  -8.5  -8.2
100  -8.5  -8.3
⋮  ⋮  ⋮

REFERENCES

[1] R. Fourer and S. Mehrotra, Solving symmetric indefinite systems in an interior-point method for linear programming, Mathematical Programming, 62 (1993), pp. 15–39.

[2] G. H. Golub and C. F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, 2nd ed., 1989.

[3] N. J. Higham, Accuracy and Stability of Numerical Algorithms (provisional title), 1994. In preparation.

[4] J. Ji, F. A. Potra, and S. Huang, A predictor-corrector method for linear complementarity problems with polynomial complexity and superlinear convergence, Technical Report 18, Department of Mathematics, University of Iowa, Iowa City, Iowa, August 1991.

[5] M. Kojima, Y. Kurita, and S. Mizuno, Large–step interior point algorithms for linear complementarity problems, SIAM Journal on Optimization, 3 (1993), pp. 398–412.

[6] M. Kojima, S. Mizuno, and A. Yoshise, An O(\sqrt{n}L) iteration potential reduction algorithm for linear complementarity problems, Mathematical Programming, 50 (1991), pp. 331–342.

[7] I. J. Lustig, R. E. Marsten, and D. F. Shanno, Computational experience with a primaldual interior point method for linear programming, Linear Algebra and Its Applications, 152 (1991), pp. 191–222.

[8] , Computational experience with a globally convergent primal-dual predictor-corrector algorithm for linear programming, Technical Report SOR 92–10, Program in Statistics and Operations Research, Princeton University, Princeton, N. J., 1992.

[9] S. Mehrotra, On the implementation of a primal-dual interior point method, SIAM Journal on Optimization, 2 (1992), pp. 575–601.

[10] R. D. C. Monteiro and I. Adler, Interior path-following primal-dual algorithms. part II: Convex quadratic programming, Mathematical Programming, 44 (1989), pp. 43–66.

[11] D. B. Ponceleón, Barrier methods for large-scale quadratic programming, PhD thesis, Stanford University, 1990.

[12] F. A. Potra, An o(nl) infeasible-interior-point algorithm for LCP with quadratic convergence, Report on Computational Mathematics 50, Department of Mathematics, University of Iowa, Iowa City, Iowa, January 1994.

[13] R. J. Vanderbei, LOQO User’s Manual, Technical Report SOR 92–5, Program in Statistics and Operations Research, Princeton University, Princeton, N. J., 1992.

[14] S. J. Wright, A path-following interior-point algorithm for linear and quadratic optimization problems, Preprint MCS–P401–1293, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill., December 1993.

[15] , An infeasible-interior-point algorithm for linear complementarity problems, Mathematical Programming, 67 (1994), pp. 29–52.

[16] , Stability of linear algebra computations in interior-point methods for linear programming, Preprint MCS–P446–0694, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill., June 1994.

[17] X. Xu, P. Hung, and Y. Ye, A simplified homogeneous and self-dual linear programming algorithm and its implementation. Manuscript, September 1993.

[18] Y. Ye and K. Anstreicher, On quadratic and O(\sqrt{n}L) convergence of a predictor-corrector algorithm for LCP, Mathematical Programming, Series A, 62 (1993), pp. 537–551.

[19] Y. Zhang, On the convergence of a class of infeasible-interior-point methods for the horizontal linear complementarity problem, SIAM Journal on Optimization, 4 (1994), pp. 208–227.