A Modified Cholesky Algorithm Based on a Symmetric Indefinite Factorization.
Given a symmetric and not necessarily positive definite matrix A, a modified Cholesky algorithm computes a Cholesky factorization permutation matrix and E is a perturbation chosen to make A+E positive definite. The aims include producing a small-normed E and making A+E reasonably well conditioned. Modified Cholesky factorizations are widely used in optimization. We propose a new modified Cholesky algorithm based on a symmetric indefinite factorization computed using a new pivoting strategy of Ashcraft, Grimes, and Lewis. We analyze the effectiveness of the algorithm, both in theory and practice, showing that the algorithm is competitive with the existing algorithms of Gill, Murray, and Wright and Schnabel and Eskow. Attractive features of the new algorithm include easy-to-interpret inequalities that explain the extent to which it satisfies its design goals, and the fact that it can be implemented in terms of existing software.
Introduction
. Modified Cholesky factorization is a widely used technique in
optimization; it is used for dealing with indefinite Hessians in Newton methods [11],
[21] and for computing positive definite preconditioners [6], [20]. Given a symmetric
matrix A, a modified Cholesky algorithm produces a symmetric perturbation E such
that A + E is positive definite, along with a Cholesky (or LDL T ) factorization of
A+E. The objectives of a modified Cholesky algorithm can be stated as follows [21].
O1. If A is "su#ciently positive definite" then E should be zero.
O2. If A is indefinite, #E# should not be much larger than
positive definite }
for some appropriate norm.
O3. The matrix should be reasonably well conditioned.
O4. The cost of the algorithm should be the same as the cost of standard Cholesky
factorization to highest order terms.
Two existing modified Cholesky algorithms are one by Gill, Murray, and Wright
[11, section 4.4.2.2], which is a refinement of an earlier algorithm of Gill and Murray
[10], and an algorithm by Schnabel and Eskow [21].
The purpose of this work is to propose an alternative modified Cholesky algorithm
that has some advantages over the existing algorithms. In outline, our approach is to
compute a symmetric indefinite factorization
# Received by the editors April 26, 1996; accepted for publication (in revised form) by P. Gill June
4, 1997; published electronically July 17, 1998. The research of the second author was supported by
Engineering and Physical Sciences Research Council grant GR/H/94528.
http://www.siam.org/journals/simax/19-4/30289.html
Department of Mathematics, University of Manchester, Manchester, M13 9PL, England
(chengsh@ma.man.ac.uk, na.nhigham@na-net.ornl.gov).
where P is a permutation matrix, L is unit lower triangular, and D is block diagonal
with diagonal blocks of dimension 1 or 2, and to provide the factorization
where F is chosen so that D + F (and hence also A + E) is positive definite. 1 This
approach is not new; it was suggested by Mor-e and Sorensen [19] for use with factorizations
computed with the Bunch-Kaufman [3] and Bunch-Parlett [4] pivoting
strategies. However, for neither of these pivoting strategies are all the conditions
satisfied, as is recognized in [19]. The Bunch-Parlett pivoting strategy
requires O(n 3 ) comparisons for an n-n matrix, so condition (O4) does not hold. For
the Bunch-Kaufman strategy, which requires only O(n 2 ) comparisons, it is di#cult
to satisfy conditions (O1)-(O3), as we explain in section 3.
We use a new pivoting strategy for the symmetric indefinite factorization devised
by Ashcraft, Grimes, and Lewis [2], for which conditions (O1)-(O3) are satisfied
to within factors depending only on n and for which the cost of the pivot searches
is usually negligible. We describe this so-called bounded Bunch-Kaufman (BBK)
pivoting strategy and its properties in the next section.
There are two reasons why our algorithm might be preferred to those of Gill, Mur-
ray, and Wright and of Schnabel and Eskow (henceforth denoted the GMW algorithm
and the SE algorithm, respectively). The first is a pragmatic one: we can make use of
any available implementation of the symmetric indefinite factorization with the BBK
pivoting strategy, needing to add just a small amount of post-processing code to form
the modified Cholesky factorization. In particular, we can use the e#cient implementations
for both dense and sparse matrices written by Ashcraft, Grimes, and Lewis [2],
which make extensive use of levels 2 and 3 BLAS for e#ciency on high-performance
machines. In contrast, in coding the GMW and SE algorithms one must either begin
from scratch or make nontrivial changes to an existing Cholesky factorization code.
The second attraction of our approach is that we have a priori bounds that explain
the extent to which conditions (O1)-(O3) are satisfied-essentially, if L is well
conditioned then an excellent modified Cholesky factorization is guaranteed. For the
GMW and SE algorithms it is di#cult to describe under what circumstances the
algorithms can be guaranteed to perform well.
2. Pivoting strategies. We are interested in symmetric indefinite factorizations
computed in the following way. If the symmetric matrix A # R n-n is nonzero,
we can find a permutation # and an integer so that
with Having chosen such a # we can factorize
I n-s
# .
Strictly, (1.2) is not a Cholesky factorization, since we allow D + F to have 2 - 2 diagonal
blocks, but since any such blocks are positive definite it seems reasonable to use the term "modified
MODIFIED CHOLESKY ALGORITHM 1099
This process is repeated recursively on the (n - s) - (n - s) Schur complement
yielding the factorization (1.1) on completion. This factorization costs n 3 /3 operations
(the same cost as Cholesky factorization of a positive definite matrix) plus the cost
of determining the permutations #.
The Bunch-Parlett pivoting strategy [4] searches the whole submatrix S at each
stage, requiring a total of O(n 3 ) comparisons, and it yields a matrix L whose maximum
element is bounded by 2.781. The Bunch-Kaufman pivoting strategy [3], which is used
with the symmetric indefinite factorization in both LAPACK [1] and LINPACK [7],
searches at most two columns of S at each stage, so it requires only O(n 2 ) comparisons
in total. The Bunch-Kaufman pivoting strategy yields a backward stable factorization
[16], but #L# is unbounded, even relative to #A# , which makes this pivoting
strategy unsuitable for use in a modified Cholesky algorithm, for reasons explained in
section 3.
To describe the BBK pivoting strategy [2] it su#ces to describe the pivot choice
for the first stage of the factorization.
Algorithm BBK (BBK pivoting strategy). This algorithm determines the pivot
for the first stage of the symmetric indefinite factorization applied to a symmetric
matrix A # R n-n .
maximum magnitude of any subdiagonal entry in column 1.
there is nothing to do on this stage of the factorization.
if |a 11 | # 1
use a 11 as a 1 - 1 pivot
else
repeat
r := row index of first (subdiagonal) entry of maximum magnitude
in column i.
# r := maximum magnitude of any o#-diagonal entry in column r.
if |a rr | # r
use a rr as a 1 - 1 pivot rows and columns
1 and r).
else if #
use # a ii a ri
a ri a rr
# as a 2 - 2 pivot
columns 1 and i, and 2 and r).
else
until a pivot is chosen
The repeat loop in Algorithm BBK searches for an o#-diagonal element a ri that
is simultaneously the largest in magnitude in the rth row and the ith column, and
it uses this element to build a 2 - 2 pivot; the search terminates prematurely if a
suitable
The following properties noted in [2] are readily verified, using the property that
any 2 - 2 pivot satisfies
# a ii a ri
a ri a rr
# .
1. Every entry of L is bounded by max{1/(1 - #), 1/# 2.78.
2. Every 2 - 2 pivot block D ii satisfies # 2 (D ii
3. The growth factor for the factorization, defined in the same way as for Gaussian
elimination, is bounded in the same way as for the Bunch-Kaufman
pivoting strategy, namely, by (1
Since the value of # i increases strictly from one pivot step to the next, the search
in Algorithm BBK takes at most n steps. The cost of the searching is intermediate
between the cost for the Bunch-Kaufman strategy and that for the Bunch-Parlett
strategy. Matrices are known for which the entire remaining submatrix must be
searched at each step, in which case the cost is the same as for the Bunch-Parlett
strategy. However, Ashcraft, Grimes, and Lewis [2] found in their numerical experiments
that on average less than 2.5k comparisons were required to find a pivot from a
k -k submatrix, and they give a probabilistic analysis which shows that the expected
number of comparisons is less than ek # 2.718k for matrices with independently
distributed random elements. Therefore we regard the symmetric indefinite factorization
with the BBK pivoting strategy as being of similar cost to Cholesky factorization,
while recognizing that in certain rare cases the searching overhead may increase the
operation count by about 50%.
The symmetric indefinite factorization with the BBK pivoting strategy is backward
stable; the same rounding error analysis as for the Bunch-Kaufman pivoting
strategy is applicable [2], [16].
The modified Cholesky algorithm of the next section and the corresponding analysis
are not tied exclusively to the BBK pivoting strategy. We could use instead the
"fast Bunch-Parlett" pivoting strategy from [2], which appears to be more e#cient
than the BBK strategy when both are implemented in block form [2]. We mention
in passing that a block implementation of the SE algorithm has been developed by
Dayd-e [5]. Alternatively, we could use one of the pivoting strategies from [8], [9].
3. The modified Cholesky algorithm. We begin by defining the distance
from a symmetric matrix A # R n-n to the symmetric matrices with minimum eigenvalue
# min at least #, where # 0:
The distances in the 2- and Frobenius norms, and perturbations that achieve them,
are easily evaluated (cf. [12, Thms. 2.1, 3.1]).
Theorem 3.1. Let the symmetric matrix A # R n-n have the spectral decomposition
for the Frobenius norm,
and there is a unique optimal perturbation in (3.1), given by
MODIFIED CHOLESKY ALGORITHM 1101
For the 2-norm,
and an optimal perturbation is #)I. The Frobenius norm perturbation
(3.2) is also optimal in the 2-norm.
Our modified Cholesky algorithm has a parameter # 0 and it attempts to
produce the perturbation (3.2).
Algorithm MC (modified Cholesky factorization). Given a symmetric matrix
A # R n-n and a parameter # 0 this algorithm computes a permutation matrix P ,
a unit lower triangular matrix L, and a block diagonal matrix D with diagonal blocks
of dimension 1 or 2 such that
and A+E is symmetric positive definite (or symmetric positive semidefinite if
The algorithm attempts to ensure that if # min (A) < # then # min E) #.
1. Compute the symmetric indefinite factorization PAP
using the
pivoting strategy.
2. Let
D, where #
D is the minimum Frobenius norm perturbation
that achieves # min (
D) # (thus #
D ii is the
minimum Frobenius norm perturbation that achieves # min (
D ii +#
To what extent does Algorithm MC achieve the objectives (O1)-(O4) listed in
section 1? Objective (O4) is clearly satisfied, provided that the pivoting strategy
does not require a large amount of searching, since the cost of step 2 is negligible. For
objectives (O1)-(O3) to be satisfied we need the eigenvalues of A to be reasonably
well approximated by those of
D. For the Bunch-Kaufman pivoting strategy the
elements of L are unbounded and the eigenvalues of
D can di#er greatly from those
of A (subject to A and
D having the same inertia), as is easily shown by example.
This is the essential reason why the Bunch-Kaufman pivoting strategy is unsuitable
for use in a modified Cholesky algorithm.
To investigate objectives (O1)-(O3) we will make use of a theorem of Ostrowski
[18, p. 224]. Here, the eigenvalues of a symmetric n-n matrix are ordered # n #
Theorem 3.2 (Ostrowski). Let M # R n-n be symmetric and S # R n-n
nonsingular. Then for each
Assuming first that # min (A) > 0 and applying the theorem with
D and
D).
Now E will be zero if # min (
D) #, which is certainly true if
Next, we assume that # min (A) is negative and apply Theorems 3.1 and 3.2 to
obtain
Using Theorem 3.2 again, with (3.4), yields
D)
A final invocation of Theorem 3.2 gives
E) # min (LL T )# min (
D) # min (LL T )#
and
D)
# .
Hence
E)
# .
We can now assess how well objectives (O1)-(O3) are satisfied. To satisfy objective
(O1) we would like E to be zero when # min (A) #, and to satisfy (O2) we would
like #E# 2 to be not much larger than # min (A) when A is not positive definite.
The su#cient condition (3.3) for E to be zero and inequality (3.5) show that these
conditions do hold modulo factors # max,min (LL T ). Inequality (3.6) bounds # 2 (A+E)
with the expected reciprocal dependence on #, again with terms # max,min (LL T ). The
conclusion is that the modified Cholesky algorithm is guaranteed to perform well if
are not too far from 1.
Note that, since L is unit lower triangular, e T
which implies that
1. For the BBK pivoting strategy we have
Furthermore,
using a bound from [15, Thm. 8.13 and Prob. 8.5]. These upper bounds are approximately
attainable, but in practice are rarely approached. In particular, the upper
bound of (3.8) can be approached only in the unlikely event that most of the subdiagonal
elements of L are negative and of near maximal magnitude. Note that each
causes a subdiagonal element l i+1,i to be zero and so further reduces the
likelihood of #L being large.
In the analysis above we have exploited the fact that the extent to which the
eigenvalues of A and
D agree can be bounded in terms of the condition of L. If L is
well conditioned then the singular values of A are close to the moduli of the eigenvalues
of
D. We are currently exploring the application of this fact to the computation of
rank-revealing factorizations.
MODIFIED CHOLESKY ALGORITHM 1103
4. Comparison with the GMW and SE algorithms. The GMW and SE
algorithms both carry out the steps of a Cholesky factorization of a symmetric matrix
A # R n-n , increasing the diagonal entries as necessary in order to ensure that
negative pivots are avoided. (Actually, the GMW algorithm works with an LDL T
factorization, where D is diagonal, but the di#erence is irrelevant to our discussion.)
Hence both algorithms produce Cholesky factors of P T E)P with a diagonal E.
From Theorem 3.1 we note that the "optimal" perturbation in objective (O2) of section
1 is, in general, full for the Frobenius norm and can be taken to be diagonal for
the 2-norm (but is generally not unique). There seems to be no particular advantage
to making a diagonal perturbation to A. Our algorithm perturbs the whole matrix,
in general.
By construction, the GMW and SE algorithms make perturbations E to A that
are bounded a priori by functions of n and #A# only. The GMW algorithm produces
a perturbation E for which
where # 0 is a tolerance,
and u is the unit roundo# [11, p. 110]. For the SE algorithm the perturbation is
bounded in terms of a certain eigenvalue bound # obtained by applying Gershgorin's
theorem:
(#),
where # is a tolerance, suggested in [21] to be chosen as . The quantity #
satisfies # n(#), so (4.2) is a smaller bound than (4.1) by about a factor n.
The bounds (4.1) and (4.2) can be compared with (3.5) for Algorithm MC. The
bound (3.5) has the advantage of directly comparing the perturbation made by Algorithm
MC with the optimal one, as defined by (3.1) and evaluated in Theorem 3.1, and
it is potentially a much smaller bound than (4.1) and (4.2) if |# min (A)| # max (A)|
and # 2 (LL T ) is not too large. On the other hand, the bound (3.5) can be much larger
than (4.1) and (4.2) if # 2 (LL T ) is large.
All three algorithms satisfy objective (O1) of not modifying a su#ciently positive
definite matrix, though for the GMW and SE algorithms no condition analogous
to (3.3) that quantifies "su#ciently" in terms of # min (A) is available. Bounds for
that are exponential in n hold for the GMW and SE algorithms [21]. The
same is true for Algorithm MC: see (3.6)-(3.8).
To summarize, in terms of the objectives of section 1 for a modified Cholesky algo-
rithm, Algorithm MC is theoretically competitive with the GMW and SE algorithms,
with the weakness that if # 2 (LL T ) is large then the bound on #E# 2 is weak.
When applied to an indefinite matrix, the GMW and SE algorithms provide information
that enables a direction of negative curvature of the matrix to be produced;
these directions are required in certain algorithms for unconstrained optimization in
order to move away from nonminimizing stationary points. For an indefinite matrix,
Algorithm MC provides immediate access to a direction of negative curvature from the
1104 SHEUNG HUN CHENG AND NICHOLAS J. HIGHAM
computed in step 1, and because #(L) is bounded, this direction
satisfies conditions required for convergence theory [19].
Finally, we consider the behavior of the algorithms in the presence of rounding
errors. Algorithm MC is backward stable because the underlying factorization is [2]:
barring large element growth in the symmetric indefinite factorization with the BBK
pivoting strategy, the algorithm produces LDL T factors not of P but of
a constant. Although no
comments on numerical stability are given in [11] and [21], a simple argument shows
that the GMW and SE algorithms are backward stable. Apply either algorithm to
A, obtaining the Cholesky factorization P apply the same
algorithm to P (A+E)P T : it will not need to modify P (A+E)P T , so it will return the
same computed R factor. But since no modification was required, the algorithm must
have carried out a standard Cholesky factorization. Since Cholesky factorization is
a backward stable process, the modified Cholesky algorithm must itself be backward
stable.
5. Numerical experiments. We have experimented with Matlab implementations
of Algorithm MC and the GMW and SE algorithms. The M-file for the GMW
algorithm was provided by M. Wright and sets the tolerance (which is the
value of Matlab's variable eps). The M-file for the SE algorithm was provided by
E. Eskow and sets the tolerance In Algorithm MC we set # u#A# .
The aims of the experiments are as follows: to see how well the Frobenius norm of
the perturbation E produced by Algorithm MC approximates the distance -F (A, #)
defined in (3.1), and to compare the norms of the perturbations E and the condition
numbers of A produced by the three algorithms. We measure the perturbations
E by the ratios
-F (A, #)
|# min (A)|
which di#er only in their normalization and the choice of norm. Algorithm MC
attempts to make r F close to 1. The quantity r 2 is used by Schnabel and Eskow to
compare the performance of the GMW and SE algorithms; since E is diagonal for
these algorithms, r 2 compares the amount added to the diagonal with the minimum
diagonal perturbation that makes the perturbed matrix positive semidefinite.
First, we note that the experiments of Schnabel and Eskow [21] show that the SE
algorithm can produce a substantially smaller value of r 2 than the GMW algorithm.
Schnabel and Eskow also identified a 4 - 4 matrix for which the GMW algorithm
significantly outperforms the SE algorithm:
1538.3 284.9 -2706.6
We give results for this matrix in Table 5.1; they show that Algorithm MC can also
significantly outperform the SE algorithm.
We ran a set of tests similar to those of Schnabel and Eskow [21]. The matrices
A are of the form with the eigenvalues # i from one
MODIFIED CHOLESKY ALGORITHM 1105

Table
Measures of E for 4 - 4 matrix (5.1).
MC GMW SE
matrix
a) n=25, eig. range[-1,10000]
matrix
c) n=25, eig. range[-1,1]
matrix
eig. range[-10000,-1]
matrix
eig. range[-1,10000]
matrix
d) n=25, eig. range[-1,1]
f) n=25, eig. range[-10000,-1]
Fig. 5.1. Measures of E for
-, MC -.
of three random uniform distributions: [-1, 10 4 ], [-1, 1], and [-10 4 , -1]. For the first
range, one eigenvalue is generated from the range [-1, 0) to ensure that A has at least
one negative eigenvalue. The matrix Q is a random orthogonal matrix from the Haar
distribution, generated using the routine qmult from the Test Matrix Toolbox [14],
which implements an algorithm of Stewart [22]. For each eigenvalue distribution we
generated di#erent matrices, each corresponding to a fresh sample of # and of
Q. We took 100. The ratios r F and r 2 are plotted in Figures 5.1-5.3.

Figure

5.4 plots the condition numbers # 2 the condition numbers
show a very similar behavior. Table 5.2 reports the number
of comparisons used by the BBK pivoting strategy on these matrices for each n; the
maximum number of comparisons is less than n 2 in each case.
In

Figure

5.5 we report results for three nonrandom matrices from the Test Matrix
1106 SHEUNG HUN CHENG AND NICHOLAS J. HIGHAM
matrix
a) n=50, eig. range[-1,10000]
matrix
c) n=50, eig. range[-1,1]
eig. range[-10000,-1]
matrix
eig. range[-1,10000]
matrix
d) n=50, eig. range[-1,1]
matrix
f) n=50, eig. range[-10000,-1]
Fig. 5.2. Measures of E for
-, MC -.

Table
Number of comparisons for BBK pivoting strategy.
n:
mean 343.9 1432.8 5998.4
Toolbox. Clement is a tridiagonal matrix with eigenvalues plus and minus the numbers
n- 1, n- 3, n- 5, . , (1 or 0). Dingdong is the symmetric n-n Hankel matrix with
(i, whose eigenvalues cluster around #/2 and -#/2.
Ipjfact is the Hankel matrix with (i,
Our conclusions from the experiments are as follows.
1. None of the three algorithms is uniformly better than the others in terms of
producing a small perturbation E, whichever measure r F or r 2 is used. All
three algorithms can produce values of r F and r 2 significantly greater than
1, depending on the problem.
2. Algorithm MC often achieves its aim of producing r F # 1. It produced r F of
for the eigenvalue distribution [-1, for each n, and the values
of # 2 (LL T ) (not shown here) were approximately 100r F in each such case.
MODIFIED CHOLESKY ALGORITHM 1107
matrix
a) n=100, eig. range[-1,10000]
matrix
c) n=100, eig. range[-1,1]
eig. range[-10000,-1]
matrix
eig. range[-1,10000]
matrix
d) n=100, eig. range[-1,1]
matrix
f) n=100, eig. range[-10000,-1]
Fig. 5.3. Measures of E for
-, MC -.
However, often r F was of order 1 when #
large value of # 2 (LL T ) is only a necessary condition, not a su#cient one, for
poor performance of Algorithm MC; in other words, the bounds of section 3
can be weak.
3. The condition numbers # 2 (A+E) vary greatly among the algorithms. Our experience
is that for # u#A# Algorithm MC fairly consistently produces
condition numbers of order 100/ # u; the condition number is, as predicted by
(3.6), much smaller for the random matrices with eigenvalues on the range
because the algorithm attempts to perturb all the eigenvalues to
#. The condition numbers produced by the GMW and SE algorithms vary
greatly with the type of matrix.
The fact that r F is close to 1 for the random matrices with eigenvalues in the
range [-10 4 , -1] for Algorithm MC is easily explained. Let A be negative definite.
Then Algorithm MC computes P Hence
matrix
a) n=25, eig. range[-1,10000]
matrix
eig. range[-1,1]
matrix
c) n=25, eig. range[-10000,-1]
Fig. 5.4. Condition numbers # 2 E) for
GMW -, SE -, MC -.
using (3.7), so r F can exceed 1 only by a tiny amount for Algorithm MC applied to
a negative definite matrix, irrespective of # 2 (LL T ).
6. Concluding remarks. Algorithm MC, based on the symmetric indefinite
factorization with the bounded Bunch-Kaufman pivoting strategy, merits consideration
as an alternative to the algorithms of Gill, Murray, and Wright and Schnabel and
Eskow. The results in section 5 suggest that the new algorithm is competitive with
the GMW and SE algorithms in terms of the objectives (O1)-(O4) listed in section 1.
Algorithm MC has the advantages that the extent to which it satisfies the objectives
is neatly, although not sharply, described by the bounds of section 3 and that it can be
implemented by augmenting existing software with just a small amount of additional
code.
Since all three modified Cholesky algorithms can "fail," that is, they can produce
unacceptably large perturbations, it is natural to ask how failure can be detected and
what should be done about it. The GMW and SE algorithms produce their (diagonal)
MODIFIED CHOLESKY ALGORITHM 1109
a) clement(n,1)
c) dingdong(n)
d) dingdong(n)
f) ipjfact(n,1)
Fig. 5.5. Measures of E for three nonrandom matrices. Key: GMW -, SE -, MC -.
perturbations explicitly, so it is trivial to evaluate their norms. For Algorithm MC,
the perturbation to A is (see which would require
operations to form explicitly. However, we can estimate #E# using the norm
estimator from [13] (which is implemented in LAPACK). The estimator requires the
formation of products Ex for certain vectors x, and these can be computed in O(n 2 )
operations; the estimate produced is a lower bound that is nearly always within a
factor 3 of the true norm. For all three algorithms, then, we can inexpensively test
whether the perturbation produced is acceptably small. Unfortunately, for none of
the algorithms is there an obvious way to improve a modified Cholesky factorization
that makes too big a perturbation; whether improvement is possible, preferably
cheaply, is an open question. Of course one can always resort to computing an optimal
perturbation by computing the eigensystem of A and using the formulae in
Theorem 3.1.
We note that we have approached the problem of modified Cholesky factorization
from a purely linear algebra perspective. An important test of a modified Cholesky
algorithm is to evaluate it in an optimization code on representative problems, as
was done by Schlick [20] for the GMW and SE algorithms. This we plan to do for
Algorithm MC in future work.
Finally, we mention that a generalization of the modified Cholesky problem motivated
by constrained optimization is analyzed in detail in [17].



--R


Accurate symmetric indefinite linear equation solvers
Some stable methods for calculating inertia and solving symmetric linear systems
Direct methods for solving symmetric indefinite systems of linear equations


Society for Industrial and Applied Mathematics
The factorization of sparse symmetric indefinite matrices
Direct solution of sets of linear equations whose matrix is sparse


Computing a nearest symmetric positive semidefinite matrix
FORTRAN codes for estimating the one-norm of a real or complex matrix
The Test Matrix Toolbox for Matlab (Version 3.0)
Accuracy and Stability of Numerical Algorithms
Stability of the diagonal pivoting method with partial pivoting
Modifying the inertia of matrices arising in optimiza- tion
Matrix Analysis

Modified Cholesky factorizations for sparse preconditioners
A new modified Cholesky factorization

--TR
modified Cholesky factorization;newton's method;symmetric indefinite factorization;optimization