2.1 Sparse matrix data structures

The simplest sparse matrix data structure is a list of the non-zero entries in arbitrary order, also called the triplet form. This is easy to generate but hard to use in most sparse direct methods, so the format is often used in an interface to a package but not in its internal representation. This data structure can be easily converted into compressed-column form in linear time via a bucket sort. In this format, each column is represented as a list of values and their corresponding row indices. To create this structure, the first pass counts the number of entries in each column of the matrix, and the column pointer array is constructed as the cumulative sum of the column counts. The entries are placed in their appropriate columns in a second pass. In the compressed-column form, an m $\times$ n sparse matrix that can contain up to nzmax entries is represented with an integer array p of length n+1, an integer array i of length nzmax, and a real array a of length nzmax. Row indices of entries in column j are stored in i[p[j]] to i[p[j+1]-1], and the numerical values are stored in the same locations in a. In zero-based form (where rows and columns start at zero) the first entry p[0] is always zero, and p[n] is the number of entries in the matrix. An example of a matrix and its zero-based compressed-column form is given by

(2.1)$$\begin{eqnarray}A=\left[\begin{array}{@{}cccc@{}}4.5 & 0 & 3.2 & 0\\ 3.1 & 2.9 & 0 & 0.9\\ 0 & 1.7 & 3.0 & 0\\ 3.5 & 0.4 & 0 & 1.0\\ \end{array}\right]\!.\end{eqnarray}$$

Exact numerical cancellation is rare, and most algorithms ignore it. An entry in the data structure that is computed but found to be numerically zero is still called a ‘non-zero’, by convention. Leaving these entries in the matrix leads to much simpler algorithms and more elegant graph-theoretical statements about the algorithms.

Accessing a column of this data structure is very fast, but extracting a given row is very costly. Algorithms that operate on such a data structure must be designed accordingly. Likewise, modifying the non-zero pattern of a compressed-column matrix is not trivial. Deleting or adding single entry takes $O(|A|)$ time, if gaps are not tolerated between columns. MATLAB uses the compressed-column data structure for its sparse matrices, with the extra constraint that no numerically zero entries are stored Gilbert, Moler and Schreiber (Reference Gilbert, Moler and Schreiber1992). McNamee (Reference Mcnamee1971, Reference Mcnamee1983a, Reference Mcnamee1983b) provided the first published software for basic sparse matrix data structures and operations (multiply, add, and transpose). Gustavson (Reference Gustavson, Rose and Willoughby1972) summarizes a range of methods and data structures, including how this data structure can be modified dynamically to handle fill-in during factorization.

Finite element methods generate a matrix as a collection of elements, or dense submatrices. A solver dedicated to solving finite element problems will often accept its input in this form. Each element requires a non-zero pattern of the rows and columns it affects. The complete matrix is a summation of the elements, and two or more elements may contribute to the same matrix entry. Thus, it is a common practice for any data structure that any duplicate entries be summed, which can be done in linear time.

2.2 Matrix–vector multiplication

One of the simplest sparse matrix algorithms is matrix–vector multiplication, $z=Ax+y$, where $y$ and $x$ are dense vectors and $A$ is sparse. If $A$ is split into $n$ column vectors, the result $z=Ax+y$ is

$$\begin{eqnarray}z=\left[\begin{array}{@{}ccc@{}}A_{\ast 1} & \cdots \, & A_{\ast n}\end{array}\right]\left[\begin{array}{@{}c@{}}x_{1}\\ \vdots \\ x_{n}\end{array}\right]+y.\end{eqnarray}$$

Allowing the result to overwrite the input vector $y$, the $j$th iteration computes $y=y+A_{\ast j}x_{j}$. The 1-based pseudocode for computing $y=Ax+y$ is given below.

Algorithm 2.1 illustrates the use of gather/scatter operations, which are very common in sparse direct methods. The innermost loop has the form y[i[k]] = y[i[k]] + a[k]*x[j], which requires both a gather and a scatter. A gather operation is a subscripted subscript appearing on the right-hand side of an assignment statement, of the form s = y[i[k]], where y is a dense vector and s is a scalar. A scatter operation occurs when an expression y[i[k]] appears as the target of an assignment. Gather/scatter operations result in a very irregular access to memory, and are thus very costly. Memory transfers are fastest when access has spatial or temporal locality, although the computations themselves can be pipelined on modern architectures (Lewis and Simon Reference Lewis and Simon1988). One of the goals of supernodal, frontal, and multifrontal methods (Sections 9–11) is to replace most of the irregular gather/scatter operations with regular operations on dense submatrices.

Sparse matrix–vector multiplication is a common kernel in iterative methods, although they typically use a compressed-row format (or one of many other specialized data structures), since it results in a faster method because of memory traffic. It requires only a gather operation, not the gather/scatter of Algorithm 2.1.

2.3 Transpose

The algorithm for transposing a sparse matrix ($C=A^{T}$) is very similar to the method for converting a triplet form to a compressed-column form. The algorithm computes the row counts of A and the cumulative sum to obtain the row pointers. It then iterates over each non-zero entry in A, placing the entry in its appropriate row vector. Transposing a matrix twice results in sorted row indices. Gustavson (Reference Gustavson1978) describes algorithms for matrix multiply ($C=AB$) and permuted transpose ($C=(PAQ)^{T}$ where $P$ and $Q$ are permutation matrices). The latter takes linear time ($O(n+|A|)$ where $n$ is the matrix dimension and $|A|$ is the number of non-zeros).

2.4 Matrix multiplication and addition

Algorithms are tuned to their data structures. For example, if two sparse matrices $A$ and $B$ are stored in compressed-column form, then a matrix multiplication $C=AB$ where $C$ is $m\times n$, $A$ is $m\times k$, and $B$ is $k\times n$, should access $A$ and $B$ by column, and create $C$ one column at a time. If $C_{\ast j}$ and $B_{\ast j}$ denote column $j$ of $C$ and $B$, then $C_{\ast j}=AB_{\ast j}$. Splitting $A$ into its $k$ columns and $B_{\ast j}$ into its $k$ individual entries results in

(2.2)$$\begin{eqnarray}C_{\ast j}=\left[\begin{array}{@{}ccc@{}}A_{\ast 1} & \cdots \, & A_{\ast k}\end{array}\right]\left[\begin{array}{@{}c@{}}b_{1j}\\ \vdots \\ b_{kj}\end{array}\right]=\mathop{\sum }_{t=1}^{k}A_{\ast t}b_{tj}.\end{eqnarray}$$

The non-zero pattern of $C$ is given by the following theorem.

Theorem 2.1 (Gilbert Reference Gilbert1994).

The non-zero pattern of $C_{\ast j}$ is the set union of the non-zero pattern of $A_{\ast t}$ for all $t$ for which $b_{tj}$ is non-zero. If ${\mathcal{C}}_{j}$, ${\mathcal{A}}_{i}$, and ${\mathcal{B}}_{j}$ denote the set of row indices of non-zero entries in $C_{\ast j}$, $A_{\ast t}$, and $B_{\ast j}$, then (ignoring numerical cancellation)

(2.3)$$\begin{eqnarray}{\mathcal{C}}_{j}=\mathop{\bigcup }_{t\in {\mathcal{B}}_{j}}{\mathcal{A}}_{t}.\end{eqnarray}$$

A matrix multiplication algorithm must compute both the numerical values $C_{\ast j}$ and the pattern ${\mathcal{C}}_{j}$, for each column $j$. A variant of Gustavson’s algorithm below uses a temporary dense vector $x$ to construct each column, and a flag vector $w$, both of size $n$ and both initially zero (Gustavson Reference Gustavson1978). The matrix $C$ grows one entry at a time, column by column. The three ‘for each …’ loops access a single column of a sparse matrix in a compressed-column data structure.

As stated, this algorithm requires a dynamic allocation of $C$ since it starts out as empty and grows one entry at a time to a size that is not known a priori. Computing the pattern ${\mathcal{C}}$ or even just its size is actually much harder than computing the pattern or size of a sparse factorization. The latter is discussed in Section 4. Another approach computes the size of the pattern ${\mathcal{C}}$ in an initial symbolic pass, followed by the second numeric pass above, which allows for $C$ to be statically allocated. The first pass and second passes take the same time, in a big-${\mathcal{O}}$ sense. For sparse factorization, by contrast, the symbolic analysis is asymptotically faster than the numeric factorization. The time taken by Algorithm 2.2 is $O(n+f+|B|)$, where $f$ is the number of floating-point operations performed, which is typically dominated by $f$.

Matrix addition $C=A+B$ is very similar to matrix multiplication. The only difference is that two nested ‘for each’ loops above are replaced by two non-nested loops, one that accesses the $j$th column of $B$ and the other the $j$th column of $A$.

2.5 Permutations

A sparse matrix must typically be permuted either before or during its numeric factorization, either for reducing fill-in or for numerical stability, or often for both. Fill-in is the introduction of new non-zeros in the factors that do not appear in the corresponding positions in the matrix being factorized. Finding a suitable fill-reducing permutation is discussed in Section 8.

Permutations are typically represented as integer vectors. If $C=PA$ is to be computed, then the row permutation $P$ can be represented as a permutation vector p of size n. If row $i$ of $A$ becomes the $k$th row of $C$, then i=p[k]. The inverse permutation is k=invp[i]. Permuting $A$ to obtain $C$ requires the latter, since a traversal of $A$ gives a set of row indices $i$ of $A$ that must be translated to rows $k$ of $C$, via the computation k=invp[i]. Column permutations $C=AQ$ are most simply done with a permutation vector q where j=q[k] if column $j$ of $A$ corresponds to column $k$ of $C$. This allows $C$ to be constructed one column at a time, from left to right.

3.1 A dense right-hand side

There are many ways of solving $Lx=b$ but if $L$ is stored as a compressed-column sparse matrix, accessing $L$ by columns is the most natural. Assuming $L$ has unit diagonal, we obtain a simple Algorithm 3.1, below, that is very similar to matrix–vector multiplication. The vector $x$ is accessed via gather/scatter. Solving related systems, such as $L^{T}x=b$, $Ux=b$, and $U^{T}x=b$ (where $U$ is upper triangular) is similar, except that the transposed solves are best done by accessing the transpose of the matrix row-by-row, if the matrix is stored in column form.

3.2 A sparse right-hand side

When the right-hand side is sparse, the solution $x$ is also sparse, and not all columns of $L$ take part in the computation since the $j$th loop of Algorithm 3.1 can be skipped if $x_{j}=0$. Scanning all of $x$ for this condition adds $O(n)$ to the time complexity, which is not optimal. However, if the set of indices $j$ for which $x_{j}$ will be non-zero is known, denoted ${\mathcal{X}}=\{j\,|\,x_{j}\neq 0\}$, the algorithm is as follows.

Figure 3.1.

Sparse triangular solve (Davis Reference Davis2006).

The runtime is now optimal, but ${\mathcal{X}}$ must be determined first, as Gilbert and Peierls (Reference Gilbert and Peierls1988) describe. Finding ${\mathcal{X}}$ requires that we follow a set of implications. Entries in $x$ become non-zero in two places, the first and last lines of Algorithm 3.2. Neglecting numerical cancellation, rule one states that $b_{i}\neq 0\Rightarrow x_{i}\neq 0$, and rule two is $x_{j}\neq 0\wedge \exists i(l_{ij}\neq 0)\Rightarrow x_{i}\neq 0$. These two rules lead to a graph traversal problem. Consider a directed acyclic graph $G_{L}=(V,E)$ where $V=\{1\ldots n\}$ and $E=\{(j,i)\,|\,l_{ij}\neq 0\}$. Rule one marks all nodes in ${\mathcal{B}}$, the non-zero pattern of $b$. Rule two states that if a node is marked, all its neighbours become marked. These rules are illustrated in Figure 3.1. In graph terminology,

(3.1)$$\begin{eqnarray}{\mathcal{X}}=\text{Reach}_{L}({\mathcal{B}}).\end{eqnarray}$$

Computing ${\mathcal{X}}$ requires a depth-first search of the directed graph $G_{L}$, starting at nodes in ${\mathcal{B}}$. The time taken is proportional to the number of edges traversed, which is exactly equal to the floating-point operation count. A depth-first search computes ${\mathcal{X}}$ in topological order, and performing the numerical solve in that order preserves the numerical dependencies. An example of a triangular solve is shown in Figure 3.2. If ${\mathcal{B}}=\{4,6\}$, then in this example ${\mathcal{X}}=\{6,10,11,4,9,12,13,14\}$, listed in topological order as produced by the depth-first search. The sparse triangular solve forms the basis of the left-looking sparse LU method presented in Section 6.2. For the parallel triangular solve, see Section 12.2.

Figure 3.2.

Solving $Lx=b$, where $L$, $x$, and $b$ are sparse (Davis Reference Davis2006).

4.1 Elimination tree

The elimination tree (or etree) appears in many sparse factorization algorithms and in many related theorems. It guides the symbolic analysis and provides a framework for numeric factorization, both parallel and sequential.

Consider equation (4.1), applied to the leading $k\times k$ submatrix of $A$. The vector $l_{12}$ is computed with a sparse triangular solve, $L_{11}l_{12}=a_{12}$, and its transpose becomes the $k$th row of $L$. Its non-zero pattern is thus ${\mathcal{L}}_{k}=\text{Reach}_{G_{k-1}}({\mathcal{A}}_{k})$, where $G_{k-1}$ is the directed graph of $L_{11}$, ${\mathcal{L}}_{k}$ is the non-zero pattern of the $k$th row of $L$, and ${\mathcal{A}}_{k}$ is the non-zero pattern of the upper triangular part of the $k$th column of $A$.

The depth-first search of $G_{k-1}$ is sufficient for computing ${\mathcal{L}}_{k}$, but a simpler method exists, taking only $O(|{\mathcal{L}}_{k}|)$ time. It is based on a pruning of the graph $G_{k-1}$ (Figure 4.1). Computing the $k$th row of $L$ requires a triangular solve of $L_{1:k-1,1:1-k}x=b$, where the right-hand $b$ is the $k$ row of $A$. Thus, if $a_{ik}$ is non-zero, so is $b_{i}$ and $x_{i}$. This entry $x_{i}$ becomes $l_{ki}$ and thus $a_{ik}\neq 0$ implies $x_{i}=l_{ki}\neq 0$.

Figure 4.1.

Pruning the directed graph $G_{L}$ yields the elimination tree.

If there is another non-zero $l_{ji}$ with $j<k$, then there is an edge from $i$ to $j$ in the graph, and so when doing the graph traversal from node $i$, node $j$ will be reached. As a result, $x_{j}\neq 0$. This becomes the entry $l_{kj}$, and thus the existence of a pair of non-zeros $l_{ji}$ and $l_{ki}$ implies that $l_{kj}$ is non-zero (Parter Reference Parter1961). Outer-product Gaussian elimination is another way to look at this ($A=LU$): for a symmetric matrix $A$, two non-zero entries $u_{ij}$ and $l_{ki}$ cause $l_{kj}$ to become non-zero when the $i$th pivot entry is eliminated.

For the sparse triangular solve, if the graph traversal starts at node $i$, it will see both nodes $j$ and $k$. The reach of node $i$ is not changed if the edge $(i,k)$ is ignored; it can still reach $k$ from a path of length two: $i$ to $j$ to $k$. As a result, any edge $(i,k)$ corresponding to the non-zero $l_{ki}$ can be ignored in the graph search, just so long as there is another non-zero $l_{ji}$ above it in the same column. Only the first off-diagonal entry is needed, that is, the smallest $j>i$ for which $l_{ji}$ is non-zero. Pruning all but this edge results in the same reach (3.1) for the triangular solve, $L_{11}l_{12}=a_{12}$, but the resulting structure is faster to traverse. A directed acyclic graph with at most one outgoing edge per node is a tree (or forest). This is the elimination tree. It may actually be a forest, but by convention is still called the elimination tree.

In terms of the graph of $A$, if there is a path of length two between the two nodes $k$ and $j$, where the intermediate node is $i<\min (k,j)$, then this causes the edge $(k,j)$ to appear in the filled graph. Rose, Tarjan and Lueker (Reference Rose, Tarjan and Lueker1976) generalized Parter’s result in their path lemma: $l_{ij}$ is non-zero if and only if there is a path path $i{\rightsquigarrow}j$ in $A$ where all intermediate nodes are numbered less than $\min (i,j)$.

The elimination tree is based on the pattern of $L$, but it can be computed much more efficiently without finding the pattern of $L$, in time essentially linear in $|A|$ (Liu Reference Liu1986a, Schreiber Reference Schreiber1982). The algorithm relies on the quick traversal of paths in the tree as it is being constructed. A key theorem by Liu states that if $a_{ki}\neq 0$ (where $k>i$), then $i$ is a descendant of $k$ in the elimination tree (Liu Reference Liu1990). The tree is constructed by ensuring this property holds for each entry in $A$.

Let ${\mathcal{T}}$ denote the elimination tree of $L$, and let ${\mathcal{T}}_{k}$ denote the elimination tree of submatrix $L_{1\ldots k,1\ldots k}$, the first $k$ rows and columns of $L$.

The tree is constructed incrementally, for each leading submatrix of $A$. That is, the tree ${\mathcal{T}}_{k}$ for the matrix $A_{1..k,1..k}$ is constructed from ${\mathcal{T}}_{k-1}$. For each entry $a_{ki}\neq 0$, it suffices to ensure that $i$ is a descendant of $k$ in the next tree ${\mathcal{T}}_{k}$. Walking the path from $i$ to a root $t$ in ${\mathcal{T}}_{k-1}$ means that $t$ must be a child of $k$ in ${\mathcal{T}}_{k}$. This path could be traversed one node at a time, but it is much more efficient to compress the paths as they are found, to speed up subsequent traversals. Future traversals from $i$ must arrive at $k$, so as the path is traversed, each node along the path is given a shortcut to its ancestor $k$. This method is an application of the disjoint-set-union-find algorithm of Tarjan (Reference Tarjan1975) (see also Cormen et al. Reference Cormen, Leiserson and Rivest1990). Thus, its runtime is essentially linear in the number of entries in $A$ (this is called nearly $O(|A|)$). This is much faster than $O(|L|)$, and it means that the tree can be computed prior to computing the non-zero pattern of $L$, and then used to construct the pattern of $L$ or to explore its properties.

The column elimination tree is the elimination tree of $A^{T}A$, and is used in the QR and LU factorization algorithms. It can be computed without forming $A^{T}A$, also in nearly $O(|A|)$ time. Row $i$ of $A$ creates a dense submatrix, or clique, in the graph of $A^{T}A$. To speed up the construction of the tree without forming the graph of $A^{T}A$ explicitly, a sparser graph is used, in which each clique is replaced with a path amongst the nodes in the clique. The resulting graph/matrix has the same Cholesky factor and the same elimination tree as $A^{T}A$.

Once the tree is found, some algorithms and many theorems require it to be postordered. In a postordered tree, the $d$ proper descendants of any node $k$ are numbered $k-d$ to $k-1$. If the resulting node renumbering is written as a permutation matrix $P$, then the filled graph of $A$ and $PAP^{T}$ are isomorphic (Liu Reference Liu1990). In other words, the two Cholesky factorizations have the same number of non-zeros and require the same amount of work. Their elimination trees are also isomorphic. Even if the symbolic or numeric factorization does not require a postordering of the etree, doing so makes the computations more regular by placing similar submatrices close to each other, thus improving memory traffic and speeding up the work (Liu Reference Liu1987c).

Figure 4.2.

Example matrix $A$, factor $L$, and elimination tree. Fill-in entries in $L$ are shown as circled x’s (Davis Reference Davis2006).

Figure 4.2 shows a matrix $A$, its Cholesky factor $L$, and its elimination tree ${\mathcal{T}}$. Figure 4.3 illustrates the matrix $PAP^{T}$, its Cholesky factor, and its elimination tree, where $P$ is the postordering of the elimination tree in Figure 4.2.

Figure 4.3.

After elimination tree postordering (Davis Reference Davis2006).

The elimination tree is fundamental to sparse matrix computations; Liu (Reference Liu1990) describes its many uses. Jess and Kees (Reference Jess and Kees1982) applied a restricted definition of the tree for matrices $A$ that had perfect elimination orderings, and Duff and Reid (Reference Duff and Reid1982) defined a related tree for the multifrontal method (Section 11). The elimination tree was first formally defined by Schreiber (Reference Schreiber1982). Kumar, Kumar and Basu (Reference Kumar, Kumar and Basu1992) present a parallel algorithm for computing the tree, and for symbolic factorization (Section 4.3).

4.2 Row and column counts

The row/column counts for a Cholesky factorization are the number of non-zeros in each row/column of $L$. These can be found without actually computing the pattern of $L$, and in time that is nearly $O(|A|)$, by using just the elimination tree and pattern of $A$. The column count is useful for laying out a good data structure for $L$, and row counts are helpful in determining data dependencies for parallel numeric factorization.

4.2.1 Row counts and row subtrees

The $k$th row subtree, ${\mathcal{T}}^{k}$, consists of a subtree of ${\mathcal{T}}$ and defines the non-zero pattern of the $k$th row of $L$. Each node $j$ in ${\mathcal{T}}^{k}$ corresponds to a non-zero $l_{kj}$. Its leaves are a subset of the $k$th row of the lower triangular part of $A$. The root of the tree is $k$. It is actually a tree, not a forest, and its nodes are the same as those traversed in the sparse lower triangular solve for the $k$th step of the up-looking Cholesky factorization. Since the tree can be found without finding the pattern of $L$, and since the pattern of $A$ is also known, the $k$th row subtree precisely defines the non-zero pattern of the $k$th row of $L$, without the need for storing each entry explicitly.

One simple but non-optimal method for computing the row and column counts is to just traverse each subtree and determine its size using the sparse triangular solve. The size of the tree ${\mathcal{T}}^{k}$ is number of non-zeros in row $k$ of $L$ (the row count). Visiting a node $j$ in ${\mathcal{T}}^{k}$ adds one to the column count of $j$. This method is simple to understand, but it takes $O(|L|)$ time.

The time can be reduced to nearly linear in $|A|$ (Gilbert, Li, Ng and Peyton Reference Gilbert, Li, Ng and Peyton2001, Gilbert, Ng and Peyton Reference Gilbert, Ng and Peyton1994). The basic idea is to decompose each row subtree into a set of disjoint paths, each starting with a leaf node and terminating at the least common ancestor of the current leaf and the prior leaf node. The paths are not traversed one node at a time. Instead, the lengths of these paths are found via the difference in the levels of their starting and ending nodes, where the level of a node is its distance from the root. The row count algorithm exploits the fact that all subtrees are related to each other; they are all subtrees of the elimination tree. Summing up the sizes of these disjoint paths gives the row counts.

The row count algorithm postorders the elimination tree, and then finds the level and first descendant of each node, which is the lowest numbered descendant of each node. This phase takes $O(n)$ time.

The next step traverses all row subtrees, considering column (node) $j$ in each row subtree in which it appears (that is, all non-zeros $a_{ij}$ in column $j$, starting at column $j=1$ and proceeding to $j=n$). If $j$ is a leaf in the $i$th row subtree, then the least common ancestor $c$ of this node and the prior leaf seen in the $i$th row subtree define a unique disjoint path, starting at $j$ and preceding up to but not including $c$. Finding $c$ relies on another application of the disjoint-set-union-find algorithm. Once $c$ is found, the length of the path from $j$ to the child of $c$ can be found quickly by taking the difference of the levels $j$ and $c$.

However, not all entries $a_{ij}$ are leaves $j$ in ${\mathcal{T}}^{i}$. The subset of $A$ that contains just leaves of row subtrees is called the skeleton matrix. The factorization of skeleton matrix of $A$ has the same non-zero pattern as $A$ itself. The skeleton is found by looking at the first descendants of the nodes in the tree. When the entry $a_{ij}$ is considered, node $j$ is a leaf of the $i$th row subtree, ${\mathcal{T}}^{i}$, if (and only if) the first descendant of $j$ is larger than any first descendant yet seen in that subtree.

Figure 4.4 shows the postordered skeleton matrix of $A$, denoted $\widehat{A}$, its factor $L$, its elimination tree ${\mathcal{T}}$, and the row subtrees ${\mathcal{T}}^{1}$ to ${\mathcal{T}}^{11}$ (compare with Figure 4.3). Entries in the skeleton matrix are shown as dark circles. A white circle denotes an entry in $A$ that is not in the skeleton matrix $\widehat{A}$.

Figure 4.4.

Skeleton matrix, factor $L$, elimination tree, and the row subtrees (Davis Reference Davis2006).

The first disjoint path found in the $i$th row subtree goes from $j$ up to $i$. Subsequent disjoint paths are defined by the pair of nodes $p$ and $j$, where $p$ is the prior entry $a_{ip}$ in the skeleton matrix, and also the prior leaf seen in this subtree. These two nodes define the unique disjoint path from $j$ to the common ancestor $c$ of $p$ and $j$. Summing up the lengths of these disjoint paths gives the number of nodes in the $i$th row subtree, which is the number of entries in row $i$ of $L$.

Liu (Reference Liu1986a) defined the row subtree and the skeleton matrix and used these concepts to create a compact data structure that represents the pattern of the Cholesky factor $L$ row-by-row using only $O(A)$ space, holding just the skeleton matrix and the elimination tree. Bank and Smith (Reference Bank and Smith1987) present a symbolic analysis algorithm in which the $i$th iteration traverses the $i$th row subtree, one disjoint path at a time, in the same manner as the symbolic up-looking Algorithm 4.1, except that it computes the elimination tree on the fly. They do not use the phrase ‘row subtree’ nor do they use the term elimination tree. However, their function $m(i)$ is a representation of the elimination tree since $m(i)$ is identical to the parent of $i$. Their algorithm computes the row counts and the elimination tree in $O(|L|)$ time.

4.2.2 Column counts

We now consider how to determine the number of non-zeros in each column of $L$. Let ${\mathcal{L}}_{j}$ denote the non-zero pattern of the $j$th column of $L$, and let ${\mathcal{A}}_{j}$ denote the pattern of the $j$th column of the strictly lower triangular part of $A$. Let $\widehat{{\mathcal{A}}}_{j}$ denote entries in the same part of the skeleton matrix $\widehat{A}$. George and Liu (Reference George and Liu1981) show that ${\mathcal{L}}_{j}$ is the union of its children in the elimination tree, plus the original entries in $A$:

(4.3)$$\begin{eqnarray}{\mathcal{L}}_{j}={\mathcal{A}}_{j}\cup \{j\}\cup \biggl(\mathop{\bigcup }_{j=\text{parent}(c)}{\mathcal{L}}_{c}\setminus \{c\}\biggr).\end{eqnarray}$$

A key corollary by Schreiber (Reference Schreiber1982) states that the non-zero pattern of the $j$th column of $L$ is a subset of the path $j{\rightsquigarrow}r$ from $j$ to the root of the elimination tree ${\mathcal{T}}$. The column counts are the sizes of each of these sets, $|{\mathcal{L}}_{j}|$. If $j$ is a leaf in ${\mathcal{T}}$, then the column count is simple: $|{\mathcal{L}}_{j}|=|{\mathcal{A}}_{j}|+1=|\widehat{{\mathcal{A}}}_{j}|+1$. Let $e_{j}$ denote the number of children of $j$ in ${\mathcal{T}}$. When $j$ is not a leaf, the skeleton entries are leaves of their row subtrees, and do not appear in any child, and thus

(4.4)$$\begin{eqnarray}|{\mathcal{L}}_{j}|=|\widehat{{\mathcal{A}}}_{j}|-e_{j}+\left|\mathop{\bigcup }_{j=\text{parent}(c)}{\mathcal{L}}_{c}\right|.\end{eqnarray}$$

Suppose there were an efficient way to find the overlap $o_{j}$ between the children $c$ in (4.4), replacing the set union with a summation:

(4.5)$$\begin{eqnarray}|{\mathcal{L}}_{j}|=|\widehat{{\mathcal{A}}}_{j}|-e_{j}-o_{j}+\biggl(\mathop{\sum }_{j=\text{parent}(c)}|{\mathcal{L}}_{c}|\biggr).\end{eqnarray}$$

As an example, consider column $j=4$ of Figure 4.4. Its two children are ${\mathcal{L}}_{2}=\{2,4,10,11\}$ and ${\mathcal{L}}_{3}=\{3,4,11\}$. In the skeleton matrix, $\widehat{{\mathcal{A}}}_{j}$ is empty, so

$$\begin{eqnarray}{\mathcal{L}}_{4}=\widehat{{\mathcal{A}}}_{4}\cup {\mathcal{L}}_{2}\setminus \{2\}\cup {\mathcal{L}}_{3}\setminus \{3\}=\emptyset \cup \{4,10,11\}\cup \{4,11\}=\{4,10,11\}.\end{eqnarray}$$

The overlap $o_{4}=2$, because rows 4 and 11 each appear twice in the children. The number of children is $e_{4}=2$. Thus, $|{\mathcal{L}}_{4}|=0-2-2+(4+3)=3$.

The key observation is see that if $j$ has $d$ children in a row subtree, then it will be the least common ancestor of exactly $d-1$ successive pairs of leaves in that row subtree. For example, node 4 is the least common ancestor of leaves 1 and 3 in ${\mathcal{T}}^{4}$, and the least common ancestor of leaves 2 and 3 in ${\mathcal{T}}^{11}$. Thus, each time column $j$ becomes a least common ancestor of any successive pair of leaves in the row count algorithm, the overlap $o_{j}$ can be incremented. When the row count algorithm completes, the overlaps are used to recursively compute the column counts of $L$ using (4.5). As a result, the time taken to compute the row and column counts is nearly $O(|A|)$.

The column counts give a precise definition of the optimal amount of floating-point work required for any sparse Cholesky factorization (Rose Reference Rose and Read1972, Bunch Reference Bunch and Traub1973):

(4.6)$$\begin{eqnarray}\mathop{\sum }_{j=1}^{n}|{\mathcal{L}}_{j}|^{2}.\end{eqnarray}$$

Hogg and Scott (Reference Hogg and Scott2013a) extend the row/column count algorithm of Gilbert et al. (Reference Gilbert, Ng and Peyton1994) to obtain an efficient analysis phase for finite element problems, where the matrix $A$ can be viewed as a collection of cliques, each clique being a single finite element. The extension of the row/column count algorithm to LU and QR factorization is discussed in Section 7.1.

4.3 Symbolic factorization

The final step in the symbolic analysis is the symbolic factorization, which is to find the non-zero pattern ${\mathcal{L}}$ of the Cholesky factor $L$, for equation (4.3). The time taken for this step is $O(|L|)$ if each entry in $L$ is explicitly represented. George and Liu (Reference George and Liu1980c) reduce this to time proportional to the size of a compressed representation related to the supernodal and multifrontal factorizations. We will first consider the $O(|L|)$-time methods.

Both the up-looking (4.1) and left-looking factorizations provide a framework for constructing ${\mathcal{L}}$ in compressed-column form. Recall that ${\mathcal{L}}_{j}$ denotes the non-zero pattern of the $j$th column of $L$. Since the row and column counts are known, the data structure for $L$ can be statically allocated with enough space ($|{\mathcal{L}}_{j}|$) in each column $j$ to hold all its entries.

The left-looking method simply computes (4.3), one column at a time, from left to right (for $j=1$ to $n$). Since each child takes part in only a single computation of (4.3), and since the work required for the set unions is the sum of the set sizes, the total time is $O(|L|)$. This method does not produce ${\mathcal{L}}_{j}$ in sorted order, however. If this is needed by the numeric factorization, another $O(|L|)$-time bucket sort is required.

The up-looking symbolic factorization is based on the up-looking numeric factorization (4.1), in which the $k$th row is computed from a triangular solve. The symbolic method constructs ${\mathcal{L}}$ one row $k$ at a time by traversing each row subtree (Davis Reference Davis2006). Selecting an entry $a_{kj}$ in the lower triangular part of $A$, Algorithm 4.1 walks up the elimination tree until it sees a marked node, and marks all nodes along the way. It will stop at the root $k$ since this starts out as marked. Each node seen in this traversal is one node of the $k$th row subtree. As a by-product, each ${\mathcal{L}}_{j}$ is in sorted order.

Rose et al. (Reference Rose, Tarjan and Lueker1976) present the first algorithm for computing ${\mathcal{L}}$, using what could be seen as an abbreviated right-looking method. The downside of their method is that it requires a more complex dynamic data structure for ${\mathcal{L}}$ that takes more than $O(|L|)$ space. In their method, the sets ${\mathcal{L}}_{j}$ are initialized to the non-zero pattern of the lower triangular part of the matrix $A$ (that is, ${\mathcal{A}}_{j}\cup \{j\}$ from (4.3)). At each step $j$, duplicates in ${\mathcal{L}}_{j}$ are pruned, and then the set ${\mathcal{L}}_{j}$ is added to its parent in the elimination tree,

$$\begin{eqnarray}{\mathcal{L}}_{\text{parent}(j)}={\mathcal{L}}_{\text{parent}(j)}\cup ({\mathcal{L}}_{j}\setminus \{j\}).\end{eqnarray}$$

As a result, by the time the $j$th step starts, (4.3) has been computed for column $j$ since all the children of $j$ have been considered. Normally, computing the union of two sets takes time proportional to the size of the two sets, but the algorithm avoids this by allowing duplicates to temporarily reside in ${\mathcal{L}}_{\text{parent}(j)}$. Note that Rose et al. did not use the term ‘elimination tree’ since the term had not yet been introduced, but we use the term here since it provides a concise description of their method.

George and Liu (Reference George and Liu1980c) present the first linear-time symbolic factorization method using $O(|L|)$ space, based on the path lemma and the quotient graph. To describe this method, we first need to consider two sub-optimal methods: one using only the graph of $A$, and the second using the elimination graph.

Recall that the path lemma states that $l_{ij}\neq 0$ if and only there is a path in the graph of $A$ from $i$ to $j$ whose intermediate vertices are all numbered less than $\min (i,j)$ (that is, excluding $i$ and $j$ themselves). One method for constructing ${\mathcal{L}}_{j}$ is to apply the path lemma directly on $A$, and determine all nodes $i$ reachable from $j$ by traversing only nodes 1 to $j-1$ in the graph of $A$. This first alternative method based on the graph of $A$ would result in a compact data representation but would require a lot of work.

The second alternative is to mimic the outer-product factorization (4.2), by constructing a sequence of elimination graphs as in Algorithm 4.2 below. Let $G_{0}$ be the graph of $A$ (with an edge $(i,j)$ if and only if $a_{ij}\neq 0$). Let $G_{j}$ be the graph of $A_{22}$ in (4.2) after the first $k$ nodes have been eliminated. Eliminating a node $j$ adds a clique to the graph between all neighbours of $j$.

Algorithm 4.2 is not practical since it takes far too much time (the same time as the numeric factorization). However, Eisenstat, Schultz and Sherman (Reference Eisenstat, Schultz, Sherman, Bunch and Rose1976a) show how to represent the sequence of graphs in a more compact way, which George and Liu (Reference George and Liu1980c) refer to as the quotient graph. The quotient graph ${\mathcal{G}}$ is a graph that contains two kinds of nodes: those corresponding to uneliminated nodes in $G$, and those corresponding to a subset of eliminated nodes in $G$. Each eliminated node $e$ results in a clique in $G$. If the eliminated node $e$ has $t$ neighbours, this requires $O(t^{2})$ new edges in $G$. Instead, we can create a new kind of node which George and Liu call a supernode. This is not quite the same as the term supernode used elsewhere in this article, so we will use Eisenstat et al.’s (Reference Eisenstat, Schultz, Sherman, Bunch and Rose1976a) term element to avoid confusion. The element $e$ represents the clique implicitly, with only $O(t)$ edges.

Algorithm 4.3 uses the notation of Amestoy, Davis and Duff (Reference Amestoy, Davis and Duff1996a). Let ${\mathcal{A}}_{j}$ denote the (regular) nodes adjacent to $j$ in ${\mathcal{G}}$, and let ${\mathcal{E}}_{j}$ represent the elements adjacent to $j$. Finally, let ${\mathcal{L}}_{e}$ represent the regular nodes adjacent to element $e$; we use this notation because it is the same as the pattern of column $e$ of $L$. If one clique is contained within another, it can be deleted without changing the graph. As a result, no elements are adjacent to any other elements in the quotient graph.

Edges are pruned from ${\mathcal{A}}$ as Algorithm 4.3 progresses. Suppose there is an original edge $(i,k)$ between two nodes $i$ and $k$, but an element $j$ is formed that includes both of them. The two nodes are still adjacent in $G$ but this is represented by element $j$. The explicit edge $(i,k)$ is no longer needed. The quotient graph elimination process takes $O(|L|)$ time, and because edges are pruned and assuming ${\mathcal{L}}_{e}$ is deleted after it is used, Algorithm 4.3 only needs $O(|A|)$ memory. Keeping the pattern of $L$ requires $O(|L|)$ space.

Figures 4.5 and 4.6 from Amestoy et al. (Reference Amestoy, Davis and Duff1996a) present the sequence of elimination graphs, quotient graphs, and the pattern of the factors after three steps of elimination. All three represent the same thing in different forms. In the graph ${\mathcal{G}}_{2}$, element 2 represents a clique of nodes 5, 6, and 9, and thus the original edge (5,9) is redundant and has been pruned. In ${\mathcal{G}}_{5}$, element 5 is adjacent to prior elements 2 and 3, and thus the latter two elements are deleted.

Figure 4.5.

Elimination graph, quotient graph, and matrix (Amestoy et al. Reference Amestoy, Davis and Duff1996a).

Figure 4.6.

Elimination graph, quotient graph, and matrix (continued) (Amestoy et al. Reference Amestoy, Davis and Duff1996a).

Algorithm 4.3 takes $O(|L|)$ time and memory to construct the non-zero pattern of $L$. Eisenstat, Schultz and Sherman (Reference Eisenstat, Schultz, Sherman, Bunch and Rose1976b) describe how the non-zero pattern of $L$ and $U$ can be represented in less space than one integer per non-zero entry. Suppose the matrix $L$ is stored by columns. If the non-zero pattern of a column $j$ of $U$ is identical to a subsequence of a prior row, then the explicit storage of the pattern of the $j$th column can be deleted, and replaced with a reference to the subsequence in the prior column. George and Liu (Reference George and Liu1980c) exploit this property to reduce the time and memory required for symbolic factorization, and extend it by relying on the elimination tree to determine when this property occurs. Suppose a node $j$ in the etree has a single child $c$, and that ${\mathcal{L}}_{j}={\mathcal{L}}_{c}\setminus \{c\}$. This is a frequently occurring special case of (4.3). If the etree is postordered, $c=j-1$ as well, further simplifying the representation. There is no need to store the pattern of ${\mathcal{L}}_{j}$ explicitly. If the row indices in ${\mathcal{L}}_{j-1}$ appear in sorted order, then ${\mathcal{L}}_{j}$ is the same list, excluding the first entry.

This observation can also be used to speed up the up-looking, left-looking, and right-looking symbolic factorization algorithms just described above. Each of them can be implemented in time and memory proportional to the size of the representation of the pattern of $L$. The savings can be significant: for an $s\times s$ two-dimensional mesh using a nested dissection ordering, the matrix $A$ is $n\times n$ with $n=s^{2}$. The numeric factorization takes $O(s^{3})$ time and $|L|$ is $O(s^{2}\log s)$, but the compact representation of $L$ takes only $O(s^{2})$, or $O(n)$ space, and takes the same time to compute.

George, Poole and Voigt (Reference George, Poole and Voigt1978) and George and Rashwan (Reference George and Rashwan1980) extend the quotient graph model to combine nodes in the graph of $A$ for the symbolic and numeric Cholesky factorization of a block partitioned matrix, motivated by the finite element method.

Gilbert (Reference Gilbert1994) surveys the use of graph algorithms for the symbolic analysis for Cholesky, QR, LU factorization, eigenvalue problems, and matrix multiplication. Most of his paper focuses on the use of directed graphs for non-symmetric problems, but many of the results he covers are closely related to the Cholesky factorization. For example, the Cholesky factorization of $A^{T}A$ provides an upper bound for the non-zero pattern of QR and LU factorization of $A$. Pothen and Toledo (Reference Pothen, Toledo, Mehta and Sahni2004) give a more recent tutorial survey of algorithms and data structures for the symbolic analysis of both symmetric and unsymmetric factorizations. They discuss the elimination tree, skeleton graph, row/column counts, symbolic factorization (topics discussed above), and they also consider unsymmetric structures for LU and QR factorization.

Symbolic factorization algorithms take time proportional to the size of their output, which is often a condensed form of the pattern of $L$ not much bigger than $O(|A|)$ itself. As a result, parallel speedup is very difficult to achieve. The practical goal of a parallel symbolic factorization algorithm is to exploit distributed memory to solve problems too large to fit on any one processor, which all of the following parallel algorithms achieve.

The first parallel algorithm to accomplish this goal was that of George, Heath, Ng and Liu (Reference George, Heath, Ng and Liu1987, Reference George, Heath, Liu and Ng1989a), which computes ${\mathcal{L}}$ column-by-column in a distributed-memory environment. Their method assumes the elimination tree is already known, and obtains a very modest speedup. Ng (Reference Ng1993) extends this method by exploiting supernodes, which improves both sequential and parallel performance. Zmijewski and Gilbert (Reference Zmijewski and Gilbert1988) show how to compute the tree in parallel. The entries $A$ can be partitioned, and each processor computes its own elimination tree (forest, to be precise) on the entries it owns. These forests are then merged to construct the elimination tree of $A$. Once the elimination tree is known, their parallel symbolic factorization constructs the row subtrees of $L$, each one as an independent problem. They do not obtain any parallel speedup when computing the elimination tree, however. Likewise, Kumar et al. (Reference Kumar, Kumar and Basu1992) construct both the etree and ${\mathcal{L}}$ in parallel, and they obtain a modest parallel speedup in constructing the tree and in the symbolic factorization. Gilbert and Hafsteinsson (Reference Gilbert and Hafsteinsson1990) present a highly parallel shared-memory (CRCW) algorithm for finding the elimination tree and the pattern ${\mathcal{L}}$ using one processor for each entry in $L$, in $O(\log ^{2}n)$ time. They do not present an implementation, however.

Parallel symbolic factorization methods for LU factorization (Grigori, Demmel and Li Reference Grigori, Demmel and Li2007b), supernodal Cholesky (Ng Reference Ng1993), and multifrontal methods (Gupta, Karypis and Kumar Reference Gupta, Karypis and Kumar1997) are discussed in Sections 6.1, 9.1, and 11, respectively.

5.1 Envelope, skyline, and profile methods

Current algorithms exploit every zero entry in the factor, or nearly all of them, using the graph theory and algorithms discussed in Section 4. However, the earliest methods did not have access to these developments. Instead, they were only able to exploit the envelope. Using current terminology, the $k$th row subtree ${\mathcal{T}}^{k}$ describes the non-zero pattern of the $k$th row of $L$ (refer to Figure 4.4). The leaves of this tree are entries in the skeleton matrix $\widehat{A}$, and the non-zero $a_{kj}$ with the smallest column index $j$ will always be the leftmost leaf. The $k$th row subtree will be a subset of nodes $j$ to $k$. That is, any non-zeros in the $k$th row of $L$ can only appear in columns $j$ to $k$. An algorithm that stores all of the entries $j$ to $k$ is called an envelope method. These methods are also referred to as profile or skyline methods (the latter is a reference to $L^{T}$).

Jennings (Reference Jennings1966) created the first envelope method for sparse Cholesky. Felippa (Reference Felippa1975) adapted the method for finite element problems by partitioning the matrix into two sets, according to internal degrees of freedom (those appearing only with a single finite element) and external degrees of freedom. Neither considered permutations to reduce the profile. While envelope/profile/skyline factorization methods are no longer commonly used, profile reduction orderings (Section 8.2) are still an active area of research since they are very well-suited to frontal methods (Section 10).

George and Liu (Reference George and Liu1978a, Reference George and Liu1978b, Reference George and Liu1979a) went beyond just two partitions in their recursive block partitioning method for irregular finite element problems. The diagonal block of each partition is factorized via an envelope method. They also consider ordering methods to find the partitions and to reduce the profiles of the diagonal blocks. The factorizations of the off-diagonal blocks are not stored, but computed column-by-column when needed in the forward/backsolves (George Reference George1974).

Bjorstad (Reference Bjorstad1987) uses a similar partitioned strategy, also exploiting parallelism by factorizing multiple partitions in parallel (each with a sequential profile method). Updates from each partition are applied in a right-looking manner and held on disk, similar to the strategy of George and Rashwan (Reference George and Rashwan1985), discussed in Section 5.3.

5.2 Up-looking Cholesky

Unlike the envelope method described in the previous section, the up-looking Cholesky factorization method presented here can exploit every zero entry in $L$, and is asymptotically optimal in the work it performs. It computes each row of $L$, one at a time, starting with row 1 and proceeding to row $n$. The $k$th step requires a sparse triangular solve, with a sparse right-hand side, using the rows of $L$ already computed ($L_{1:k-1,1:k-1}$). It is also called the bordering method and row-Cholesky. It is not the first asymptotically optimal sparse Cholesky factorization algorithm, but it is presented first since it is closely related to presentation of the sparse triangular solve and symbolic analysis in Sections 3.2 and 4. It appears below in MATLAB.

Consider the up-looking factorization Algorithm 4.1 presented in Section 4.3. The algorithm traverses each disjoint subpath of the $k$th row subtree, ${\mathcal{T}}^{k}$, but it does not traverse them in the topological order required for the numeric factorization, since this is not required for the symbolic factorization.

A simple change to the method results in a proper topological order that satisfies the numerical dependencies for the sparse triangular solve. Consider the computation of row 11 of $L$, from Figure 4.4. Suppose the non-zeros of $A$ in columns 2, 3, 4, 8, and 10 are visited in that order (they can actually be visited in any order). The first disjoint path seen is $(2,4,10,11)$, from the first non-zero $a_{11,2}$ to the root of ${\mathcal{T}}^{11}$. These nodes are all marked. Next, the single node 3 is in a disjoint path by itself, since the path starts from the non-zero $a_{11,3}$ and halts when the marked node 4 is seen. Node 4 is considered because $a_{11,4}$ is non-zero, but node 4 is already marked so no path is formed. Node 8 gives the path $(8,9)$. Finally, node 10 is skipped because it is already marked.

The resulting disjoint paths are $(2,4,10,11)$, $(3)$, and $(8,9)$, and are found in that order. These three paths cannot be traversed in that order for the triangular solve since (for example) the non-zero $l_{4,3}$ requires node 3 to appear before node 4. However, if these three paths are reversed, as $(8,9)$, then $(3)$, and finally $(2,4,10,11)$, we obtain a non-zero pattern for the sparse triangular solve as ${\mathcal{X}}=(8,9,3,2,4,10,11)$. This ordering of ${\mathcal{X}}$ is topological. If you consider any pair of nodes $j$ and $i$ in the list ${\mathcal{X}}$, then any edge in the graph $G_{L}$ (with an edge $(j,i)$ for each $l_{ij}\neq 0$) will go from left to right in that list. Performing the triangular solve from Section 3.2 in this order satisfies all numerical dependencies.

Sparse Cholesky factorization via the up-looking method was first considered by Rose, Whitten, Sherman and Tarjan (Reference Rose, Whitten, Sherman and Tarjan1980), but they did not provide an algorithm. Liu (Reference Liu1986a) introduces the row subtrees and the skeleton matrix as a method for a compact row-by-row storage scheme for $L$. Each row subtree can be represented by only its leaves (each of which is an entry in the skeleton matrix). Liu did not describe a corresponding up-looking factorization method, however.

Bank and Smith (Reference Bank and Smith1987) describe an up-looking numeric factorization algorithm that is a companion to their up-looking symbolic factorization algorithm. They suggest two methods to handle the numerical dependencies in the triangular solve: (1) explicit sort, and (2) a traversal of the entries in each row of $A$ in reverse order, which produces a topological order. This second method is optimal, and it appears to be the first instance of a topological ordering for a sparse triangular solve. However, to save space, they do not store the non-zero pattern of $L$. Instead, they construct each row or column as needed by traversing each row subtree. This reduces space but it leads to a non-optimal amount of work for their numeric factorization since it must rely on dot products of pairs of sparse vectors for the triangular solve. The computational cost of their method is analysed by Bank and Rose (Reference Bank and Rose1990).

Liu (Reference Liu1991) was the first to implement an asymptotically optimal up-looking method, via a generalization of the envelope method. The method partitions the matrix $L$ into blocks. Each diagonal block corresponds to a chain of consecutive nodes in the elimination tree. That is, each chain consists of a sequence of $t$ columns $j,j+1,\ldots ,j+t$, where the parent of each column is the very next column in the sequence. The key observation is that the diagonal block $L_{j:j+t,j:j+t}$ has a full envelope, and the submatrix below this diagonal block, namely, $L_{j+t+1:n,j:j+t}$, also has a full envelope structure. In any given row $i$ of the subdiagonal block, if $l_{ik}$ is non-zero for some $k$ in the sequence $j,j+1,\ldots ,j+t$, then all entries in $L_{i,k:j+t}$ must also be non-zero. This is because each row subtree is composed of a set of disjoint paths, and the postordering of the elimination tree ensures that the subpaths consist of contiguous subsequences of the diagonal blocks.

Davis’s up-looking method (Reference Davis2005, Reference Davis2006) relies on a vanilla compressed-column data structure for $L$, which also results in an asymptotically optimal amount of work.

Although the left-looking, supernodal, and multifrontal Cholesky factorizations are widely used and appear more frequently in the literature, the up-looking method is also widely used in practice. In particular, MATLAB relies on CHOLMOD for   and chol(A) when A is a sparse symmetric definite matrix (Chen, Davis, Hager and Rajamanickam Reference Chen, Davis, Hager and Rajamanickam2008). CHOLMOD’s symbolic analysis computes the row and column counts (Section 4.2), which also gives the floating-point work. If the ratio of the work over the number of non-zeros in $L$ is less than 40, CHOLMOD uses the up-looking method presented here, and a left-looking supernodal method otherwise (Section 9.1). This is because the up-looking method is fast in practice for very sparse matrices, as compared to supernodal and multifrontal methods.

5.3 Left-looking Cholesky

The left-looking Cholesky factorization algorithm is widely used and has been the focus of more research articles than the up-looking method. It is also called the fan-in or backward-looking method. It forms the foundation of the left-looking supernodal Cholesky factorization (Section 9.1). The method computes $L$ one column at a time, and thus it is also called column-Cholesky. The method can be derived from the expression

(5.1)$$\begin{eqnarray}\left[\begin{array}{@{}ccc@{}}L_{11} & & \\ l_{12}^{T} & l_{22} & \\ L_{31} & l_{32} & L_{33}\end{array}\right]\left[\begin{array}{@{}ccc@{}}L_{11}^{T} & l_{12} & L_{31}^{T}\\ & l_{22} & l_{32}^{T}\\ & & L_{33}^{T}\end{array}\right]=\left[\begin{array}{@{}ccc@{}}A_{11} & a_{12} & A_{31}^{T}\\ a_{12}^{T} & a_{22} & a_{32}^{T}\\ A_{31} & a_{32} & A_{33}\end{array}\right]\!,\end{eqnarray}$$

where the middle row and column of each matrix are the $k$th row and column of $L$, $L^{T}$, and $A$, respectively. If the first $k-1$ columns of $L$ are known, then the $k$th column of $L$ can be computed as follows:

(5.2)$$\begin{eqnarray}\begin{array}{@{}l@{}}l_{22}=\sqrt{a_{22}-l_{12}^{T}l_{12}},\\ l_{32}=(a_{32}-L_{31}l_{12})/l_{22}.\end{array}\end{eqnarray}$$

These two expressions can be folded together into a single operation, a sparse matrix times sparse vector, that computes a vector $c$ of length $n-k+1$,

(5.3)$$\begin{eqnarray}c=\left[\begin{array}{@{}c@{}}c_{1}\\ c_{2}\end{array}\right]=\left[\begin{array}{@{}c@{}}a_{22}\\ a_{32}\end{array}\right]-\left[\begin{array}{@{}c@{}}l_{12}^{T}\\ L_{31}\end{array}\right]l_{12},\end{eqnarray}$$

where $c_{1}$ and $a_{22}$ are scalars, and $c_{2}$ and $a_{32}$ are vectors of size $n-k$. Computing $c$ is the bulk of the work for step $k$ of the algorithm; this is followed by

(5.4)$$\begin{eqnarray}\displaystyle l_{22} & \displaystyle =\sqrt{c_{1}}, & \displaystyle \\ \displaystyle l_{32} & \displaystyle =c_{2}/l_{22}.\end{eqnarray}$$

The MATLAB equivalent is given below.

The key observation is that this sparse matrix times sparse vector multiply (5.3) needs to be done only for columns for which $l_{12}^{T}$ is non-zero, which corresponds to a traversal of each node in the $k$th row subtree, ${\mathcal{T}}^{k}$. In MATLAB notation the computation of c(k:n) is then as follows.

In Algorithm 5.1 below, $L$ is stored in compressed-column form. It requires access to the $k$th row of $L$, or $l_{12}^{T}$ in (5.3) at the $k$ step. Accessing an arbitrary row in this data structure would take too much time, but fortunately the access to the rows of $L$ is not arbitrary. Rather, the algorithm requires access to each row in a strictly ascending order. To accomplish this, the left-looking algorithm keeps track of a pointer into each column that advances to the next row as each row is accessed (the workspace $p$).

The algorithm traverses ${\mathcal{T}}^{k}$ the same way as the $k$th iteration of Algorithm 4.1, so the details are not given here. The cmod and cdiv operations defined in Algorithm 5.1 are discussed below.

Rose et al. (Reference Rose, Whitten, Sherman and Tarjan1980) give an overview of three sparse Cholesky methods (right-looking, up-looking, and left-looking), and ordering methods. The left-looking method they consider is YSMP, by Eisenstat et al. (Reference Eisenstat, Schultz and Sherman1975, Reference Eisenstat, Gursky, Schultz and Sherman1982, Reference Eisenstat, Schultz and Sherman1981). YSMP is actually described as computing $L^{T}$ one row at a time, which is equivalent to the left-looking column-by-column method in Algorithm 5.1 above. SPARSPAK is also based on the left-looking sparse Cholesky algorithm (Chu, George, Liu and Ng Reference Chu, George, Liu and Ng1984, George and Liu Reference George and Liu1981, George and Ng Reference George and Ng1985a, George and Ng Reference George and Ng1984b).

YSMP and SPARSPAK both rely on a set of $n$ dynamic linked lists to traverse the row subtrees, one for each row. Each column $j$ is in only a single linked list at a time. In Algorithm 5.1, $p_{j}$ refers to the location in column $j$ of the next entry that will be accessed. If the entry has row index $t$, then column $j$ will reside in list $t$. When $p_{j}$ advances, column $j$ is moved to the linked list for the next row. As a result, when step $k$ commences, the $k$th linked list contains all nodes in the $k$th row subtree, ${\mathcal{T}}^{k}$.

The left-looking sparse Cholesky algorithm has been the basis of many parallel algorithms for both shared-memory and distributed-memory parallel computers. It is easier to exploit parallelism in this method as compared to the up-looking method.

The elimination tree plays a vital role in all parallel sparse direct methods. Different methods (e.g. left-looking, right-looking, supernodal, multifrontal) define their tasks in many different ways, but in all the methods the tree governs the parallelism. In some methods, each node of the tree is a single task, but more often, the work at a given node is split into multiple tasks. Tasks on independent branches of the tree can be computed in parallel.

George, Heath, Liu and Ng (Reference George, Heath, Liu and Ng1986a) present the first parallel left-looking sparse Cholesky method, using a shared-memory computer. Algorithm 5.1 contains two basic tasks, each of which operate on the granularity of one or two columns: (1) cmod($k,j$), the modification of column $k$ by column $j$, and (2) cdiv($k$), the division of column $k$ by the diagonal, $l_{kk}$. Their parallel version of Algorithm 5.1 uses the same linked list structure as their sequential method, and one independent task per column $k$. Task $k$ waits until a column $j$ appears in the $k$th link list, and then performs the cmod($k,j$) and moves the column $j$ to the next linked list (incrementing $p_{j}$). When all nodes in row $k$ have been processed, task $k$ finishes by performing the cdiv($k$) task. Only task $k$ needs write-access to column $k$. Each node of the elimination tree defines a single task. In their method, a task can overlap with other tasks below it and above it in the tree, but a column $j$ can only modify one ancestor column $k$ at a time. As an example, consider the matrix and tree in Figure 4.3. Suppose columns 2 and 4 have been computed (tasks 2 and 4 are done), and cmod(10,2) has been completed by task 10, which then places column 2 in the list for task 11. Then task 10 can do cmod(10,4) at the same time task 11 does cmod(11,2). However, tasks 10 and 11 cannot do cmod(10,2) and cmod(11,2), respectively, at the same time; task 11 must wait until task 10 finishes cmod(10,2) before it can use column 2 for cmod(11,2).

Liu (Reference Liu1986b) considers three models for sparse Cholesky: fine-, medium-, and coarse-grained, placing them in a common framework. The fine-grained model considers each floating-point operation as its own task; see for example the LU factorization method by Wing and Huang (Reference Wing and Huang1980). The coarse-grained model is exemplified by Jess and Kees (Reference Jess and Kees1982), who propose a parallel right-looking LU factorization for matrices with symmetric structure (the methods of Wing and Huang (Reference Wing and Huang1980) and Jess and Kees (Reference Jess and Kees1982) are discussed in Section 6.3, on the right-looking LU factorization). The medium-grained model introduced by Liu (Reference Liu1986b) considers each cmod and cdiv as its own task. Each node of the $n$ nodes in Liu’s graph is a cdiv, and each edge corresponds to a cmod, and thus the graph has same structure as ${\mathcal{L}}$, with $|{\mathcal{L}}|$ tasks. Task cmod($k,j$) for any non-zero $l_{kj}$ in row $k$ must precede cdiv($k$), which in turn must precede cmod($i,k$) for any non-zero $l_{ik}$ in column $k$. Using the same example of Figure 4.3, in this model cmod(10,2) and cmod(11,2) can be done in parallel. Coalescing the tasks cdiv($k$) with all tasks cmod($i,k$) for each non-zero $l_{ik}$ in column $k$ results in the coarse-grained right-looking Cholesky, whereas combining all cmod($k,j$) for each $l_{kj}$ in row $k$, with cdiv($k$), results in a parallel left-looking Cholesky method.

Once the left-looking factorization progresses to step $k$, it no longer needs the first $k-1$ rows of $L$ (the $L_{11}$) matrix. George and Rashwan (Reference George and Rashwan1985) exploit this property in their out-of-core method. They partition the matrix with incomplete nested dissection. Submatrices are factorized with a left-looking method, and then the remainder of the unfactorized matrix is updated (a right-looking phase) and written to disk. It is read back in when subsequent submatrices are factorized. Liu (Reference Liu1987a) also exploits this property in his out-of-core method. Unlike George and Rashwan’s (Reference George and Rashwan1985) method, Liu’s method is purely left-looking. The columns of $L$ are computed one at a time, and the memory space required grows as a result. If memory is exhausted, the $L_{11}$ matrix is written to disk, and the factorization continues. The $L_{21}$ matrix (rows $k$ to $n$) remains in core. To improve memory usage, the matrix is permuted via a generalization of the pebble game, applied to the elimination tree (Liu Reference Liu1987b).

George, Heath, Liu and Ng (Reference George, Heath, Liu and Ng1988a) extend their parallel left-looking method to the distributed-memory realm, where no processors share any memory and all data must be explicitly communicated by sending and receiving messages. Using the nomenclature of Ashcraft, Eisenstat, Liu and Sherman (Reference Ashcraft, Eisenstat, Liu and Sherman1990b), the method becomes right-looking, and so it is considered in Section 5.4.

If $l_{ij}$ is non-zero, then at some point column $j$ must update column $i$, via cmod($i,j$). In a left-looking method, the update from column $j$ to column $i$ (cmod($i,j$)) is done at step $i$, according to the target column. In a right-looking method, cmod($i,j$) is done at step $j$, according to the source column. The difference between left/right-looking in a parallel context is subtle, because it depends on where you are standing and which processor is being considered: the one sending an update or the one receiving it. Multiple steps can execute in parallel, although numerical dependencies must be followed. That is, if $l_{ij}$ is non-zero, then column $i$ must be finalized before cmod($i,j$) can be computed.

Ashcraft, Eisenstat and Liu (Reference Ashcraft, Eisenstat and Liu1990a) observe that the method of George et al. (Reference George, Heath, Liu and Ng1989a) sends more messages than necessary. Suppose that one processor A owns both columns $k_{1}$ and $k_{2}$, and both columns need to update a target column $i$ owned by a second processor B. In George et al.’s (Reference George, Heath, Liu and Ng1989a) method processor A sends two messages to processor B – column $k_{1}$ and column $k_{2}$ – so that processor B can compute cmod($i,k_{1}$) and cmod($i,k_{2}$). In the left-looking (fan-in) method of Ashcraft et al. (Reference Ashcraft, Eisenstat and Liu1990a), the update (5.3) for these two columns of $L$ is combined by processor A, which then sends only a single column to processor B, as an aggregated update. Constructing aggregate updates takes extra memory, and if this is not available, Eswar, Huang and Sadayappan (Reference Eswar, Huang and Sadayappan1994) describe a left-looking method that delays the construction of the update until other aggregate updates have been computed, sent, and freed.

All of the parallel methods described so far assign one or more entire columns of $L$ to a single processor, resulting in a one-dimensional assignment of tasks to processors (columns only). Schreiber (Reference Schreiber, George, Gilbert and Liu1993) shows that any one-dimensional mapping is inherently non-scalable because of communication overhead, and that a two-dimensional mapping of computations to processors is required instead (Gilbert and Schreiber Reference Gilbert and Schreiber1992). These mappings are possible in supernodal and multifrontal methods, discussed in Sections 9 and 11, and in a two-dimensional right-looking method considered in the next section.

5.4 Right-looking Cholesky

Right-looking Cholesky, also known as fan-out or submatrix-Cholesky, is based on equation (4.2). It is described as follows in MATLAB notation as chol_right.

At step $k$, the outer product L(k+1:n,k)*L(k+1:n,k)’ is subtracted from the lower right $(n-k)\times (n-k)$ submatrix. This is difficult to implement since it can take extra work to find the entries in $A$ to modify. There are many target columns, and any given target column may have many non-zeros that are not modified by the $k$th update. This is not a concern in the left- or up-looking methods, since in those methods the single target column is temporarily held in a work vector of size $n$, accessed via gather/scatter.

George et al. (Reference George, Heath, Liu and Ng1988a) present a parallel right-looking method for a distributed-memory computer. Each processor owns a set of columns. Any column $k$ with no non-zeros in row $k$ of $L$ can be processed immediately via cdiv($k$); this node $k$ is a leaf of the elimination tree. After a processor owning column $k$ does its cdiv($k$), it sends column $k$ of $L$ to any processor owning any column $i$ for which $l_{ik}$ is non-zero. If a receiving processor owns more than one such column $i$, the message is sent only once. When a processor receives a column $j$, it computes cmod($k,j$) for any column $k$ it owns. This is also a right-looking view, since a single column $j$ is applied to the submatrix of all columns owned by this processor. If given a single processor, the method is identical to the right-looking method, since one processor owns the entire matrix. When all cmods have been applied to a column $k$ that it owns, it does cdiv($k$) and sends it out, as just described. Column tasks correspond to each node of the elimination tree, and are assigned to processors in a simple wrap-around manner. That is, the leaves are all doled out to processors in a round-robin manner, followed by all nodes one level up from the levels, and so on to the root.

George et al. (Reference George, Heath, Liu and Ng1989a) extend their method to the hypercube. They change the task assignment so that whole subtrees of the elimination tree are given to individual processors. In this case, the subtree is a node $k$ and all its descendants, not to be confused with the $k$th row subtree, ${\mathcal{T}}^{k}$. They also show that performance can be improved via pipelining. When a processor receives a column $j$, it does all cmod($k,j$) for all columns $k$ that it owns. However, rather than waiting until all such cmod($k,j$)s are finished, it finalizes any column $k$ that is now ready for its cdiv($k$) and sends it out, before continuing with the rest of the cmods for this incoming column $j$. George, Liu and Ng (Reference George, Liu and Ng1989b) analyse the communication between processors in this method and show that it is asymptotically optimal for a two-dimensional mesh with the nested dissection ordering. Gao and Parlett (Reference Gao and Parlett1990) augment the analysis of George et al. (Reference George, Liu and Ng1989b), showing that not only is the total communication volume minimized, but it is also balanced across the processors.

Ordering methods such as minimum degree can generate unbalanced trees. Geist and Ng (Reference Geist and Ng1989) generalize the subtree-to-subcube task assignment of George et al.’s (Reference George, Heath, Liu and Ng1989a) method via a bin-packing to account for heavily unbalanced trees. Eswar, Huang and Sadayappan (Reference Eswar, Huang and Sadayappan1995) consider task assignments that combine multiple strategies (wrap-based and subtree-based) for mapping processors to columns.

Zhang and Elman (Reference Zhang and Elman1992) explore several shared-memory variants: the left-looking methods of George et al. (Reference George, Heath, Liu and Ng1986a) and Ashcraft et al. (Reference Ashcraft, Eisenstat and Liu1990a), and the task-scheduling method of Geist and Ng (Reference Geist and Ng1989). The latter two are distributed-memory algorithms. Zhang and Elman report that Geist and Ng’s (Reference Geist and Ng1989) method works well in a shared-memory environment.

Ashcraft et al. (Reference Ashcraft, Eisenstat, Liu and Sherman1990b) and Eswar, Sadayappan and Visvanathan (Reference Eswar, Sadayappan, Visvanathan, Özgüner and Erçal1993a) both observe that the right-looking method of Geist and Ng (Reference Geist and Ng1989) sends more messages than both the distributed multifrontal method and the left-looking method with aggregated updates.

Jess and Kees (Reference Jess and Kees1982) describe a parallel right-looking LU factorization method for matrices with symmetric non-zero pattern, and with no pivoting, so their method can also be viewed as a way to compute the right-looking Cholesky factorization. Each node $k$ of the elimination tree corresponds to cdiv($k$) followed by cmod($i,k$) for all non-zeros $l_{ik}$. The tree describes the parallelism, since nodes that do not have an ancestor/descendant relationship can be computed in parallel. This assumes that multiple updates to the same column are handled correctly. The updates cmod($i,k_{1}$) and cmod($i,k_{2}$) can be computed in any order, but a race condition could occur if two processes attempt to update the same target column $i$ at the same time. Jess and Kees (Reference Jess and Kees1982) use a critical-path scheduling strategy, with the simplifying assumption that each node takes the same amount of time to compute.

Manne and Haffsteinsson (Reference Manne and Haffsteinsson1995) describe a right-looking method for a SIMD parallel computer with a two-dimensional processor mesh. In a SIMD machine, all processors perform the same operation in lock-step fashion. With the correct mapping scheme, this strategy simplifies the synchronization between processors, so that the problem of simultaneous updates to the same entry is resolved. They rely on a two-dimensional mapping scheme that assigns all entries in a column of $L$ to one column of processors and all entries in a row of $L$ to a single row of processors. The same mapping is used for both rows and columns. Thus, $a_{ij}$ is given to the processor in row $M(i)$ and column $M(j)$ of the two-dimensional processor mesh. Since each entry is owned by a single processor, multiple updates to the same entry are done in successive steps. Each outer loop of the chol_right algorithm above is done sequentially. Within the $k$th iteration, each processor iterates over the cmod operations it must compute. Note that a single processor only does part of any given cmod, since it does not own an entire column of the matrix.

Aggregating updates and sending a single message instead of one message is a common theme for many algorithms presented here. Ashcraft et al. (Reference Ashcraft, Eisenstat and Liu1990a) use this strategy for the left-looking method. The multifrontal method handles the submatrix update in an entirely different manner by delaying the updates and aggregating them. Hulbert and Zmijewski (Reference Hulbert and Zmijewski1991) consider a non-multifrontal right-looking method that also aggregates messages. Their method is based on the hypercube algorithm of George et al. (Reference George, Heath, Liu and Ng1989a) and Geist and Ng (Reference Geist and Ng1989). After a processor computes cdiv($k$) for a column $k$ that it owns, cmod($j,k$) must be computed for each non-zero $l_{jk}$ in column $k$. Some of these will be owned by the same processor, and these computations can be done right away. Others columns are owned by other processors. The method of George et al. (Reference George, Heath, Liu and Ng1989a) and Geist and Ng (Reference Geist and Ng1989) would send column $k$ to those other processors right away. Hulbert and Zmijewski (Reference Hulbert and Zmijewski1991) take another approach. The algorithm operates in two distinct phases. In the first phase, the method places column $k$ in an outgoing message queue, waits until all local computations that can be completed have finished. The queue acquires multiple messages for any given target processor. If there are multiple messages from this processor to the same target column, these are summed into a single aggregate column update (just as in the left-looking method of Ashcraft et al. Reference Ashcraft, Eisenstat and Liu1990a). As soon as the process runs out of local work it enters the second phase, which is identical to the original method of George et al. (Reference George, Heath, Liu and Ng1989a) and Geist and Ng (Reference Geist and Ng1989). In the second phase, the outgoing message queue is no longer used, and updates are sent as soon as they are computed. Hulbert and Zmijewski (Reference Hulbert and Zmijewski1991) show that this strategy results in a significant reduction in message traffic.

Ashcraft (Reference Ashcraft, George, Gilbert and Liu1993) generalizes the notion of left/right-looking and fan-in/out. Each non-zero $l_{ij}$ defines a cmod($i,j$). Recall that in the left-looking/fan-in method, cmod($i,j$) is done at step $i$ and in a right-looking/fan-out method, cmod($i,j$) is done at step $j$. There is no reason to use all one method or the other; each cmod($i,j$) can be individually assigned to task $i$ or task $j$. In Ashcraft’s (Reference Ashcraft, George, Gilbert and Liu1993) fan-both method the processors are aligned in a two-dimensional grid. The processor grid tiles the matrix, and if a processor owns the diagonal entry $l_{kk}$, then it owns column $k$. If it also owns $l_{kj}$, then cmod($k,j$) is treated in a left-looking (fan-in) fashion, and if it owns $l_{ik}$, then cmod($j,k$) is treated in a right-looking (fan-out) fashion.

We consider parallel supernodal, frontal, and multifrontal Cholesky factorization methods in Sections 9–11.

6.1 Symbolic analysis for LU factorization

If $A$ is square and symmetric, and no numerical pivoting is required, then the non-zero pattern of $L+U$ is the same as the Cholesky factorization of $A$. This observation provides the framework for several LU factorization algorithms presented here. Other methods consider the symbolic analysis of a matrix with unsymmetric structure, but with no pivoting. If arbitrary row interchanges can occur due to partial pivoting, then LU factorization is more closely related to QR factorization, a fact that other methods rely on. Both strategies are discussed immediately below, in this section. Other methods allow for arbitrary row and column pivoting during numeric factorization, and for such methods a prior symbolic analysis is not possible at all. Many of the results in this section for sparse LU (and also QR factorization) are surveyed by Gilbert and Ng (Reference Gilbert, Ng, George, Gilbert and Liu1993) and Gilbert (Reference Gilbert1994). Pothen and Toledo (Reference Pothen, Toledo, Mehta and Sahni2004) consider both the symmetric and unsymmetric cases in their recent survey of graph models of sparse elimination.

6.1.2 Symbolic LU with pivoting and its relationship to QR factorization

Consider both $LU=PA$ and $QR=A$, where $P$ is determined by partial pivoting. George and Ng (Reference George and Ng1985b, Reference George and Ng1987) have shown that $R$ is an upper bound on the pattern of $U$. More precisely, $u_{ij}$ can be non-zero if and only if $r_{ij}\neq 0$. Gilbert and Ng (Reference Gilbert, Ng, George, Gilbert and Liu1993) and Gilbert and Grigori (Reference Gilbert and Grigori2003) strengthened this result, showing that the bound is tight if $A$ is strong-Hall. A matrix is strong-Hall if it cannot be permuted into block upper triangular form with more than one block (Section 8.7.3). This upper bound is tight in a one-at-a-time sense; for any $r_{ij}\neq 0$, there exists an assignment of numerical values to entries in the pattern of $A$ that makes $u_{ij}\neq 0$. The outline of the proof can be seen by comparing Gaussian elimination with Householder reflections. Additional details are given in the qr_right_householder function discussed in Section 7.

Both LU and QR factorization methods eliminate entries below the diagonal. For a Householder reflection, George, Liu and Ng (Reference George, Liu and Ng1988b) show that the non-zero pattern of all rows affected by the transformation take on a non-zero pattern that is the union of all of these rows. With partial pivoting and row interchanges, these rows are candidate pivot rows for LU factorization. Only one of them is selected as the pivot row. Every other candidate pivot row is modified by adding to it a scaled copy of the pivot row. Thus, an upper bound on the pivot row pattern is the union of all candidate pivot rows. This also establishes a bound on $L$, namely, the non-zero pattern of $V$, which is a matrix whose column vectors correspond to the Householder reflections used for QR factorization.

With this relationship, a symbolic QR ordering and analysis (Section 7.1) becomes one possible symbolic ordering and analysis method for LU factorization. It is also possible to statically preallocate space for $L$ and $U$. The bound can be loose, however. In particular, if the matrix is diagonally dominant, then no pivoting is needed to maintain numerical accuracy. This is called static pivoting, where all pivoting is done prior to numeric factorization. If the matrix $A$ also has a symmetric non-zero pattern (or if all entries in the pattern of $A+A^{T}$ are considered to be ‘non-zero’), then the non-zero patterns of $L$ and $U$ are identical to the pattern of the Cholesky factors $L$ and $L^{T}$, respectively, of a symmetric positive definite matrix with the same non-zero pattern as $A+A^{T}$. In this case, a symmetric fill-reducing ordering of $A+A^{T}$ is appropriate. Alternatively, the permutation matrix $Q$ can be selected to reduce the worst-case fill-in for $PAQ=LU$ for any $P$, and then the permutation $P$ can be selected solely on the basis of partial pivoting, with no regard for sparsity.

Thus, LU factorization can rely on three basic strategies for finding a fill-reducing ordering. Two of them are methods used prior to factorization: (1) a symmetric pre-ordering of $A+A^{T}$, and (2) a column pre-ordering suitable for QR factorization. These options are discussed in Section 8. The third option is to dynamically choose pivots during numeric factorization, as discussed in Section 6.4.

With the QR upper bound, LU factorization can proceed using a statically allocated memory space. This bound can be quite high, however. It is sometimes better just to make a guess at the final $|L|$ and $|U|$, or to guess that no partial pivoting will be needed and to use a symbolic Cholesky analysis to determine a guess for $|L|$ and $|U|$. Sometimes a good guess is available from the LU factorization of a similar matrix in the same application. The only penalty for making a wrong guess is that the memory space for $|L|$ or $|U|$ must be reallocated if the guess is too low, or memory may run out if the guess is too high.

In contrast to Rose and Tarjan (Reference Rose and Tarjan1978) and Eisenstat and Liu (Reference Eisenstat and Liu1992), George and Ng (Reference George and Ng1987) consider partial pivoting. They rely on their result that QR forms an upper bound for LU to create a symbolic factorization method for both QR and LU. The resulting non-zero patterns for $L$ and $U$ can accommodate any partial pivoting with row interchanges. The symbolic factorization takes $O(|L|+|U|)$ time, which is much faster than Rose and Tarjan’s (Reference Rose and Tarjan1978) method. Their method (Algorithm 6.1 below) is much like the row-merge QR factorization of George and Heath (Reference George and Heath1980), which we discuss in Section 7.2. In this algorithm, ${\mathcal{L}}_{k}$ is set of row indices that is the upper bound on the pattern of the $k$th column of $L$, ${\mathcal{U}}_{k}$ is the upper bound on the $k$th row of $U$, and ${\mathcal{A}}_{i}$ is the pattern of the $i$th row of $A$.

At the $k$th step, partial pivoting can select any row whose leftmost non-zero falls in column $k$. Thus, ${\mathcal{U}}_{k}$ is the union of all such candidate pivot rows. This pivot row causes fill-in in all other candidate rows, which becomes the upper bound ${\mathcal{L}}_{k}$. Since these rows now all have the same non-zero pattern, their patterns are discarded and replaced with ${\mathcal{U}}_{k}$ itself. Thus, untouched rows of $A$ need only be considered once, and future steps need only look at the pattern of ${\mathcal{U}}_{k}$. The next time this row (representing a set of candidate pivot rows, ${\mathcal{L}}_{k}\setminus \{k\}$) is considered is at the step $p$ corresponding to the next non-zero entry in ${\mathcal{U}}_{k}$; namely, $p=\min {\mathcal{U}}_{k}\setminus \{k\}$. This is the first off-diagonal entry in the $k$th row of $U$. At that step $p$, the pattern ${\mathcal{U}}_{k}\setminus \{k\}$ is the upper bound of one (or more) unselected candidate pivot rows. The resulting algorithm takes time proportional to the upper bound on the LU factors, $O(|L|+|U|)$.

It should be noted that a single dense row destroys sparsity, causing this upper bound to become an entirely dense matrix. Such rows can be optimistically withheld, in an ad hoc manner, and placed as the last pivot rows. In this case, arbitrary partial pivoting is no longer possible.

Assuming the matrix $A$ is strong-Hall, the column elimination tree for the LU factorization of $A$ is given by the expression $p=\min {\mathcal{U}}_{k}\setminus \{k\}$ in Algorithm 6.1 above, where $p$ is the parent of $k$, the first off-diagonal entry in the $k$th row of $U$. It is identical to the elimination tree of $A^{T}A$ in this case. If $A$ is not strong-Hall, the tree is not the same, and it is referred to as the row-merge tree instead. Grigori, Cosnard and Ng (Reference Grigori, Cosnard and Ng2007a) provide a general characterization of this row-merge tree and its properties, and describe an efficient algorithm for computing it. Grigori, Gilbert and Cosnard (Reference Grigori, Gilbert and Cosnard2009) consider cases where numerical cancellation can result in sparser LU factors, and they show that if numerical cancellation is ignored the row-merge tree provides a tight bound on the structure of $L$ and $U$.

Algorithm 6.1 provides an upper bound on the QR factorization of $A$, as discussed in Section 7.1, which also gives an example of a matrix and its QR and LU factorizations as computed by this algorithm.

Gilbert and Liu (Reference Gilbert and Liu1993) generalize the elimination tree for a symmetric matrix to a pair of directed acyclic graphs (elimination DAGs, or EDAGs) for the LU factorization of an unsymmetric matrix (without pivoting). In LU factorization, the $k$th row of $L$ can be constructed via a sparse triangular solve using the first $k-1$ columns of $U$, and the $k$th column of $U$ arises from a triangular solve with the first $k-1$ rows of $L$. This comes from an unsymmetric analogue of the up-looking Cholesky factorization, namely,

(6.1)$$\begin{eqnarray}\left[\begin{array}{@{}cc@{}}L_{11} & \\ l_{21} & 1\end{array}\right]\left[\begin{array}{@{}cc@{}}U_{11} & u_{12}\\ & u_{22}\end{array}\right]=\left[\begin{array}{@{}cc@{}}A_{11} & a_{12}\\ a_{21} & a_{22}\end{array}\right]\!,\end{eqnarray}$$

where all leading submatrices are $(n-1)\times (n-1)$. Then the LU factorization can be computed with (1) $L_{11}U_{11}=A_{11}$ (a recursive LU factorization of the leading submatrix of $A$), (2) $L_{11}u_{12}=a_{12}$ (a sparse triangular solve for $u_{12}$), (3) $U_{11}^{T}l_{21}^{T}=a_{21}^{T}$ (a sparse triangular solve for $l_{21}$), and a dot product to compute $u_{22}$. As a result, the non-zero pattern of the $k$th row and $L$ and $k$th column of $U$ can be found as the reach in the acyclic graphs of $U$ and $L$, respectively, using equation (3.1). Basing a symbolic factorization on this strategy would result in a method taking the same time as the numeric factorization. Gilbert and Liu (Reference Gilbert and Liu1993) show how the graphs can be pruned via transitive reduction, giving the EDAGs (one for $L$ and another for $U$). For a symmetric matrix, the EDAGs are the same as the elimination tree of $A$. The transitive reduction of a graph preserves all paths (if $i{\rightsquigarrow}j$ is a path in a DAG, then a path still exists in the transitively reduced DAG). Gilbert and Liu (Reference Gilbert and Liu1993) characterize the row and column structures of $L$ and $U$ with these EDAGs, and present a left-looking symbolic factorization algorithm that constructs ${\mathcal{L}}$ and ${\mathcal{U}}$, building the EDAGs one node at a time. The EDAGs are more highly pruned graphs than the symmetric reductions of Eisenstat and Liu (Reference Eisenstat and Liu1992), but they take more time to compute.

Eisenstat and Liu (Reference Eisenstat and Liu2005a) provide a single tree that takes the place of the two EDAGs of Gilbert and Liu (Reference Gilbert and Liu1993). The symmetric elimination tree (a forest, actually) is given by edges in the Cholesky factor $L$: $i$ is the parent of $k$ if $i$ is the first off-diagonal non-zero in column $k$ of $L$ ($\min \{i>k|l_{ik}\neq 0\}$). By contrast, Eisenstat and Liu’s single tree for LU (also technically a forest) is defined in terms of paths in $L$ and $U$ instead of edges in $L$. In this tree, $i$ is the parent of $k$ if $i>k$ is the smallest node for which there is a path in $L$ form $i$ to $k$, and also a path from $k$ back to $i$ in the graph of $U$. Analogous to the $k$th row subtree for a symmetric matrix, they characterize the non-zero patterns of the rows of $L$ and columns of $U$ in terms of sub-forests of this path-based elimination tree. In a sequel to this paper (Eisenstat and Liu Reference Eisenstat and Liu2008), they present an algorithm for constructing this path-based tree/forest, and show how it characterizes the graph of $L+U$ in a recursive bordered block triangular form.

Gilbert et al. (Reference Gilbert, Li, Ng and Peyton2001) describe an algorithm that computes the row and column counts for sparse QR and LU factorization, as an extension of their prior work (Gilbert et al. Reference Gilbert, Ng and Peyton1994). Details are given in Section 4.2.

Grigori et al. (Reference Grigori, Demmel and Li2007b) present a parallel symbolic LU factorization method, based on a left-looking approach of Gilbert and Liu (Reference Gilbert and Liu1993), which computes the $k$th column of $L$ and the $k$th row of $U$ at the $k$th step. Their parallel algorithm generates a symbolic structure for $L$ and $U$ that can accommodate arbitrary partial pivoting. It starts with a vertex separator of the graph of $A+A^{T}$ to determine what parts of the symbolic factorization can be done in parallel. Vertex separators are an important tool for sparse direct methods and form the basis of the nested dissection ordering method discussed in Section 8.6.

6.2 Left-looking LU

The left-looking LU factorization algorithm computes $L$ and $U$ one column at a time. At the $k$th step, it accesses columns 1 to $k-1$ of $L$ and column $k$ of $A$. If partial pivoting is ignored, it can be derived from the following 3 $\times$ 3 block matrix expression, which is very similar to (5.1) for the left-looking Cholesky factorization algorithm. The matrix $L$ is assumed to have a unit diagonal:

(6.2)$$\begin{eqnarray}\left[\begin{array}{@{}ccc@{}}L_{11} & & \\ l_{21} & 1 & \\ L_{31} & l_{32} & L_{33}\end{array}\right]\left[\begin{array}{@{}ccc@{}}U_{11} & u_{12} & U_{13}\\ & u_{22} & u_{23}\\ & & U_{33}\end{array}\right]=\left[\begin{array}{@{}ccc@{}}A_{11} & a_{12} & A_{13}\\ a_{21} & a_{22} & a_{23}\\ A_{31} & a_{32} & A_{33}\end{array}\right]\!.\end{eqnarray}$$

The middle row and column of each matrix is the $k$th row and column of $L$, $U$, and $A$, respectively. If the first $k-1$ columns of $L$ and $U$ are known, three equations can be used to derive the $k$th columns of $L$ and $U$: $L_{11}u_{12}=a_{12}$ is a triangular system that can be solved for $u_{12}$ (the $k$th column of $U$), $l_{21}u_{12}+u_{22}=a_{22}$ can be solved for the pivot entry $u_{22}$, and $L_{31}u_{12}+l_{32}u_{22}=a_{32}$ can then be solved for $l_{32}$ (the $k$th column of $L$). However, these three equations can be rearranged so that nearly all of them are given by the solution to a single triangular system:

(6.3)$$\begin{eqnarray}\left[\begin{array}{@{}ccc@{}}L_{11} & & \\ l_{21} & 1 & \\ L_{31} & 0 & I\end{array}\right]\left[\begin{array}{@{}c@{}}x_{1}\\ x_{2}\\ x_{3}\end{array}\right]=\left[\begin{array}{@{}c@{}}a_{12}\\ a_{22}\\ a_{32}\end{array}\right]\!.\end{eqnarray}$$

The solution to this system gives $u_{12}=x_{1}$, $u_{22}=x_{2}$, and $l_{32}=x_{3}/u_{22}$. Partial pivoting with row interchanges is simple with this method. Once $x$ is found, entries in rows $k$ to $n$ can be searched for the entry with the largest magnitude. Permutation of the indices of $L$ is delayed until the end of the factorization. The non-zero patterns of the candidate pivot rows are not available (this would require a right-looking method) and thus the pivot row cannot be chosen for its sparsity. Fill-reducing orderings must instead be applied to the columns of $A$ only, as discussed in Section 8.5.

Relying on their optimal sparse triangular solve (Section 3.2), Gilbert and Peierls (Reference Gilbert and Peierls1988) show that their left-looking method takes time proportional to the number of floating-point operations. Other LU factorization methods can be faster in practice, but no other method provides this guarantee. This may seem like an obvious goal, but it can be quite difficult to achieve; it can take more time to search for entries and modify the data structure for $L$ and $U$ than the floating-point work to compute them. This method was the first sparse LU factorization in MATLAB (Gilbert et al. Reference Gilbert, Moler and Schreiber1992). Recent versions of MATLAB no longer use it for , but it is still relied upon for the [L,U,P]=lu(A) syntax when A is sparse. When a column ordering is required, [L,U,P,Q]=lu(A) relies on UMFPACK instead (Davis and Duff Reference Davis and Duff1997, Reference Davis and Duff1999, Davis Reference Davis2004a, Reference Davis2004b, discussed in Section 11.4).

The earliest left-looking method by Sato and Tinney (Reference Sato and Tinney1963) did not accommodate any pivoting. Dembart and Neves (Reference Dembart, Neves and Anderson1977) show how the left-looking method can be implemented on a vector machine with hardware gather/scatter operations. Their method does take time proportional to the floating-point work, but only because they rely on a precomputed sparsity pattern of $L$ and $U$. This restriction was lifted by Gilbert and Peierls’ (Reference Gilbert and Peierls1988) method. The first method to allow for partial pivoting was NSPIV by Sherman (Reference Sherman1978a, Reference Sherman1978b). It relied on a dynamic data structure with set unions performed as a merge operation, and as a result it could take more time than the floating-point work required. Sadayappan and Visvanathan (Reference Sadayappan and Visvanathan1988, Reference Sadayappan and Visvanathan1989) consider a parallel left-looking LU factorization method for circuit simulation matrices that does not allow for pivoting.

The majority of the methods described here (Sato and Tinney Reference Sato and Tinney1963, Sherman Reference Sherman1978a, Sherman Reference Sherman1978b, Sadayappan and Visvanathan Reference Sadayappan and Visvanathan1988, Eisenstat and Liu Reference Eisenstat and Liu1993a) actually store $L$ and $U$ by rows and compute them one row at a time, but this is identical to the left-looking method applied to $A^{T}$.

Eisenstat and Liu (Reference Eisenstat and Liu1993a) reduce the work that Gilbert and Peierls’ (Reference Gilbert and Peierls1988) method requires for computing the reach in the graph of $L$ when finding the pattern of $x$ (the $k$th column of $L$ and $U$) in the sparse triangular solve, ${\mathcal{X}}=\text{Reach}_{L}({\mathcal{B}})$. They observe that the reach of a node in the graph is unaffected if edges are pruned. For each $k$, let $i$ be the smallest index such that both $l_{ik}$ and $u_{ki}$ are non-zero. This entry forms a symmetric pair, or FSNZ($k$) (Eisenstat and Liu Reference Eisenstat and Liu1992). Any entries below this in $L$ can be pruned from the graph of $L$, and the reach is unaffected. In the best case when the non-zero pattern of $L$ and $U$ is symmetric, this pruning results in the elimination tree, and the time to find the pattern ${\mathcal{X}}$ reduced to $O(|{\mathcal{X}}|)$.

Davis (Reference Davis2006) provides an implementation of the left-looking method of Gilbert and Peierls (Reference Gilbert and Peierls1988) in the CSparse package. It also forms the basis of KLU (Davis and Palamadai Natarajan Reference Davis and Palamadai Natarajan2010), a solver targeted for circuit simulation matrices. These matrices are too sparse for methods based on dense submatrices (supernodal, frontal, and multifrontal) to be efficient.

Gustavson, Liniger and Willoughby (Reference Gustavson, Liniger and Willoughby1970) and Hachtel, Brayton and Gustavson (Reference Hachtel, Brayton and Gustavson1971) present an alternative method for sparse LU factorization. Their symbolic analysis produces not only the non-zero patterns of the $L$ and $U$, but also a loop-free code, with a sequence of operations that factorizes the matrix and is specific to its particular non-zero pattern. Norin and Pottle (Reference Norin and Pottle1971) consider fill-reducing orderings for this method. The method of generating loop-free code requires significant memory for the compiled code (proportional to the number of floating-point operations), which Gay (Reference Gay1991) shows is not required for obtaining good performance.

Chen, Wang and Yang (Reference Chen, Wang and Yang2013) present a multicore algorithm, NICSLU, based on the left-looking sparse LU. It accommodates partial pivoting during numerical factorization, and relies on the column elimination tree discussed in Section 6.1 for its parallel scheduling. Each task consists of one node in this tree, corresponding to the computation of a single column of $L$ and $U$. The first phase handles nodes towards the bottom of this tree, one level at a time. In the second phase, each task updates its column with any prior columns that affect it and which have already completed. Prior columns not yet finished are skipped in the first pass of this task, and then handled in a second pass after they are complete.

6.3 Right-looking LU factorization with a static data structure

Gaussian elimination is a right-looking variant of LU factorization. At each step, an outer product of the pivot column and the pivot row is subtracted from the lower right submatrix of $A$. After the $k$th step, the lower right submatrix $A^{[k]}$ is a Schur complement of the upper left $k\times k$ submatrix, also called the active submatrix. Numerical pivoting is typically essential, but ignoring it for the moment simplifies the derivation. The derivation of the method starts with an equation very similar to (4.2) for the right-looking Cholesky factorization,

(6.4)$$\begin{eqnarray}\left[\begin{array}{@{}cc@{}}l_{11} & \\ l_{21} & L_{22}\end{array}\right]\left[\begin{array}{@{}cc@{}}u_{11} & u_{12}\\ & U_{22}\end{array}\right]=\left[\begin{array}{@{}cc@{}}a_{11} & a_{12}\\ a_{21} & A_{22}\end{array}\right]\!,\end{eqnarray}$$

where $l_{11}=1$ is a scalar, and all three matrices are square and partitioned identically. Other choices for $l_{11}$ are possible; this choice leads to a unit lower triangular $L$ and the four equations

(6.5)$$\begin{eqnarray}\displaystyle u_{11} & = & \displaystyle a_{11},\end{eqnarray}$$

(6.6)$$\begin{eqnarray}\displaystyle u_{12} & = & \displaystyle a_{12},\end{eqnarray}$$

(6.7)$$\begin{eqnarray}\displaystyle l_{21}u_{11} & = & \displaystyle a_{21},\end{eqnarray}$$

(6.8)$$\begin{eqnarray}\displaystyle l_{21}u_{12}+L_{22}U_{22} & = & \displaystyle A_{22}.\end{eqnarray}$$

Each equation is solved in turn, and can be expressed in MATLAB notation as the lu_right function below, where after the $k$ step, A(k+1:n,k+1:n) holds the $k$th Schur complement, $A^{[k]}$.

One advantage of the right-looking method over left-looking sparse LU factorization is that it can select a sparse pivot row and pivot column. The left-looking method does not keep track of the non-zero pattern of the $A^{[k]}$ submatrix, and thus cannot determine the number of non-zeros in its pivot rows or columns. With pivoting of both rows and columns for the dual purposes of maintaining sparsity and ensuring numerical accuracy, the resulting factorization is $LU=PAQ$, where $Q$ is the column permutation. The disadvantage of the right-looking method is that it is significantly more difficult to implement, particularly when pivoting is done during numerical factorization. This variant is discussed in Section 6.4.

If the pivoting is determined prior to numeric factorization, however, then a simpler static data structure can be used for $L$ and $U$. Eisenstat, Schultz and Sherman (Reference Eisenstat, Schultz, Sherman, Duff and Stewart1979) describe such a method that can be considered as a precursor to the multifrontal method (Section 11). It is a numeric form of their symbolic method (Eisenstat et al. Reference Eisenstat, Schultz, Sherman, Bunch and Rose1976a), which represents the symbolic graph elimination via a set of elements, or cliques that form during factorization. The method computes each update in a dense submatrix (much like the multifrontal method). To save space, it then discards the rows and columns of $L$ and $U$ inside these elements, except for the top-level separator. The rows and columns of $L$ and $U$ are then recalculated when needed.

With a static data structure, the right-looking method is more amenable to a parallel implementation, as compared to a dynamic data structure. Several of the earliest methods rely on pre-pivoting to enhance parallelism in a right-looking LU factorization by finding independent sets of pivots (Huang and Wing Reference Huang and Wing1979, Wing and Huang Reference Wing and Huang1980, Jess and Kees Reference Jess and Kees1982, Srinivas Reference Srinivas1983); none of these methods modify the pivot ordering during numeric factorization. In (6.4), the $a_{11}$ scalar becomes an $s\times s$ diagonal matrix. This ordering strategy is discussed in Section 8.8, while we discuss the numeric factorization here.

Huang and Wing (Reference Huang and Wing1979) and Wing and Huang (Reference Wing and Huang1980) analyse the data dependencies in a fine-grained parallel method where each floating-point operation is a single task, either a division by the pivot,

$$\begin{eqnarray}l_{ik}=\frac{a_{ik}^{[k]}}{u_{kk}},\end{eqnarray}$$

or an update to compute one entry in the Schur complement,

$$\begin{eqnarray}a_{ij}^{[k]}=a_{ij}^{[k-1]}-l_{ik}u_{kj}.\end{eqnarray}$$

Each computation is placed in a directed acyclic graph, where the edges represent the data dependencies in these two kinds of computations, and a scheduling method is presented based on a local greedy heuristic. Note that several of these papers include the word ‘optimal’ in their title. This is an incorrect use of the term, since finding an optimal pivot ordering and an optimal schedule are NP-hard problems. Srinivas (Reference Srinivas1983) refines Wing and Huang’s (Reference Wing and Huang1980) scheduling method to reduce the number of parallel steps required to factorize the matrix.

Jess and Kees (Reference Jess and Kees1982) introduced the term elimination tree for their parallel right-looking LU factorization algorithm, which assumes a symmetric non-zero pattern. Their definition of the tree was limited to a filled-in graph of $L+U$; this was later generalized by Schreiber (Reference Schreiber1982), who defined the elimination tree in the way it is currently used (Section 4.1). Jess and Kees used the tree for a coarse-grained parallel algorithm, where each node is a single step $k$ in the lu_right function above. Two nodes $k_{1}$ and $k_{2}$ can be executed in parallel if they do not share an ancestor/descendant relationship in the tree, and Jess and Kees define a method for scheduling the $n$ tasks. They note that two independent tasks can still require synchronization, however, since both can update the same set of entries in the lower right submatrix.

George and Ng (Reference George and Ng1985b) present a completely different approach for right-looking LU factorization by defining a static data structure for $L$ and $U$ that allows for arbitrary row interchanges to accommodate numerical pivoting during factorization. This is based on a symbolic QR factorization, and they show that the non-zero pattern of $R$ for the QR factorization ($QR=A$) is an upper bound on the pattern of $U$ and $L^{T}$. In a companion paper, George et al. (Reference George, Liu and Ng1988b) refine this method with a more compact data structure. They demonstrate and exploit the fact that the pattern of $L$ is contained in the non-zero pattern of the $n$ Householder vectors (concatenated to form the matrix $V$) for the QR factorization. Their numeric factorization is much like the symbolic right-looking LU factorization method of Algorithm 6.1 in Section 6.1. At each step, they maintain a set of candidate pivot rows for each step $k$, analogous to the sets ${\mathcal{S}}_{1..n}$ in the symbolic LU factorization Algorithm 6.1. In that algorithm, the set of candidate pivot rows is replaced by a single representative row, since they all have the same non-zero pattern. In the numeric factorization, however, each row must be included. Instead of ${\mathcal{S}}_{k}$, the numeric factorization uses ${\mathcal{Z}}_{k}$ as the set of candidate pivot rows for step $k$. Initially, ${\mathcal{Z}}_{k}$ holds all original rows of $A$ whose leftmost non-zero entry resides in column $k$. At the $k$th step, one pivot row is chosen from ${\mathcal{Z}}_{k}$, and the remainder are added into the parent set ${\mathcal{Z}}_{p}$, where $p$ is the parent of $k$ in the column elimination tree.

George and Ng (Reference George and Ng1990) construct a parallel algorithm for their method, suitable for a shared-memory system. In contrast to the other papers on parallel algorithms discussed in this section, they describe an actual implementation. Their method relies on the fact that the sets ${\mathcal{Z}}_{1:n}$ are always disjoint. Step $k$ of lu_right corresponds to node $k$ in the column elimination tree. If two nodes $a$ and $b$ do not have an ancestor/descendant relationship in the tree, and if all their descendants have been computed, then these two steps have no candidate pivot rows in common (${\mathcal{Z}}_{a}$ and ${\mathcal{Z}}_{b}$ are disjoint). Each step selects one of their candidate pivot rows and uses it to update the remaining rows in their own set. This work for nodes $a$ and $b$ can be done in parallel.

6.4 Right-looking LU with dynamic pivoting

Unlike the left-looking method, the right-looking LU factorization method has the Schur complement available to it, in the active submatrix $A^{[k]}$. It can thus select pivots from the entire active submatrix during factorization, using both sparsity and numerical criteria. For dense matrices, the complete pivoting strategy selects the entry in the active submatrix with the largest magnitude. The search for this pivot requires $O((n-k)^{2})$ work at the $k$th step, the same work as the subsequent numerical update.

Sparse row-oriented QR factorization with Givens rotations

George and Heath (Reference George and Heath1980) present the first row-oriented QR factorization, based on Givens rotations, and it can be simply described (hereafter referred to as the row-Givens method). Row-Givens starts with a symbolic analysis of the normal equations, since under a few simplifying assumptions, the non-zero pattern of $R$ is given by the Cholesky factor of $A^{T}A$. All of the methods for Cholesky symbolic analysis in Section 4 can thus be used on $A^{T}A$, and the pattern ${\mathcal{R}}$ of $R$ is known prior to numeric factorization.

The numeric factorization starts with an $R$ that is all zero, but with a known final non-zero pattern. It is stored row-by-row in a static data structure. The method selects a row of $A$ and finds its leftmost non-zero; suppose this entry is in column $k$. This entry is annihilated with a Givens rotation (a 2 $\times$ 2 orthogonal matrix $G$) that uses the $k$th row of $R$ to annihilate the $k$ entry of the selected row of $A$. Both rows are modified. The next leftmost entry in this row of $A$ is then found, and the process continues until the row of $A$ has been annihilated to zero. If the diagonal $r_{kk}$ is zero, the Givens rotation becomes a swap. This happens if the $k$th row of $R$ is all zero, and thus the process stops early.

The MATLAB rowgivens script below illustrates the row-Givens sparse QR factorization. To keep it simple, it does not consider the pattern of $R$ and stores R in dense format. However, it does consider the sparsity of the rows of A. The function implicitly stops early if it encounters a zero on the diagonal of R. The matrix $Q$ could be computed, or kept implicitly as a sequence of Givens rotations, but the script simply discards it.

Dense column-oriented QR factorization with Householder reflections

A Householder reflection is a symmetric orthogonal matrix of the form $H=I-{\it\beta}vv^{T}$, where ${\it\beta}$ is a scalar and $v$ is a column vector. The vector $v$ and scalar ${\it\beta}$ can be chosen based on a vector $x$ so that $Hx$ is zero except for the first entry $(Hx)_{1}=\pm \Vert x\Vert _{2}$. Computing the matrix–vector product $Hx$ takes only $O(n)$ work. The non-zero patterns of $v$ and $x$ are the same. The MATLAB script qr_right_householder uses Householder reflections in a right-looking manner to compute the QR factorization. It represents $Q$ implicitly as a sequence of Householder reflections (V and Beta) which could be discarded if they are not needed. Householder reflections can also be used for a left-looking QR factorization, assuming they are all kept.