Stability of Augmented System Factorizations in Interior-Point Methods

Stephen Wright^{∗}

July 7, 1995

PREPRINT MCS-P446-0694, MATHEMATICS AND COMPUTER SCIENCE DIVISION, ARGONNE NATIONAL LABORATORY

Abstract

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.

1 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 ∞ 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 [4]). 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.

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

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 surprising. 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 breakdown 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 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 [8]; the review paper of Todd [16] contains a wealth of historical information on potential-reduction methods. Early developments in path-following methods are surveyed by Gonzaga [6], while Mizuno, Todd, and Ye [14] describe an important variant of these methods that does not require the iterates to stay within a cramped neighborhood of the central path. Zhang [24] extended the path-following approach further, allowing the iterates to be infeasible while retaining global convergence and polynomial complexity; see also Wright [20]. 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 [9] generated search directions of the type described in this paper. They compute the maximum step α^{∗} that could be taken along this direction without violating the positivity bounds, then set the actual step length to .995 α^{∗}. Mehrotra’s [13] 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 [13], Fourer and Mehrotra [4], Lustig, Marsten, and Shanno [11], Vanderbei [17], and Xu, Hung, and Ye [22] all implement Mehrotra’s predictor-corrector strategy. These newer codes continue to use step lengths based on α^{∗}; 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ón [15] and Wright [21]. Ponceleón [15] showed that these systems are not too sensitive to structured perturbations from a certain class provided that the underlying optimization problem is well conditioned. Wright [21] analyzed Gaussian elimination in the context of interior-point algorithms for linear complementarity problems.

Simultaneously with the original version of this paper, and independently, Forsgren, Gill, and Shinnerl [3] 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 [18] gives an illuminating discussion of the augmented system in other contexts besides 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 [2] 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:

min c^{T} x,  Ax = b,  x ≥ 0,  (1)

where x ∈ IR^{n} and b ∈ IR^{m}. The dual of (1) is

max b^{T} λ,  A^{T} λ + s = c,  s ≥ 0,  (2)

where s ∈ IR^{n} and λ ∈ IR^{m}. A vector triple (λ^{∗}, x^{∗}, s^{∗}) is a primal-dual solution if x^{∗} is feasible for (1), (λ^{∗}, s^{∗}) is feasible for (2), and s^{∗} and x^{∗} are complementary; that is,

x^{∗}^{T} s^{∗} = c^{T} x^{∗} − b^{T} λ^{∗} = 0.  (3)

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 (x, s) > 0. Search directions are found by applying a modification of Newton’s method to the following system of nonlinear equations:

Ax − b = 0,  A^{T} λ + s − c = 0,  XSe = 0,  (4)

where X = diag(x_{1}, x_{2}, …, x_{n}) and S = diag(s_{1}, s_{2}, …, s_{n}). Specifically, the search direction (Δλ, Δx, Δs) satisfies the linear equations

[ \begin{matrix}0 & A^{T} & I \\ A & 0 & 0 \\ S & 0 & X \end{matrix} ] [ \begin{matrix}Δx \\ Δλ \\ Δs \end{matrix} ] = [ \begin{matrix}−A^{T} λ − s + c \\ b − Ax \\ −XSe + σμe \end{matrix} ] ,  (5)

where σ ∈ [0, 1] is known as the centering parameter and the important quantity μ is defined by

μ = x^{T} s/n.

The step length α along the search direction is determined by various factors; minimally, the updated x and s components are required to stay strictly positive:

(x, s) + α(Δx, Δs) > 0.  (6)

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 α can be quite close to one without violating the property (6), when the search direction (Δλ, Δx, Δs) is an exact solution of (5). If perturbations caused by roundoff are present in the critical components of (Δλ, Δx, Δs), the requirement (6) can severely curtail the allowable step length and slow the convergence. Hence, it is important that the critical components of (Δλ, Δx, Δs) 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 +, −, ×, /, and fl(z) denotes the floating-point approximation of any real number z, we have

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

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(·) in a slightly unconventional way. When ξ and η are two nonnegative numbers, we write ξ = O(η) if there is a positive constant C (not too large) such that ξ ≤ Cη. We say that a matrix or vector is O(η) if its norm is O(η). We say that ξ = Ω(η) if ξ = O(η) and η = O(ξ).

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(·).

If G is a matrix, G_{·j} denotes its j-th column, while G_{i·} denotes the i-th row. The matrix whose elements are |G_{ij}| is denoted by |G|.

We use ∥ · ∥ to denote any one of the equivalent matrix norms ∥ · ∥_{1}, ∥ · ∥_{2}, or ∥ · ∥_{∞}. 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

κ(G) = \frac{svmax (G)}{svmin (G)}.

If G is square and nonsingular, this definition coincides with the usual definition

κ(G) = ∥G∥_{2}∥G^{−1}∥_{2}.

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 Ax = b, A^{T} λ + s = c, (x, s) ≥ 0. Feasibility 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 (x, s) > 0. Then there exists a solution (Δλ, Δx, Δs) to (5).

Proof. The proof follows from Section 6 of Wright [20]. See, in particular, Lemma 6.2, Theorem 6.3, and the remarks in the last two paragraphs of [20]. Note that A need not have full rank for Theorem 3.1 to hold.

The set of basic indices B ⊂ {1, 2, …, n} can be defined as

B = {i | s_{i}^{∗} = 0 for all (λ^{∗}, x^{∗}, s^{∗}) ∈ S},  (8)

while the nonbasic set N is

N = {i | x_{i}^{∗} = 0 for all (λ^{∗}, x^{∗}, s^{∗}) ∈ S}.  (9)

It is well known that B and N form a partition of {1, 2, …, n} and that there is at least one solution (λ^{∗}, x^{∗}, s^{∗}) that is strictly complementary, that is, x^{∗} +s^{∗} > 0 (Goldman and Tucker [5]). The cardinality of B is denoted by |B|. By partitioning the columns of A according to B and N, we define

B = [A_{·j}]_{j∈B},  N = [A_{·j}]_{j∈N} ,  (10)

so that B is m × |B| and N is m × |N |. We say that the linear program is nondegenerate if |B| = m 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:

x_{i} = Ω(1) (i ∈ B),  s_{i} = Ω(1) (i ∈ N),  (11a)
x_{i} = Ω(μ) (i ∈ N),  s_{i} = Ω(μ) (i ∈ B).  (11b)

In addition, the infeasibility is O(μ); that is,

b − Ax = O(μ),  c − A^{T} λ − s = O(μ).  (12)

Assumption 1 is not very strong. Güler and Ye [7] study algorithms in which all iterates are strictly feasible; that is

Ax = b,  A^{T} λ + s = c,  (x, s) > 0.  (13)

In fact they require that x and s be slightly separated from the boundary of the positive orthant, in the sense that

x_{i}s_{i} ≥ γμ,  i = 1, 2, …, n,  (14)

for some constant γ ∈ (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 [19] satisfies Assumption 1. So does the algorithm in [20], provided that the sequence or iterates (x, s) is bounded. Implemented algorithms such as those of Vanderbei [17], Lustig, Marsten, and Shanno [9, 10], and Xu, Hung, and Ye [22] 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

H = [ \begin{matrix}H_{11} & H_{12} \\ H_{21} & H_{22} \end{matrix} ] ,

where H_{11} and H_{22} are also square. Suppose that H_{11} and H_{22}−H_{21}H_{11}^{−1}H_{12} are nonsingular. Then H is nonsingular, and

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

4 Exact and Approximate Search Directions

By defining r_{b} = Ax − b and r_{c} = A^{T} λ + s − c in (5), we obtain

[ \begin{matrix}0 & A & 0 \\ A^{T} & 0 & I \\ 0 & S & X \end{matrix} ] [ \begin{matrix}Δλ \\ Δx \\ Δs \end{matrix} ] = [ \begin{matrix}−r_{b} \\ −r_{c} \\ −XSe + σμe \end{matrix} ] .  (15)

By eliminating Δs from this system, we obtain the augmented system formulation:

[ \begin{matrix}0 & A \\ A^{T} & −X^{−1}S \end{matrix} ] [ \begin{matrix}Δλ \\ Δx \end{matrix} ] = [ \begin{matrix}−r_{b} \\ −r_{c} + s − σμX^{−1}e \end{matrix} ]  (16a)
Δs = −s + σμX^{−1}e − X^{−1}SΔx.  (16b)

In Wright [21], 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 [21] are relevant to the present case of (16), as we discuss later.

Potential difficulties with the formulation (16) arise from two sources — possible rankdeficiency in certain submatrices of A, and the fact that some diagonal elements of X^{−1}S and S^{−1}X approach zero while others approach +∞. 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 interiorpoint 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 α^{∗} to denote the largest number in [0, 1] such that

(x + αΔx, s + αΔs) ≥ 0 for all α ∈ [0, α^{∗}];  (17a)
(x + αΔx)^{T} (s + αΔs)  is decreasing for α ∈ [0, α^{∗}].  (17b)

Theorem 4.1 Suppose that Assumption 1 holds. Let (Δλ, Δx, Δs) be the exact solution of (5) (equivalently, (16)), and let (Δ̂λ, Δ̂x, Δ̂s) be an approximation to this step. Suppose that the centering parameter σ in (5) lies in the range [0, 1/2] and that the following conditions hold:

(Δx, Δs) = O(μ),  (18a)
(Δs_{N}, Δx_{B}) − (Δ̂s_{N}, Δ̂x_{B}) = O(u),  (18b)
(Δs_{B}, Δx_{N}) − (Δ̂s_{B}, Δ̂x_{N}) = O(μu).  (18c)

Define α^{∗} as in (17), and suppose α̂^{∗} is obtained by replacing (Δx, Δs) with (Δ̂x, Δ̂s) in (17). Then for all μ sufficiently small, we have

1 − α^{∗} = O(μ),  (19)
α̂^{∗} = α^{∗} + O(u) = 1 + O(μ) + O(u),  (20)

and

(x + α̂^{∗}Δ̂x)^{T} (s + α̂^{∗}Δ̂s)/n = σO(μ) + O(μ(μ + u)).  (21)

Proof. From (11a) and (18a), we have

s_{N} + αΔs_{N} > 0,  x_{B} + αΔx_{B} > 0,  for all α ∈ [0, 1],

so these components do not restrict the value of α^{∗}. Since u is much smaller than 1, we use (18b) as well to deduce that

s_{N} + αΔ̂s_{N} ≥ s_{N} + αΔs_{N} + α(Δ̂s_{N} − Δs_{N}) > 0, for all α ∈ [0, 1].

Similarly, we can show that x_{B} + αΔ̂x_{B} > 0 for all α ∈ [0, 1].

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 α^{∗} or α̂^{∗}. For the exact direction, we have from (5), (18a), and σ ∈ [0, 1/2] that

\frac{d}{dα} (x + αΔx)^{T} (s + αΔs) = x^{T} Δs + s^{T} Δx + 2αΔx^{T} Δs
≤ −(1 − σ)nμ + 2∥Δx∥∥Δs∥
≤ −nμ/2 + O(μ^{2}),

for all α ∈ [0, 1]. Hence, for μ sufficiently small, the duality gap is decreasing over [0, 1]. For the approximate direction (Δ̂x, Δ̂s), 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

\frac{d}{dα} (x + αΔ̂x)^{T} (s + αΔ̂s) ≤ −nμ/2 + O(μu + μ^{2}).

Again, we find that the duality gap is decreasing over the whole interval α ∈ [0, 1].

Hence, the only condition that can bound α^{∗} and α̂^{∗} away from 1 is (17a), and then only for the N -components of x and the B-components of s. In fact, α^{∗} satisfies

\frac{1}{α^{∗}} = max (1, max_{i∈B} −\frac{Δs_{i}}{s_{i}} , max_{i∈N} −\frac{Δx_{i}}{x_{i}} ) .  (22)

From (5), we have x_{i}Δs_{i} + s_{i}Δx_{i} = −x_{i}s_{i} + σμ. Hence, since x_{i}s_{i} = Ω(μ) from (11), we have

−\frac{Δs_{i}}{s_{i}} = 1 + \frac{Δx_{i}}{x_{i}} − σ \frac{μ}{x_{i}s_{i}} < 1 + \frac{Δx_{i}}{x_{i}} .

For i ∈ B we have from (11a) and (18a) that |Δx_{i}/x_{i}| = O(μ) and therefore

max_{i∈B} −\frac{Δs_{i}}{s_{i}} ≤ 1 + O(μ).

An identical argument can be used for the other term in (22), so we have

\frac{1}{α^{∗}} ≤ max(1, 1 + O(μ)) ⇒ 1 − α^{∗} = O(μ),

proving (19).

For the maximum step length α̂^{∗} along the approximate direction (Δ̂x, Δ̂s), we have from (18c) and (11b) that

\frac{Δ̂s_{i}}{s_{i}} − \frac{Δs_{i}}{s_{i}} = \frac{O(μu)}{Ω(μ)} = O(u), (i ∈ B), \frac{Δ̂x_{i}}{x_{i}} − \frac{Δx_{i}}{x_{i}} = O(u), (i ∈ N).

Hence, from (22), we have

\frac{1}{α̂^{∗}} = max (1, max_{i∈B} −\frac{Δ̂s_{i}}{s_{i}} , max_{i∈N} −\frac{Δ̂x_{i}}{x_{i}} ) = \frac{1}{α^{∗}} + O(u).  (23)

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

(x + αΔ̂x)^{T} (s + αΔ̂s)
= [x + αΔx + α(Δ̂x − Δx)]^{T} [s + αΔs + α(Δ̂s − Δs)]
≤ nμ(1 − α(1 − σ)) + O(μu) + O(μu^{2}),  (24)

where we have used Assumption 1 and (18) to estimate the remainder terms. Finally, we obtain (21) by substituting α = α̂^{∗} = 1 + O(μ + u) into (24).

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 (Δ̂λ, Δ̂x, Δ̂s) 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

LDL^{T} = P T P^{T} ,  (25)

where P is a permutation matrix, L is unit lower triangular, and D is a block-diagonal matrix with 1 × 1 and 2 × 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 ẑ 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 of

[ \begin{matrix}0 & B & N \\ B^{T} & 0 & 0 \\ N^{T} & 0 & Λ \end{matrix} ] + O(μ + u),  (26)

where - μ > 0 and u ≥ 0 are small;

- Λ is diagonal with all diagonal elements of magnitude Ω(μ^{−1});

- B = Ω(1) and κ(B) = O(1); and

- N = O(1). We call T a degenerate canonical matrix if it has the form

[ \begin{matrix}0 & 0 \\ 0 & Λ \end{matrix} ] + O(μ + u),  (27)

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 | N], respectively, and n̄ = m + n 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 Δ is a matrix of the same dimension as T such that

|Δ| ≤ Δ_{u} + |T |δ_{u},  (28)

where δ_{u} and the elements of Δ_{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,

χ_{i} = max{|T_{ij}| | j = 1, 2, …, n̄, j ≠ i}.  (29)

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 ẑ that satisfies

(T + Δ)ẑ = d̂,  (30)

where Δ is a canonical error matrix associated with T and d̂− d = 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_{N}^{−1}e may give rise to errors of similar magnitude.

Theorem 5.1 Suppose that Assumption 1 holds and that the problem is nondegenerate, that is, |B| = m, with κ(B) moderate. Suppose that the procedure for solving (16) satisfies Condition 1, and denote the approximate solution to (16a) by (Δ̂λ, Δ̂x). Then for all sufficiently small μ, we have

(Δλ, Δx, Δs) = O(μ)  (31)

and

(Δλ − Δ̂λ, Δx_{B} − Δ̂x_{B}) = O(u),  Δx_{N} − Δ̂x_{N} = O(μu).  (32)

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 of

[ \begin{matrix}0 & B & N & 0 & 0 \\ B^{T} & 0 & 0 & I & 0 \\ N^{T} & 0 & 0 & 0 & I \\ 0 & S_{B} & 0 & X_{B} & 0 \\ 0 & 0 & S_{N} & 0 & X_{N} \end{matrix} ] .  (33)

Because of (11), the diagonal elements in X_{B} and S_{N} are Ω(1), while the matrices S_{B} and X_{N} 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

(Δλ, Δx, Δs) = O(∥r_{b}∥ + ∥r_{c}∥ + ∥XSe − σμe∥),

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:

[ \begin{matrix}0 & B & N \\ B^{T} & −X_{B}^{−1}S_{B} \\ N^{T} & & −X_{N}^{−1}S_{N} \end{matrix} ] [ \begin{matrix}Δλ \\ Δx_{B} \\ Δx_{N} \end{matrix} ] = [ \begin{matrix}−r_{b} \\ −(r_{c})_{B} + s_{B} − σμX_{B}^{−1}e \\ −(r_{c})_{N} + s_{N} − σμX_{N}^{−1}e \end{matrix} ] .  (34)

From (11), we have for sufficiently small μ that the diagonals in X_{B}^{−1}S_{B} are Ω(μ) while the diagonals of X_{N}^{−1}S_{N} are Ω(μ^{−1}), so this coefficient matrix is canonical.

By defining

M_{B} = [ \begin{matrix}0 & B \\ B^{T} & −X_{B}^{−1}S_{B} \end{matrix} ] , M_{N} = [ \begin{matrix}N \\ 0 \end{matrix} ] ,

Λ = −X_{N}^{−1}S_{N},  z_{N} = Δx_{N},  z_{B} = [ \begin{matrix}Δλ \\ Δx_{B} \end{matrix} ] ,

d_{B} = [ \begin{matrix}−r_{b} \\ −(r_{c})_{B} + s_{B} − σμX_{B}^{−1}e \end{matrix} ] , d_{N} = −(r_{c})_{N} + s_{N} − σμX_{N}^{−1}e,

we can restate the system as

[ \begin{matrix}M_{B} & M_{N} \\ M_{N}^{T} & Λ \end{matrix} ] [ \begin{matrix}z_{B} \\ z_{N} \end{matrix} ] = [ \begin{matrix}d_{B} \\ d_{N} \end{matrix} ] .  (35)

From our assumption on B, we have M_{B} = O(1) and M_{B}^{−1} = O(1).

Because of Condition 1, the computed solution ẑ of (35) satisfies

([ \begin{matrix}M_{B} & M_{N} \\ M_{N}^{T} & Λ \end{matrix} ] + Δ) [ ẑ_{B} ; ẑ_{N} ] = [ d̂_{B} ; d̂_{N} ] ,  (36)

where d̂− d = O(u) and the canonical error matrix Δ satisfies

|Δ| ≤ O(u) + [ \begin{matrix}|M_{B}| & |M_{N}| \\ |M_{N}|^{T} & |Λ| \end{matrix} ] O(u) = O(u) + [ \begin{matrix}0 & 0 \\ 0 & |Λ| \end{matrix} ] O(u).

By combining this estimate with (35) and (36), we obtain

([ \begin{matrix}M_{B} & M_{N} \\ M_{N}^{T} & Λ \end{matrix} ] + Δ) [ \begin{matrix}ẑ_{B} − z_{B} \\ ẑ_{N} − z_{N} \end{matrix} ] = −Δ [ \begin{matrix}z_{B} \\ z_{N} \end{matrix} ] + [ \begin{matrix}d̂_{B} − d_{B} \\ d̂_{N} − d_{N} \end{matrix} ] .  (37)

Since z = O(μ) from (31), we have

|Δ [ \begin{matrix}z_{B} \\ z_{N} \end{matrix} ]| ≤ (O(u) + [ \begin{matrix}0 & 0 \\ 0 & |Λ| \end{matrix} ] O(u)) [ \begin{matrix}O(μ) \\ O(μ) \end{matrix} ] ≤ [ \begin{matrix}O(μu) \\ O(u) \end{matrix} ] ,  (38)

so when we add the effect of d̂− 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

H_{11} = M_{B} + O(u),
H_{12} = M_{N} + O(u) = O(1),
H_{21} = M_{N}^{T} + O(u) = O(1),
H_{22} = O(u) + Λ(I + O(u)) = Λ(I + O(u)).

Lemma 3.2 yields the following estimates:

(H^{−1})_{22} = (H_{22} − H_{21}H_{11}^{−1}H_{12})^{−1} = Λ^{−1}(I + O(u + μ)) = O(μ),
(H^{−1})_{12} = O(μ),  (H^{−1})_{21} = O(μ),
(H^{−1})_{11} = M_{B}^{−1} + O(μ + u).

By combining these observations with (38), we obtain

[ \begin{matrix}ẑ_{B} − z_{B} \\ ẑ_{N} − z_{N} \end{matrix} ] = [ \begin{matrix}(H^{−1})_{11} & (H^{−1})_{12} \\ (H^{−1})_{21} & (H^{−1})_{22} \end{matrix} ] O(u) = [ \begin{matrix}O(u) \\ O(μu) \end{matrix} ] ,

giving (32).

Next, we examine the accuracy of Δ̂s, which is calculated by substituting Δ̂λ and Δ̂x into (16b).

Theorem 5.2 Suppose that the assumptions of Theorem 5.1 are satisfied and that Δ̂s is evaluated in floating-point arithmetic from the formula (16b), with Δ̂x replacing Δx. We then have

Δs_{B} − Δ̂s_{B} = O(μu),  (39a)
Δs_{N} − Δ̂s_{N} = O(u).  (39b)

Proof. Standard roundoff error analysis applied to (16b) shows that

Δ̂s = −s + σμX^{−1}e − X^{−1}SΔ̂x + [|s| + σμ|X^{−1}|e + |X^{−1}S||Δ̂x|] O(u).  (40)

By differencing (16b) and (40), we obtain

|Δs − Δ̂s| ≤ |X^{−1}S||Δx − Δ̂x| + [|s| + σμ|X^{−1}|e + |X^{−1}S||Δ̂x|] O(u).  (41)

If i ∈ B, we have from (11) that

|X_{B}^{−1}S_{B}| = O(μ),  |s_{B}| = O(μ),  σμ|X_{B}^{−1}|e = O(μ).

By combining these estimates with (32) and (41), we obtain the desired result (39a). For i ∈ N, we have from (11) again that

|X_{N}^{−1}S_{N}| = O(μ^{−1}),  |s_{N}| = O(1),  σμ|X_{N}^{−1}|e = O(1),

while from (31) and (32) we have

Δx_{N} − Δ̂x_{N} = O(μu),  Δ̂x_{N} = O(μ).

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 nondegenerate, that is |B| = m with κ(B) moderate. Suppose that the procedure for solving (16) satisfies Condition 1. If the approximate step is computed with σ ∈ [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 [1] is as follows.

Choose δ ∈ (0, 1); find r such that χ_{1} = |T_{r1}|; if χ_{1} > 0

if |T_{11}| ≥ δχ_{1}

1 × 1 pivot, P_{1} = I

else find χ_{r}; if χ_{r}|T_{11}| ≥ δχ_{1}^{2}

1 × 1 pivot, P_{1} = I

elseif |T_{rr}| ≥ δχ_{r}

1 × 1 pivot; choose P_{1} so that (P_{1}T P_{1}^{T})_{11} = T_{rr}

else

2 × 2 pivot; choose P_{1} so that (P_{1}T P_{1}^{T})_{21} = T_{r1} end

end end.

If we denote the 1 × 1 or 2 × 2 pivot block by E and write

P_{1}T P_{1}^{T} = [ \begin{matrix}E & C^{T} \\ C & T̂ \end{matrix} ] ,  (42)

the first step of the factorization yields

P_{1}T P_{1}^{T} = \left[{I \atop CE^{−1} I}\right] [ E ; T̄ ] \left[{I E^{−1}C^{T} \atop I}\right] ,  (43)

where T̄ = T̂ − CE^{−1}C. The algorithm continues by applying this procedure to T̄. Note that the χ_{i} are generally changed by each stage of the factorization. The submatrix CE^{−1} contains the subdiagonals in the first one or two columns of the L factor.

Bunch and Kaufman [1] show that for the particular choice δ = (1 + \sqrt{17})/8, we have

max_{i,j} |T̄_{ij}| ≤ (2.57) max_{i,j} |T_{ij}|,  (44)

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 × 1 block, chosen from among the diagonal elements of Λ; or

(ii) a 2 × 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 |C| = O(1), while

(i) |E| = O(μ^{−1}) and |E^{−1}| = O(μ) if E is a 1 × 1 pivot; and

(ii) |E| = O(1) and |E^{−1}| = O(1) if E is a 2 × 2 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 × 1 block, chosen from any of the diagonals (large or small); or

(ii) a 2 × 2 block, in which all the elements are O(μ + u);

(b) The matrix remaining after the elimination is canonical (not necessarily degenerate), and the absolute change to the remaining matrix is O(μ + u).

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 × 1 pivot, the dimensions of B are unchanged, while a 2 × 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 Λ + O(μ + 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

L̂D̂L̂^{T} = P T P^{T} + P Δ̄P^{T} ,  (45)

where Δ̄ is a canonical error matrix associated with T.

Proof. We prove the result by an induction argument on the dimension n̄ = m + n 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 m ≤ n.

For n̄ = 1, we must have m = 0 and so trivially P = 1, L̂ = 1, D̂ = T_{11}. Therefore (45) holds with Δ̄ = 0.

For n̄ = 2, we have two cases m = 0 and m = 1. For m = 0, there are two elements of magnitude Ω(μ^{−1}) on the diagonal, while the off-diagonals are O(μ + u). Hence, a 1 × 1 pivot is chosen. If there is no pivoting, the first step of elimination yields

L̂_{21} = T_{21}/T_{11} + |T_{21}/T_{11}|O(u),
D̂_{11} = T_{11},
D̂_{22} = T_{22} − T_{21}^{2}/T_{11} + (|T_{22}| + |T_{21}^{2}/T_{11}|)O(u).

Since L̂ has unit diagonals, we obtain by expanding the factors that

L̂D̂L̂^{T} = T + [ \begin{matrix}0 & |T_{21}|O(u) \\ |T_{21}|O(u) & |T_{21}^{2}/T_{11}|O(u) + |T_{22}|O(u) \end{matrix} ] = T + Δ̄,

where

|Δ̄| ≤ |T |O(u) + O(u),

so Δ̄ is a canonical error matrix associated with T. The same logic applies if pivoting occurs.

In the remaining case m = 1, the pivot is 2 × 2, we have L̂ = I, P = I, and D̂ = T, and (45) holds trivially with Δ̄ = 0.

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 L̂_{1} and D̂_{1}, where

L̂_{1} = [ \begin{matrix}I & 0 \\ CE^{−1} + Δ_{L} & I \end{matrix} ] ,  D̂_{1} = [ \begin{matrix}E & 0 \\ 0 & T̄ + Δ_{D} \end{matrix} ] ,  (46)

where

|Δ_{L}| ≤ |C||E^{−1}|O(u), |Δ_{D}| ≤ |T̂|O(u) + |C||E^{−1}||C|^{T} O(u) = |T̂|O(u) + O(u).

Note that Δ_{D} is a canonical error matrix corresponding to T̂. 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

L̂_{2}D̂_{2}L̂_{2}^{T} = P_{2}(T̄ + Δ_{D})P_{2}^{T} + P_{2}Δ̄_{2}P_{2}^{T}  (47)

for some permutation matrix P_{2} and some canonical error matrix Δ̄_{2} corresponding to (T̄ + Δ_{D}). We compose the overall factors of T as follows:

L̂ = [ \begin{matrix}I & 0 \\ P_{2}(CE^{−1} + Δ_{L}) & L̂_{2} \end{matrix} ] ,  D̂ = [ \begin{matrix}E & 0 \\ 0 & D̂_{2} \end{matrix} ] ,  P = [ \begin{matrix}I & 0 \\ 0 & P_{2} \end{matrix} ] P_{1}.

Now,

L̂D̂L̂^{T} = [ \begin{matrix}E & (C + Δ_{2})^{T} P_{2}^{T} \\ P_{2}(C + Δ_{2}) & L̂_{2}D̂_{2}L̂_{2}^{T} + P_{2}CE^{−1}C^{T} P_{2}^{T} + P_{2}Δ_{1}P_{2}^{T} \end{matrix} ] ,  (48)

where

Δ_{1} = Δ_{L}C^{T} + CΔ_{L}^{T} + Δ_{L}EΔ_{L}^{T}, Δ_{2} = Δ_{L}E,

and so

|Δ_{1}| ≤ |C||E^{−1}||C|^{T} O(u) = O(u), |Δ_{2}| ≤ |C||E^{−1}||E|O(u) = O(u).

By substituting (47) and (46) into (48), we obtain

L̂D̂L̂^{T} = [ \begin{matrix}E & (C + Δ_{2})^{T} P_{2}^{T} \\ P_{2}(C + Δ_{2}) & P_{2} [T̄ + Δ_{D} + CE^{−1}C^{T} + Δ_{1} + Δ̄_{2}] P_{2}^{T} \end{matrix} ]

= [ \begin{matrix}E & (C + Δ_{2})^{T} P_{2}^{T} \\ P_{2}(C + Δ_{2}) & P_{2} [T̂ + Δ_{D} + Δ_{1} + Δ̄_{2}] P_{2}^{T} \end{matrix} ]
= P T P^{T} + P Δ̄P^{T} ,

where

Δ̄ = P_{1}^{T} [ \begin{matrix}0 & Δ_{2}^{T} \\ Δ_{2} & Δ_{D} + Δ_{1} + Δ̄_{2} \end{matrix} ] P_{1}.

Since |Δ_{1}| = O(u), |Δ_{2}| = O(u), and Δ_{D} and Δ̄_{2} are canonical error matrices corresponding to T̂, we have

|Δ̄| ≤ P_{1}^{T} [ \begin{matrix}0 & |Δ_{2}|^{T} \\ |Δ_{2}| & |Δ_{D}| + |Δ_{1}| + |Δ̄_{2}| \end{matrix} ] P_{1}
≤ O(u) + P_{1}^{T} [ \begin{matrix}0 & 0 \\ 0 & |T̂| \end{matrix} ] P_{1}O(u)
≤ O(u) + |T |O(u).

Hence, Δ̄ is a canonical error matrix corresponding to T.

We complete the proof by noting that Theorem 6.1 can be applied to the remaining matrix, 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 ẑ is found by performing two vector permutations with P, triangular substitutions with L̂ and L̂^{T}, and a blockwise inversion of D̂. The 2 × 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 × 2 linear system Ey = g in which E is symmetric with

|E_{11}| ≤ δ|E_{12}|,  |E_{11}||E_{22}| ≤ δ^{2}|E_{12}|^{2},  (49)

for some δ ∈ (0, 1). Then if we compute the solution by applying Gaussian elimination to the permuted system

[ \begin{matrix}E_{12} & E_{22} \\ E_{11} & E_{12} \end{matrix} ] [ \begin{matrix}y_{1} \\ y_{2} \end{matrix} ] = [ \begin{matrix}g_{2} \\ g_{1} \end{matrix} ] ,  (50)

then the computed solution ŷ satisfies

(E + Δ_{E})ŷ = g,

where

|Δ_{E}| ≤ |E|O(u).  (51)

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 ẑ to the system L̂D̂L̂^{T} z = P d satisfies

(L̂D̂L̂^{T} + P Δ̂P^{T})ẑ = P d,  (52)

where Δ̂ is a canonical error matrix associated with T.

Proof. From standard results for triangular substitution, the computed solution of L̂z_{a} = P d satisfies

(L̂ + Δ̂_{L1})ẑ_{a} = P d, |Δ̂_{L1}| ≤ |L̂|O(u).

A similar result holds for triangular substitution with the transpose L̂^{T}.

For solution of D̂z_{b} = ẑ_{a}, we note that D̂ is block-diagonal with 1 × 1 and 2 × 2 blocks. For the 2 × 2 pivot blocks that arise in the Bunch-Kaufman procedure, the assumptions of Lemma 6.4 hold, so the computed solution ŷ of a 2 × 2 subsystem Ey = g satisfies

(E + Δ_{E})ŷ = g,  |Δ_{E}| = |E|O(u).  (53)

When E is a 1 × 1 block, the estimate (53) holds trivially. Hence, the computed solution ẑ_{b} of D̂z_{b} = ẑ_{a} satisfies

(D̂ + Δ̂_{D})ẑ_{b} = ẑ_{a}, |Δ̂_{D}| ≤ |D̂|O(u).

By combining the error expressions for the three component systems, we find that our computed solution ẑ satisfies

(L̂ + Δ̂_{L1})(D̂ + Δ̂_{D})(L̂ + Δ̂_{L2})^{T} ẑ = P d.

Multiplying the matrix products, we find that (52) is satisfied with

P |Δ̂|P^{T} ≤ |L̂||D̂||L̂|^{T} O(u) + O(μu + u^{2}).

From our earlier discussions on the composition of L̂ and D̂, it is easy to see that the absolute matrix product |L̂||D̂||L̂|^{T} contains all O(1) elements, except for the large diagonals, which occur in the same positions as in P T P^{T}. Hence P Δ̂P^{T} is a canonical error matrix corresponding to P T P^{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 ẑ to T z = d satisfies

(L̂D̂L̂^{T} + P Δ̂P^{T})ẑ = d̂,

Substituting from (45), we obtain

(P T P^{T} + P Δ̄P^{T} + P Δ̂P^{T})ẑ = P d̂,

so Condition 1 follows when we set Δ = Δ̄ + Δ̂.

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 (Δ̂λ, Δ̂x, Δ̂s) 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 σ is small or zero, a large reduction in the duality gap μ can be expected. In the extreme case σ = 0 (the “affinescaling” choice), linear convergence with a rate constant of O(u) can be attained if the actual step length is close to α̂^{∗}. Most practical algorithms choose the step length to be a fixed multiple — typically .95 or .9995 — of α̂^{∗}, 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 [21, Section 4], for instance, extra restrictions are applied to α to ensure that (12) and (14) continue to hold at the next iterate. These restrictions may result in α being much smaller than one. This case is analyzed in [21, Section 4], so we do not repeat it here.

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 δ ∈ (0, 1), χ_{off} = |T_{rs}| = max_{i≠j} |T_{ij}|, χ_{diag} = |T_{pp}| = max_{i} |T_{ii}|; if χ_{diag} ≥ δχ_{off}

s = 1 and choose P_{1} so that (P_{1}T P_{1}^{T})_{11} = T_{pp}

else

s = 2 and choose P_{1} so that (P_{1}T P_{1}^{T})_{21} = 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 δ = (1 + \sqrt{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 × 1 pivots;

2. 2 × 2 pivots of the type described in Theorem 6.1(a) are chosen;

3. When no more 2 × 2 pivots like this are available and the remaining matrix contains only elements of size O(μ + u), a combination of small 1 × 1 and 2 × 2 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 χ_{off} = O(1) and χ_{diag} = Ω(μ^{−1}). The pivot element will therefore be one of the large diagonals corresponding to Λ. The remaining matrix is updated by subtracting CE^{−1}C, where clearly C = O(1) and E^{−1} = 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

[ \begin{matrix}0 & B \\ B^{T} & 0 \end{matrix} ] + O(μ + u).  (54)

Stage 2 now begins. If B is not vacuous, we have χ_{off} = O(1) and χ_{diag} = O(μ + u). In fact, by the assumption B = Ω(1), we have χ_{off} = Ω(1), and the element T_{rs} that achieves the maximum comes from B. The 2 × 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 × 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 |rows(B) − columns(B)|, and all its elements have size O(μ+u). In stage 3, both 1×1 and 2×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 [4]) 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 [4], with a slight modification noted below.

For each index i = 1, 2, …, n̄ 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

n̂_{ij} = min(n_{i} + n_{j} − 4, n̄ − 2).

A 2 × 2 pivot block

E = [ \begin{matrix}T_{ii} & T_{ij} \\ T_{ij} & T_{jj} \end{matrix} ]  (55)

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

oxo: (n_{i} − 1)(n_{j} − 1),
tile: (n_{i} − 1)(n̂_{ij} + 1) if T_{ii} = 0, (n_{j} − 1)(n̂_{ij} + 1) if T_{jj} = 0,
full: n̂_{ij}^{2},

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 δ ∈ (0, 1) and the off-diagonal norms χ_{i} defined in (29). Any 1 × 1 pivot must satisfy

|T_{ii}^{−1}|χ_{i} ≤ 2/δ,  (56)

while a 2 × 2 pivot (55) must have

|[ \begin{matrix}T_{ii} & T_{ij} \\ T_{ij} & T_{jj} \end{matrix} ]^{−1}| [ \begin{matrix}χ_{i} \\ χ_{j} \end{matrix} ] ≤ [ \begin{matrix}1/δ \\ 1/δ \end{matrix} ] .  (57)

The pivot selection procedure is as follows. for r = 1, 2, …

for i with n_{i} = r consider T_{ii} with degree r; if any of these elements satisfy (56)

accept as a 1 × 1 pivot and exit; else label it as unstable;

end

for unstable pivots T_{ii} from the previous loop

consider 2 × 2 pivots involving T_{ii}, with costs at most (r − 1)^{2}, (r − 1)(2r − 3), and (2r − 4)^{2} for oxo, tile, and full pivots, respectively;

if any of these blocks satisfy (57)

accept as a 2 × 2 pivot and exit;

end 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 × 1 pivots. In the description of [4], the right-hand side is 1/δ rather than 2/δ. With the original choice, the sparse Bunch-Parlett algorithm applied to a degenerate canonical matrix could allow another type of pivot: a 2 × 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 — |B| ≠ m — 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 × 1 or 2 × 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_{N}^{−1}S_{N} are used as 1 × 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

[ \begin{matrix}NX_{N}S_{N}^{−1}N^{T} & B \\ B^{T} & −X_{B}^{−1}S_{B} \end{matrix} ] [ \begin{matrix}Δλ \\ Δx_{B} \end{matrix} ] = [ \begin{matrix}−r_{b} + NS_{N}^{−1}X_{N}[−(r_{c})_{N} + s_{N} − σμX_{N}^{−1}e] \\ −(r_{c})_{B} + s_{B} − σμX_{B}^{−1}e \end{matrix} ] ,

(58)

Δx_{N} = X_{N}^{−1}S_{N} [(r_{c})_{N} − s_{N} + σμX_{N}^{−1}e + N^{T} Δλ] .  (59)

The coefficient matrix in (58) is an O(μ) perturbation of the matrix

[ \begin{matrix}0 & B \\ B^{T} & 0 \end{matrix} ] .  (60)

Since B is well conditioned by Definition 2, the matrix in (60) has 2 min(|B|, m) nonzero singular values of magnitude Ω(1). In the nondegenerate case, (60) is well conditioned. Otherwise, it has |m − |B|| zero singular values. When |B| < m, the null space of (60) is spanned by

[ Z̄ ; 0 ] ,  (61)

where Z̄ is an m × (m − |B|) matrix of full rank such that B^{T} Z̄ = 0. When |B| > m, the null space of (60) is spanned by the matrix

[ 0 ; Ẑ ] ,  (62)

where Ẑ 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 |B| < m, the perturbations occur mostly in the range space of the matrix (61), that is, in the components of Δλ. Similarly, when |B| > m, the perturbations occur in the range space of the matrix (61), that is, in the components of Δx_{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 |B| > 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 Δx_{B}. The large relative errors in Δx_{B} induce large relative errors in Δs_{B} through the formula (16b). The step length to the boundary α^{∗} may therefore be sharply curtailed because of the nonnegativity requirements (17a). In the case |B| < m, the large relative errors in Δλ induce errors in Δx_{N} through the formula (59), while in turn induce large relative errors in Δs_{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 r = 0 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 [4] 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 [23]) or crossover to the simplex method (Meggido [12]) 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 infeasible-interior-point path-following algorithm described in Wright [20]. In exact arithmetic, this algorithm achieves superlinear convergence because it eventually always takes affine-scaling steps (σ = 0 in (5)) with step length α 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_{B}^{∗}, and s_{N}^{∗}. The matrix A is dense and random, with elements defined by

A_{1j} = τ10^{6τ−3},  j = 1, …, n,
A_{ij} = (τ − .5)10^{6τ−3},  i = 2, …, m, j = 1, …, n,

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

N = {1, 2, …, |N |}, B = {1, 2, …, n}\N .

We choose a particular solution (λ^{∗}, x^{∗}, s^{∗}) by setting

λ^{∗} = e,  s_{B}^{∗} = 0,  x_{N}^{∗} = 0,
s_{i}^{∗} = 10^{4τ−2}, i ∈ N , x_{i}^{∗} = 10^{3τ−1}, i ∈ B,

where each τ is as before. The vectors b and c are determined by the choices of A and (λ^{∗}, x^{∗}, s^{∗}).

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 arithmetic, giving u ≈ 10^{−14} on the SPARC-5 on which these results were obtained.

We report on problems with m = 6, n = 12. (In problems smaller than this, exactly zero pivots often occur in degenerate cases, leading to breakdown.) Termination occurs when μ ≤ 10^{−30} — an artificially stringent criterion, chosen to give us a clear look at asymptotic effects.

The first result is for a nondegenerate problem, for which |B| = m = 6. Table 1 shows the sizes of μ and ∥r∥ on each iterate. For the reasons that we outlined immediately following Condition 1, ∥r∥ 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^{−10}. 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 17, when μ ≈ 10^{−7}. As expected, there are six 1 × 1 pivots of magnitude Ω(μ^{−1}), and six 2 × 2 pivots in which the diagonals are tiny and the off-diagonals are Ω(1). The same structure is present in D at every iteration after iteration 15.

Our second example is for a dual degenerate problem with |B| = 6 > m. 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 Table 1: Nondegenerate problem: m = 6, n = 12

k log_{10} μ_{k} log_{10} ∥r^{k}∥_{1} Affine Step?
1  5.4  3.1
2  4.7  2.3
3  4.3  1.6
4  3.8  0.8
5  3.1  -12.0
⋮  ⋮  ⋮
15  -3.2  -14.0
16  -4.6  -13.7  ∗
17  -7.2  -14.4  ∗
18  -12.3  -14.1  ∗
19  -22.1  -13.8  ∗
20  -33.3  -14.2 termination

Table 2: The D factor at iteration 17 of the nondegenerate test problem (∗ = magnitude less than 10^{−6})

Row/Column  Pivot Block
1,2 ∗ .94(1)
.94(1)  ∗
3,4 ∗ -.91(2)
-.91(2)  ∗
5  .26(7)
6  .30(11)
7  .33(10)
8  .47(7)
9,10  -.30(-5) .71(2)
.71(2)  ∗
11,12  -.27(-3) -.15(2)
-.15(2)  ∗
13,14 ∗ -.31(0)
-.31(0) -.49(-5)
15,16 ∗ .16(0)
.16(0)  ∗
17  .27(4)
18  .32(6)

Table 3: Dual degenerate problem: m = 6, n = 12, |B| = 8

k log_{10} μ_{k} log_{10} ∥r^{k}∥_{1} Affine Step?
1  5.4  3.1
⋮  ⋮  ⋮
19  -6.0  -13.8
20  -9.8  -14.1  ∗
21  -13.6  -14.2  ∗
22  -14.8  -13.8  ∗
23  -15.4  -13.2  ∗
⋮  ⋮  ⋮
99  -17.5  -13.5
100  -17.5  -13.4
⋮  ⋮  ⋮

that there are exactly min(m, |B|) of the stable 2 × 2 pivots with an off-diagonal from B, and |N | = n − |B| of the large 1 × 1 pivots. Together, these stable pivots account for

2 min(m, |B|) + |N | = n + m − |m − |B||  (63)

stages of the factorization, so unstable pivots are used on the remaining submatrix whose dimension is |m − |B||. In Table 4, we see that the last two 1 × 1 pivots are unstable, as expected. As we described in the first part of Section 9, the errors in Δ̂x_{B} and Δ̂s_{B} are preventing further progress. On iteration 100, the computed affine step has ∥Δ̂x_{B}∥_{∞} = .17(6), while its exact counterpart would have ∥Δx_{B}∥_{∞} = O(μ). By comparing components of Δ̂s_{B} with s_{B}, we find that the step to the boundary is sharply curtailed by the restriction s_{B} + αΔ̂s_{B} ≥ 0 (cf. (23)). The remaining components of the step do not contain deleterious errors; we have

∥Δ̂x_{N}∥_{∞} = .59(−18),  ∥Δ̂λ∥_{∞} = .66(−14),  ∥Δ̂s_{N}∥_{∞} = .11(−12).

Finally, we consider a primal degenerate problem with |B| = 4 < 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 × 1 blocks, which again matches the prediction (63). As discussed in Section 9, the deleterious errors occur in the subvector Δ̂λ, so errors are induced in Δ̂s_{N} and Δ̂x_{N} through formulas 59 and (16b). On iteration 100, we have ∥Δ̂λ∥_{∞} = .32(5) and ∥Δ̂s_{N}∥_{∞} = .30(7) for the affine scaling step. The components Δ̂x_{B} and Δ̂s_{B} are not affected; their ∞-norms are .17(−18) and .51(−12), respectively.

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:

Table 4: The D factor at iteration 17 of the degenerate test problem with m = 6, n = 12, |B| = 8 (∗ = magnitude less than 10^{−6})

Row/Column  Pivot Block
1,2 ∗ .95(1)
.95(1)  ∗
3,4 ∗ -.92(2)
-.92(2)  ∗
5,6 ∗ .26(2)
.26(2)  ∗
7  .86(23)
8  .85(18)
9  .55(20)
10  .29(17)
11,12 ∗ .71(2)
.71(2)  ∗
13,14 ∗ -.30(0)
-.30(0)  ∗
15,16 ∗ .15(0)
.15(0)  ∗
17  .20(-13)
18  -.60(-19)

Table 5: Primal degenerate problem: m = 6, n = 12, |B| = 4

k log_{10} μ_{k} log_{10} ∥r^{k}∥_{1} Affine Step?
1  5.4  3.1
⋮  ⋮  ⋮
15  -5.3  -13.9
16  -8.8  -13.7  ∗
17  -13.7  -14.2  ∗
18  -14.0  -11.6  ∗
⋮  ⋮  ⋮
99  -17.6  -13.9
100  -17.6  -14.0
⋮  ⋮  ⋮

Table 6: The D factor at iteration 17 of the degenerate test problem with m = 6, n = 12, |B| = 8 (∗ =magnitude less than 10^{−6})

Row/Column  Pivot Block
1,2 ∗ .95(1)
.95(1)  ∗
3  .49(23)
4  .53(19)
5  .58(19)
6  .27(20)
7  .53(9)
8  .12(21)
9,10 ∗ .71(2)
.71(2)  ∗
11,12 ∗ -.15(2)
-.15(2)  ∗
13  .25(18)
14  .76(17)
15,16 ∗ -.16(1)
-.16(1)  ∗
17  -.15(-8)
18  -.52(-18)

(iii) The pivot is 1 × 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

χ_{1} ≤ \frac{|T_{11}|}{δ}  or  χ_{1} ≤ \sqrt{\frac{χ_{r}|T_{11}|}{δ}}.  (64)

Now, since χ_{r} is the maximum off-diagonal in some column of (26), we have χ_{r} = O(1), while since T_{11} comes from either the (1, 1) or (2, 2) block of (26), we have |T_{11}| = O(μ + u). Since δ ∈ (0, 1) is fixed, we have from (64) that

χ_{1} = O(μ^{1/2} + u^{1/2}).  (65)

Since χ_{1} is the magnitude of the largest off-diagonal in some row/column of (26), we have that χ_{1} is the ∞-norm of some row or column of B. But (65) is incompatible with B = O(1) and κ(B) = O(1). Hence |T_{11}| 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 × 2 and involves at least one element from Λ. Since all the off-diagonals in (26) are O(1), the quantities χ_{i}, i = 1, 2, …, n̄ are all O(1). A 2 × 2 pivot with diagonal elements T_{11} and T_{rr} must have

|T_{11}| ≤ δχ_{1}, |T_{rr}| ≤ δχ_{r},

which implies that T_{11} and T_{rr} are both O(1). Since all the diagonals of Λ are Ω(μ^{−1}), they cannot be candidates for T_{11} and T_{rr}.

(v) The pivot is 2 × 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 |T_{1r}| = χ_{1} — is O(μ + u). Since T_{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 = Ω(1) and κ(B) = O(1).

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

[ \frac{1}{P̃} ] [ \begin{matrix}(Λ + O(μ + u))_{ii} & N_{·i}^{T} & 0 & 0 \\ N_{·i} & 0 & B & Ñ \\ 0 & B^{T} & 0 & 0 \\ 0 & Ñ^{T} & 0 & Λ̃ \end{matrix} ] [ \frac{1}{P̃^{T}} ] + O(μ + u),

where P̃ is some permutation matrix, N_{·i} denotes the i-th column of N, Ñ is obtained from N by removing N_{·i}, and Λ̃ is obtained from Λ by removing its i-th row and column. Since |(Λ + O(μ + u))_{ii}^{−1}| = O(μ), the submatrix that remains after elimination is

P̃ [ \begin{matrix}0 & B & Ñ \\ B^{T} & 0 & 0 \\ Ñ^{T} & 0 & Λ̃ \end{matrix} ] P̃^{T} − P̃ [ \begin{matrix}N_{·i} \\ 0 \\ 0 \end{matrix} ] Λ_{ii}^{−1} [ N_{·i}^{T} 0 0 ] P̃^{T} + O(μ + u)

= P̃ [ \begin{matrix}0 & B & Ñ \\ B^{T} & 0 & 0 \\ Ñ^{T} & 0 & Λ̃ \end{matrix} ] P̃^{T} + O(μ + u).  (66)

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

[ \frac{I}{P̂} ] [ \begin{matrix}0 & B_{ij} & 0 & B_{i·;j} & N_{i·} \\ B_{ij} & 0 & B_{·j;i}^{T} & 0 & 0 \\ 0 & B_{·j;i} & 0 & B̂ & N̂ \\ B_{i·;j}^{T} & 0 & B̂^{T} & 0 & 0 \\ N_{i·}^{T} & 0 & N̂^{T} & 0 & Λ \end{matrix} ] [ \frac{I}{P̂^{T}} ] + O(μ + u),

where - P̂ is some permutation matrix;

- N_{i·} is the i-th row of N;

- N̂ is N with N_{i·} removed;

- B_{i·;j} is the i-th row of B, with its j-th element removed;

- B_{·j;i} is the j-th column of B, with its i-th element removed;

- B̂ is B with 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 |B_{ij}| = Ω(1). Denoting the pivot block by E, we have

E = B_{ij} [ \begin{matrix}0 & 1 \\ 1 & 0 \end{matrix} ] + O(μ + u),  E^{−1} = \frac{1}{B_{ij}} [ \begin{matrix}0 & 1 \\ 1 & 0 \end{matrix} ] + O(μ + u).  (67)

Therefore the elimination step yields the remaining matrix

P̂ [ \begin{matrix}0 & B̂ & N̂ \\ B̂^{T} & 0 & 0 \\ N̂^{T} & 0 & Λ \end{matrix} ] P̂^{T} − \frac{1}{B_{ij}} P̂ [ \begin{matrix}0 & B_{·j;i} \\ B_{i·;j}^{T} & 0 \\ N_{i·}^{T} & 0 \end{matrix} ] [ \begin{matrix}0 & 1 \\ 1 & 0 \end{matrix} ] [ \begin{matrix}0 & B_{i·;j} & N_{i·} \\ B_{·j;i}^{T} & 0 & 0 \end{matrix} ] P̂^{T} + O(μ + u)

= P̂ [ \begin{matrix}0 & B̂ & N̂ \\ B̂^{T} & 0 & 0 \\ N̂^{T} & 0 & Λ \end{matrix} ] P̂^{T} + O(μ + u),  (68)

where

B̄ = B̂ − \frac{1}{B_{ij}} B_{·j;i}B_{i·;j},  N̄ = N̂ − \frac{1}{B_{ij}} B_{·j;i}N_{i·}.

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 κ(B̄) = O(1).

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 C = O(1) whether the pivot block is 1 × 1 or 2 × 2. For 1 × 1 pivots, we have |E| = Ω(μ^{−1}) and |E^{−1}| = Ω(μ). For 2 × 2 pivots, we have from (67) and |B_{ij}| = Ω(1) that |E| = O(1) and |E^{−1}| = O(1).

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 × 2 and contains at least one element from Λ. In a degenerate canonical matrix, we have χ_{i} = O(μ + u), i = 1, 2, …, n̄. A 2 × 2 pivot with diagonal elements T_{11} and T_{rr} must have

|T_{11}| ≤ δχ_{1}, |T_{rr}| ≤ δχ_{r},

which implies that both diagonals are O(μ + u), so neither element can come from Λ.

In the case of either a 1 × 1 or 2 × 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 |C| = O(μ + u) and |E| = Ω(μ^{−1}). Hence, the update to the remaining submatrix is bounded by

|C||E^{−1}||C|^{T} = O(μ(μ + u)^{2}),

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:

[ \begin{matrix}1 & 0 \\ E_{11}/E_{12} + δ_{1} & 1 \end{matrix} ] ,  [ \begin{matrix}E_{12} & E_{22} \\ 0 & E_{12} − E_{11}E_{22}/E_{12} + δ_{2} \end{matrix} ] ,  (69)

where

δ_{1} = |\frac{E_{11}}{E_{12}} | O(u), δ_{2} = |E_{12}|O(u) + |E_{11}E_{22}/E_{12}| O(u).

It is well known that for triangular substitution applied to any triangular system Uz = h, the computed solution ẑ satisfies (U + Δ_{U})ẑ = h, where |E_{U}| = |U|O(u). By applying this observation to each of the matrices in (69), we find that the computed solution ŷ of (50) satisfies

[ \begin{matrix}1 & 0 \\ E_{11}/E_{12} + δ_{3} & 1 \end{matrix} ] [ \begin{matrix}E_{12} + δ_{4} & E_{22} + δ_{5} \\ 0 & E_{12} − E_{11}E_{22}/E_{12} + δ_{6} \end{matrix} ] [ ŷ_{1} ; ŷ_{2} ] = [ \begin{matrix}g_{2} \\ g_{1} \end{matrix} ] , (70)

where

δ_{3} = δ_{1} + |E_{11}/E_{12}| O(u) = |E_{11}/E_{12}| O(u),
δ_{4} = |E_{12}|O(u),
δ_{5} = |E_{22}|O(u),
δ_{6} = δ_{2} + (|E_{12}| + |E_{11}E_{22}/E_{12}|) O(u) = |E_{12}|O(u) + |E_{11}E_{22}/E_{12}| O(u).

By multiplying out the coefficient matrix in (70), we obtain

[ \begin{matrix}E_{12} + δ_{4} & E_{22} + δ_{5} \\ E_{11} + δ_{7} & E_{12} + δ_{8} \end{matrix} ] ,  (71)

where

δ_{7} = |E_{12}|δ_{3} + |E_{11}/E_{12}| δ_{4} = |E_{11}|O(u),
δ_{8} = |E_{11}/E_{12}| δ_{5} + |E_{22}|δ_{3} + δ_{+} (|E_{12}| + |E_{11}E_{22}/E_{12}|) O(u)
= |E_{11}|O(u) + |E_{11}E_{22}/E_{12}| O(u) + |E_{12}|O(u).
= |E_{12}|O(u).

(The last equality follows from (49).) Hence, (71) can be written as

[ \begin{matrix}0 & 1 \\ 1 & 0 \end{matrix} ] (E + Δ_{E}),

where Δ_{E} 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 × 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 χ_{i} = O(μ + u), which implies that one of the rows or columns of B is O(μ + u). However, this estimate is incompatible with B = Ω(1) and κ(B) = O(1), so this kind of pivot cannot occur.

(iv) The pivot is 2 × 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 × 1 pivot at an earlier point in the algorithm. But if it was considered, it would have been accepted, since

|T_{ii}^{−1}|χ_{i} = O(μ)O(1) = O(μ) ≤ 2/δ

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

T_{ii} = Ω(μ^{−1}),  T_{jj} = O(μ + u),  T_{ij} = O(1),

and so

|[ \begin{matrix}T_{ii} & T_{ij} \\ T_{ij} & T_{jj} \end{matrix} ]^{−1}| = \frac{1}{|T_{ii}T_{jj} − T_{ij}^{2}|} |[ \begin{matrix}T_{jj} & −T_{ij} \\ −T_{ij} & T_{ii} \end{matrix} ]| .

|[ \begin{matrix}T_{jj} & −T_{ij} \\ −T_{ij} & T_{ii} \end{matrix} ]| [ \begin{matrix}χ_{i} \\ χ_{j} \end{matrix} ] ≤ |T_{ii}T_{jj} − T_{ij}^{2}| [ \begin{matrix}1/δ \\ 1/δ \end{matrix} ] = [ \begin{matrix}O(1) \\ O(1) \end{matrix} ] .

From the second row of this inequality, we have

χ_{j} ≤ \frac{1}{|T_{ii}|}O(1) = O(μ).

But χ_{j} is the ∞-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 × 2, and the pivot block E is drawn either entirely from the (1, 1) block 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

|[ \begin{matrix}T_{jj} & −T_{ij} \\ −T_{ij} & T_{ii} \end{matrix} ]| [ \begin{matrix}χ_{i} \\ χ_{j} \end{matrix} ] ≤ |T_{ii}T_{jj} − T_{ij}^{2}|O(1).

Taking the second row of this relation, we obtain

|T_{ij}|χ_{i} + |T_{ii}|χ_{j} ≤ |T_{ii}T_{jj} − T_{ij}^{2}|O(1) ≤ (|T_{ii}T_{jj}| + |T_{ij}|^{2}) O(1),  (72)

where, by definition, χ_{i} and χ_{j} are both nonnegative. Consider two cases. When |T_{ij}|^{2} ≥ |T_{ii}T_{jj}| we have from (72) that

|T_{ij}|χ_{i} ≤ |T_{ij}|^{2}O(1) ⟹ χ_{i} = O(|T_{ij}|) = O(μ + u).

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 |T_{ij}|^{2} < |T_{ii}T_{jj}|, we have

|T_{ii}|χ_{j} ≤ |T_{ii}T_{jj}|O(1) ⟹ χ_{j} = O(|T_{jj}|) = O(μ + u),

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 × 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 × 2 block in which both diagonals are from Λ, then one of them (T_{ii}, say) must have been considered previously as a 1 × 1 pivot. But if it was considered, it would have been accepted, since

|T_{ii}^{−1}|χ_{i} = O(μ)O(μ + u) ≤ 2/δ.

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 × 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 × 1 pivot and rejected. Because of (56), T_{ii} must satisfy

|T_{ii}| < \frac{δ}{2}χ_{i}.  (73)

On the other hand, since the 2 × 2 pivot is accepted, we must have

|[ \begin{matrix}T_{jj} & −T_{ij} \\ −T_{ij} & T_{ii} \end{matrix} ]| [ \begin{matrix}χ_{i} \\ χ_{j} \end{matrix} ] ≤ |T_{ii}T_{jj} − T_{ij}^{2}| [ \begin{matrix}1/δ \\ 1/δ \end{matrix} ] .  (74)

Consider first the case of T_{ij}^{2} ≥ |T_{ii}T_{jj}|. Then from the first block row in (74), this inequality implies that

|T_{jj}|χ_{i} ≤ |T_{ii}T_{jj} − T_{ij}^{2}|\frac{1}{δ} ≤ 2T_{ij}^{2} \frac{1}{δ}.

Since |T_{ij}| ≤ χ_{i}, we have

|T_{jj}| ≤ 2|T_{ij}|\frac{1}{δ} = O(μ + u),

which contradicts our assumption that T_{jj} has size Ω(μ^{−1}). The remaining case has T_{ij}^{2} < |T_{ii}T_{jj}|. From (74) and (73), we have

|T_{jj}|χ_{i} ≤ 2|T_{ii}T_{jj}|\frac{1}{δ} < 2|T_{jj}|\frac{1}{δ} \frac{δ}{2}χ_{i} = |T_{jj}|χ_{i},

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

|E^{−1}C^{T} | ≤ |E^{−1}||C^{T} | = O(1/δ) = O(1).

Hence, the update matrix CE^{−1}C^{T} is bounded as follows:

|CE^{−1}C^{T} | = ∥C∥O(1) = O(μ + u),

giving the result.

References

[1] J. Bunch and L. Kaufman, Some stable methods for calculating inertia and solving symmetric linear systems, Mathematics of Computation, 31 (1977), pp. 163–179.

[2] I. S. Duff, The solution of augmented systems, Technical Report RAL–93–084, Rutherford Appleton Laboratory, Oxon, U. K., November 1993.

[3] A. Forsgren, P. Gill, and J. Shinnerl, Stability of symmetric ill-conditioned systems arising in interior methods for constrained optimization, Report TRITA-MAT-1994-24, Royal Institute of Technology, June 1994.

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

[5] A. J. Goldman and A. W. Tucker, Theory of linear programming, in Linear Equalities and Related Systems, H. W. Kuhn and A. W. Tucker, eds., Primceton University Press, Princeton, N. J., 1956, pp. 53–97.

[6] C. Gonzaga, Path-following methods in linear programming, SIAM Review, 34 (1991), pp. 167–224.

[7] O. Güler and Y. Ye, Convergence behavior of interior-point algorithms, Mathematical Programming, (1993), pp. 215–228.

[8] 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.

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

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

[11] , Interior-point methods for linear programming: Computational state of the art, ORSA Journal on Computing, 6 (1994), pp. 1–14.

[12] N. Megiddo, On finding primal- and dual-optimal bases, ORSA Journal on Computing, 3 (1991), pp. 63–65.

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

[14] S. Mizuno, M. Todd, and Y. Ye, On adaptive step primal-dual interior-point algorithms for linear programming, Mathematics of Operations Research, 18 (1993), pp. 964– 981.

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

[16] M. Todd, Potential reduction methods in mathematical programming, Technical Report 1112, Cornell University, Ithaca, N.Y., 14853-3801, 1995.

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

[18] S. Vavasis, Stable numerical algorithms for equilibrium systems, SIAM Journal on Matrix Analysis and Applications, 15 (1994), pp. 1108–1131.

[19] S. J. Wright, A path-following infeasible-interior-point algorithm for linear complementarity problems, Optimization Methods and Software, 2 (1993), pp. 79–106.

[20] , 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. (to appear in Annals of Operations Research).

[21] , Stability of linear equations solvers in interior-point methods, Preprint MCS– P400–1293, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill., December 1993. (to appear in SIAM Journal on Matrix Analysis and Applications).

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

[23] Y. Ye, On the finite convergence of interior-point algorithms for linear programming, Tech. Rep. 91–5, Department of Management Sciences, University of Iowa, Iowa City, February 1991.

[24] 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.