Multiple-Rank Modifications of a Sparse Cholesky Factorization.
Given a sparse symmetric positive definite matrix $\mathbf{AA}\tr$ and an associated sparse Cholesky factorization $\mathbf{LDL}\tr$ or $\mathbf{LL}\tr$, we develop sparse techniques for updating the factorization after either adding a collection of columns to A or deleting a collection of columns from A. Our techniques are based on an analysis and manipulation of the underlying graph structure, using the framework developed in an earlier paper on rank-1 modifications [T. A. Davis and W. W. Hager, SIAM J. Matrix Anal. Appl., 20 (1999), pp. 606--627]. Computationally, the multiple-rank update has better memory traffic and executes much faster than an equivalent series of rank-1 updates since the multiple-rank update makes one pass through L computing the new entries, while a series of rank-1 updates requires multiple passes through L.
Introduction
. This paper presents a method for evaluating a multiple rank
update or downdate of the sparse Cholesky factorization LDL T or LL T of the matrix
AA T , where A is m by n. More precisely, given an m r matrix W, we evaluate the
Cholesky factorization of AA T either  is +1 (corresponding to an
update) and W is arbitrary, or  is 1 (corresponding to a downdate) and W consists
of columns of A. Both AA T and AA T +WW T must be positive denite. It follows
that in the case of an update, and n r  m in the case of a downdate.
One approach to the multiple rank update is to express it as a series of rank-1
updates and use the theory developed in [10] for updating a sparse factorization after
a rank-1 change. This approach, however, requires multiple passes through L as it is
updated after each rank-1 change. In this paper, we develop a sparse factorization
algorithm that makes only one pass through L.
For a dense Cholesky factorization, a one-pass algorithm to update a factorization
is obtained from Method C1 in [18] by making all the changes associated with one
column of L before moving to the next column, as is done in the following algorithm
that overwrites L and D with the new factors of AA T
performs
oating-point operations.
Algorithm 1 (Dense rank-r update/downdate).
to r do
end for
do
to r do
This work was supported by the National Science Foundation.
y davis@cise.u
.edu/~davis, PO Box 116120, Department of Computer
and Information Science and Engineering, University of Florida, Gainesville, FL 32611-6120. Phone
(352) 392-1481. Fax (352) 392-1220. TR-99-006 (June 1999, revised Sept. 2000)
z hager@math.u
.edu/~hager, PO Box 118105, Department of Mathe-
matics, University of Florida, Gainesville, FL 32611-8105. Phone (352) 392-0281. Fax (352) 392-8357.
A. DAVIS AND WILLIAM W. HAGER
end for
do
to r do
l
end for
end for
end for
We develop a sparse version of this algorithm that only accesses and modies those
entries in L and D which can change. For the theory in our rank-1 paper [10]
shows that those columns which can change correspond to the nodes in an elimination
tree on a path starting from the node k associated with the rst nonzero element w k1
in W. For r > 1 we show that the columns of L which can change correspond
to the nodes in a subtree of the elimination tree, and we express this subtree as a
modication of the elimination tree of AA T . Also, we show that with a reordering
of the columns of W, it can be arranged so that in the inner loop where elements in
row p of W are updated, the elements that change are adjacent to each other. The
sparse techniques that we develop lead to sequential access of matrix elements and to
e-cient computer memory tra-c. These techniques to modify a sparse factorization
have many applications including the Linear Program Dual Active Set Algorithm
least-squares problems in statistics, the analysis of electrical circuits
and power systems, structural mechanics, sensitivity analysis in linear programming,
boundary condition changes in partial dierential equations, domain decomposition
methods, and boundary element methods (see [19]).
Section 2 describes our notation. In section 3, we present an algorithm for computing
the symbolic factorization of AA T using multisets, which determines the location
of nonzero entries in L. Sections 4 and 5 describe our multiple rank symbolic update
and downdate algorithms for nding the nonzero pattern of the new factors. Section 6
describes our algorithm for computing the new numerical values of L and D, for either
an update or downdate. Our experimental results are presented in Section 7.
2. Notation and background. Given the location of the nonzero elements of
AA T , we can perform a symbolic factorization (this terminology is introduced by
George and Liu in [15]) of the matrix to predict the location of the nonzero elements
of the Cholesky factor L. In actuality, some of these predicted nonzeros may be
zero due to numerical cancellation during the factorization process. The statement
will mean that l ij is symbolically nonzero. The main diagonals of L and
D are always nonzero since the matrices that we factor are positive denite (see [26,
p. 253]). The nonzero pattern of column j of L is denoted
while L denotes the collection of patterns:
Similarly, A j denotes the nonzero pattern of column j of A,
while A is the collection of patterns:
The elimination tree can be dened in terms of a parent map  (see [22]). For any
node j, (j) is the row index of the rst nonzero element in column j of L beneath
the diagonal element:
where \min X" denotes the smallest element of
i:
Our convention is that the min of the empty set is zero. Note that j < (j) except
in the case where the diagonal element in column j is the only nonzero element. The
children of node j is the set of nodes whose parent is j:
The ancestors of a node j, denoted P(j), is the set of successive parents:
for each j, the ancestor sequence is nite. The sequence of nodes
j, (j), ((j)),   , forming P(j), is called the path from j to the associated tree
root, the nal node on the path. The collection of paths leading to a root form an
elimination tree. The set of all trees is the elimination forest. Typically, there is a
single tree whose root is m, however, if column j of L has only one nonzero element,
the diagonal element, then j will be the root of a separate tree.
The number of elements (or size) of a set X is denoted jX j, while jAj or jLj denote
the sum of the sizes of the sets they contain.
3. Symbolic factorization. For a matrix of the form AA T , the pattern L j of
column j is the union of the patterns of each column of L whose parent is j and each
column of A whose smallest row index of its nonzero entries is j (see [16, 22]):
min Ak=j
To modify (3.1) during an update or downdate, without recomputing it from
scratch, we need to keep track of how each entry i entered into L j [10]. For example,
if (c) changes, we may need to remove a term L c n fcg. We cannot simply perform
a set subtraction, since we may remove entries that appear in other terms. To keep
track of how entries enter and leave the set L j , we maintain a multiset associated
with column j. It has the form
4 TIMOTHY A. DAVIS AND WILLIAM W. HAGER
where the multiplicity m(i; j) is the number of children of j that contain row index i
in their pattern plus the number of columns of A whose smallest entry is j and that
contain row index i. Equivalently, for i 6= j,
For we increment the above equation by one to ensure that the diagonal entries
never disappear during a downdate. The set L j is obtained from L ]
by removing the
multiplicities.
We dene the addition of a multiset X ] and a set Y in the following way:
where
Similarly, the subtraction of a set Y from a multiset X ] is dened by
where
The multiset subtraction of Y from X ] undoes a prior addition. That is, for any
multiset X ] and any set Y , we have
In contrast ((X [ Y) n Y) is equal to X if and only if X and Y are disjoint sets.
Using multiset addition instead of set union, (3.1) leads to the following algorithm
for computing the symbolic factorization of AA T .
Algorithm 2 (Symbolic factorization of AA T , using multisets).
do
for each c such that do
end for
for each k where min A do
end for
end for
4. Multiple rank symbolic update. We consider how the pattern L changes
when AA T is replaced by AA T +WW T . Since
we can in essence augment A by W in order to evaluate the new pattern of column
in L. According to (3.1), the new pattern L j of column j of L after the update is
min Ak=j
A
where W i is the pattern of column i in W. Throughout, we put a bar over a matrix
or a set to denote its new value after the update or downdate.
In the following theorem, we consider a column j of the matrix L, and how its
pattern is modied by the sets W i . Let L ]
j denote the multiset for column j after the
rank-r update or downdate has been applied.
Theorem 4.1. To compute the new multiset L ]
j and perform
the following modications:
Case A: For each i such that to the pattern for column
Case B: For each c such that
(c is a child of j in both the old and new elimination tree).
Case C: For each c such that
(c is a child of j in the new tree, but not the old one).
Case D: For each c such that
(c is a child of j in the old tree, but not the new one).
Proof. Cases A{D account for all the adjustments we need to make in L j in order
to obtain L j . These adjustments are deduced from a comparison of (3.1) with (4.1).
In case A, we simply add in the W i multisets of (4.1) that do not appear in (3.1). In
case B, node c is a child of node j both before and after the update. In this case, we
must adjust for the deviation between L c and L c . By [10, Prop. 3.2], after a rank-1
update, L c  L c . If w i denotes the i-th column of W, then
Hence, updating AA T by WW T is equivalent to r successive rank-1 updates of AA T .
By repeated application of [10, Prop. 3.2], L c  L c after a rank-r update of AA T . It
6 TIMOTHY A. DAVIS AND WILLIAM W. HAGER
follows that L c and L c deviate from each other by the set L c n L c . Consequently, in
case B we simply add in L c n L c .
In case C, node c is a child of j in the new elimination tree, but not in the old
tree. In this case we need to add in the entire set L c n fcg since the corresponding
term does not appear in (3.1). Similarly, in case D, node c is a child of j in the old
elimination tree, but not in the new tree. In this case, the entire set L c nfcg should be
deleted. The case where c is not a child of j in either the old or the new elimination
tree does not result in any adjustment since the corresponding L c term is absent from
both (3.1) and (4.1).
An algorithm for updating a Cholesky factorization that is based only on this
theorem would have to visit all nodes j from 1 to m, and consider all possible children
c < j. On the other hand, not all nodes j from 1 to m need to be considered since not
all columns of L change when AA T is modied. In [10, Thm. 4.1] we show that for
the nodes whose patterns can change are contained in P(k 1 ) where we dene
. For a rank-r update, let P (i) be the ancestor map associated with the
elimination tree for the Cholesky factorization of the matrix
Again, by [10, Thm. 4.1], the nodes whose patterns can change during the rank-r
update are contained in the union of the patterns Although we
could evaluate each i, it is di-cult to do this e-ciently since we need to
perform a series of rank-1 updates and evaluate the ancestor map after each of these.
On the other hand, by [10, Prop. 3.1] and [10, Prop. 3.2], P (i) (j)  P (i+1) (j) for each
and j, from which it follows that P (i) Consequently, the
nodes whose patterns change during a rank-r update are contained in the set
1ir
Theorem 4.2, below, shows that any node in T is also contained in one or more
of the sets P (i) (k i ). From this it follows that the nodes in T are precisely those nodes
for which entries in the associated columns of L can change during a rank-r update.
Before presenting the theorem, we illustrate this with a simple example shown in

Figure

4.1. The left of Figure 4.1 shows the sparsity pattern of original matrix AA T ,
its Cholesky factor L, and the corresponding elimination tree. The nonzero pattern
of the rst column of W is 2g. If performed as a single rank-1 update, this
causes a modication of columns 1, 2, 6, and 8 of L. The corresponding nodes in the
original tree are encircled; these nodes form the path P (1) 8g from node
1 to the root (node 8) in the second tree. The middle of Figure 4.1 shows the matrix
after this rank-1 update, and its factor and elimination tree. The entries in the second
1 that dier from the original matrix AA T are shown as small
pluses. The second column of W has the nonzero pattern W 7g. As a rank-1
update, this aects columns P (2) of L. These columns
form a single path in the nal elimination tree shown in the right of the gure.
For the rst rank-1 update, the set of columns that actually change are P (1)
8g. This is a subset of the path in the nal tree. If we
use P(1) to guide the work associated with column 1 of W, we visit all the columns
after second update
after first update
elimination tree
Elimination tree
After first update
Elimination tree
Original factor L Factor after second update
Factor after first update
Original matrix
After second update1
A T
A T6743
A T
Original
A
Fig. 4.1. Example rank-2 update
that need to be modied, plus column 7. Node 7 is in the set of nodes P(3) aected
by the second rank-1 update, however, as shown in the following theorem.
Theorem 4.2. Each of the paths contained in T and conversely, if
contained in P (i)
Proof. Before the theorem, we observe that each of the paths contained
in T . Now suppose that some node j lies in the tree T . We need to prove that it is
contained in P (i) s be the largest integer such that P(k s ) contains
j and let c be any child of j in T . If c lies on the path P(k i ) for some i, then j lies
on the path P(k i ) since j is the parent of c. Since j does not lie on the path P(k i )
for any i > s, it follows that c does not lie on the path P(k i ) for any i > s. Applying
this same argument recursively, we conclude that none of the nodes on the subtree
of T rooted at j lie on the path P(k i ) for any i > s. Let T j denote the subtree of T
rooted at j. Since contained in P(k i ) for each i, none of the nodes of T j
lie on any of the paths Thm. 4.1], the patterns of all nodes
outside the path are unchanged for each i. Let L (i)
c be the pattern of column
c in the Cholesky factorization of (4.2). Since any node c contained in T j does not lie
8 TIMOTHY A. DAVIS AND WILLIAM W. HAGER
(d,c)
(b,e)
e
f
c
d
a
Fig. 4.2. Example rank-8 symbolic update and subtree T
on any of the paths
c for all i, l  s. Since k s is a node
of T j , the path P must include j.

Figure

4.2 depicts a subtree T for an example rank-8 update. The subtree consists
of all those nodes and edges in one or more of the paths P(k 1
These paths form a subtree, and not a general graph, since they are all paths from
an initial node to the root of the elimination tree of the matrix L. The subtree T
might actually be a forest, if L has an elimination forest rather than an elimination
tree. The rst nonzero positions in w 1 through w 8 correspond to nodes k 1 through
k 8 . For this example node k 4 happens to lie on the path P (1) (k 1 ). Nodes at which
paths rst intersect are shown as smaller circles, and are labeled a through f . Other
nodes along the paths are not shown. Each curved arrow denotes a single subpath.
For example, the arrow from nodes b to e denotes the subpath from b to e in P(b).
This subpath is denoted as P(b; e) in Figure 4.2.
The following algorithm computes the rank-r symbolic update. It keeps track of
an array of m \path-queues," one for each column of L. Each queue contains a set
of path-markers in the range 1 to r, which denote which of the paths P(k 1 ) through
next. If two paths have merged, only one of the paths
needs to be considered (we arbitrarily select the higher-numbered path to represent
the merged paths). This set of path-queues requires O(m Removing and
inserting a path-marker in a path-queue takes O(1) time. The only outputs of the
algorithm are the new pattern of L and its elimination tree, namely, L ]
and (j) for
all columns are aected by the rank-r update. We dene L
and node j not in T .
Case C will occur for c and j prior to visiting column (c), since
We thus place c in the lost-child-queue of column (c) when encountering case C
for nodes c and j. When the algorithm visits node (c), its lost-child-queue will
contain all those nodes for which case D holds. This set of lost-child-queues is not
the same as the set of path-queues (although there is exactly one lost-child-queue and
one path-queue for each column j of L).
Algorithm 3 (Symbolic rank-r update, add new matrix W).
Find the starting nodes of each path
to r do
place path-marker i in path-queue of column k i
end for
Consider all columns corresponding to nodes in the paths P(k 1
to m do
if path-queue of column j is non-empty do
for each path-marker i on path-queue of column j do
Let c be the prior column on this path (if any), where
do
Case A: j is the rst node on the path P(k i ), no prior c
else if
Case B: c is an old child of j, possibly changed
else
Case C: c is a new child of j and a lost child of (c)
place c in lost-child-queue of column (c)
endif
end for
Case D: consider each lost child of j
for each c in lost-child-queue of column j do
end for
Move up one step in the path(s)
Let i be the largest path-marker in path-queue of column j
Place path-marker i in path-queue of column (j)
if path-queue of column j non-empty
end for
The optimal time for a general rank-r update is
A. DAVIS AND WILLIAM W. HAGER
The actual time taken by Algorithm 3 only slightly higher, namely,
because of the O(m) book-keeping required for the path-queues. In most practical
cases, the O(m) term will not be the dominant term in the run time.
Algorithm 3 can be used to compute an entire symbolic factorization. We start
by factorizing the identity matrix I = II T into LDL III. In this case, we have
j. The initial elimination tree is a forest of m nodes and no
edges. We can now determine the symbolic factorization of I +AA T using the rank-r
update algorithm above, with m. This matrix has identical symbolic
factors as AA T . Case A will apply for each column in A, corresponding to the
min Ak=j
term in (3.1). Since (c) = 0 for each c, cases B and D will not apply. At column j,
case C will apply for all children in the elimination tree, corresponding to the
term in (3.1). Since duplicate paths are discarded when they merge, we modify
each column j once, for each child c in the elimination tree. This is the same work
performed by the symbolic factorization algorithm, Algorithm 2, which is O(jLj).
Hence, Algorithm 3 is equivalent to Algorithm 2 when we apply it to the update
I +AA T . Its run time is optimal in this case.
5. Multiple rank symbolic downdate. The downdate algorithm is analogous.
The downdated matrix is AA T WW T where W is a subset of the columns of A.
In a downdate, P(k)  P(k), and thus rather than following the paths P(k i ), we
follow the paths P(k i ). Entries are dropped during a downdate, and thus L j  L j
and (j)  (j). We start with L ]
j and make the following changes.
Case A: If then the pattern W i is removed from
column j,
Case B: If then c is a child of j in both the
old and new tree. We need to remove from L ]
entries in the old pattern L c
but not in the new pattern L c ,
Case C: If for some node c, then c is a child of j in the old
elimination tree, but not the new tree. We compute
MULTIPLE-RANK MODIFICATIONS 11
Case D: If for some node c, then c is a child of j in the new
tree, but not the old one. We compute
Case C will occur for c and j prior to visiting column (c), since
We thus place c in the new-child-queue of (c) when encountering case C for nodes c
and j. When the algorithm visits node (c), its new-child-queue will contain all those
nodes for which case D holds.
Algorithm 4 (Symbolic rank-r downdate, remove matrix W).
Find the starting nodes of each path
to r do
place path-marker i in path-queue of column k i
end for
Consider all columns corresponding to nodes in the paths P(k 1
to m do
if path-queue of column j is non-empty do
for each path-marker i on path-queue of column j do
Let c be the prior column on this path (if any), where
do
Case A: j is the rst node on the path P(k i ), no prior c
else if
Case B: c is an old child of j, possibly changed
else
Case C: c is a lost child of j and a new child of (c)
place c in new-child-queue of column (c)
endif
end for
Case D: consider each new child of j
for each c in new-child-queue of j do
end for
Move up one step in the path(s)
Let i be the largest path-marker in path-queue of column j
Place path-marker i in path-queue of column (j)
if path-queue of column j non-empty
end for
A. DAVIS AND WILLIAM W. HAGER
The time taken by Algorithm 4 is
which slightly higher than the optimal time,
In most practical cases, the O(m) term in the asymptotic run time for Algorithm 4
will not be the dominant term.
6. Multiple rank numerical update and downdate. The following numerical
rank-r update/downdate algorithm, Algorithm 5, overwrites L and D with the
updated or downdated factors. The algorithm is based on Algorithm 1, the one-pass
version of Method C1 in [18] presented in Section 1. The algorithm is used after
the symbolic update algorithm (Algorithm 3) has found the subtree T corresponding
to the nodes whose patterns can change, or after the symbolic downdate algorithm
(Algorithm 4) has found T . Since the columns of the matrix W can be reordered
without aecting the product WW T , we reorder the columns of W using a depth-rst
search [6] of T (or T ) so that as we march through the tree, consecutive columns
of W are utilized in the computations. This reordering improves the numerical up-
date/downdate algorithm by placing all columns of W that aect any given subpath
next to each other, eliminating an indexing operation. Reordering the columns of
a sparse matrix prior to Cholesky factorization is very common [3, 22, 23, 25]. It
improves data locality and simplies the algorithm, just as it does for reordering W
in a multiple rank update/downdate. The depth rst ordering of the tree changes as
the elimination tree changes, so columns of W must be ordered for each update or
downdate.
To illustrate this reordering, consider the subtree T in Figure 4.2 for a rank-8
update. If the depth-rst-search algorithm visits child subtrees from left to right, the
resulting reordering is as shown in Figure 6.1. Each subpath in Figure 6.1 is labeled
with the range of columns of W that aect that subpath, and with the order in which
the subpath is processed by Algorithm 5. Consider the path from node c to e. In

Figure

4.2, the columns of L corresponding to nodes on this subpath are updated by
columns 2, 8, 3, and 5 of W, in that order. In the reordered subtree (Figure 6.1), the
columns on this subpath are updated by columns 5 through 8 of the reordered W.
Algorithm 5 (Sparse numeric rank-r modication, add WW T ).
The columns of W have been reordered.
to r do
end for
for each subpath in depth-rst-search order in T
Let c 1 through c 2 be the columns of W that aect this subpath
for each column j in the subpath do
do
e
f
c
d
a
13th
2nd
1st
6th
3rd 4th
7th
8th
9th
10th
11th
12th
3 47-885th
Fig. 6.1. Example rank-8 update after depth-rst-search reordering
end for
for all
do
l
end for
end for
end for
end for
The time taken by r rank-1 updates [10] is
O@ r
where L (i)
j is the pattern of column j after the i-th rank-1 update. This time is
asymptotically optimal. A single rank-r update cannot determine the paths
but uses P(k i ) instead. Thus, the time taken by Algorithm 5 for a rank-r update is
O@ r
This is slightly higher than (6.1), because
14 TIMOTHY A. DAVIS AND WILLIAM W. HAGER

Table
Dense matrix performance for 64-by-64 matrices and 64-by-1 vectors
operation M
ops
DGEMM (matrix-matrix multiply) 171.6
DGEMV (matrix-vector multiply) 130.0
DTRSV (solve
DAXPY (the vector computation
DDOT (the dot product
the i-th column of W does not necessarily aect all of the columns
in the path P(k i ). If w i does not aect column j, then w ji and
will both be zero in
the inner loop in Algorithm 5. An example of this occurs in Figure 4.1, where column
1 of W does not aect column 7 of L. We could check this condition, and reduce the
asymptotic run time to
O@ r
In practice, however, we found that the paths dier much.
Including this test did not improve the overall performance of our algorithm. The
time taken by Algorithm 5 for a rank-r downdate is similar, namely,
O@ r
The numerical algorithm for updating and downdating LL T is essentially the
same as that for LDL T [4, 24]; the only dierence is a diagonal scaling. For either
LL T or LDL T , the symbolic algorithms are identical.
7. Experimental results. To test our methods, we selected the same experiment
as in our earlier paper on the single-rank update and downdate [10], which
mimics the behavior of the Linear Programming Dual Active Set Algorithm [20]. The
rst
consists of 5446 columns from a larger 6071-
arising in an airline scheduling problem
(DFL001) [13]. The 5446 columns correspond to the optimal solution of the linear
programming problem. Starting with an initial LDL T factorization of the matrix
T , we added columns from B (corresponding to an update) until we
obtained the factors of 10 6 I +BB T . We then removed columns in a rst-in-rst-out
order (corresponding to a downdate) until we obtained the original factors. The LP
DASA algorithm would not perform this much work (6784 updates and 6784 down-
dates) to solve this linear programming problem.
Our experiment took place on a Sun Ultra Enterprise running the Solaris 2.6
operating system, with eight 248 Mhz UltraSparc-II processors (only one processor
was used) and 2GB of main memory. The dense matrix performance in millions of
oating-point operations per second (M
ops) of the BLAS [12] is shown in Table 7.1.
All results presented below are for our own codes (except for colmmd, spooles, and
the BLAS) written in the C programming language and using double precision
oating
point arithmetic.
We rst permuted the rows of B to preserve sparsity in the Cholesky factors of
BB T . This can be done e-ciently with colamd [7, 8, 9, 21], which is based on an

Table
Average update and downdate performance results
ops
r in seconds
update downdate update downdate
9
14
approximate minimum degree ordering algorithm [1]. However, to keep our results
consistent with our prior rank-1 update/downdate paper [10], we used the same permutation
as in those experiments (from colmmd [17]). Both colamd and Matlab's
colmmd compute the ordering without forming BB T explicitly. A symbolic factorization
of BB T nds the nonzero counts of each column of the factors. This step takes
an amount of space this is proportional to the number of nonzero entries in B. It
gives us the size of a static data structure to hold the factors during the updating
and downdating process. The numerical factorization of BB T is not required. A
second symbolic factorization nds the rst nonzero pattern L. An initial numerical
factorization computes the rst factors L and D. We used our own non-supernodal
factorization code (similar to SPARSPAK [5, 15]), since the update/downdate algorithms
do not use supernodes. A supernodal factorization code such as spooles [3] or
a multifrontal method [2, 14] can get better performance. The factorization method
used has no impact on the performance of the update and downdate algorithms.
We ran dierent experiments, each one using a dierent rank-r update and
downdate, where r varied from 1 to 16. After each rank-r update, we solved the
sparse linear system LDL T using a dense right-hand side b. To compare the
performance of a rank-1 update with a rank-r update (r > 1), we divided the run time
of the rank-r update by r. This gives us a normalized time for a single rank-1 update.
The average time and M
ops rate for a normalized rank-1 update and downdate for
the entire experiment is shown in Table 7.2. The time for the update, downdate, or
solve increases as the factors become denser, but the performance in terms of M
ops
is fairly constant for all three operations. The rst rank-16 update when the factor
L is sparsest takes 0.47 seconds (0.0294 seconds normalized) and runs at 65.5 M
ops
compared to 65.1 M
ops in

Table

7.2 for the average speed of all the rank-16 updates.
The performance of each step is summarized in Table 7.3. A rank-5 update takes
about the same time as using the updated factors to solve the sparse linear system
even though the rank-5 update performs 2.6 times the work.
The work, in terms of
oating-point operations, varies only slightly as r changes.
With rank-1 updates, the total work for all the updates is 17.293 billion
A. DAVIS AND WILLIAM W. HAGER

Table
Dense matrix performance for 64-by-64 matrices and 64-by-1 vectors
Operation Time (sec) M
ops Notes
colamd ordering 0.45 -
Symbolic factorization (of BB T
Symbolic factorization for rst L 0.46 - 831 thousand nonzeros
Numeric factorization for rst L (our code) 20.07 24.0
Numeric factorization for rst L (spooles) 18.10 26.6
Numeric factorization of BB T (our code) 61.04 18.5 not required
Numeric factorization of BB T (spooles) 17.80 63.3 not required
Average rank-16 update 0.63 65.1 compare with rank-1
Average rank-5 update 0.25 51.0 compare with solve step
Average rank-1 update 0.084 30.3
Average solve LDL T
point operations, or 2.55 million per rank-1 update. With rank-16 updates (the
worst case), the total work increases to 17.318 billion
oating-point operations. The
downdates take a total of 17.679 billion
oating-point operations (2.61 million
per rank-1 downdate), while the rank-16 downdates take a total of 17.691 billion
operations. This conrms the near-optimal operation count of the multiple rank
update/downdate, as compared to the optimal rank-1 update/downdate.
Solving when L is sparse and b is dense, and computing the sparse LDL T
factorization using a non-supernodal method, both give a rather poor computation-
to-memory-reference ratio of only 2/3. We tried the same loop unrolling technique
used in our update/downdate code for our sparse solve and sparse LDL T factorization
codes, but this resulted in no improvement in performance.
A sparse rank-r update or downdate can be implemented in a one-pass algorithm
that has much better memory tra-c than that of a series of r rank-1 modications. In
our numerical experimentation with the DFL001 linear programming test problem, the
rank-r modication was more than twice as fast as r rank-1 modications for r  11.
The superior performance of the multiple rank algorithm can be explained using the
computation-to-memory-reference ratio. If c in Algorithm 5 (a subpath aected
by only one column of W), it can be shown that this ratio is about 4/5 when L j is
large. The ratio when c aected by 16 columns of W) is
about 64/35 when L j is large. Hence, going from a rank-1 to a rank-16 update
improves the computation-to-memory-reference ratio by a factor of about 2.3 when
column j of L has many nonzeros. By comparison, the level-1 BLAS routines for
dense matrix computations (vector computations such as DAXPY and DDOT) [11]
have computation-to-memory-reference ratios between 2/3 and 1. The level-2 BLAS
(DGEMV and DTRSV, for example) have a ratio of 2.
8.

Summary

. Because of improved memory locality, our multiple-rank sparse
update/downdate method is over twice as fast as our prior rank-1 update/downdate
method. The performance of our new method (65.1 M
ops for a sparse rank-16
update) compares favorably with both the dense matrix performance (81.5 M
ops to
solve the dense system and the sparse matrix performance (18.0 M
ops to
solve the sparse system and an observed peak numerical factorization of 63.3
ops in spooles) on the computer used in our experiments. Although not strictly
optimal, the multiple-rank update/downdate method has nearly the same operation
count as the rank-1 update/downdate method, which has an optimal operation count.
MULTIPLE-RANK MODIFICATIONS 17



--R

An approximate minimum degree ordering algorithm
Vectorization of a multiprocessor multifrontal code
SPOOLES: an object-oriented sparse matrix library
A Cholesky up- and downdating algorithm for systolic and SIMD architectures
SPARSPAK: Waterloo sparse matrix package
Introduction to Algorithms
A column approximate minimum degree ordering algorithm
A column approximate minimum degree ordering algorithm

Modifying a sparse Cholesky factorization
Philadelphia: SIAM Publications
A set of level-3 basic linear algebra subprograms
Distribution of mathematical software via electronic mail
The multifrontal solution of inde
Computer Solution of Large Sparse Positive De
A data structure for sparse QR and LU factorizations
Sparse matrices in MATLAB: design and implementation
Methods for modifying matrix factorizations
Updating the inverse of a matrix

An approximate minimum degree column ordering algorithm
The role of elimination trees in sparse factorization
A supernodal Cholesky factorization algorithm for shared-memory multiprocessors


New York
--TR

--CTR
W. Hager, The Dual Active Set Algorithm and Its Application to Linear Programming, Computational Optimization and Applications, v.21 n.3, p.263-275, March 2002
Ove Edlund, A software package for sparse orthogonal factorization and updating, ACM Transactions on Mathematical Software (TOMS), v.28 n.4, p.448-482, December 2002
Matine Bergounioux , Karl Kunisch, Primal-Dual Strategy for State-Constrained Optimal Control Problems, Computational Optimization and Applications, v.22 n.2, p.193-224, July 2002
Nicholas I. M. Gould , Jennifer A. Scott , Yifan Hu, A numerical evaluation of sparse direct solvers for the solution of large sparse symmetric linear systems of equations, ACM Transactions on Mathematical Software (TOMS), v.33 n.2, p.10-es, June 2007
numerical linear algebra;matrix updates;cholesky factorization;sparse matrices;mathematical software;direct methods