Modifying a Sparse Cholesky Factorization.
Given a sparse symmetric positive definite matrix ${\bf AA}^{\sf T}$ and an associated sparse Cholesky factorization ${\bf LDL}^{\sf T}$ or ${\bf LL}^{\sf T}$, we develop sparse techniques for obtaining the new factorization associated with either adding a column to ${\bf A}$ or deleting a column from ${\bf A}$. Our techniques are based on an analysis and manipulation of the underlying graph structure and on ideas of Gill et al.\ [ Math. Comp., 28 (1974), pp. 505--535] for modifying a dense Cholesky factorization. We show that our methods extend to the general case where an arbitrary sparse symmetric positive definite matrix is modified. Our methods are optimal in the sense that they take time proportional to the number of nonzero entries in ${\bf L}$ and ${\bf D}$ that change.
Introduction
This paper presents a method for updating and downdating the sparse Cholesky
factorization LL T of the matrix AA T , where A is m by n. More precisely, we evaluate
the Cholesky factorization of AA T either oe is +1 (corresponding to
an update) and w is arbitrary, or oe is \Gamma1 (corresponding to a downdate) and w is a
column of A. Both AA T and AA T must be symmetric and positive definite.
The techniques we develop for the matrix AA T can be extended to determine the
effects on the Cholesky factors of a general symmetric positive definite matrix M of
any symmetric change of the form M preserves positive definiteness.
Our methods are optimal in the sense that they take time proportional to the number
of nonzero entries in L that change.
There are many applications of the techniques presented in this paper. In
the Linear Program Dual Active Set Algorithm (LP DASA) [19], the A matrix
corresponds to the basic variables in the current basis of the linear program, and
in successive iterations, we bring variables in and out of the basis, leading to changes
of the form AA T Other application areas where the techniques developed
in this paper are applicable include least-squares problems in statistics, the analysis
of electrical circuits, structural mechanics, sensitivity analysis in linear programming,
boundary condition changes in partial differential equations, domain decomposition
methods, and boundary element methods. For a discussion of these application areas
and others, see [18].
Section 2 introduces our notation. For an introduction to sparse matrix
techniques, see [4, 8]. In x3 we discuss the structure of the nonzero elements in
the Cholesky factorization of AA T and in x4, we discuss the structure associated
with the Cholesky factors of AA T . The symbolic update and downdate
methods provide the framework for our sparse version of Method C5 of Gill, Golub,
Murray, and Saunders [15] for modifying a dense Cholesky factorization. We discuss
our sparse algorithm in x5. Section 6 presents the general algorithm for modifying
the sparse Cholesky factorization for any sparse symmetric positive definite matrix.
Implementation details for the modification algorithm associated with AA T
are summarized in x7, while the results of a numerical experiment with a large
optimization problem from Netlib [3] are presented in x8. Section 9 concludes with a
discussion of future work.
Notation
Throughout the paper, matrices are capital bold letters like A or L, while vectors are
lower case bold letters like x or v. Sets and multisets are in calligraphic style like A,
L, or P. Scalars are either lower case Greek letters, or italic style like oe, k, or m.
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 [8]) 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 "l ij 6= 0" will mean that l ij is symbolically
nonzero. The diagonal of L is always nonzero since the matrices that we factor are
positive definite (see [24, 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 defined in terms of a parent map - (see [20]), which
gives the row index associated with a given node j of the first nonzero element in
column j of L beneath the diagonal element:
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
of the parent map is the children multifunction. That is, the children of
node k is the set defined by
The ancestors of a node j, denoted P(j), is the set of successive parents:
Here the powers of a map are defined in the usual way: - 0 is the identity while - i
is the i-fold composition of - with itself. The sequence of nodes j, -(j),
forming P(k), is called the path from j to the associated tree root. The
collection of paths leading to a root form an elimination tree, and the set of all trees
is the elimination forest. Typically, there is a single tree whose root is m, however,
if column j of AA T 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. Define the directed graph G(M) of an m
by m matrix M with nonzero pattern M as the vertex and edge sets:
is the vertex set and is the directed
edge set.
3 Symbolic factorization
Any approach for generating the pattern set L is called
23]. The symbolic factorization of a matrix of the form AA T is given in Algorithm 1
(see [9, 20]).
Algorithm 1 (Symbolic factorization of AA T )
do
min A k =j
A kA
end for
Algorithm 1 basically says that the pattern of column j of L can be expressed as
the union of the patterns of each column of L whose parent is j and the patterns of the
columns of A whose first nonzero element is j. The elimination tree, connecting each
child to its parent, is easily formed during the symbolic factorization. Algorithm 1
can be done in O(jLj
1 Asymptotic complexity notation is defined in [2]. We write there exist positive
constants c and n 0
such that 0 - f(n) - cg(n) for all n ? n 0
Observe that the pattern of the parent of node j contains all entries in the pattern
of column j except j itself [22]. That is,
Proceeding by induction, if k is an ancestor of j, then
This leads to the following relation between L j and the path P(j). The first part of
this proposition, and its proof, is given in [22]. Our proof differs slightly from the one
in [22]. We include it here since the same proof technique is exploited later.
Proposition 3.1 For each j, we have L j ' P(j); furthermore, for each k and j 2
Proof. Obviously, be any given element of L j with i 6= j. Since
we see that the following relation holds for
Now suppose that (2) holds for some integer l - 0, and let k denote - l (j). By (1)
and the fact that k - i, we have which implies that
Hence, either or (2) holds with l replaced by l + 1. Since (2) is violated
for l sufficiently large, we conclude that there exists an l for which
Consequently, each element of L j is contained in P(j), we have
is an ancestor of k, and P(j) ' P(k).
Since we have already shown that L j ' P(j), the proof is complete.
As we will see, the symbolic factorization of AA T can be obtained by
updating the symbolic factorization of AA T using an algorithm that has the same
structure as that of Algorithm 1. The new pattern L j is equal to the old pattern L j
union entries that arise from new children and from the pattern of the new column
(We put a bar over a matrix or a set or a multiset to denote its value after the
update or downdate is complete.)
Downdating is not as easy. Once a set union has been computed, it cannot be
undone without knowledge of how entries entered the set. We can keep track of
this information by storing the elements of L as multisets rather than as sets. The
multiset associated with column j has the form
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,
With this definition, we can undo a set union by subtracting multiplicities.
We now define some operations involving multisets. First, if X ] is a multiset
consisting of pairs (i; m(i)) where m(i) is the multiplicity associated with i, then X
is the set obtained by removing the multiplicities. In other words, the multiset X
and the associated base set X satisfy the relation:
We define 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 defined 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.
Algorithm 2 below performs a symbolic factorization of AA T , with each set union
operation replaced by a multiset addition. This algorithm is identical to Algorithm 1
except for the bookkeeping associated with multiplicities.
Algorithm 2 (Symbolic factorization of AA T , using multisets)
do
for each c 2 - \Gamma1 (j) do
end for
for each k where min A do
end for
end for
We conclude this section with a result concerning the relation between the patterns
of AA T and the patterns of AA T +ww T .
Proposition 3.2 Let C and D be the patterns associated with the symmetric positive
definite matrices C and D respectively. Neglecting numerical cancellation, C j ' D j
for each j implies that (L C for each j, where L C and LD are the patterns
associated with the Cholesky factors of C and D respectively.
Proof. In [8, 21] it is shown that an edge (i; j) is contained in the graph of the
Cholesky factor of a symmetric positive definite matrix C if and only if there is a path
from i to j in the graph of C with each intermediate vertex of the path between 1 and
min fi; jg. If C j ' D j for each j, then the paths associated with the graph of C are
a subset of the paths associated with the graph of D. It follows that (L C
for each j.
Ignoring numerical cancellation, the edges in the graph of AA T are a subset of the
edges in the graph of AA T +ww T . By Proposition 3.2, we conclude that the edges in
the graphs of the associated Cholesky factors satisfy the same inclusion. As a result,
if the columns of A and the vectors w used in the update AA T +ww T are all chosen
from the columns of some fixed matrix B, and if a fill-reducing permutation P can be
found for which the Cholesky factors of PBB T P T are sparse ([1] for example), then
by Proposition 3.2, the Cholesky factors of PAA T P T and of P(AA T +ww T )P T will
be at least as sparse as those of PBB T P T .
4 Modifying the symbolic factors
Let A be the modified version of A. Again, we put a bar over a matrix or a set or a
multiset to denote its value after the update or downdate is complete. In an update
A is obtained from A by appending the column w on the right, while in a downdate,
A is obtained from A by deleting the column w from A. Hence, we have
where oe is either +1 and w is the last column of A (update) or oe is \Gamma1 and w is a
column of A (downdate). Since A and A differ by at most a single column, it follows
from Proposition 3.2 that L j ' L j for each j during an update, while L j ' L j during
a downdate. Moreover, the multisets associated with the Cholesky factor of either the
updated or downdated matrix have the structure described in the following theorem:
Theorem 4.1 Let k be the index associated with the first nonzero component of w.
For an update, P(k) ' P(k) and L ]
the complement of P(k).
That is, L
i for all i except when i is k or one of the new ancestors of k. For
a downdate, P(k) ' P(k) and L ]
. That is, L
except when i is k or one of the old ancestors of k.
Proof. To begin, let us consider an update. We will show that each element of P(k)
is a member of P(k) as well. Clearly, k lies in both P(k) and P(k). Proceeding by
induction, suppose that - 0 (k),
We need to show that to complete the induction.
P(k), the induction step is complete, and - j+1
If -(l) 6= -(l), then by Proposition 3.2, -(l) ! -(l), and the following relation
holds for
Now suppose that (3) holds for some integer p - 1, and let q denote - p (l). By
Proposition 3.2, -(l) 2 L l ' L l , and combining this with (3),
It follows from (1) that -(l) 2 L q for (l). By the definition of the parent,
Hence, either - p+1 holds with p replaced by p + 1. Since (3) is
violated for p sufficiently large, we conclude that there exists an integer p such that
-(l), from which it follows that
P(k), the induction step is complete and P(k) ' P(k).
Suppose that l 2 P(k) c . It is now important to recall that k is the index of the
first nonzero component of w, the vector appearing in the update. Observe that l
cannot equal k since l 2 P(k) c and k 2 P(k). The proof that L ]
l is by induction
on the depth d defined by
If then l has no children, and the child loop of Algorithm 2 will be skipped
when either L ]
l or L ]
l are evaluated. And since l 6= k, the pattern associated with w
cannot be added into L ]
l . Hence, when
l is trivial. Now,
assuming that for some p - 0, we have L
l whenever l 2 P(k) c and d(l) - p,
let us suppose that
by the induction assumption. And since l 6= k, the pattern of w is not
added to L ]
l . Consequently, when Algorithm 2 is executed, we have L
l , which
completes the induction step.
Now consider the downdate part of the theorem. Rearranging the downdate
relation AA T
Hence, in a downdate, we can think of A as the updated version of A. Consequently,
the second part of the theorem follows directly from the first part.
4.1 Symbolic update algorithm
We now present an algorithm for evaluating the new pattern L associated with an
update. Based on Theorem 4.1, the only sets L j that change are those associated
with P(k) where k is the index of the first nonzero component of w. Referring to
Algorithm 2, we can set marching up the path
from j, and evaluate all the changes induced by the additional column in A. In order
to do the bookkeeping, there are at most four cases to consider:
Case 1. At the start of the new path, we need to add the pattern for w to L ]
.
Case 2. c 2 P(k), In this case, c is a child
of j in both the new and the old elimination tree. Since the pattern L c may
differ from L c , we need to add the difference to L ]
j . Since j has a unique child
on the path P(k), there is at most one node c that satisfies these conditions.
Also note that if
then by Theorem 4.1, L and hence, this node does not lead to an
adjustment to L ]
in Algorithm 2.
Case 3. In this case, c is a child of j in the
new elimination tree, but not in the old tree, and the entire set L c should be
added to L ]
since it was not included in L ]
. By Theorem 4.1,
all
and hence, c 62 P(k) c , or equivalently, c 2 P(k). Again, since each node on the
path P(k) from k has only one child on the path, there is at most one node c
satisfying these conditions, and it lies on the path P(k).
Case 4. In this case, c is a child of j
in the old elimination tree, but not in the new tree, and the set L c should be
subtracted from L ]
since it was previously added to L ]
each , the fact that implies that c 2 P(k). In the
algorithm that follows, we refer to nodes c that satisfy these conditions as lost
children.
In each of the cases above, every node c that led to adjustments in the pattern was
located on the path P(k). In the detailed algorithm that appears below, we simply
march up the path from k to the root making the adjustments enumerated above.
Algorithm 3 (Symbolic update, add new column w)
Case 1: first node in the path
while j 6= 0 do
Case 2: c is an old child of j, possibly changed
else
Case 3: c is a new child of j and a lost child of -(c)
place c in lost-child-queue of -(c)
Case 4: consider each lost child of j
for each c in lost-child-queue of j do
end for
while
The time taken by this algorithm is given by the following lemma.
Lemma 4.2 The time to execute Algorithm 3 is bounded above by a constant times
the number of entries associated with patterns for nodes on the new path P(k). That
is, the time is
OB @
Proof. In Algorithm 3, we simply march up the path P(k) making adjustments to
j as we proceed. At each node j, we take time proportional to jL j j, plus the time
taken to process the children of j. Each node is visited as a child c at most twice,
since it falls into one or two of the four cases enumerated above (a node c can be a
new child of one node, and a lost child of another). If this work (proportional to jL c j
or jL c j) is accounted to step c instead of j, the time to make the adjustment to the
pattern is bounded above by a constant times either jL j j or jL j j. Since jL
by Theorem 4.1, the proof is complete.
In practice, we can reduce the execution time for Algorithm 3 by skipping over
any nodes for which L ]
. That is, if the current node j has no lost children, and
if its child c falls under case 2 with L
j and the update
can be skipped. The execution time for this modified algorithm, which is the one we
implement, is
4.2 Symbolic downdate algorithm
Let us consider the removal of a column w from A and let k be the index of the first
nonzero entry in w. The symbolic downdate algorithm is analogous to the symbolic
update algorithm, but the roles of P(k) and P(k) are interchanged in accordance
with Theorem 4.1. Instead of adding entries to L ]
j , we subtract entries, instead of
lost child queues, we have new child queues, instead of walking up the path P(k), we
walk up the path P(k) ' P(k).
Algorithm 4 (Symbolic downdate, remove column w)
Case 1: first node in the path
while j 6= 0 do
Case 2: c is an old child of j, possibly changed
else
Case 3: c is a lost child of j and a new child of -(c)
place c in new-child-queue of -(c)
Case 4: consider each new child of j
for each c in new-child-queue of j do
end for
while
Algorithm 4
Similar to Algorithm 3, the execution time obeys the following estimate:
Lemma 4.3 The time to execute Algorithm 4 is bounded above by a constant times
the number of entries associated with patterns for nodes on the old path P(k). That
is, the time is
Again, we achieve in practice some speedup when we check whether or not L ]
changes for any given node j 2 P(k). The time taken by this modified Algorithm 4
is
5 The numerical factors
When we add or delete a column in A, we update or downdate the symbolic
factorization in order to determine the location in the Cholesky factor of either new
nonzero entries or nonzero entries that are now zero. Knowing the location of the
nonzero entries, we can update the numerical value of these entries. We first consider
the case when A and L are dense and draw on the ideas of [15]. Then we show how
the method extends to the sparse case.
5.1 Dense matrices
Our algorithm to implement the numerical update and downdate is based on the
Method C5 in [15] for dense matrices. To summarize their approach, we start by
writing
I
is positive on account
of our assumption that both AA T and AA T are positive definite. That is, since
the matrices AA T
and I are congruent, and by
Sylvester's law of inertia (see [24]), they have the same number of positive, negative,
and zero eigenvalues. The eigenvalues of I are one, with multiplicity
(corresponding to the eigenvectors orthogonal to v) and 1 (corresponding
to the eigenvector v). Since AA T
is positive definite, its eigenvalues are all positive,
which implies that 1
Combining (4) and (5), we have
A sequence of Givens rotations, the product being denoted G 1 , is chosen so that
1 where e 1 is the vector with every entry zero except for the first entry.
The matrix has a lower Hessenberg structure. A second sequence of
Givens rotations, the product being denoted G 2 , is chosen so that (I
is a lower triangular matrix, denoted G, of the following form:
where ffi and fl are vectors computed in Algorithm 5 below. The diagonal elements of
G are given by g while the elements below the diagonal are
Algorithm 5 (Compute G, dense case)
else
find G 1 , to zero out all but first entry of v:
to 1 do
end for
find G 2 , and combine to obtain G:
end for
Algorithm 5
In Algorithm 6 below, we evaluate the new Cholesky factor
forming G explicitly [15] by taking advantage of the special structure of G. The
product is computed column by column, moving from right to left. In practice, L can
be overwritten with L.
Algorithm 6 (Compute
to 1 do
to m do
end for
end for
Algorithm 6
Observe that about 1/3 of the multiplies can be eliminated if L is stored as a
product ~
LD, where D is a diagonal matrix. The components of ffi can be absorbed
into D, with fl adjusted accordingly, and ffi in Algorithm 6 can be eliminated. We did
not exploit this simplification for the numerical results reported in x8.
5.2 The sparsity pattern of v
In the sparse case, sparse, and its nonzero pattern is crucial. Since the
elements of G satisfy for each i ? j, we conclude that only the columns
of L associated with the nonzero components of v enter into the computation of L.
The nonzero pattern of v can be found using the following lemma.
Lemma 5.1 The nodes reachable from any given node k by path(s) in the directed
graph G(L T ) coincide with the path P(k).
Proof. If P(k) has a single element, the lemma holds. Proceeding by induction,
suppose that the lemma holds for all k for which jP(k)j - j. Now, if P(k) has
elements, then by the induction hypothesis, the nodes reachable from -(k) by path(s)
in the directed graph G(L T ) coincide with the path P(-(k)). The nodes reachable
in one step from k consist of the elements of L k . By Proposition 3.1, each of the
elements of L k is contained in the path P(k). If
By the induction hypothesis, the nodes reachable from i coincide with P(i) ' P(k).
The nodes reachable from k consist of fkg union the nodes reachable from L k . Since
that the nodes reachable from k are contained in P(k). On
the other hand, for each p, the element of L T in row - p (k) and column - p+1 (k) is
nonzero. Hence, all the elements of P(k) are reachable from k. Since the nodes in
coincide with the the nodes reachable from k by path(s) in the directed graph
G(L T ), the induction step is complete.
Theorem 5.2 During symbolic downdate AA T
(where w is a column
of A), the nonzero pattern of equal to the path P(k) in the (old)
elimination tree of L where
Proof. Let 0g. Theorem 5.1 of Gilbert [10, 11, 14] states that
the nonzero pattern of v is the set of nodes reachable from the nodes in W by paths
in the directed graph G(L T ). By Algorithm 1, W ' L k . Hence, each element of W
is reachable from k by a path of length one, and the nodes reachable from W are a
subset of the nodes reachable from k. Conversely, since k 2 W, the nodes reachable
from k are a subset of the nodes reachable from W. Combining these inclusions, the
nodes reachable from k and from W are the same, and by Lemma 5.1, the nodes
reachable from k coincide with the path P(k).
Corollary 5.3 During symbolic update AA T
the nonzero pattern
of equal to the path P(k) in the (new) elimination tree of L where k is
defined in (6).
Proof. Since , we can view L as the Cholesky factor for the
downdate Hence, we can apply Theorem 5.2, in effect replacing P by
P.
5.3 Sparse matrices
In the dense algorithm presented at the start of this section, we write
To be specific, let us consider the case where (7) corresponds to
an update. The nonzero elements of I lie along the diagonal or at
the intersection of rows i 2 P(k) and columns j 2 P(k), where k is defined in (6).
In essence, we can extract the submatrix of I corresponding to these rows
and columns, we can apply Algorithm 5 to this dense submatrix, we can modify the
submatrix of L consisting of rows and columns associated with P(k), and then we
can reinsert this dense submatrix in the m by m sparse matrix L. At the algebraic
level, we can think of this process in the following way: If P is the permutation matrix
with the property that the columns associated with P(k) are moved to the first jP(k)j
columns by multiplication on the right, then we have
The matrix P T diagonal with the identity in the lower
right corner and a dense upper left corner V,
!/
I
The elements of L 21 correspond to the elements l ij of L for which i 2 P(k) c and
follows that i 2 L c
Note that in general
which is why Algorithm 6 is column-oriented. When (8) is multiplied by the
Givens rotations computed by Algorithm 5, we obtain
!/
I
We apply P to the left and P T to the right in (9) to obtain L. Neither L 12 nor L 22
are modified.
When (7) corresponds to a downdate, the discussion is the same as for an
update except that P(k) replaces P(k), and the role of L j and L j are interchanged.
In summary, the sparse update or downdate of a Cholesky factorization can be
accomplished by first evaluating nonzero entries are contained
in the path P(k) or P(k) where k is defined in (6). We apply the dense Algorithm 5
to the vector of (symbolically) nonzero entries of v. Then we update the entries in L
in the rows and columns associated with indices in P(k) or P(k) using Algorithm 6.
As the following result indicates, Algorithm 6 does not have to be applied to all
of L 11 because of its specific structure.
Proposition 5.4 During symbolic downdate, L 11 is a lower triangular matrix with
a dense profile. That is, for each row i ? 1, there is an integer
Proof. The rows and columns columns of L 11 correspond to the nodes in the path
P(k). We assume that the nodes have been relabeled so that
It follows that [L 11
by (1), if Consequently, [L
Letting p i be the smallest such index p, the proof is complete.
Corollary 5.5 During symbolic update L 11 is a lower triangular matrix with a dense
profile. That is, for each row i ? 1, there is an integer
Proof. This follows immediately from Proposition 5.4, replacing P(k) with P(k).
As a consequence of Proposition 5.4, Corollary 5.5, and Proposition 3.2, we can
skip over any index i in Algorithm 6 for which Algorithm 6 is applied
to the sparse L, as opposed to the dense submatrix, the indices j and i take values
in P(k) and L j respectively for a downdate, while they take values in P(k) and L j ,
respectively, for an update.
6 Arbitrary symbolic and numerical factors
The methods we have developed for computing the modification to the Cholesky
factors of AA T corresponding to the addition or deletion of columns in A can be
used to determine the effect on the Cholesky factors of a general symmetric positive
definite matrix M of any symmetric change of the form M preserves
positive definiteness. We briefly describe how Algorithms 1 through 6 are modified
for the general case.
Let M j denote the nonzero pattern of the lower triangular part of M:
The symbolic factorization of M [5, 6, 7, 8, 23] is obtained by replacing the union of
A k terms in Algorithm 1 with the set M j . With this change, L j of Algorithm 1 is
given by
This leads to a change in Algorithm 2 for computing the multiplicities. The
multiplicity of an index i in L j becomes
The for loop involving the A k terms in Algorithm 2 is replaced by the single statement
More precisely, we have:
for each k where min A do
end for
Entries are removed or added symbolically from AA T by the deletion or addition
of columns of A, and numerical cancellation is ignored. Numerical cancellation of
entries in M should not be ignored, however, because this is the only way that entries
can be dropped from M. When numerical cancellation is taken into account, neither
of the inclusions M may hold. We resolve this problem by using
a symbolic modification scheme with two steps: a symbolic update phase in which
new nonzero entries in M are taken into account, followed by a separate
downdate phase to handle entries that become numerically zero. Since each
modification step now involves an update phase followed by a downdate phase, we
attach (in this section) a overbar to quantities associated with the update and an
underbar to quantities associated with the downdate.
Let W be the nonzero pattern of w, namely 0g. In the first
phase, entries from W are symbolically added to M j for each j 2 W. That
In the second symbolic phase, entries from W are symbolically deleted for each j 2 W:
In practice, we need to introduce a drop tolerance t and replace the equality
in (10) by the inequality t. For a general matrix, the analogue of
Theorem 4.1 is the following:
Theorem 6.1 If ff is the first index for which M ff 6= M ff , then P(ff) ' P(ff)
and L ]
i for all i 2 P(ff) c . If fi is the first index for which M fi 6= M fi , then
In evaluating the modification in the symbolic factorization associated with
we start at the first index ff where M ff 6= M ff and we march up
the path P(ff) making changes to L ]
j . In the second phase, we start at
the first index where M fi 6= M fi , and we march up the path P(fi) making changes to
. The analogue of Algorithm 3 in the general case only differs in the
starting index (now ff) and in the addition of the sets M j n M j in each pass through
the j-loop:
Algorithm 7a (Symbolic update phase, general matrix)
Case 1: first node in the path
while j 6= 0 do
Case 2: c is an old child of j, possibly changed
else
Case 3: c is a new child of j and a lost child of -(c)
place c in lost-child-queue of -(c)
Case 4: consider each lost child of j
for each c in lost-child-queue of j do
end for
while
Algorithm 7a
Similarly, the analogue of Algorithm 4 in the general case only differs in the
starting index (now fi), in the subtraction of the sets and M j n M j in each pass
through the j-loop.
Algorithm 7b (Symbolic downdate phase, general matrix)
Case 1: first node in the path
while j 6= 0 do
Case 2: c is an old child of j, possibly changed
else
Case 3: c is a lost child of j and a new child of -(c)
place c in new-child-queue of -(c)
Case 4: consider each new child of j
for each c in new-child-queue of j do
end for
while
Algorithm 7b
Algorithms 5 and 6 are completely unchanged in the general case. They can be
applied after the completion of Algorithm 7b so that we know the location of new
nonzero entries in the Cholesky factor. They process the submatrix associated with
rows and columns in P(k) where k is the index of the first nonzero element of w.
When M has the form AA T and when M is gotten by either adding or deleting a
column in A, then assuming no numerical cancellations, Algorithm 7b can be skipped
when we add a column to A since M for each j. Similarly, when a column
is removed from A, Algorithm 7a can be skipped since M for each j. Hence,
when Algorithm 7a followed by Algorithm 7b are applied to a matrix of the form
AA T , only Algorithm 7a takes effect during a update while only Algorithm 7b takes
effect during a downdate. Thus the approach we have presented in this section for an
arbitrary symmetric positive definite matrix generalizes the earlier approach where
we focus on matrices of the form AA T .
7 Implementation issues
In this section, we discuss implementation issues in the context of updates and
downdates to a matrix of the form AA T . A similar discussion applies to a general
symmetric positive definite matrix. We assume that the columns of the matrix A
in the product AA T are all chosen from among the columns of some fixed matrix
B. The update and downdate algorithms can be used to compute the modification
to a Cholesky factor corresponding to additions or deletions to the columns of A.
The Cholesky factorization of the initial AA T (before any columns are added or
deleted) is often preceded by a fill reducing permutation P of the rows and columns.
For example, we could compute a permutation to reduce the fill for BB T since the
Cholesky factors of AA T will be at least as sparse as those of BB T by Proposition 3.2,
regardless of how the columns of A are chosen from the columns of B. Based on the
number of nonzeros in each column of the Cholesky factors of BB T , we could allocate
a static storage structure that will always contain the Cholesky factors of each AA T .
On the other hand, this could lead to wasted space if the number of nonzeros in the
Cholesky factors of AA T are far less than the number of nonzeros in the Cholesky
factors of BB T . Alternatively, we could store the Cholesky factor of the current AA T
in a smaller space, and reallocate storage during the updates and downdates, based
on the changes in the nonzero patterns.
The time for the initial Cholesky factorization of AA T is given by (see [8]) the
following
which is O(m 3 ) if L is dense. Each update and downdate proceeds by first finding
the new symbolic factors and the required path (P(k) for updating, and P(k) for
downdating), using Algorithm 3 or 4 (modified to skip in constant time any column
j that does not change). Algorithms 5 and 6 are then applied to the columns in the
path. The lower triangular solve of the system Algorithm 5 can be done
in time
during downdating, and
OB @
during updating [14]. The remainder of Algorithm 5 takes time O(jP(k)j) during
downdating, and O(jP(k)j) during updating. Our sparse form of Algorithm 6
discussed in x5.3, which computes the product
during downdating, and
OB @
during updating. Since L j ' L j during an update, it follows that the asymptotic
time for the entire downdate or update, both symbolic and numeric, is equal to the
asymptotic time of Algorithm 6. This can be much less than the O(m 2
by Algorithms 5 and 6 in the dense case.
8 Experimental results
We have developed Matlab codes to experiment with all the algorithms presented in
this paper, including the algorithms of x6 for a general symmetric, positive definite
matrix. In this section, we present the results of a numerical experiment with a large
sparse optimization problem from Netlib [3]. The computer used for this experiment
was a Model 170 UltraSparc, equipped with 256MB of memory, and with Matlab
Version 4.2c.
8.1 Experimental design
We selected an optimization problem from airline scheduling (DFL001). Its constraint
matrix B is 6071-by-12,230 with 35,632 nonzeros. The matrix BB T has 37,923
nonzeros in its strictly lower triangular part. Its Cholesky factor LB has 1.49 million
nonzeros (with a fill-minimizing permutation PB of the rows of B, described below)
and requires 1.12 billion floating-point operations and 115 seconds to compute (the
LP Dual Active Set Algorithm does not require this matrix, however, as noted earlier,
this is an upper bound on the number of nonzeros that can occur during the execution
of the LP DASA). This high level of fill-in in LB is the result of the highly irregular
nonzero pattern of B. The basis matrix A 0 corresponding to an optimal solution of
the linear programming problem has 5,446 columns.
We wrote a set of Matlab scripts that implements our complete Cholesky
update/downdate algorithm, discussed in x7. We first found PB , using 101 trials
of Matlab's column multiple minimum degree ordering algorithm (colmmd [12]), 100
of them with a different random permutation of the rows of B. We then took the
best permutation found. With this permutation (PB ), the factor L of A 0 A T
0 has 831
thousand nonzeros, and took 481 million floating-point operations and 51 seconds to
compute (using Matlab's chol). Following the method used in LP DASA, we added
\Gamma12 to the diagonal to ensure positive definiteness. We used the same permutation
PB for the entire experiment. The initial symbolic factorization took 15 seconds
(Algorithm 2). It is this matrix and its factor that are required by the LP DASA.
We did not use Matlab's sparse matrix data structure, since Matlab removes
explicit zeros. Changing the nonzero pattern by a single entry can cause Matlab to
make a new copy of the entire matrix. This would defeat the asymptotic performance
of our algorithms. Instead, the column-oriented data structure we use for L, L ] , and
L consists of three arrays of length jL B j, an array of length m that contains indices
to the first entry in each column, and an array of length m holding the number of
nonzeros in each column. The columns are allocated so that each column can hold as
many nonzeros as the corresponding column of LB , without reallocation.
Starting with the optimal basis of 5,446 columns, we added one column at a time
until the basis included all 12,230 columns, and then removed them one at a time
to obtain the optimal basis again. The total time and work required to modify the
factors was 76,598 seconds and 61.5 billion floating point operations. This divides
into
(a) 441 seconds for bookkeeping to keep track of the basis set,
(b) 6,051 seconds for the symbolic updates and downdates (Algorithms 3 and 4),
(c) 21,167 seconds to solve
(d) 4,977 seconds for the remainder of Algorithm 5, and
43,962 seconds for Algorithm 6.
Algorithm 6 clearly dominates the modification algorithm. The symbolic updates and
downdates took very little time, even though Algorithms 3 and 4 are not well suited
to implementation in Matlab. By comparison, using the Cholesky factorization to
solve a linear system LL T with a dense right-hand-side b (using our column-
oriented data structure for L) at each step took a total of 93,418 seconds and 67.7
billion floating point operations (each solve takes O(jLj) time). Note that this is
much higher than the time taken to solve Algorithm 5, because v and w
are sparse. The time taken for the entire update/downdate computation would be
much smaller if our code was written in a compiled language. Solving one system
with a dense right hand side (using the factorization of the optimal basis)
takes 5.53 seconds using our column-oriented data structure, 1.26 seconds using a
Matlab sparse matrix for L, and 0.215 seconds in Fortran 77.
8.2 Numerical accuracy
In order to measure the error in the computed Cholesky factorization, we evaluated
the difference kAA
is the computed Cholesky factorization.
For the airline scheduling matrix of x8, - L has up to 1.49 million nonzeros and it is
impractical to compute the product -
after each update. To obtain a quick and
accurate estimate for kEk 1 , where
, we applied the strategy presented
in [16] (see [17, p. 139] for a symbolic statement of the algorithm) to estimate the
1-norm of the inverse of a matrix. That is, we used a gradient accent approach to
compute a local maximum for the following problem:
L is used multiple times in the following algorithm, we copied our data structure
for -
L into a Matlab sparse matrix.
Algorithm 8 (Estimate 1-norm of an m by m matrix E)
(ae is the current estimate for kEk)
while do
do
end for
while
Algorithm 8
To improve the accuracy of the 1-norm estimate, we used Algorithm 8 three times.
In the second and third trials, a different starting vector x was used as described in
[16]. Observe that Algorithm 8 only makes use of the product between the matrix
E and a vector. This feature is important in the context of sparse matrices since E
contains the term -
and it is impractical to compute the product -
, but it is
practical to multiply -
by a vector. For the airline scheduling matrix of x8 the
values for kEk 1 initially, at step 6,784, and at the end were 2:49
and 1:54\Theta10 \Gamma10 respectively. The estimates obtained using Algorithm 8 were identical
at the same three steps. On the other hand, the times to compute kEk 1 at the initial
step and at step 6,784 were 137.9 and 279.5 seconds, while the times for three trials of
Algorithm 8 were 8.7 and 14.7 seconds, respectively (excluding the time to construct
the Matlab sparse matrix for -
L). Our methods were quite accurate for this problem.
After 6,784 updates and 6,784 downdates, or 13,568 changes in A, the 1-norm of E
increased by only a factor 618! Figure 1 shows the estimated value of kEk 1 computed
every steps using Algorithm 8. The 1-norm of the matrix AA T increases from
458.0 initially to 1107.0 at iteration 6,784, then returns to 458.0 at iteration 13,568.
Hence, the product of the computed Cholesky factors agrees with the product AA T
to about 15 significant digits initially, while the products agree to about 12 significant
digits after 13,568 modifications of A.
step
norm,
as
estimated
by
Algorithm
Figure

1: Estimated 1-norm of error in Cholesky factorization
8.3 Alternative permutations
Our methods are optimal in the sense that they take time proportional to the number
of nonzero entries in L that change at each step. However, they are not optimal with
respect to fill-in, since we assume a single initial permutation, and no subsequent
permutations. A fill-reducing ordering of BB T might not be the best ordering to use
for all basis matrices. A simple pathological example is the m-by-n matrix B, where
and the nonzero pattern of each of the column of B is a unique
pair of integers from the set f1; 2; mg. In this case, every element of BB T is
nonzero, while the nonzero pattern of AA T is arbitrary. As the matrix A changes, it
might be advantageous to compute a fill-reducing ordering of AA T if the size of its
factors grow "too large." A refactorization with the new permutation would then be
required.
We found a fill-reducing permutation PA of the optimal basis A 0 A T
(again, the
best of 101 trials of colmmd). This results in a factor L with 381 thousand nonzeros,
requiring only 169 million floating point operations to compute. This is significantly
less than the number of nonzeros (831 thousand) and floating point operations (481
million) associated with the fill reducing permutation for BB T . We also computed
an ordering of AA T at each step, using colmmd just once, and then computed the
number of nonzeros in the factor if we were to factorize A using this permutation
Although it only takes about 1 second to compute the ordering [12] and symbolic
step
in
Nonzeros in L for three different permutations

Figure

2: Nonzeros in L using three different permutations
factorization [13], it is not practical to use 101 random trials at each step.

Figure

2 depicts the nonzero counts of L for these three different permutations,
at each of the 13,568 steps. The fixed permutation PB results in the smooth curve
starting at 831 thousand and peaking at 1.49 million. The fixed permutation PA
results in a number of nonzeros in L that starts at 381 thousand and rises quickly,
leaving the figure at step 1,206 and peaking at 7.4 million in the middle. It surpasses
PB at step 267. Using a permutation, P s , computed at each step s, gives the
erratic line in the figure, starting at 390 thousand and peaking at 1.9 million in the
middle. These results indicate that it might be advantageous to start with the fixed
permutation PA , use it for 267 steps, and then refactorize with the permutation P s
computed at step 267. This results in a new factor with only 463 thousand nonzeros.
Near the center of the figure, however, the basis A includes most of the columns in
B, and in this case the PB permutation should be used.
9

Summary

We have presented a new method for updating and downdating the factorization LL T
of a sparse symmetric positive definite matrix AA T . Our experimental results show
that the method should be fast and accurate in practice. Extensions to an arbitrary
sparse symmetric positive definite matrix, M, have been discussed. We mention here
an additional extension to our work that would be useful.
We do not make use of the supernodal form of the factorization, nor do we use
the related compressed pattern of L [8]. Any method can be used for the numerical
factorization of the first basis matrix, of course, but the factor would then be copied
into a simple column-oriented data structure. Keeping the supernodal form has the
potential of reducing time taken by the symbolic factorization (Algorithm 2), the
symbolic update and downdate (Algorithms 3 and 4), and the numerical update and
downdate (dense matrix kernels could be used in Algorithms 5 and 6). However,
supernodes would merge during update, and split during downdate, complicating the
supernodal form of the factorization.



--R

An approximate minimum degree ordering algorithm
Introduction to Algorithms
Distribution of mathematical software via electronic mail
Direct Methods for Sparse Matrices
Yale sparse matrix package
The design of a user interface for a sparse matrix package


A data structure for sparse QR and LU factorizations
Predicting structure in sparse matrix computations

Sparse matrices in MATLAB: design and implementation
An efficient algorithm to compute row and column counts for sparse Cholesky factorization
Sparse partial pivoting in time proportional to arithmetic operations
Methods for modifying matrix factorizations



Active set strategies in the LP dual active set algorithm
The role of elimination trees in sparse factorization
Algorithmic aspects of vertex elimination on graphs
A new implementation of sparse Gaussian elimination
On the efficient solution of sparse systems of linear and nonlinear equations
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
Frank Dellaert , Michael Kaess, Square Root SAM: Simultaneous Localization and Mapping via Square Root                 Information Smoothing, International Journal of Robotics Research, v.25 n.12, p.1181-1203, December  2006
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
Olga Sorkine , Daniel Cohen-Or , Dror Irony , Sivan Toledo, Geometry-Aware Bases for Shape Approximation, IEEE Transactions on Visualization and Computer Graphics, v.11 n.2, p.171-180, March 2005
numerical linear algebra;matrix updates;cholesky factorization;sparse matrices;mathematical software;direct methods