Modified Cholesky Factorizations in Interior-Point Algorithms for Linear Programming.
We investigate a modified Cholesky algorithm typical of those used in most 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 algorithms (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 alternative factorizations; and they perform robustly on most practical problems, yielding good interior-point steps even when the coefficient matrix of the main linear system to be solved for the step components is ill conditioned. We investigate this surprisingly robust performance by using analytical tools from matrix perturbation theory and error analysis, illustrating our results with computational experiments. Finally, we point out the potential limitations of this approach.
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
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'esz'aros [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
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 Newton-
logarithmic 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 infor-
mation, 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 interior-point
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
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 oe i (A). We use oe 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 . The ith column of M is denoted by M \Deltai ,
and the column submatrix consisting of columns j 2 J is denoted by M \DeltaJ .
Unit roundoff error is denoted by u. Higham [4, Chapter 1] defines u implicitly
by the statement that when ff and i are any two floating-point numbers, op denotes
+, \Gamma, \Theta, and =, and fl(\Delta) denotes the floating-point representation of a real number,
we have
For any positive integer m with mu ! 1, we define
(1)
(see Higham [4, Lemma 2.1]).
The notation k \Delta k denotes the Euclidean vector norm k \Delta k 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 jAj. We use 1 to denote the vector
Finally, we mention the parameter ffl that defines the pivot threshold in the modified
algorithm. A second quantity -ffl, which is related to ffl by
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:
subject to
(2)
. The dual of (2) is
subject to A T -
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)
(4b)
(4c)
(4d)
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 f1; ng can be partitioned into two sets B and N such that for all
primal-dual solutions
x
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,
In general, the search direction (\Deltax; \Delta-; \Deltas) is obtained from the following linear
\Delta-
\Deltas5 =4 \Gammar c
where the coefficient matrix is the Jacobian of (6) and the right-hand side components
r b and r c are defined by
In a pure Newton (affine-scaling) method, the remaining right-hand side component
r xs is defined by
and, in this case, we denote the solution of (7) by (\Deltax aff ; \Delta- aff ; \Deltas aff ). In a path-following
method, we have
where - is the duality gap defined by
is a centering parameter. In the "Mehrotra predictor-corrector" al-
gorithm, which is used as the basis of many practical codes, the search direction is
calculated by setting
where \DeltaX aff and \DeltaS aff are the diagonal matrices formed from the affine-scaling step
components \Deltax aff and \Deltas 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 i in (12).
Once a search direction has been determined, the primal-dual algorithm takes a
step of the form
where ff is chosen to maintain strict positivity of the x and s components; that is,
In most codes, ff is chosen to be some fraction of the step-to-boundary ff max defined
as
A typical strategy is to set
as the interior-point method approaches the solution
set.
By applying block elimination to (7) and using the notation
we obtain the following equivalent system:
(16a)
(16c)
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 \Delta- is obtained by triangular
substitution with the computed factor. The remaining direction components
are recovered from (16b) and (16c). This technique yields steps (\Deltax; \Delta-; \Deltas) that
are useful search directions for the interior-point algorithm, even when the matrix
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 \Delta- from (16a),
we are not interested so much in the error in \Delta- itself as in the effect of this error
on the remaining step components \Deltas and \Deltax that are recovered from (16b) and
(16c), respectively. If the relative errors in these components are large, the positivity
requirement may cause the step length ff 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 ill-conditioned
matrices for which small or negative pivots may arise during the factorization

Algorithm modchol accepts an m \Theta m symmetric positive definite matrix M as
input, together with a small positive user-defined parameter ffl, 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 ae
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 ffl with
if M (i\Gamma1)
(* skip this elimination step *)
im
.
else
(* perform the usual Cholesky elimination step *)
~
~
L ki .
The ith column of ~
L is zero for each i 2 J ; that is, ~
If we denote
and denote the complement of J in f1;
J , it follows from (17) that
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
which we denote for convenience by ~
M . That is, we have
~
By partitioning this equation into its J and -
J components and using ~
(19), we obtain
~
(21a)
~
The exact Cholesky factor L (whose existence is guaranteed by the assumed positive
definiteness of M ) satisfies
Given the linear system
where M is the matrix factored by modchol, the exact solution obviously satisfies
The approximate solution ~
z is chosen so that the partial vector ~ z -
J solves the reduced
system M -
z -
J , 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 substitu-
tions; that is,
~
z -
~
z
Note that z = ~
z when on the other hand, the difference between
~
z and z can be large in a relative sense. We have
z -
z -
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
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 ~
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 columns of M (i)
is symmetric positive definite, for all 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 positive definite. Suppose that the desired property
holds for M (i\Gamma1) . If i 2 J , then the lower right (m \Gamma i) \Theta (m \Gamma i) submatrix of
M (i) is identical to the lower right (m \Gamma i) \Theta (m \Gamma i) submatrix of M (i\Gamma1) , which
is positive definite by assumption. Otherwise, if
is obtained by
applying one step of Cholesky reduction to M (i\Gamma1) . It is known that the remaining
submatrix resulting from this operation is positive definite; hence, the lower right
in question is positive definite, and the desired property
holds.
The second claim follows immediately from the fact that M ii -
and the fact that the diagonal elements cannot increase during Algorithm
modchol.
Lemma 3.2. For each i 2 J , we have
Therefore,
Proof. From Lemma 3.1, we have (M (i\Gamma1)
l;l for each
ffl. Since the diagonals of each submatrix
M (i\Gamma1) are bounded by 1, we have M (i\Gamma1)
i;l
l;l
Hence, we have
thereby proving the first claim. By (18), we have
thereby proving (27).
In the case in which all the small pivots appear in the bottom right corner of the
matrix (that is, index p), the estimate (27) can be
improved to
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 ~
Theorem 3.3. For the exact solution z and approximate solution ~
z defined in
(24) and (25), respectively, we have that
Proof. From (24) together with (21), we have
JJ z J
~
~
JJ
z J
while from (25), we have
~
z -
J ~ z J
J ~ z:
By combining these two relations, we obtain
~
Since ~
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 ffl. Accordingly,
we introduce a little more notation. The eigenvalues of M are denoted by oe 2
We define the diagonal matrix \Sigma by
It follows that there exists an orthogonal matrix Q such that
Because the largest diagonal in M is 1, we have by elementary analysis that
In the subsequent analysis, we assume that there is an integer p with 1 -
such that
ffl ffl is small relative to oe 2
there is a significant gap in the spectrum of M between oe 2
p and
p+1 .
(We will be more specific about these two assumptions presently.) By partitioning
the spectrum at the gap, we obtain
From (33), Q can be partitioned accordingly to obtain
are the singular values of L. In
fact, we must have
orthogonal matrix U , where \Sigma and Q are defined as above.
We use ~
to denote the eigenvalues of the perturbed matrix ~
M .
It follows immediately from (20) that the singular values of ~
are ~ oe i ,
The rank of ~
J j, because ~
J is lower triangular with nonzero diagonals while
~
Therefore, we have
~
As in (36), there are orthogonal m \Theta m matrices ~
U and ~
Q such that
~
U ~
where
~
It is an immediate consequence of an eigenvalue perturbation result of Stewart and
Sun [10, Corollary IV.4.13] and Lemma 3.2 that
The main assumption of this section is that j -
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 ffl. 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 ffl.
Lemma 3.4. If -
Proof. If
we have from (37) and (39) that
contradicting our assumption that -ffl 1=2 ! oe 2
.
However, the conditions on ffl, oe p , and oe p+1 needed to prove the other half of
the result-j -
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 -
and that the values oe p and oe p+1 from
(31) and ffl from Lemma 3.2 satisfy the conditions
(40a)
(40b)
Then there is a p \Theta p symmetric positive definite matrix ~
and an orthonormal m \Theta p
matrix ~
~
~
~
(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 ~
use (33) and partition as in (35) to obtain
We now make the following identifications with the quantities in the cited result:
~
where sep(\Delta; \Delta) is the minimum distance between the spectra of its two arguments.
From the given result, there is a matrix P of dimension (m \Gamma p) \Theta p such that the
matrix ~
defined by
~
is an invariant subspace for ~
~
Moreover, the representation of ~
M with respect to ~
~
The bound (42) follows from (44), (45), and kQ 2 It follows immediately from
the first equality in (46) that ~
is symmetric, and we have
verifying the inequality (43). This inequality implies that the smallest singular value
of ~
is no smaller than oe 2
is symmetric positive definite.
The cited result states further that the matrix ~
is orthogonal to
~
defines an invariant subspace for ~
M . In fact, we have

for some (m \Gamma p) \Theta (m \Gamma p) symmetric matrix -
. Since ~
and ~
both have rank b,
we must have -
0, so we have
~
~
~
Hence, (41) is also satisfied, and the proof is complete.
Combining (40b) with (39), we obtain
Another quantity that enters into our error bounds is the norm of ~
J . We denote
J jth singular value of ~
J . (The lower bound of 1 in
simplifies our analysis.) Note from (21a) that
1 , we have from (34) and (49) that
Under the assumption j -
the nonzero part of ~
L-the submatrix ~
J -has
full rank p and singular values ~ oe
oe p . Since ~
J differs from ~
J in the presence
of the additional rows ~
J , we have
and therefore
The additional rows ~
J can have nontrivial magnitude relative to ~
J , so ~ oe p may
be significantly larger than - \Gamma1 . However, ~ oe p cannot be too large, since from (39),
(40b), and (34), we have that
~
For the purposes of our analysis, we make the assumption that - ~
p is moderate in
size. Specifically, we assume that
Because of (39) and (40b), we have ~
implies that
and, in addition,
3:
We can now prove the main result of this section.
Theorem 3.6. Suppose that j -
m, that the conditions (40) hold, and
that the estimate (51) is satisfied. We then have
Proof. From (36), we have
since U is orthogonal. Now from the partition (35), and using the fact that kQ 2
(unless of course \Sigma 2 and Q 2 are vacuous), we obtain
~
~
The first term in this expression is easiest to bound. From (35), we have k\Sigma \Gamma1
. Applying the relations (41), (20), (38), (47), (29), (27), and (48), respectively,
we obtain
~
We therefore have
~
The second and third terms in (53) require a bound on k~z \Gamma zk. From (30) and the
fact that ~ z
~
and therefore
kz J k:
From (47), we have k ~
while from (27), we have
Substituting these estimates into (56) and using (52), we obtain
Finally, using ~
z together with - 1, (34), and (57), we obtain
Turning specifically to the second term in (53), we have from (34), Lemma 3.5,
(47), and (40) that
~
oe 22-ffl 1=2
By combining this bound with (58) and k\Sigma \Gamma1
, we obtain
~
For the third term in (53), we have from k\Sigma 2
The result of the theorem is obtained by substituting (54), (59), and (60) into
(53).
Note that if
m), we have ~
so the conclusion of
Theorem 3.6 holds trivially in this case as well is we define oe
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:
z -
while the subvector ~ z J is fixed at zero. In this section, we use - z to denote the finite
precision analog of ~
z. We examine errors in -
z due to
ffl roundoff error in Algorithm modchol,
ffl error arising during the triangular substitutions in (61), and
ffl 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 r -
J in the system (61) is replaced by r -
J .
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 M -
J .
Before stating the main results, we introduce two more assumptions. The first
concerns the relative sizes of - and u, specifically,
where fl m+1 is defined as in Section 1. Since - ? 1 and m - 1, it follows immediately
that
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\Gamma1)
does not affect whether the threshold criterion M (i\Gamma1)
ii - fi ffl passes or fails for each
i. This assumption concerns the relative sizes of u and ffl, and it requires some ex-
planation. 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 fi ffl to fall comfortably inside the "gap" in both settings. In finite precision,
the size of rounding error introduced into M (i\Gamma1)
ii by earlier steps of Algorithm mod-
chol is comparable to fiu. (Each time M ii is updated by the algorithm, a positive
number no larger than itself is subtracted from it. Since jM ii j - fi, the floating-point
error introduced here is bounded by fiu.) We want these errors to be smaller than
the threshold fi ffl, 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:
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
J to yield an approximate solution - z. Suppose, too, that (62) holds and that
roundoff error does not affect the composition of J . We then have
where z is the exact solution from (23).
Proof. Algorithm modchol operates as a standard Cholesky factorization on the
J , so we can apply a standard perturbation theorem to bound the error
in the subvector - z -
J . From Higham [4, Theorem 10.4], we find that -
z -
J satisfies
where
Comparing (66) with (61), we find that
z -
Manipulating in the usual way, we obtain
z -
It follows immediately from (67) that
Combining (50), (62), and (63), we obtain
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
Finally, we bound k~z -
J k in terms of kzk. From (34) and (57), we have
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
J to yield an approximate solution - z. Suppose, too, that (62) holds and that
roundoff errors do not affect the composition of J . Finally, suppose that either
m, the conditions (40) hold, and the estimate (51) is satisfied.
We then have
ae
oe
Proof. When the result is immediate from Lemma 4.1 and ~
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 kz J k - kzk and that, from (34), we have
zk:
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 ffl and
unit roundoff u on the quality of the computed search direction ( c
\Deltax; c
\Delta-; c
\Deltas). 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 ( c
\Deltax; c
\Delta-; c
\Deltas) 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 ffl.
In this section, we revert to an informal style of analysis, using order notation to
hide constants of moderate size. Thus if j and i are two positive numbers, we write
if the ratio j=i is not too large. Similarly, we
O(j). Conventionally, order notation is used only when j and i 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
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
We assume the coefficient matrix A to be well conditioned; that is, oe 1 (A) and
are both
\Omega\Gammah/1 We assume further that the submatrix A \DeltaB of columns A \Deltai ,
well conditioned. It follows from this assumption together with the estimate
(73) that the matrix A \DeltaB D 2
\DeltaB has full rank min(jBj; m). In fact, since A \DeltaB is well
conditioned, all nonzero singular values of A \DeltaB D 2
\DeltaB
it follows from (15) and (73) that A \DeltaN D 2
so we conclude that
(74a)
Since the largest diagonal element of AD 2 A T is
scaled coefficient
matrix for (16a) is
For consistency with Section 3, the singular values of the matrix in (75) are denoted
by oe 2
. From this definition together with (74) and (75), we deduce that
(76a)
Recalling our notation p of Section 3, we have in this case that
The exact Cholesky factor L (see Sections 3 and
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 c
\Delta- aff
for (16a). The remaining step components
c
\Deltas aff
and c
\Deltax 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 (\Deltax aff ; \Delta- aff ; \Deltas aff ). Standard results for
methods (see, for example, [15, Theorem 7.5]), together with
the conditions (72), imply that
(This estimate holds only when - falls below a data-dependent threshold ffl(A; b; c)
defined by Wright [15, Chapter 3].) From (16b) and (72), we have
so it follows from our assumptions about the well conditioning of A that
\Delta-
We can be more specific about the sizes of the critical components \Deltax aff
and \Deltas aff
we multiply the third block row in (7) by (XS) \Gamma1 and use the
definition (9), we obtain
\Deltax aff
\Deltas aff
Therefore, from (72) and (78), we have for i 2 N that
\Deltax aff
and therefore, using (72) again, we have
\Deltax aff
In a similar way, we obtain
\Deltas aff
From the estimates (78), (80), and (81), we can show that a near-unit step can
be taken along the direction (\Deltax aff ; \Delta- aff ; \Deltas aff ) without violating positivity of the
x and s components. Substituting (\Deltax; \Delta-;
have
To verify this estimate, suppose that s i
(81), we have
so it follows from (72) that
For the corresponding component x i , we have from (72) and (78) that x i
and \Deltax aff
O(-). Hence, for all - sufficiently small and all ff 2 [0; 1], we have
logic can be applied to the remaining indices i 2 N , thereby
completing our verification of (82).
Returning to the computed affine-scaling step ( c
\Deltax aff
\Delta- aff
\Deltas aff
), we now apply
Theorem 4.2 after checking that its assumptions of are satisfied for small enough
- and reasonable values of u and ffl. For double-precision computations, we have
. 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 ffl. The assumption that Algorithm
modchol correctly identifies the numerical rank (that is,
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 ffl - u (64).
A good choice for ffl-one that satisfies the assumptions just mentioned while
keeping the bound (71) as small as possible-is therefore
For generality, we continue to use ffl and -
ffl 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 oe 1 =oe O(1), while from (76b), we have oe O(-). The general estimate
while the definition of fl m+1 gives the estimate fl
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
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
u, and similarly for r b . Because of the estimate (73), the errors
in r c are magnified to (- \Gamma1 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 \Gamma1 r xs
are less significant; they lead to additional terms of sizes O(u) and O(- \Gamma1 u 2 ). In
summary, the total right-hand-side evaluation error is O(- \Gamma1 u). Hence, after scaling
by the factor ae defined in (75), we have
where e is the error vector of Section 4.
Substituting the estimates (76), (79), and (84) into (71), we have
\Delta- aff
If
(a reasonable estimate when the Cholesky factorization correctly identifies the numerical
rank and A \DeltaB is well conditioned), the error bound above simplifies to
\Delta- aff
From (77) we have that
ae
for some orthogonal matrix Q. Since orthogonal transformations do not affect the
Euclidean norm of a vector, we can substitute ae 1=2 DA T for L T in (86) and use (75)
to write
\Delta-
aff
\Delta-
aff
Note too that from (58), (65), (79), and (84), we have
\Delta- aff
\Delta- aff
\Delta- aff
\Delta- aff
where ~
\Delta- aff
is the approximate solution that would be obtained by Algorithm mod-
chol if it was used to solve (16a) in exact arithmetic.
Next, we examine the effect of the error in c
\Delta- aff
and the evaluation error in the
right-hand side of (16b) on the calculated step c
\Deltas aff
. From (79) and (88), we have
that
\Delta-
aff
\Delta-
aff
Hence, taking into account the O(u) evaluation error in the term r c , we have immediately
from (16b) that
\Deltas
\Deltas aff
\Delta- aff
Clearly, for the "large" components of s-namely, the i 2 N components-errors
of this magnitude do not affect the step length ff max to the boundary defined in (14).
However, for the critical components i 2 B, the estimate (90) is not good enough to
guarantee that ff max is close to 1. (Repeating the argument that follows (82), we find
only that Fortunately, a refined estimate of the error in the B
components is available. As in (90), we have
\Deltas
\Deltas aff
\Delta- aff
where from (87) we have
\Delta-
aff
From (73), we have D B, so from (91) we obtain
c
\Deltas aff
As in the discussion following (82), we find that s i
\Deltas aff
possible only if
This estimate suggests that near-unit steps can be taken, at least in the c
\Deltas aff
com-
ponents, provided that - is significantly larger that u. When all bets are
off!
Finally, we estimate the errors in the computed version of \Deltax aff (obtained from
(16c)) and estimate their effect on the ff max . Again, we consider the components
separately.
For B, the O(-u) evaluation error in (r xs ) i is magnified by the term s \Gamma1
replacement of \Deltas aff by c
\Deltas aff
yields an additional error of size
O(-ffl which is also magnified by
arithmetic errors are less significant. In summary, we find that
c
\Deltax aff
By the usual reasoning, we find that x i
\Deltax aff
satisfying (94).
For the O(-u) evaluation error in (r xs ) i is not magnified appreciably
by s
, while from (90), the O(- + u) error in \Deltas aff is actually diminished after
multiplication by s
We find that
c
\Deltax aff
Hence, we can have x
\Deltax aff
From (94) and (97), we conclude that the value of ff max defined by (14), with the
calculated direction ( c
\Deltax aff
\Delta- aff
\Deltas aff
replacing the exact search direction, satisfies
the estimate
Note from (89), (90), and (96) that, in an absolute sense, the errors in c
\Delta- aff
c
\Deltas aff
, and c
\Deltax aff
are small. By contrast, the O(- \Gamma1 u) term in (95) implies
that the errors in c
\Deltax
aff
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 suggest strongly that the
algorithm should be terminated when
-
When - reaches this threshold, all three terms in the estimate (98) are in balance.
Below this threshold, the O(- \Gamma1 u) term in c
\Deltax 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:
where the default value of tol is 10 \Gamma8 . (Note that 10 \Gamma8 - 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 i. We can write the search direction as
where (\Deltax cc ; \Delta- cc ; \Deltas cc ) is the "corrector-centering" step component that satisfies
the following linear system:4 0 A T I
\Delta- cc
\Deltas cc5 =4
Block elimination on this system yields the following special case of (16a):
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
A result of Stewart [9] and Todd [11] states that the norm k(AD 2 A T )
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
From (72), we have kX while from (78), it follows that k\DeltaX aff \DeltaS aff
O(- 2 ). Hence, we have
A typical heuristic for choosing the centering parameter i is to set
where - aff is the value of - that results from a full step-to-boundary ff max along the
affine-scaling direction. If the search direction is exact, we have -
this heuristic yields Use of the calculated direction ( c
\Deltax aff
\Delta- aff
\Deltas aff
together with the estimate (98) leads us to expect - case too,
provided that - u 1=2 . Hence, we have from (101) that k\Delta- cc
from (100) and (79), we have
where \Delta- is the - component of the Mehrotra search direction.
We also can apply the Stewart-Todd result to formulae for \Deltax cc and \Deltas cc to show
that k(\Deltax cc ; \Deltas cc Therefore, we have
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 ff max in terms of this calculated
\Deltax; c
\Delta-; c
\Deltas), 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 fi, and a small pivot check immediately before
the point at which the computation L
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
12. The matrix A is fully dense, with elements
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
A primal solution x   is constructed with
x
where - is randomly drawn from the uniform distribution on [0; 1]. We choose the
dual solution -   to be the vector (1; fix an optimal dual slack vector
s   to be
s
where - is random as above. Finally, we set
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 k c
\Deltax aff
k1 , k c
\Delta- aff
k1 ,
and k c
\Deltas aff
k1 alongside the duality measure - and residual norm k(r b ; r c )k 1 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 ffl was set to
\Gamma12 , 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 jJ j, the base-10 logarithms of -, k(r
the affine-scaling step norms mentioned above. The step-to-boundary ff 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 making the linear program nondegenerate
and the primal-dual solution unique. It is clear that c
\Delta- aff
and c
\Deltas aff
satisfy the
estimates (88) and (90), respectively, even when the algorithm is continues past the
point of normal termination. The component c
\Deltax aff
, on the other hand, clearly shows
the influence of the O(- \Gamma1 u) error term in (95) when - becomes comparable to or
smaller than u. Note, too, that the error in c
\Deltax 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 ff 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 c
\Deltax aff
and r b make further
progress impossible.

Table

2 shows the interesting case in which we choose 4, so that the co-efficient
matrix in (16a) has four singular values of
and two of
-). 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
encounters small pivots. Apart from this feature, the behavior is the same
as in

Table

1, with errors in c
\Deltax 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
so that (28) rather than the general estimate (27) applies to the
perturbation matrix E. In this case, we can replace -ffl 1=2 by -ffl in estimates of Section
5 such as (93), (95), and (98).

Table

3 illustrates another case in which 4, with the added complication
that A is rank deficient. (We forced rank deficiency by setting A
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 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 c
\Deltax 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 \Deltax 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 c
\Deltax aff
on the residual r b at the
following iterate is obvious.



--R

Technical Report OTC 96/01
Stability of symmetric ill-conditioned systems arising in interior methods for constrained optimization
Computer Solution of Large Sparse Positive Definite Systems
Accuracy and Stability of Numerical Algorithms
Solving Least Squares Problems
Computational experience with a primal-dual interior point method for linear programming

Block sparse Cholesky algorithms on advanced uniprocessor com- puters
On scaled projections and pseudoinverses
Matrix Perturbation Theory
A Dantzig-Wolfe-like variant of Karmarkar's interior-point linear programming algorithm
Some properties of the Hessian of the logarithmic barrier function



A simplified homogeneous and self-dual linear programming algorithm and its implementation
Solving large-scale linear programs by interior-point methods under the MATLAB enviroment
--TR

--CTR
Francis R. Bach , Michael I. Jordan, Kernel independent component analysis, The Journal of Machine Learning Research, 3, p.1-48, 3/1/2003
error analysis;interior-point algorithms and software;cholesky factorization;matrix perturbations