MODIFIED CHOLESKY FACTORIZATIONS IN INTERIOR-POINT ALGORITHMS FOR LINEAR PROGRAMMING

STEPHEN WRIGHT^{∗}

PREPRINT ANL/MCS-P600-0596, MAY, 1996, MATHEMATICS AND COMPUTER SCIENCE DIVISION, AR-GONNE NATIONAL LABORATORY, ARGONNE, IL 60439, USA.

Abstract. We investigate a modified Cholesky algorithm similar to those used in current interior-point codes for linear programming. Cholesky-based interior-point codes are popular for three reasons: their implementation requires only minimal changes to standard sparse Cholesky codes (allowing us to take full advantage of software written by specialists in that area); they tend to be more efficient than competing approaches that use different factorizations; and they perform robustly on most practical problems, yielding good interior-point steps even when the coefficient matrix is ill conditioned. We explain the surprisingly good performance of the Cholesky-based approach by using analytical tools from matrix perturbation theory and error analysis, illustrating our results with computational experiments. Finally, we point out the limitations of this approach.

Key words. Interior-point algorithms and software, Cholesky factorization, Matrix perturbations, Error analysis.

1. Introduction. Most interior-point codes for linear programming share a common feature: their major computational operation—solution of a large linear system of equations—is performed by a direct sparse Cholesky algorithm. In this algorithm, row and column orderings are determined a priori by well-known heuristics (minimum degree and enhancements, minimum local fill, nested dissection) that are based solely on the sparsity pattern and not on the numerical values of the nonzero elements. The ordering phase is followed by a symbolic factorization phase, in which the nonzero structure of the Cholesky factor is determined and storage is allocated. Finally, a numerical factorization phase fills in the numerical values of the lower triangular Cholesky factor. In interior-point codes, the first two phases usually are performed just once, during either the first interior-point iteration or computation of a starting point.

In the interior-point context, the unadorned Cholesky algorithm can run into difficulties because of extreme ill conditioning. Some of the diagonal pivots encountered during the numerical factorization phase can be zero or negative, causing the standard Cholesky procedure to break down. Instead of crashing, most codes apply a “patch” to the algorithm to handle such pivots. The offending pivot element is sometimes replaced by a huge number, as in LIPSOL [17] or PCx [1]. In other codes such as IPMOS [16], the pivot is replaced by a moderate number, but the corresponding right-hand side element is set to zero, as are the off-diagonal elements in the corresponding column of the Cholesky factor. The first practical interior-point code, OB1 [6], explicitly zeroes the components of the solution vector that correspond to small pivots. All these strategies are essentially equivalent to the algorithm we describe in this paper. To date, there has been little investigation of them from a numerical analysis viewpoint.

The “patches” described above have the advantage that they can be implemented by changing just a few lines in general sparse Cholesky codes. It is therefore possible to take advantage of the long-term development effort that has gone into designing such codes and their underlying algorithms. The recent codes LIPSOL [17] and PCx

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

[1] make explicit use of the freely available sparse Cholesky code of Ng and Peyton [8]. Other codes either modify the well-known SPARSPAK routines of George and Liu [3] or include customized linear algebra routines that implement well established algorithmic ideas. (At least one author has experimented with modifications to the standard heuristics: Mészáros [7] describes an inexact version of the minimum local fill ordering.)

One possible remedy for small pivots is diagonal pivoting. At each iteration, a “large” diagonal element is selected from the unreduced portion of the matrix and moved to the pivot position by symmetric row and column pivoting. The algorithm is terminated when none of the remaining diagonal elements is sufficiently large, and an approximate solution is computed with the partial factors. (See Higham [4, Chapter 10] for details and error analysis.) This strategy is not particularly appealing in the context of interior-point linear programming codes because of the loss of efficiency due to shifting of data during the numerical factorization. Moreover, there is little incentive to test this strategy because the simple patches described above perform so well in practice.

In this article, we use standard results from numerical analysis to explain the good performance of these patching strategies on the vast majority of problems. We also gain some insight into their limitations and into how and why they fail.

Our error analysis for the modified Cholesky algorithm is rigorous, with explicitly stated assumptions and precise bounds (see Sections 3 and 4). We revert, however, to a more informal style when applying these results to the interior-point context (Section 5). The reason is pure pragmatism. A fully rigorous analysis would be impossibly complex, notationally speaking, and unduly pessimistic. The informal analysis yields adequate insight into the typical performance of the algorithm, as our computational results in Section 6 demonstrate.

A number of other papers on linear algebra operations in barrier and interior-point methods have appeared in recent years. Wright [12] has considered the Newtonlogarithmic barrier method for general constrained optimization, in which the linear system to be solved for the Newton step is positive semidefinite and ill conditioned during later iterations. She uses a Cholesky factorization with diagonal pivoting to identify the subspace spanned by the active constraint Jacobian. From this information, an accurate solution of the Newton equations can be obtained, in which the components of the step in both the range space of the active constraint Jacobian and the null space of its transpose are well resolved. Our analysis has a similar flavor to Wright’s, but the application is somewhat different. The unknowns in our linear system are the unconstrained dual variables rather than the primals and, since this problem is linear, we have little interest in resolving the component of the step in the near-null space of the coefficient matrix. We focus too on Cholesky algorithms that perform no pivoting during the numerical factorization, reflecting computational practice in the current generation of interior-point linear programming codes.

In an earlier paper [14], we considered the stability of algorithms for the symmetric indefinite form of the step equations at each iteration of a interior-point method for linear programming. We showed that, despite the ill-conditioning of the coefficient matrix, the steps obtained by this approach are good search directions for the interiorpoint method. Forsgren, Gill, and Shinnerl [2] perform a similar analysis in the context of logarithmic barrier methods.

The remainder of this paper is organized as follows. In Section 2, we introduce primal-dual interior-point methods and derive the linear equations to be solved at each iteration of these methods. Section 3 introduces Algorithm modchol, the modified Cholesky procedure, and examines the accuracy of the solution obtained with this factorization, under certain assumptions on the eigenvalues of the factored matrix. In Section 4, we account for the effect of finite-precision floating-point arithmetic on solution accuracy. We return to the interior-point application in Section 5, showing that Algorithm modchol yields good steps for these methods until the duality gap becomes very small, even if the linear program is primal or dual degenerate. The analytical results are verified by computational experiments with an interior-point code using Algorithm modchol, which are reported in Section 6.

Notation. We summarize here the notation used in the remainder of the paper. The ith singular value of a matrix A is denoted by σ_{i}(A). We use σ_{i} alone to

denote the ith singular value of the exact Cholesky factor L in Section 3.

For any matrix M and index steps I and J, M_{IJ} denotes the submatrix formed by the elements M_{ij}, for i ∈ I and j ∈ J. The ith column of M is denoted by M_{·i}, and the column submatrix consisting of columns j ∈ J is denoted by M_{·J}.

Unit roundoff error is denoted by u. Higham [4, Chapter 1] defines u implicitly by the statement that when α and ζ are any two floating-point numbers, op denotes +, −, ×, and /, and fl(·) denotes the floating-point representation of a real number, we have

fl(α op ζ) = (α op ζ)(1 + δ) for some δ satisfying |δ| ≤ u.

For any positive integer m with mu < 1, we define

γ_{m} = \frac{mu}{1 − mu} (1)

(see Higham [4, Lemma 2.1]).

The notation ∥ · ∥ denotes the Euclidean vector norm ∥ · ∥_{2} and also its induced matrix norm, unless otherwise noted. For any matrix A, the matrix consisting of the absolute values of each element is denoted by |A|. We use 1 to denote the vector (1, 1, ⋯, 1)^{T}.

Finally, we mention the parameter ϵ that defines the pivot threshold in the modified Cholesky algorithm. A second quantity ϵ̄, which is related to ϵ by

ϵ̄ =^{def} 2m^{2}ϵ,

appears frequently in the analysis because the incorporation of the scaling term 2m^{2} saves notational clutter.

2. Primal-Dual Algorithms for Linear Programming. We consider the linear programming problem in standard form:

(2)  min c^{T} x subject to Ax = b,  x ≥ 0,

where x ∈ IR^{n}, c ∈ IR^{n}, A ∈ IR^{m×n}, and b ∈ IR^{m}. The dual of (2) is

(3)  max b^{T} π subject to A^{T} π + s = c,  s ≥ 0,

where s ∈ IR^{n} and π ∈ IR^{m}. We assume throughout the paper that A has full row rank, so that m ≤ n. The Karush-Kuhn-Tucker (KKT) conditions, which identify a vector triple (x, π, s) as a primal-dual solution for (2), (3), can be stated as follows:

(4a)  A^{T} π + s = c,
(4b)  Ax = b,
(4c)  x_{i}s_{i} = 0, i = 1, 2, …, n,
(4d)  (x, s) ≥ 0.

We assume throughout the paper that a primal-dual solution exists. We make no assumptions about uniqueness or nondegeneracy; our analysis in Section 5 continues to hold when the problem (2) is primal or dual degenerate. It is well known that the index set {1, 2, …, n} can be partitioned into two sets B and N such that for all primal-dual solutions (x^{∗}, π^{∗}, s^{∗}) we have

(5)  x_{i}^{∗} = 0 for all i ∈ N ,  s_{i}^{∗} = 0 for all i ∈ B.

Primal-dual interior-point algorithms generate a sequence of iterates (x, π, s) that satisfy the strict inequality (x, s) > 0. They find search directions by applying a modification of Newton’s method to the system of nonlinear equations formed by the first three KKT conditions (4a),(4b),(4c), namely,

Ax − b = 0,  A^{T} π + s − c = 0,  XS1 = 0, (6)

where X = diag(x_{1}, x_{2}, …, x_{n}), S = diag(s_{1}, s_{2}, …, s_{n}), and 1 = (1, 1, …, 1)^{T}. In general, the search direction (Δx, Δπ, Δs) is obtained from the following linear system:

[ \begin{matrix}0 & A^{T} & I \\ A & 0 & 0 \\ S & 0 & X \end{matrix} ] [ \begin{matrix}Δx \\ Δπ \\ Δs \end{matrix} ] = [ \begin{matrix}−r_{c} \\ −r_{b} \\ −r_{xs} \end{matrix} ] , (7)

where the coefficient matrix is the Jacobian of (6) and the right-hand side components r_{b} and r_{c} are defined by

r_{b} = Ax − b, r_{c} = A^{T} π + s − c. (8)

In a pure Newton (affine-scaling) method, the remaining right-hand side component r_{xs} is defined by

r_{xs} = XS1, (9)

and, in this case, we denote the solution of (7) by (Δx^{aff}, Δπ^{aff}, Δs^{aff}). In a pathfollowing method, we have

r_{xs} = XS1 − ζμ1, (10)

where μ is the duality gap defined by

μ = x^{T} s/n, (11)

and ζ ∈ [0, 1] is a centering parameter. In the “Mehrotra predictor-corrector” algorithm, which is used as the basis of many practical codes, the search direction is calculated by setting

r_{xs} = XS1 + ΔX^{aff}ΔS^{aff}1 − ζμ1, (12) where ΔX^{aff} and ΔS^{aff} are the diagonal matrices formed from the affine-scaling step components Δx^{aff} and Δs^{aff}. Hence, Mehrotra’s method requires the solution of two linear systems at each iteration—the affine scaling system (7), (8), (9), and the search direction system (7), (8), (12). A heuristic based on the effectiveness of the affine scaling direction is used to determine the value of ζ in (12).

Once a search direction has been determined, the primal-dual algorithm takes a step of the form

(x, π, s) + α(Δx, Δπ, Δs),

where α is chosen to maintain strict positivity of the x and s components; that is,

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

In most codes, α is chosen to be some fraction of the step-to-boundary α_{max} defined as

(14) α_{max} = sup_{α∈[0,1]}{α | (x, s) + α(Δx, Δs) ≥ 0}.

A typical strategy is to set

α = ηα_{max},

where η ∈ [.9, 1.0) approaches 1 as the interior-point method approaches the solution set.

By applying block elimination to (7) and using the notation

D^{2} = S^{−1}X, (15)

we obtain the following equivalent system:

(16a)  AD^{2}A^{T} Δπ = −r_{b} + AD^{2}(r_{c} − X^{−1}r_{xs}),
(16b)  Δs = −r_{c} − A^{T} Δπ,
(16c)  Δx = −S^{−1}(r_{xs} + XΔs).

In many codes, the solution is obtained from just this formulation. A sparse Cholesky factorization, modified to handle small pivots, is applied to the symmetric positive definite coefficient matrix AD^{2}A^{T} in (16a) and the solution Δπ is obtained by triangular substitution with the computed factor. The remaining direction components are recovered from (16b) and (16c). This technique yields steps (Δx, Δπ, Δs) that are useful search directions for the interior-point algorithm, even when the matrix AD^{2}A^{T} is ill conditioned, as often happens during later iterations. This observation is somewhat surprising, since a naive application of error analysis results would suggest that the combination of ill-conditioning and roundoff would corrupt the direction hopelessly. The results of Sections 3, 4, and 5 provide an explanation for this phenomenon.

The following observation is crucial to our analysis: In computing Δπ from (16a), we are not interested so much in the error in Δπ itself as in the effect of this error on the remaining step components Δs and Δx that are recovered from (16b) and (16c), respectively. If the relative errors in these components are large, the positivity requirement (13) may cause the step length α to be significantly shortened, thereby curtailing the algorithm’s progress. We return to this issue in Section 5, after describing and analyzing the modified Cholesky algorithm in Sections 3 and 4.

3. A Modified Cholesky Algorithm. In this section, we describe and analyze Algorithm modchol, a modified Cholesky algorithm designed to handle illconditioned matrices for which small or negative pivots may arise during the factorization.

Algorithm modchol accepts an m × m symmetric positive definite matrix M as input, together with a small positive user-defined parameter ϵ, which defines a threshold of acceptability for the pivot elements. If a candidate pivot element is smaller than this threshold, the algorithm simply skips a step of factorization. Algorithm modchol outputs an approximate lower triangular factor L̃ and an index set J ⊂ {1, 2, …, m} containing the indices of the skipped pivots. In the following specification, we use M^{(i)} to denote the unfactored part of M that remains after i steps of the algorithm.

Algorithm modchol

Given ϵ with 0 < ϵ ≪ 1; Set M^{(0)} ← M; L̃ ← 0; J ← ∅; β = max_{i=1,2,...,m} M_{ii}; for i = 1, 2, …, m

if M_{ii}^{(i−1)} ≤ βϵ (* skip this elimination step *) Set J ← J ∪ {i} and

E^{(i)} = [ \begin{matrix}0 & 0 & ⋯ & ⋯ & 0 \\ 0 & M_{ii}^{(i−1)} & ⋯ & ⋯ & M_{im}^{(i−1)} \\ ⋮ & ⋮ & 0 & & 0 \\ ⋮ & ⋮ & & ... & ⋮ \\ 0 & M_{mi}^{(i−1)} & 0 & ⋯ & 0 \end{matrix} ] , M^{(i)} = M^{(i−1)} − E^{(i)}; (17)

else (* perform the usual Cholesky elimination step *)

L̃_{ii} ← \sqrt{M_{ii}^{(i−1)}}; M^{(i)} ← 0 for j = i + 1, i + 2, …, m

L̃_{ji} = M_{ij}^{(i−1)}/L̃_{ii} ;

for j = i + 1, i + 2, …, m

for k = i + 1, i + 2, …, m
M_{jk}^{(i)} ← M_{jk}^{(i−1)} − L̃_{ji}L̃_{ki}.

The ith column of L̃ is zero for each i ∈ J; that is, L̃_{·J} = 0. If we denote

E = ∑_{i∈J} E^{(i)} (18)

and denote the complement of J in {1, 2, …, m} by J̄, it follows from (17) that

EJ̄J̄ = 0. (19)

That is, the row or column index of each nonzero element in E must lie in J. It follows from the algorithm that L̃ is the exact Cholesky factor of the perturbed matrix M −E, which we denote for convenience by M̃. That is, we have

L̃L̃^{T} = M̃ = M − E. (20)

By partitioning this equation into its J and J̄ components and using L̃_{·J} = 0 and (19), we obtain

(21a)  MJ̄J̄ = L̃_{J̄·}L̃_{J̄·}^{T} + EJ̄J̄ = L̃J̄J̄L̃_{J̄J̄}^{T}
(21b)  M_{J̄J} = L̃_{J̄·}L̃_{J ·}^{T} + E_{J̄J} = L̃J̄J̄L̃_{J J̄}^{T} + E_{J̄J} .

The exact Cholesky factor L (whose existence is guaranteed by the assumed positive definiteness of M) satisfies

LL^{T} = M. (22)

Given the linear system

Mz = r, (23)

where M is the matrix factored by modchol, the exact solution obviously satisfies

z = M^{−1}r = L^{−T} L^{−1}r. (24)

The approximate solution z̃ is chosen so that the partial vector z̃J̄ solves the reduced system MJ̄J̄z̃J̄ = rJ̄, while the complementary subvector z_{J} is set to zero. From (21a), we see that z̃J̄ can be calculated by performing a pair of triangular substitutions; that is,

z̃J̄ = L̃_{J̄J̄}^{−T}L̃_{J̄J̄}^{−1}rJ̄,  z̃_{J} = 0. (25)

Note that z = z̃ when J = ∅. When J ≠ 0, on the other hand, the difference between z̃ and z can be large in a relative sense. We have

∥z − z̃∥ = ∥[ \begin{matrix}z_{J} − 0 \\ zJ̄ − z̃J̄ \end{matrix} ]∥ ≥ ∥z_{J} ∥,

and there is no reason to expect z_{J} to be small with respect to the full vector z. We can show, however, that the difference between L^{T} z and L^{T} z̃ is relatively small under certain assumptions; this result is the culmination of the analysis of this section (Theorem 3.6). As we see in Section 5, this difference determines the usefulness of the computed solution of (16) as a search direction for the interior-point algorithm.

To simplify the analysis, we assume implicitly throughout the paper that

β = 1. (26)

A trivial scaling, which affects neither the algorithm nor its analysis, can always be applied to the symmetric positive definite matrix M to yield (26).

We start with a sequence of three results that lead to a bound on the difference between L̃^{T} z and L̃^{T} z̃. These results require few assumptions on the matrix M and are relatively simple to prove.

Lemma 3.1. The submatrix formed by the last m − i rows and columns of M^{(i)} is symmetric positive definite, for all i = 0, 1, …, m − 1. Moreover, the diagonal elements of all these submatrices are bounded by 1.

Proof. This observation follows by a simple inductive argument. By assumption, the starting matrix M^{(0)} = M is positive definite. Suppose that the desired property holds for M^{(i−1)}. If i ∈ J, then the lower right (m − i) × (m − i) submatrix of M^{(i)} is identical to the lower right (m − i) × (m − i) submatrix of M^{(i−1)}, which is positive definite by assumption. Otherwise, if i ∉ J, then M^{(i)} is obtained by applying one step of Cholesky reduction to M^{(i−1)}. It is known that the remaining submatrix resulting from this operation is positive definite; hence, the lower right (m − i) × (m − i) submatrix in question is positive definite, and the desired property holds.

The second claim follows immediately from the fact that M_{ii} ≤ β = 1, i = 1, 2, …, m and the fact that the diagonal elements cannot increase during Algorithm modchol.

Lemma 3.2. For each i ∈ J, we have

∥E^{(i)}∥_{2} ≤ ∥E^{(i)}∥_{F} ≤ (2mϵ)^{1/2}.

Therefore,

∥E∥_{2} ≤ ∥E∥_{F} ≤ ϵ̄^{1/2}, (27)

where ϵ̄ = 2m^{2}ϵ.

Proof. From Lemma 3.1, we have (M_{i,l}^{(i−1)})^{2} ≤ M_{i,i}^{(i−1)}M_{l,l}^{(i−1)} for each l = i +

1, …, m. Suppose i ∈ J, so that M_{i,i}^{(i−1)} ≤ ϵ. Since the diagonals of each submatrix M^{(i−1)} are bounded by 1, we have M_{l,l}^{(i−1)} ≤ 1 and therefore

|M_{i,l}^{(i−1)}| ≤ (M_{i,i}^{(i−1)}M_{l,l}^{(i−1)})^{1/2} ≤ ϵ^{1/2}, l = i + 1, …, m.

Hence, we have

∥E^{(i)}∥_{2}^{2} ≤ ∥E^{(i)}∥_{F}^{2} ≤ (M_{i,i}^{(i−1)})^{2} + 2 ∑_{l=i+1}^{m}(M_{i,l}^{(i−1)})^{2} ≤ ϵ^{2} + 2(m − i)ϵ ≤ 2mϵ,

thereby proving the first claim. By (18), we have

∥E∥_{F}^{2} = ∑_{i∈J} ∥E^{(i)}∥_{F}^{2} ≤ |J |2mϵ ≤ 2m^{2}ϵ,

thereby proving (27).

In the case in which all the small pivots appear in the bottom right corner of the matrix (that is, J = {p + 1, p + 2, …, m} for some index p), the estimate (27) can be improved to

∥E∥_{2} ≤ ∥E∥_{F} ≤ m^{2}ϵ = .5ϵ̄, (28)

This stronger estimate applies in most instances of the interior-point application of Section 5.

We are now able to derive an estimate of the difference between L̃^{T} z and L̃^{T} z̃. Theorem 3.3. For the exact solution z and approximate solution z̃ defined in

(24) and (25), respectively, we have that

∥L̃^{T} [z − z̃]∥ = ∥L̃_{J̄J̄}^{−1}E_{J̄J} z_{J} ∥ ≤ ∥L̃_{J̄J̄}^{−1}∥ ∥E_{J̄J} ∥ ∥z_{J} ∥. (29)

Proof. From (24) together with (21), we have

rJ̄ = MJ̄J̄zJ̄ + M_{J̄J} z_{J}
= L̃J̄J̄L̃_{J̄J̄}^{T}zJ̄ + [L̃J̄J̄L̃_{J J̄}^{T} + E_{J̄J} ] z_{J}
= L̃J̄J̄L̃_{·J̄}^{T}z + E_{J̄J} z_{J} ,

while from (25), we have

rJ̄ = L̃J̄J̄L̃_{J̄J̄}^{T}z̃J̄ = L̃J̄J̄ [L̃_{J̄J̄}^{T}z̃J̄ + L̃_{J J̄}z̃_{J} ] = L̃J̄J̄L̃_{·J̄}^{T}z̃.

By combining these two relations, we obtain

L̃_{·J̄}^{T}[z − z̃] = −L̃_{J̄J̄}^{−1}E_{J̄J} z_{J} . (30)

Since L̃_{·J} = 0, the result follows immediately.

The remaining analysis of this section requires some additional assumptions on the distribution of the singular values of M and on the parameter ϵ. Accordingly, we introduce a little more notation. The eigenvalues of M are denoted by σ_{i}^{2}, i = 1, 2, …, m, where

σ_{1}^{2} ≥ σ_{2}^{2} ≥ ⋯ ≥ σ_{m}^{2} > 0. (31)

We define the diagonal matrix Σ by

Σ = diag(σ_{1}, σ_{2}, …, σ_{m}). (32)

It follows that there exists an orthogonal matrix Q such that

M = QΣ^{2}Q^{T} . (33)

Because the largest diagonal in M is 1, we have by elementary analysis that

1 ≤ σ_{1}^{2} ≤ m. (34)

In the subsequent analysis, we assume that there is an integer p with 1 ≤ p ≤ m such that

• ϵ is small relative to σ_{p}^{2}; and

• if p < m, there is a significant gap in the spectrum of M between σ_{p}^{2} and

σ_{p+1}^{2}. (We will be more specific about these two assumptions presently.) By partitioning the spectrum at the gap, we obtain

(35)  Σ_{1} = diag(σ_{1}, σ_{2}, …, σ_{p}),  Σ_{2} = diag(σ_{p+1}, σ_{p+2}, …, σ_{m}).

From (33), Q can be partitioned accordingly to obtain

Q = [Q_{1} | Q_{2}], M = Q_{1}Σ_{1}^{2}Q_{1}^{T} + Q_{2}Σ_{2}^{2}Q_{2}^{T} .

Since M = LL^{T}, it follows that σ_{i}, i = 1, 2, …, m are the singular values of L. In fact, we must have

L^{T} = UΣQ^{T} (36) for some m × m orthogonal matrix U, where Σ and Q are defined as above.

We use σ̃_{i}^{2}, i = 1, 2, …, m to denote the eigenvalues of the perturbed matrix M̃. It follows immediately from (20) that the singular values of L̃ are σ̃_{i}, i = 1, 2, …, m. The rank of L̃ is |J̄|, because L̃J̄J̄ is lower triangular with nonzero diagonals while L̃_{·J} = 0. Therefore, we have

σ̃_{|J̄|} > σ̃_{|J̄|+1} = ⋯ = σ̃_{m} = 0. (37)

As in (36), there are orthogonal m × m matrices Ũ and Q̃ such that

L̃^{T} = ŨΣ̃Q̃^{T} , (38)

where

Σ̃ = diag(σ̃_{1}, σ̃_{2}, …, σ̃_{m}).

It is an immediate consequence of an eigenvalue perturbation result of Stewart and Sun [10, Corollary IV.4.13] and Lemma 3.2 that

{∑_{i=1}^{m}[σ_{i}^{2} − σ̃_{i}^{2}]^{2}}^{1/2} ≤ ∥E∥_{F} = ϵ̄^{1/2}. (39)

The main assumption of this section is that |J̄| = p; that is, Algorithm modchol correctly identifies the numerical rank of the matrix M. One might expect that we should not have to assume this equality at all—that it should follow from the spectrum gap and from a judicious choice of ϵ. Practical experience supports this expectation; the algorithm has little trouble determining the numerical rank on the vast majority of problems. In fact, part of the result—the bound |J̄| ≥ p—follows from a minimal assumption on ϵ.

Lemma 3.4. If ϵ̄^{1/2} < σ_{p}^{2}, we have |J̄| ≥ p.

Proof. If |J̄| < p, we have from (37) and (39) that

σ_{p}^{2} ≤ σ_{|J̄|+1}^{2} = |σ_{|J̄|+1}^{2} − σ̃_{|J̄|+1}^{2}| ≤ ϵ̄^{1/2},

contradicting our assumption that ϵ̄^{1/2} < σ_{p}^{2}.

However, the conditions on ϵ, σ_{p}, and σ_{p+1} needed to prove the other half of the result—|J̄| ≤ p—are too rigorous to be useful. This is a consequence of the fact that poorly conditioned triangular matrices need not have particularly small diagonal elements (see Lawson and Hanson [5, p. 31] for the classic example of this phenomenon).

Our next result concerns perturbation of the subspace spanned by Q_{1}, which is the invariant subspace of “large” eigenvalues of M.

Lemma 3.5. Suppose that |J̄| = p < m and that the values σ_{p} and σ_{p+1} from (31) and ϵ from Lemma 3.2 satisfy the conditions

(40a)  \frac{σ_{p+1}^{2}}{σ_{p}^{2}} ≤ .1,
(40b)  σ_{p}^{2} − σ_{p+1}^{2} > 5ϵ̄^{1/2}.

Then there is a p × p symmetric positive definite matrix Λ̃ and an orthonormal m × p matrix Q̃_{1} such that

(41)  M̃ = Q̃_{1}Λ̃Q̃_{1}^{T} ,
(42)  ∥Q̃_{1} − Q_{1}∥ ≤ \frac{2ϵ̄^{1/2}}{σ_{p}^{2} − σ_{p+1}^{2} − 2ϵ̄^{1/2}} ,
(43)  ∥Λ̃ − Σ_{1}^{2}∥ ≤ 2ϵ̄^{1/2}.

(The constants used in (40a) and in similar expressions should not be taken too seriously. We assign them specific values only to avoid an excess of notation.)

Proof. The result is a straightforward consequence of Theorem V.2.8 of Stewart and Sun [10, p. 238]. Since M̃ = M −E, we use (33) and partition as in (35) to obtain

Q^{T} M̃Q = Q^{T} MQ − Q^{T} EQ = [ \begin{matrix}Σ_{1}^{2} & 0 \\ 0 & Σ_{2}^{2} \end{matrix} ] − [ \begin{matrix}F_{11} & F_{12} \\ F_{12}^{T} & F_{22} \end{matrix} ] .

We now make the following identifications with the quantities in the cited result:

γ̃ = ∥F_{12}^{T}∥ ≤ ∥F ∥ = ∥E∥ ≤ ϵ̄^{1/2},  η̃ = ∥F_{12}∥ ≤ ϵ̄^{1/2},
δ̃ = sep(Σ_{1}^{2}, Σ_{2}^{2}) − ∥F_{11}∥ − ∥F_{22}∥ ≥ σ_{p}^{2} − σ_{p+1}^{2} − 2ϵ̄^{1/2} > 2ϵ̄^{1/2},

where sep(·, ·) is the minimum distance between the spectra of its two arguments. From the given result, there is a matrix P of dimension (m − p) × p such that the matrix Q̃_{1} defined by

Q̃_{1} = Q_{1} + Q_{2}P (44)

is an invariant subspace for M̃, where

∥P ∥ ≤ \frac{γ̃}{δ̃} ≤ \frac{2ϵ̄^{1/2}}{σ_{p}^{2} − σ_{p+1}^{2} − 2ϵ̄^{1/2}} < 1. (45)

Moreover, the representation of M̃ with respect to Q̃_{1} is

Q̃_{1}^{T} M̃Q̃_{1} = Λ̃ = Σ_{1}^{2} + F_{11} + F_{12}P. (46)

The bound (42) follows from (44), (45), and ∥Q_{2}∥ = 1. It follows immediately from the first equality in (46) that Λ̃ is symmetric, and we have

∥Λ̃ − Σ_{1}^{2}∥ ≤ ∥F_{11}∥ + ∥F_{12}∥∥P ∥ ≤ 2ϵ̄^{1/2},

verifying the inequality (43). This inequality implies that the smallest singular value of Λ̃ is no smaller than σ_{p}^{2} − 2ϵ̄^{1/2} > 0, so Λ̃ is symmetric positive definite.

The cited result states further that the matrix Q̃_{2} = Q_{2} − Q_{1}P^{T} is orthogonal to Q̃_{1} and also defines an invariant subspace for M̃. In fact, we have

[Q̃_{1} | Q̃_{2}]^{T} M̃[Q̃_{1} | Q̃_{2}] = [ \begin{matrix}Λ̃ & 0 \\ 0 & Λ̂ \end{matrix} ] ,

for some (m − p) × (m − p) symmetric matrix Λ̂. Since Λ̃ and M̃ both have rank b, we must have Λ̂ = 0, so we have

M̃ = [Q̃_{1} | Q̃_{2}] [ \begin{matrix}Λ̃ & 0 \\ 0 & 0 \end{matrix} ] [Q̃_{1} | Q̃_{2}]^{T} = Q̃_{1}Λ̃Q̃_{1}^{T} .

Hence, (41) is also satisfied, and the proof is complete.

Combining (40b) with (39), we obtain

σ̃_{1}^{2} ≤ σ_{1}^{2} + ϵ̄^{1/2} < σ_{1}^{2} + .2σ_{p}^{2} < 2σ_{1}^{2}. (47)

Another quantity that enters into our error bounds is the norm of L̃_{J̄J̄}^{−1}. We denote

(48) τ =^{def} max(∥L̃_{J̄J̄}^{−1}∥, 1) = max (σ_{|J̄|}(L̃J̄J̄)^{−1}, 1) ,

where σ_{|J̄|}(L̃J̄J̄) denotes the |J̄|th singular value of L̃J̄J̄. (The lower bound of 1 in (48) simplifies our analysis.) Note from (21a) that

∥M_{J̄J̄}^{−1}∥ = ∥L̃_{J̄J̄}^{−1}∥^{2} ≤ τ^{2}. (49)

Since ∥MJ̄J̄∥ ≤ ∥M∥ ≤ σ_{1}^{2}, we have from (34) and (49) that

κ(MJ̄J̄) ≤ σ_{1}^{2}τ^{2} ≤ mτ^{2}. (50)

Under the assumption |J̄| = p, the nonzero part of L̃—the submatrix L̃_{·J̄}—has full rank p and singular values σ̃_{1}, …, σ̃_{p}. Since L̃_{·J̄} differs from L̃J̄J̄ in the presence of the additional rows L̃_{J J̄}, we have

σ_{|J̄|}(L̃J̄J̄) = σ_{p}(L̃J̄J̄) ≤ σ̃_{p},

and therefore

τσ̃_{p} ≥ 1.

The additional rows L̃_{J J̄} can have nontrivial magnitude relative to L̃J̄J̄, so σ̃_{p} may be significantly larger than τ^{−1}. However, σ̃_{p} cannot be too large, since from (39), (40b), and (34), we have that

σ̃_{p}^{2} ≤ σ_{p}^{2} + ϵ̄^{1/2} ≤ 1.2σ_{p}^{2} ≤ 1.2σ_{1}^{2} ≤ 1.2m.

For the purposes of our analysis, we make the assumption that τσ̃_{p}^{2} is moderate in size. Specifically, we assume that

τσ̃_{p}^{2} ≤ 10. (51)

Because of (39) and (40b), we have σ̃_{p}^{2} ≥ σ_{p}^{2} − ϵ̄^{1/2} ≥ .8σ_{p}^{2}, so (51) implies that

τσ_{p}^{2} ≤ 1.25τσ̃_{p}^{2} ≤ 13.

and, in addition,

τϵ̄^{1/2} ≤ .2τσ_{p}^{2} ≤ .25τσ̃_{p}^{2} ≤ 3. (52)

We can now prove the main result of this section.

Theorem 3.6. Suppose that |J̄| = p < m, that the conditions (40) hold, and that the estimate (51) is satisfied. We then have

∥L^{T} (z̃ − z)∥ ≤ (56\frac{σ_{1}^{3}}{σ_{p}^{3}} ϵ̄^{1/2} + 6σ_{1}σ_{p+1}) τ∥z_{J} ∥.

Proof. From (36), we have

∥L^{T} (z̃ − z)∥ = ∥UΣQ^{T} (z̃ − z)∥ = ∥ΣQ^{T} (z̃ − z)∥,

since U is orthogonal. Now from the partition (35), and using the fact that ∥Q_{2}∥ = 1 (unless of course Σ_{2} and Q_{2} are vacuous), we obtain

∥L^{T} (z̃ − z)∥

≤ ∥Σ_{1}Q_{1}^{T} (z̃ − z)∥ + ∥Σ_{2}∥ ∥z̃ − z∥
≤ ∥Σ_{1}^{−1}∥ ∥Σ_{1}^{2}Q_{1}^{T} (z̃ − z)∥ + ∥Σ_{2}∥ ∥z̃ − z∥
(53) ≤ ∥Σ_{1}^{−1}∥ ∥Λ̃Q̃_{1}^{T} (z̃ − z)∥ + ∥Σ_{1}^{−1}∥ ∥Λ̃Q̃_{1}^{T} − Σ_{1}^{2}Q_{1}^{T} ∥ ∥z̃ − z∥ + ∥Σ_{2}∥ ∥z̃ − z∥.

The first term in this expression is easiest to bound. From (35), we have ∥Σ_{1}^{−1}∥ = σ_{p}^{−1}. Applying the relations (41), (20), (38), (47), (29), (27), and (48), respectively, we obtain

∥Λ̃Q̃_{1}^{T} (z̃ − z)∥ = ∥M̃(z̃ − z)∥
= ∥L̃L̃^{T} (z̃ − z)∥
≤ σ̃_{1}∥L̃^{T} (z̃ − z)∥
≤ 2σ_{1}∥L̃^{T} (z̃ − z)∥
≤ 2σ_{1}∥L̃_{J̄J̄}^{−1}∥ ∥E_{J̄J} ∥ ∥z_{J} ∥
≤ 2σ_{1}τϵ̄^{1/2}∥z_{J} ∥.

We therefore have

∥Σ_{1}^{−1}∥ ∥Λ̃Q̃_{1}^{T} (z̃ − z)∥ ≤ 2\frac{σ_{1}}{σ_{p}} τϵ̄^{1/2}∥z_{J} ∥. (54)

The second and third terms in (53) require a bound on ∥z̃ − z∥. From (30) and the fact that z̃_{J} = 0, we have L̃_{J̄J̄}^{T}(z̃J̄ − zJ̄) = L̃_{J J̄}^{T}z_{J} + L̃_{J̄J̄}^{−1}E_{J̄J} z_{J} , (55) and therefore ∥z̃J̄ − zJ̄∥ ≤ ∥L̃_{J̄J̄}^{−T}∥ (∥L̃_{J J̄}∥ + ∥L̃_{J̄J̄}^{−1}∥ ∥E_{J̄J} ∥) ∥z_{J} ∥. (56) From (47), we have ∥L̃_{J J̄}∥ ≤ σ̃_{1} ≤ 2σ_{1}, while from (27), we have ∥E_{J̄J} ∥ ≤ ϵ̄^{1/2}. Substituting these estimates into (56) and using (52), we obtain ∥z̃J̄ − zJ̄∥ ≤ τ(2σ_{1} + τϵ̄^{1/2})∥z_{J} ∥ ≤ (2σ_{1}τ + 3τ)∥z_{J} ∥. (57) Finally, using z̃_{J} = 0 together with τ ≥ 1, (34), and (57), we obtain ∥z̃ − z∥ ≤ ∥z̃J̄ − zJ̄∥ + ∥z_{J} ∥ ≤ (2σ_{1}τ + 3τ + 1)∥z_{J} ∥ ≤ 6σ_{1}τ∥z_{J} ∥. (58)

Turning specifically to the second term in (53), we have from (34), Lemma 3.5, (47), and (40) that

∥Λ̃Q̃_{1}^{T} − Σ_{1}^{2}Q_{1}^{T} ∥ ≤ ∥Σ_{1}^{2} − Λ̃∥ ∥Q_{1}∥ + ∥Λ̃∥ ∥Q_{1} − Q̃_{1}∥

≤ 2ϵ̄^{1/2} + σ̃_{1}^{2} \frac{2ϵ̄^{1/2}}{σ_{p}^{2} − σ_{p+1}^{2} − 2ϵ̄^{1/2}}

≤ 2ϵ̄^{1/2} [1 + \frac{1.2σ_{1}^{2}}{σ_{p}^{2} − σ_{p+1}^{2} − 2ϵ̄^{1/2}} ]

≤ 2ϵ̄^{1/2} \frac{2.2σ_{1}^{2}}{σ_{p}^{2} − σ_{p+1}^{2} − 2ϵ̄^{1/2}}

≤ 4.5ϵ̄^{1/2} \frac{σ_{1}^{2}}{σ_{p}^{2}} [1 − \frac{σ_{p+1}^{2}}{σ_{p}^{2}} − 2\frac{ϵ̄^{1/2}}{σ_{p}^{2}} ]^{−1}

≤ 4.5ϵ̄^{1/2} \frac{σ_{1}^{2}}{σ_{p}^{2}} [1 − .1 − .4]^{−1}
= 9ϵ̄^{1/2} \frac{σ_{1}^{2}}{σ_{p}^{2}} .

By combining this bound with (58) and ∥Σ_{1}^{−1}∥ = σ_{p}^{−1}, we obtain

∥Σ_{1}^{−1}∥ ∥Λ̃Q̃_{1}^{T} − Σ_{1}^{2}Q_{1}^{T} ∥ ∥z̃ − z∥ ≤ \frac{1}{σ_{p}} (9ϵ̄^{1/2} \frac{σ_{1}^{2}}{σ_{p}^{2}} ) (6σ_{1}τ∥z_{J} ∥)
(59)  = 54\frac{σ_{1}^{3}}{σ_{p}^{3}} τϵ̄^{1/2}∥z_{J} ∥.

For the third term in (53), we have from ∥Σ_{2}∥ = σ_{p+1} that

∥Σ_{2}∥ ∥z̃ − z∥ ≤ 6σ_{1}σ_{p+1}τ∥z_{J} ∥. (60)

The result of the theorem is obtained by substituting (54), (59), and (60) into (53).

Note that if J = ∅ (that is, |J̄| = m), we have z̃ = z, so the conclusion of Theorem 3.6 holds trivially in this case as well is we define σ_{m+1} = 0.

4. The Effect of Finite Precision Computations. In the analysis of the preceding section, we assumed for simplicity that all arithmetic was exact. In this section, we take account of the roundoff errors that are introduced when the approximate solution z̃ is calculated in a finite-precision environment.

Our analysis above focused on the approximate solution z̃ obtained from (25), where the subvector z̃_{J} satisfies the following system:

MJ̄J̄z̃J̄ = L̃J̄J̄L̃_{J̄J̄}^{T}z̃J̄ = rJ̄, (61)

while the subvector z̃_{J} is fixed at zero. In this section, we use ẑ to denote the finite precision analog of z̃. We examine errors in ẑ due to

• roundoff error in Algorithm modchol,

• error arising during the triangular substitutions in (61), and

• evaluation error in the right-hand side r. As we see in Section 5, evaluation error in the right-hand side is a significant feature of the application to interior-point codes. We denote this error by e, so that the right-hand side rJ̄ in the system (61) is replaced by rJ̄ + eJ̄.

Fortunately, our results follow in a straightforward way from existing results for the Cholesky factorization, since a close inspection of Algorithm modchol shows that it simply performs a standard Cholesky factorization on the submatrix MJ̄J̄.

Before stating the main results, we introduce two more assumptions. The first concerns the relative sizes of τ and u, specifically,

τ^{2}γ_{m+1} ≤ \frac{1}{5m^{2}} , (62)

where γ_{m+1} is defined as in Section 1. Since τ > 1 and m ≥ 1, it follows immediately that

mγ_{m+1} ≤ .2. (63)

The second assumption is that finite precision does not affect cutoff decisions in Algorithm modchol. That is, the presence of roundoff error in each submatrix M^{(i−1)} does not affect whether the threshold criterion M_{ii}^{(i−1)} ≤ βϵ passes or fails for each i. This assumption concerns the relative sizes of u and ϵ, and it requires some explanation. We cannot expect to take care of the “borderline cases” in which some candidate pivots fall just to one side or the other of the threshold. Rather, we want the cases in which there is a clear distinction between small and large pivots in exact arithmetic to retain this distinction in finite precision arithmetic, and we want the threshold βϵ to fall comfortably inside the “gap” in both settings. In finite precision, the size of rounding error introduced into M_{ii}^{(i−1)} by earlier steps of Algorithm modchol is comparable to βu. (Each time M_{ii} is updated by the algorithm, a positive number no larger than itself is subtracted from it. Since |M_{ii}| ≤ β, the floating-point error introduced here is bounded by βu.) We want these errors to be smaller than the threshold βϵ, so that pivots that are tiny in exact arithmetic do not exceed the threshold in finite precision. Hence, we can state this assumption roughly as follows:

ϵ ≥ u. (64)

The following lemma accounts for the effects of finite precision on the approximate solution z̃ obtained from Algorithm modchol and (25).

Lemma 4.1. Suppose that Algorithm modchol and the triangular substitutions in (61) are performed in finite-precision arithmetic with perturbed right-hand side rJ̄ + eJ̄ to yield an approximate solution ẑ. Suppose, too, that (62) holds and that roundoff error does not affect the composition of J. We then have

∥z̃ − ẑ∥ ≤ 30m^{5/2}γ_{m+1}τ^{3}∥z∥ + 2τ^{2}∥eJ̄∥, (65)

where z is the exact solution from (23).

Proof. Algorithm modchol operates as a standard Cholesky factorization on the submatrix MJ̄J̄, so we can apply a standard perturbation theorem to bound the error in the subvector ẑJ̄. From Higham [4, Theorem 10.4], we find that ẑJ̄ satisfies

(MJ̄J̄ + E_{J̄J̄}^{u})ẑJ̄ = rJ̄ + eJ̄, (66)

where

∥E_{J̄J̄}^{u}∥ ≤ \frac{2mγ_{m+1}}{1 − mγ_{m+1}} ∥MJ̄J̄∥. (67)

Comparing (66) with (61), we find that

MJ̄J̄(z̃J̄ − ẑJ̄) = E_{J̄J̄}^{u}ẑJ̄ − eJ̄.

Manipulating in the usual way, we obtain

∥z̃J̄ − ẑJ̄∥ ≤ \frac{1}{1 − ∥M_{J̄J̄}^{−1}∥ ∥E_{J̄J̄}^{u}∥}∥M_{J̄J̄}^{−1}∥ (∥E_{J̄J̄}^{u}∥ ∥z̃J̄∥ + ∥eJ̄∥) . (68)

It follows immediately from (67) that

∥M_{J̄J̄}^{−1}∥ ∥E_{J̄J̄}^{u}∥ ≤ κ(MJ̄J̄)\frac{2mγ_{m+1}}{1 − mγ_{m+1}} , (69)

Combining (50), (62), and (63), we obtain

κ(MJ̄J̄)\frac{2mγ_{m+1}}{1 − mγ_{m+1}} ≤ \frac{mτ^{2}(2mγ_{m+1})}{.8} ≤ .5,

so that the denominator in (68) is bounded below by .5. Hence, by substitution into (68), using (34), (49), (69), and (63), we have that

∥z̃J̄ − ẑJ̄∥ ≤ 2∥M_{J̄J̄}^{−1}∥ (∥E_{J̄J̄}^{u}∥ ∥z̃J̄∥ + ∥eJ̄∥)
≤ 2τ^{2} (m\frac{2mγ_{m+1}}{1 − mγ_{m+1}} ∥z̃J̄∥ + ∥eJ̄∥)
(70)  ≤ 5m^{2}γ_{m+1}τ^{2}∥z̃J̄∥ + 2τ^{2}∥eJ̄∥.

Finally, we bound ∥z̃J̄∥ in terms of ∥z∥. From (34) and (57), we have

∥z̃J̄∥ ≤ ∥zJ̄∥ + ∥z̃J̄ − zJ̄∥ ≤ ∥zJ̄∥ + (2σ_{1}τ + 3τ)∥z_{J} ∥ ≤ 6σ_{1}τ∥z∥ ≤ 6m^{1/2}τ∥z∥.

By combining this bound with (70), we obtain the result.

The major results of Sections 3 and 4 can be summarized in the following theorem. Theorem 4.2. Suppose that Algorithm modchol and the triangular substitutions

in (61) are performed in finite-precision arithmetic with perturbed right-hand side rJ̄ + eJ̄ to yield an approximate solution ẑ. Suppose, too, that (62) holds and that roundoff errors do not affect the composition of J. Finally, suppose that either

- J = ∅; or

- |J̄| = p < m, the conditions (40) hold, and the estimate (51) is satisfied. We then have

∥L^{T} (ẑ − z)∥ (71)

≤ \{56\frac{σ_{1}^{3}}{σ_{p}^{3}} ϵ̄^{1/2} + 6σ_{1}σ_{p+1} + 30m^{3}γ_{m+1}τ^{2}\} τ∥z∥ + 2m^{1/2}τ^{2}∥eJ̄∥.

Proof. When J = ∅, the result is immediate from Lemma 4.1 and z̃ = z. For the remaining case, we obtain (71) by combining the results of Theorem 3.6 and Lemma 4.1. We need note only that ∥z_{J} ∥ ≤ ∥z∥ and that, from (34), we have

∥L^{T} (z̃ − ẑ)∥ ≤ ∥L∥ ∥z̃ − ẑ∥ = σ_{1}∥z̃ − ẑ∥ ≤ m^{1/2}∥z̃ − ẑ∥.

5. Application to the Interior-Point Algorithm. In this section, we return to the motivating application: primal-dual interior-point software for linear programming and, in particular, the linear system (16) that is solved at each iteration. We apply the main result—Theorem 4.2—and examine the effect of the parameter ϵ and unit roundoff u on the quality of the computed search direction (Δ̂x, Δ̂π, Δ̂s). Our focus is on the later iterations of the interior-point algorithm, during which μ is small and the ill-conditioning of AD^{2}A^{T} can become acute. Our results show how and why errors arise in (Δ̂x, Δ̂π, Δ̂s) and what effect these errors have on the step length, the convergence of the algorithm, and the accuracy that can be attained by this algorithm. They also suggest an appropriate size for the parameter ϵ.

In this section, we revert to an informal style of analysis, using order notation to hide constants of moderate size. Thus if η and ζ are two positive numbers, we write η = O(ζ) if the ratio η/ζ is not too large. Similarly, we write η = Ω(ζ) if η = O(ζ) and ζ = O(η). Conventionally, order notation is used only when η and ζ are quantities that approach zero in the limit of the algorithm in question. Here, however, we use it in connection with the unit roundoff u, which is small but fixed. This slight abuse of notation results in a much clearer insight into the behavior of Algorithm modchol in the interior-point context.

In the next subsection, we look closely at the affine-scaling step, for which r_{xs} is defined by (9). This step is important because it closely approximates the steps taken by most rapidly converging algorithms during their final iterations. Subsection 5.2 shows that the steps calculated during the final stages of Mehrotra’s predictor corrector algorithm (and therefore by most interior-point codes) have essentially the same properties as affine-scaling steps.

5.1. Affine-Scaling Steps. We start by estimating the sizes of the various constituents of the equations (16)—the residuals r_{b} and r_{c}, the B and N components of x, s, and the diagonal matrix D. In standard infeasible-interior-point algorithms (see, for example, Wright [15, Chapter 6]), we have

∥r_{b}∥ = O(μ), ∥r_{c}∥ = O(μ),

x_{i} = Ω(1) (i ∈ B), x_{i} = Ω(μ) (i ∈ N), (72)

s_{i} = Ω(μ) (i ∈ B), s_{i} = Ω(1) (i ∈ N).

These estimates are also observed to hold in practice on the majority of problems for values of μ greater than u^{1/2}. An immediate consequence of these estimates and the definition (15) is that

D_{ii}^{2} = Ω(μ^{−1}) (i ∈ B), D_{ii}^{2} = Ω(μ) (i ∈ N). (73)

We assume the coefficient matrix A to be well conditioned; that is, σ_{1}(A) and σ_{m}(A) are both Ω(1). We assume further that the submatrix A_{·B} of columns A_{·i}, i ∈ B, is well conditioned. It follows from this assumption together with the estimate (73) that the matrix A_{·B}D_{BB}^{2}A_{·B}^{T} has full rank min(|B|, m). In fact, since A_{·B} is well conditioned, all nonzero singular values of A_{·B}D_{BB}^{2}A_{·B}^{T} are Ω(μ^{−1}) in size. Likewise, it follows from (15) and (73) that A_{·N} D_{N N}^{2} A_{·N}^{T} = O(μ), so we conclude that

(74a)  σ_{i}(AD^{2}A^{T}) = Ω(μ^{−1}),  i = 1, 2, …, min(m, |B|),
(74b)  σ_{i}(AD^{2}A^{T}) = O(μ),  i = min(m, |B|) + 1, …, m.

Since the largest diagonal element of AD^{2}A^{T} is also Ω(μ^{−1}), the scaled coefficient matrix for (16a) is

ρAD^{2}A^{T} , where ρ = Ω(μ). (75)

For consistency with Section 3, the singular values of the matrix in (75) are denoted by σ_{i}^{2}. From this definition together with (74) and (75), we deduce that

(76a)  σ_{i}^{2} = Ω(1),  i = 1, 2, …, min(m, |B|),
(76b)  σ_{i}^{2} = O(μ^{2}),  i = min(m, |B|) + 1, …, m.

Recalling our notation p of Section 3, we have in this case that

p = min(m, |B|).

The exact Cholesky factor L (see Sections 3 and 4) satisfies

LL^{T} = ρAD^{2}A^{T} . (77)

Suppose now that Algorithm modchol is used to compute the solution of (16a), where the right-hand-side component r_{xs} is set to its affine-scaling value XS1. This process result in a computed solution Δ̂π^{aff} for (16a). The remaining step components

Δ̂s^{aff} and Δ̂x^{aff} are obtained by substitution into (16b) and (16c), respectively, again in finite-precision arithmetic. Our main tool for analyzing the errors in the computed step is Theorem 4.2.

Consider the exact affine scaling step (Δx^{aff}, Δπ^{aff}, Δs^{aff}). Standard results for infeasible-interior-point methods (see, for example, [15, Theorem 7.5]), together with the conditions (72), imply that

∥(Δx^{aff}, Δs^{aff})∥ = O(μ). (78)

(This estimate holds only when μ falls below a data-dependent threshold ϵ(A, b, c) defined by Wright [15, Chapter 3].) From (16b) and (72), we have

(AA^{T})Δπ = A^{T} (−r_{c} − Δs) = O(μ),

so it follows from our assumptions about the well conditioning of A that

Δπ^{aff} = O(μ). (79)

We can be more specific about the sizes of the critical components Δx_{i}^{aff}, i ∈ N and Δs_{i}^{aff}, i ∈ B. If we multiply the third block row in (7) by (XS)^{−1} and use the definition (9), we obtain

\frac{Δx_{i}^{aff}}{x_{i}} + \frac{Δs_{i}^{aff}}{s_{i}} = −1, i = 1, 2, …, n.

Therefore, from (72) and (78), we have for i ∈ N that

\frac{Δx_{i}^{aff}}{x_{i}} = −1 + \frac{O(μ)}{Ω(1)} = −1 + O(μ),

and therefore, using (72) again, we have

Δx_{i}^{aff} = −x_{i} + O(μ^{2}), i ∈ N . (80)

In a similar way, we obtain

Δs_{i}^{aff} = −s_{i} + O(μ^{2}), i ∈ B. (81)

From the estimates (78), (80), and (81), we can show that a near-unit step can be taken along the direction (Δx^{aff}, Δπ^{aff}, Δs^{aff}) without violating positivity of the x and s components. Substituting (Δx, Δπ, Δs) = (Δx^{aff}, Δπ^{aff}, Δs^{aff}) in (14), we have

1 − α_{max} = O(μ). (82)

To verify this estimate, suppose that s_{i} + αΔs_{i}^{aff} = 0 for some index i ∈ B. From (81), we have

s_{i}(1 − α) + O(μ^{2}) = 0,

so it follows from (72) that

1 − α = O(μ^{2})/s_{i} = O(μ).

For the corresponding component x_{i}, we have from (72) and (78) that x_{i} = Ω(1) and Δx_{i}^{aff} = O(μ). Hence, for all μ sufficiently small and all α ∈ [0, 1], we have x_{i} + αΔx_{i}^{aff} > 0. Similar logic can be applied to the remaining indices i ∈ N, thereby completing our verification of (82).

Returning to the computed affine-scaling step (Δ̂x^{aff}, Δ̂π^{aff}, Δ̂s^{aff}), we now apply Theorem 4.2 after checking that its assumptions of are satisfied for small enough μ and reasonable values of u and ϵ. For double-precision computations, we have u ≈ 10^{−14}. Hence, since A is well conditioned, we can expect the condition (62) to hold in all nonpathological circumstances. Because of (76), our assumption (40a) on the singular value distribution clearly holds for all sufficiently small μ. The condition (40b) is satisfied for any reasonable choice of ϵ. The assumption that Algorithm modchol correctly identifies the numerical rank (that is, |J̄| = p) is, as we discussed in Section 3, difficult to guarantee, but it was observed to hold on all problems that we tested. The assumption that rounding errors do not interfere with the makeup of the small pivot index set J is likewise impossible to verify rigorously; but, as discussed in Section 4, it can reasonably be expected to hold when ϵ ≥ u (64).

A good choice for ϵ—one that satisfies the assumptions just mentioned while keeping the bound (71) as small as possible—is therefore

ϵ = u. (83)

For generality, we continue to use ϵ and ϵ̄ in the analysis that follows, substituting the specific value (83) only at the end.

Having verified that we can reasonably expect Theorem 4.2 to hold for the system (16a), we now estimate the quantities on the right-hand side of (71). From (76a), we have σ_{1}/σ_{p} = O(1), while from (76b), we have σ_{p+1} = O(μ). The general estimate (34) yields σ_{1} = O(1), while the definition of γ_{m+1} gives the estimate γ_{m+1} = O(u).

We need to account, too, for the errors incurred in evaluating the right-hand side of (16a). The floating-point error in forming r_{xs} = XS1 is only O(μu) in magnitude, since just a single floating-point multiplication is needed to calculate each component x_{i}s_{i} of this vector, and each such element is O(μ) (see (72)). The residuals r_{b} and r_{c} have magnitude O(μ) in exact arithmetic (see (72)), but they are calculated as differences of O(1) quantities and so contain evaluation error of absolute magnitude O(u). Specifically, componentwise errors in the computed version of r_{c} are bounded by (|A|^{T} |π| + |s| + |c|) u, and similarly for r_{b}. Because of the estimate (73), the errors in r_{c} are magnified to (μ^{−1}u) when we multiply by AD^{2} in (16a). In fact, this term is the dominant one in the total right-hand-side evaluation error. The errors that occur when we perform floating-point addition of the terms r_{b}, AD^{2}r_{c}, and AS^{−1}r_{xs} are less significant; they lead to additional terms of sizes O(u) and O(μ^{−1}u^{2}). In summary, the total right-hand-side evaluation error is O(μ^{−1}u). Hence, after scaling by the factor ρ defined in (75), we have

∥e∥ = O(u), (84)

where e is the error vector of Section 4.

Substituting the estimates (76), (79), and (84) into (71), we have

∥L^{T} (Δ̂π^{aff} − Δπ^{aff})∥ ≤ {O(ϵ̄^{1/2}) + O(μ) + τ^{2}O(u)} τO(μ) + τ^{2}O(u).

If

τ = O(1) (85)

(a reasonable estimate when the Cholesky factorization correctly identifies the numerical rank and A_{·B} is well conditioned), the error bound above simplifies to

∥L^{T} (Δ̂π^{aff} − Δπ^{aff})∥ ≤ O(ϵ̄^{1/2}μ + μ^{2} + u). (86)

From (77) we have that

ρ^{1/2}DA^{T} = QL^{T} ,

for some orthogonal matrix Q. Since orthogonal transformations do not affect the Euclidean norm of a vector, we can substitute ρ^{1/2}DA^{T} for L^{T} in (86) and use (75) to write

∥DA^{T} (Δ̂π^{aff} − Δπ^{aff})∥ = ρ^{−1/2}∥L^{T} (Δ̂π^{aff} − Δπ^{aff})∥
(87)  ≤ O(ϵ̄^{1/2}μ^{1/2} + μ^{3/2} + μ^{−1/2}u).

Note too that from (58), (65), (79), and (84), we have

(88) ∥Δ̂π^{aff}∥ ≤ ∥Δ̂π^{aff} − Δ̃π^{aff}∥ + ∥Δ̃π^{aff} − Δπ^{aff}∥ + ∥Δπ^{aff}∥ = O(μ + u),

where Δ̃π^{aff} is the approximate solution that would be obtained by Algorithm modchol if it was used to solve (16a) in exact arithmetic.

Next, we examine the effect of the error in Δ̂π^{aff} and the evaluation error in the right-hand side of (16b) on the calculated step Δ̂s^{aff}. From (79) and (88), we have that

(89) ∥Δπ^{aff} − Δ̂π^{aff}∥ ≤ ∥Δπ^{aff}∥ + ∥Δ̂π^{aff}∥ = O(μ + u).

Hence, taking into account the O(u) evaluation error in the term r_{c}, we have immediately from (16b) that

(90) Δs^{aff} − Δ̂s^{aff} = O(u) − A^{T} (Δπ^{aff} − Δ̂π^{aff}) = O(μ + u).

Clearly, for the “large” components of s—namely, the i ∈ N components—errors of this magnitude do not affect the step length α_{max} to the boundary defined in (14). However, for the critical components i ∈ B, the estimate (90) is not good enough to guarantee that α_{max} is close to 1. (Repeating the argument that follows (82), we find only that 1 − α_{max} = O(1).) Fortunately, a refined estimate of the error in the B components is available. As in (90), we have

(91) Δs^{aff} − Δ̂s^{aff} = −A^{T} (Δπ^{aff} − Δ̂π^{aff}) + O(u) = D^{−1}v + O(u),

where from (87) we have

v = DA^{T} (Δ̂π^{aff} − Δπ^{aff}) = O(ϵ̄^{1/2}μ^{1/2} + μ^{3/2} + μ^{−1/2}u). (92)

From (73), we have D_{ii} = Ω(μ^{−1/2}) for i ∈ B, so from (91) we obtain

Δ̂s_{i}^{aff} − Δs_{i}^{aff} = O(ϵ̄^{1/2}μ + μ^{2} + u), i ∈ B. (93)

As in the discussion following (82), we find that s_{i} + αΔ̂s_{i}^{aff} = 0 is possible only if

1 − α = O(ϵ̄^{1/2} + μ + μ^{−1}u). (94)

This estimate suggests that near-unit steps can be taken, at least in the Δ̂s^{aff} components, provided that μ is significantly larger that u. When μ = O(u), all bets are off!

Finally, we estimate the errors in the computed version of Δx^{aff} (obtained from (16c)) and estimate their effect on the α_{max}. Again, we consider the components i ∈ B and i ∈ N separately.

For i ∈ B, the O(μu) evaluation error in (r_{xs})_{i} is magnified by the term s_{i}^{−1} =

Ω(μ^{−1}). From (93), replacement of Δs^{aff} by Δ̂s^{aff} yields an additional error of size O(ϵ̄^{1/2}μ) + O(μ^{2}) + O(u), which is also magnified by the Ω(μ^{−1}) factor. The other arithmetic errors are less significant. In summary, we find that

Δ̂x_{i}^{aff} − Δx_{i}^{aff} = O(ϵ̄^{1/2} + μ + μ^{−1}u), i ∈ B. (95)

By the usual reasoning, we find that x_{i} + αΔ̂x_{i}^{aff} = 0 is possible for i ∈ B only for α satisfying (94).

For i ∈ N, the O(μu) evaluation error in (r_{xs})_{i} is not magnified appreciably by s_{i}^{−1}, while from (90), the O(μ + u) error in Δs^{aff} is actually diminished after multiplication by s_{i}^{−1}x_{i} = O(μ). We find that

Δ̂x_{i}^{aff} − Δx_{i}^{aff} = O(μu + μ^{2}), i ∈ N . (96)

Hence, we can have x_{i} + αΔ̂x_{i}^{aff} = 0 for i ∈ N only if

|1 − α| = O(u + μ). (97)

From (94) and (97), we conclude that the value of α_{max} defined by (14), with the

calculated direction (Δ̂x^{aff}, Δ̂π^{aff}, Δ̂s^{aff}) replacing the exact search direction, satisfies the estimate

1 − α_{max} = O(ϵ̄^{1/2} + μ + μ^{−1}u). (98)

Note from (89), (90), and (96) that, in an absolute sense, the errors in Δ̂π^{aff},

Δ̂s^{aff}, and Δ̂x_{i}^{aff}, i ∈ N are small. By contrast, the O(μ^{−1}u) term in (95) implies

that the errors in Δ̂x_{i}^{aff}, i ∈ B, may become large as μ ↓ 0. These large errors may in turn cause the residuals r_{b} to grow as μ ↓ 0. These expectations are confirmed by the computational experiments of Section 6.

The estimate (98) and the parameter choice ϵ = u (83) suggest strongly that the algorithm should be terminated when

μ ≤ u^{1/2}. (99)

When μ reaches this threshold, all three terms in the estimate (98) are in balance.

Below this threshold, the O(μ^{−1}u) term in Δ̂x_{i}^{aff} may cause r_{b} to grow, making further reduction of μ counterproductive. The convergence tolerances used by most interior-point codes—arrived at by practical experience rather than any theoretical considerations—are similar to (99). The code PCx is typical. It declares optimality if the following three conditions are satisfied:

\frac{∥r_{b}∥}{1 + ∥b∥} ≤ tol, \frac{∥r_{c}∥}{1 + ∥c∥} ≤ tol, \frac{|c^{T} x − b^{T} π|}{1 + |c^{T} x|} ≤ tol,

where the default value of tol is 10^{−8}. (Note that 10^{−8} ≈ u^{1/2} in double precision arithmetic on most machines.)

5.2. Mehrotra Predictor-Corrector Steps. Having analyzed the affine-scaling search direction and its calculated approximation, we turn our attention briefly to the search direction used by Mehrotra’s predictor-corrector algorithm. As mentioned in Section 2, these steps are obtained by setting r_{xs} as in (12), for some heuristic choice of the centering parameter ζ. We can write the search direction as

(100) (Δx, Δπ, Δs) = (Δx^{aff}, Δπ^{aff}, Δs^{aff}) + (Δx^{cc}, Δπ^{cc}, Δs^{cc}),

where (Δx^{cc}, Δπ^{cc}, Δs^{cc}) is the “corrector-centering” step component that satisfies the following linear system:

[ \begin{matrix}0 & A^{T} & I \\ A & 0 & 0 \\ S & 0 & X \end{matrix} ] [ \begin{matrix}Δx^{cc} \\ Δπ^{cc} \\ Δs^{cc} \end{matrix} ] = [ \begin{matrix}0 \\ 0 \\ ζμ1 − ΔX^{aff}ΔS^{aff}1 \end{matrix} ] .

Block elimination on this system yields the following special case of (16a):

AD^{2}A^{T} Δπ^{cc} = AD^{2} [X^{−1}(ζμ1 − ΔX^{aff}ΔS^{aff}1)] .

Since we assume full rank of A, and since the diagonal elements of D are all strictly positive, the coefficient matrix is invertible, and we have

∥Δπ^{cc}∥ ≤ ∥(AD^{2}A^{T})^{−1}AD^{2}∥ ∥X^{−1}∥ ∥ζμ1 − ΔX^{aff}ΔS^{aff}1∥.

A result of Stewart [9] and Todd [11] states that the norm ∥(AD^{2}A^{T})^{−1}AD^{2}∥ is bounded independently of D over the set of all positive definite diagonal matrices D (and therefore independently of x and s with (x, s) > 0). Therefore, we have

∥Δπ^{cc}∥ = O(∥X^{−1}∥)O(∥ζμ1 − ΔX^{aff}ΔS^{aff}1∥).

From (72), we have ∥X^{−1}∥ = O(μ^{−1}), while from (78), it follows that ∥ΔX^{aff}ΔS^{aff}1∥ = O(μ^{2}). Hence, we have

∥Δπ^{cc}∥ = O(ζ + μ). (101)

A typical heuristic for choosing the centering parameter ζ is to set

ζ = (μ_{aff}/μ)^{3},

where μ_{aff} is the value of μ that results from a full step-to-boundary α_{max} along the affine-scaling direction. If the search direction is exact, we have μ_{aff} = O(μ^{2}), so

this heuristic yields ζ = O(μ^{3}). Use of the calculated direction (Δ̂x^{aff}, Δ̂π^{aff}, Δ̂s^{aff}) together with the estimate (98) leads us to expect μ_{aff} = O(μ^{2}) in this case too, provided that μ ≥ u^{1/2}. Hence, we have from (101) that ∥Δπ^{cc}∥ = O(μ), and so, from (100) and (79), we have

∥Δπ∥ = O(μ), (102)

where Δπ is the π component of the Mehrotra search direction.

We also can apply the Stewart-Todd result to formulae for Δx^{cc} and Δs^{cc} to show that ∥(Δx^{cc}, Δs^{cc})∥ = O(μ). Therefore, we have

∥(Δx, Δs)∥ = O(μ), (103)

corresponding to (78).

Because of the estimates (102) and (103), the analysis of the preceding subsection can be applied without modification to the calculated version of the search direction (100). In particular, if we redefine the step-to-boundary α_{max} in terms of this calculated step (Δ̂x, Δ̂π, Δ̂s), we find that the estimate (98) still applies. We conclude that near-unit steps can still be taken along this direction provided that μ ≥ u^{1/2}.

6. Implementation and Computational Results. Most interior-point codes use modified Cholesky algorithms with essentially the same properties as Algorithm modchol. They differ slightly, however, in the implementation. The IPMOS code of Xu, Hung, and Ye [16] replaces small pivot elements by 1 and fills out the corresponding column of the Cholesky factor with zeros and also inserts a zero in the right-hand side. The criterion for identifying a small pivot is not explained in the reference [16], but otherwise this strategy is equivalent to Algorithm modchol. Zhang’s LIPSOL code [17] and the PCx code of Czyzyk, Mehrotra, and Wright [1] replace small pivots by a huge number—10^{128}—but otherwise leave the Cholesky algorithm unchanged. The net effect is, however, almost equivalent to Algorithm modchol and the triangular substitution procedure (25). The advantage of this approach is that it involves minimal changes to a standard sparse Cholesky code. We need only add a loop to calculate the largest diagonal element β, and a small pivot check immediately before the point at which the computation L_{ii} = \sqrt{M_{ii}} is performed.

To test that the analysis of this paper was reflected in practical computations, we coded a primal-dual algorithm that used Algorithm modchol in conjunction with the formulation (16). The code was used to solve some small random linear programs in which the amount of degeneracy—the composition of index sets B and N —was carefully controlled. At each iterate, we monitored various quantities and compared them against the estimates of Section 5.

The linear programming test problems were posed in standard form (2) with m = 6 and n = 12. The matrix A is fully dense, with elements (τ_{1} − .5)10^{6(τ_{2}−.5)}, where τ_{1} and τ_{2} are random variables drawn from a uniform distribution on the interval [0, 1]. (Of course, the values of τ_{1} and τ_{2} are different for each element of the matrix.) We can reasonably expect this matrix A to satisfy the well-conditioning assumptions of Section 5. The user specifies the number of indices to appear in B, and we set

|N | = n − |B|,  N = {1, 2, ⋯, |N |},  B = {|N | + 1, ⋯, n}.

A primal solution x^{∗} is constructed with

x_{i}^{∗} = 0 (i = 1, 2, ⋯, |N |), x_{i}^{∗} = 10^{3(τ−.5)} (i = |N | + 1, ⋯, n),

where τ is randomly drawn from the uniform distribution on [0, 1]. We choose the dual solution π^{∗} to be the vector (1, 1, ⋯, 1)^{T}, and fix an optimal dual slack vector s^{∗} to be

s_{i}^{∗} = 10^{4(τ−.5)} (i = 1, 2, ⋯, |N |), s_{i}^{∗} = 0 (i = |N | + 1, ⋯, n),

where τ is random as above. Finally, we set b = Ax^{∗} and c = A^{T} π^{∗} + s^{∗}.

The code was an implementation of the infeasible-interior-point algorithm described by Wright [13]. The details of this algorithm are unimportant; we need note only that its iterates satisfy the estimates (72) in exact arithmetic and that the algorithm takes steps along the affine scaling direction during its later iterations. At each iteration of the algorithm, we calculated the affine scaling direction (whether or not it

was actually used as a search direction) and printed the norms ∥Δ̂x^{aff}∥_{∞}, ∥Δ̂π^{aff}∥_{∞},

and ∥Δ̂s^{aff}∥_{∞} alongside the duality measure μ and residual norm ∥(r_{b}, r_{c})∥_{∞} for the current point. We also kept track of the number of small pivots encountered during the factorization, that is, the number of elements in J. The parameter ϵ was set to 10^{−12}, which is about 100u on the SPARCstation 5 that was used for the experiments. The results were not particularly sensitive to this parameter.

Results are shown in Tables 1–4. For each iteration of the algorithm, these tables list the number of small pivots |J |, the base-10 logarithms of μ, ∥(r_{b}, r_{c})∥_{∞}, and the affine-scaling step norms mentioned above. The step-to-boundary α_{max} along the calculated affine-scaling direction is also tabulated. A horizontal line in each table indicates the iterate at which termination occurs according to the criterion (99).

In Table 1 we chose |B| = m = 6, making the linear program nondegenerate

and the primal-dual solution unique. It is clear that Δ̂π^{aff} and Δ̂s^{aff} satisfy the estimates (88) and (90), respectively, even when the algorithm is continues past the

point of normal termination. The component Δ̂x^{aff}, on the other hand, clearly shows the influence of the O(μ^{−1}u) error term in (95) when μ becomes comparable to or

smaller than u. Note, too, that the error in Δ̂x^{aff} is transmitted to the residual r_{b} on succeeding iterations but that this effect does not become destructive until μ is much smaller than its normal termination threshold. The values of α_{max} are also consistent with the estimate (98). This step length approaches 1 until the normal

point of termination is reached, after which the errors in Δ̂x^{aff} and r_{b} make further progress impossible.

Table 2 shows the interesting case in which we choose |B| = 4, so that the coefficient matrix in (16a) has four singular values of magnitude Ω(μ^{−1}) and two of magnitude Ω(μ). The second column shows that Algorithm modchol correctly identifies the numerical rank during the last few iterations and that the interior-point algorithm continues to generate useful steps and to make good progress even after modchol encounters small pivots. Apart from this feature, the behavior is the same as in Table 1, with errors in Δ̂x^{aff} causing the interior-point algorithm to behave poorly when it is permitted to run past its normal point of termination. We noted that for all iterations, the “small” pivots were at the bottom right corner of the Cholesky matrix, so that (28) rather than the general estimate (27) applies to the perturbation matrix E. In this case, we can replace ϵ̄^{1/2} by ϵ̄ in estimates of Section 5 such as (93), (95), and (98).

Table 3 illustrates another case in which |B| = 4, with the added complication that A is rank deficient. (We forced rank deficiency by setting A_{1j} = 0 and A_{2j} = 0 for j = 1, 2, ⋯, n − 1, so that the first and second rows each contain a single nonzero in their last column.) The (2, 2) pivot is skipped at every invocation of Algorithm modchol. As μ becomes small, the final pivot is skipped as well, and the numerical rank is correctly determined. Since the small pivots are not localized in the bottom right corner, the special bound (28) does not apply, so we cannot strengthen the bounds on the step components as in the previous paragraph. The computational behavior is qualitatively the same as in Tables 1 and 2.

Table 4 illustrates a problem for which |B| = 8. Here, the coefficient matrices retain full numerical rank at all iterates, and the behavior is similar to that reported in Table 1. One point of difference is that the errors in Δ̂x^{aff}, which start to increase after iteration 19, do not have an immediate effect on the residual r_{b}. The reason is simply that this particular interior-point algorithm chose to take a path-following step at iterations 21 and 22 rather than the affine scaling step, and the Δx components were calculated accurately in the path following step. An affine-scaling step is, however, taken at iteration 28, and the effect of the error in Δ̂x^{aff} on the residual r_{b} at the following iterate is obvious.

REFERENCES

[1] J. Czyzyk, S. Mehrotra, and S. J. Wright, PCx User Guide, Technical Report OTC 96/01, Optimization Technology Center, Argonne National Laboratory and Northwestern University, May 1996.

[2] A. Forsgren, P. Gill, and J. Shinnerl, Stability of symmetric ill-conditioned systems arising in interior methods for constrained optimization, SIAM Journal on Matrix Analysis and Applications, 17 (1996), pp. 187–211.

[3] A. George and J. W.-H. Liu, Computer Solution of Large Sparse Positive Definite Systems, Prentice-Hall, 1981.

[4] N. J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM Publications, Philadelphia, 1996.

[5] C. L. Lawson and R. J. Hanson, Solving Least Squares Problems, Prentice-Hall, Englewood Cliffs, NJ, 1974.

Table 1

Affine scaling step characteristics for a problem with m = 6, n = 12, |B| = 6. ∥ · ∥ = ∥ · ∥_{∞}, and the horizontal line represents the normal point of termination.

Small  log  log  log  log
Iteration Pivots log μ ∥(r_{b}, r_{c})∥ ∥Δ̂x^{aff}∥ ∥Δ̂π^{aff}∥ ∥Δ̂s^{aff}∥  α_{max}

⋮

12  0  -0.6  -11.1  -0.1  -0.6  0.6 .26426
13  0  -1.4  -10.7  0.4  -1.1  0.1 .77520
14  0  -2.1  -10.7  1.2  -2.3  -1.1 .39373
15  0  -3.3  -10.4  -0.3  -1.3  -0.1 .62276
16  0  -4.8  -8.1  -1.1  -5.2  -3.9 .99697
17  0  -7.2  -10.5  -3.5  -8.3  -7.1 .99999
18  0  -12.0  -12.2  -8.2  -14.0  -12.5 >.99999
19  0  -21.0  -12.0  -3.6  -14.9  -13.9 .99975
20  0  -24.2  -4.6  -1.4  -15.0  -13.9 .93989
21  0  -26.2  -1.5  1.4  -15.3  -14.5 .06843

⋮

Table 2

Affine scaling step characteristics for a problem with m = 6, n = 12, |B| = 4. ∥ · ∥ = ∥ · ∥_{∞}, and the horizontal line represents the normal point of termination.

Small  log  log  log  log
Iteration Pivots log μ ∥(r_{b}, r_{c})∥ ∥Δ̂x^{aff}∥ ∥Δ̂π^{aff}∥ ∥Δ̂s^{aff}∥  α_{max}

⋮

12  0  -0.6  -12.0  0.1  -1.3  0.7 .95133
13  0  -1.9  -11.4  -1.5  -0.2  1.8 .51719
14  0  -2.4  -9.5  -1.8  -0.9  1.0 .90453
15  1  -3.4  -9.3  -2.7  -5.5  -3.5 .98770
16  2  -5.2  -9.1  -4.4  -7.2  -5.2 .99977
17  2  -8.5  -11.1  -7.7  -10.5  -8.5 >.99999
18  2  -14.4  -13.2  -12.5  -15.9  -13.8 >.99999
19  2  -25.1  -12.3  -2.1  -15.9  -13.7 >.99999
20  2  -30.4  -1.8  5.0  -15.9  -14.5 .00016
21  3  -30.4  2.3  10.6  -16.1  -13.5 <.00001

⋮

Table 3

Affine scaling step characteristics for a problem with m = 6, n = 12, |B| = 4, in which A is rank deficient. ∥ · ∥ = ∥ · ∥_{∞}, and the horizontal line represents the normal point of termination.

Small  log  log  log  log
Iteration Pivots log μ ∥(r_{b}, r_{c})∥ ∥Δ̂x^{aff}∥ ∥Δ̂π^{aff}∥ ∥Δ̂s^{aff}∥  α_{max}

⋮

12  1  -1.4  -11.3  0.1  0.5  1.0 .78614
13  1  -2.1  -10.5  -1.9  1.7  2.0 .17726
14  1  -2.7  -9.2  -0.8  1.2  1.3 .41306
15  1  -2.9  -9.1  0.1  0.5  0.6 .00442
16  1  -3.2  -8.9  -0.6  0.2  0.4 .78585
17  1  -3.9  -8.8  -1.1  -2.0  -1.7 .93466
18  1  -4.8  -9.7  -2.0  -2.0  -1.7  99179
19  2  -6.2  -10.9  -3.4  -6.0  -5.5  99970
20  2  -8.6  -10.1  -5.8  -8.3  -7.9 >.99999
21  2  -12.7  -10.9  -10.2  -12.7  -12.0 >.99999
22  2  -20.2  -11.5  -4.1  -12.3  -12.5 .99988
23  2  -21.7  -4.6  -2.6  -12.9  -12.2 >.99999
24  2  -27.3  -3.0  2.9  -12.5  -12.3 .00711

⋮

Table 4

Affine scaling step characteristics for a problem with m = 6, n = 12, |B| = 8. ∥ · ∥ = ∥ · ∥_{∞}, and the horizontal line represents the normal point of termination.

Small  log  log  log  log
Iteration Pivots log μ ∥(r_{b}, r_{c})∥ ∥Δ̂x^{aff}∥ ∥Δ̂π^{aff}∥ ∥Δ̂s^{aff}∥ α_{max}

⋮

12  0  -0.1  -5.7  3.5  -3.6  -2.3 .42851
13  0  -0.9  -9.5  2.1  -2.4  -1.4 .60234
14  0  -1.3  -9.9  1.8  -2.8  -1.7 .38898
15  0  -2.0  -10.6  1.8  -1.8  -0.6 .30608
16  0  -2.2  -10.6  1.2  -1.8  -0.6 .37400
17  0  -3.3  -10.5  -1.0  -3.9  -2.7 .68815
18  0  -4.2  -11.3  -0.3  -4.1  -2.9 .99691
19  0  -6.7  -9.7  -3.0  -7.5  -6.3 .99998
20  0  -11.3  -10.6  0.6  -12.4  -11.2 .98674
21  0  -13.2  -10.4  2.3  -14.4  -13.1 .19743
22  0  -13.3  -10.5  2.3  -14.6  -13.5 .04721

⋮

28  0  -18.8  −∞  -2.2  -18.5  -14.4 .99999
29  0  -23.7  -5.5  1.1  -15.0  -13.7 .59156

⋮

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

[7] C. Mészáros, The inexact minimum local fill-in ordering algorithm, Working Paper 95–7, Computer and Automation Institute, Hungarian Academy of Sciences, Budapest, 1995.

[8] E. Ng and B. W. Peyton, Block sparse Cholesky algorithms on advanced uniprocessor computers, SIAM Journal on Scientific Computing, 14 (1993), pp. 1034–1056.

[9] G. W. Stewart, On scaled projections and pseudoinverses, Linear Algebra and Its Applications, 112 (1989), pp. 189–193.

[10] G. W. Stewart and J. Sun, Matrix Perturbation Theory, Computer Science and Scientific Computing, Academic Press, New York, 1990.

[11] M. J. Todd, A Dantzig-Wolfe-like variant of Karmarkar’s interior-point linear programming algorithm, Operations Research, 38 (1990), pp. 1006–1018.

[12] M. H. Wright, Some properties of the Hessian of the logarithmic barrier function, Mathematical Programming, 67 (1994), pp. 265–295.

[13] 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. To appear in Annals of Operations Research.

[14] , Stability of augmented system factorizations in interior-point methods, Preprint MCS– P446–0694, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill., June 1994. Revised, July 1995. To appear in SIAM Journal of Matrix Analysis and Applications.

[15] , Primal-Dual Interior-Point Methods, SIAM, Philadelphia, Pa, 1996. To appear.

[16] X. Xu, P. Hung, and Y. Ye, A simplified homogeneous and self-dual linear programming algorithm and its implementation. To appear in Annals of Operations Research, September 1993.

[17] Y. Zhang, Solving large-scale linear programs by interior-point methods under the MATLAB enviroment, Technical Report TR96-01, Department of Mathematics and Statistics, University of Maryland Baltimore County, Baltimore, Md, 1996.