Big Ideas in Applied Math: Sparse Matrices

Sparse matrices are an indispensable tool for anyone in computational science. I expect there are a very large number of simulation programs written in scientific research across the country which could be faster by ten to a hundred fold at least just by using sparse matrices! In this post, we’ll give a brief overview what a sparse matrix is and how we can use them to solve problems fast.

A matrix is sparse if most of its entries are zero. There is no precise threshold for what “most” means; Kolda suggests that a matrix have at least 90% of its entries be zero for it to be considered sparse. The number of nonzero entries in a sparse matrix $A$ is denoted by $\operatorname{nnz}(A)$ . A matrix that is not sparse is said to be dense.

Sparse matrices are truly everywhere. They occur in finite difference, finite element, and finite volume discretizations of partial differential equations. They occur in power systems. They occur in signal processing. They occur in social networks. They occur in intermediate stages in computations with dense rank-structured matrices. They occur in data analysis (along with their higher-order tensor cousins).

Why are sparse matrices so common? In a word, locality. If the $ij$ th entry $a_{ij}$ of a matrix $A$ is nonzero, then this means that row $i$ and column $j$ are related in some way to each other according to the the matrix $A$ . In many situations, a “thing” is only related to a handful of other “things”; in heat diffusion, for example, the temperature at a point may only depend on the temperatures of nearby points. Thus, if such a locality assumption holds, every row will only have a small number of nonzero entries and the matrix overall will be sparse.

Storing and Multiplying Sparse Matrices

A sparse matrix can be stored efficiently by only storing its nonzero entries, along with the row and column in which these entries occur. By doing this, a sparse matrix can be stored in $\mathcal{O}(\operatorname{nnz}(A))$ space rather than the standard $\mathcal{O}(N^2)$ for an $N\times N$ matrix $A$ .¹ For the efficiency of many algorithms, it will be very beneficial to store the entries row-by-row or column-by-column using compressed sparse row and column (CSR and CSC) formats; most established scientific programming software environments support sparse matrices stored in one or both of these formats. For efficiency, it is best to enumerate all of the nonzero entries for the entire sparse matrix and then form the sparse matrix using a compressed format all at once. Adding additional entries one at a time to a sparse matrix in a compressed format requires reshuffling the entire data structure for each new nonzero entry.

There exist straightforward algorithms to multiply a sparse matrix $A$ stored in a compressed format with a vector $x$ to compute the product $b = Ax$ . Initialize the vector $b$ to zero and iterate over the nonzero entries $a_{ij}$ of $A$ , each time adding $a_{ij}x_j$ to $b_i$ . It is easy to see this algorithm runs in $\mathcal{O}(\operatorname{nnz}(A))$ time.² The fact that sparse matrix-vector products can be computed quickly makes so-called Krylov subspace iterative methods popular for solving linear algebraic problems involving sparse matrices, as these techniques only interact with the matrix $A$ by computing matrix-vector products $x \mapsto Ax$ (or matrix-tranpose-vector products $y \mapsto A^\top y$ ).

Lest the reader think that every operation with a sparse matrix is necessarily fast, the product of two sparse matrices $A$ and $B$ need not be sparse and the time complexity need not be $\mathcal{O}(\operatorname{nnz}(A) + \operatorname{nnz}(B))$ . A counterexample is

(1) $\begin{equation*} A = B^\top = \begin{bmatrix} a_1 & 0 & \cdots & 0 \\ a_2 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & 0 \\ a_n & 0 & \cdots & 0 \end{bmatrix} \in \mathbb{R}^n \end{equation*}$

for $a_1,\ldots,a_n \ne 0$ . We have that $\operatorname{nnz}(A) = \operatorname{nnz}(B) = N$ but

(2) $\begin{equation*} AB = \begin{bmatrix} a_1^2 & a_1 a_2 & \cdots & a_1a_N \\ a_2a_1 & a_2^2 & \cdots & a_2a_N \\ \vdots & \vdots & \ddots & \vdots \\ a_Na_1 & a_Na_2 & \cdots & a_N^2 \end{bmatrix} \end{equation*}$

which has $\operatorname{nnz}(AB) = N^2$ nonzero elements and requires $\mathcal{O}(N^2)$ operations to compute. However, if one does the multiplication in the other order, one has $\operatorname{nnz}(BA) = 1$ and the multiplication can be done in $\mathcal{O}(N)$ operations. Thus, some sparse matrices can be multiplied fast and others can’t. This phenomena of different speeds for different sparse matrices is very much also true for solving sparse linear systems of equations.

Solving Sparse Linear Systems

The question of how to solve a sparse system of linear equations $Ax = b$ where $A$ is sparse is a very deep problems with fascinating connections to graph theory. For this article, we shall concern ourselves with so-called sparse direct methods, which solve $Ax = b$ by means of computing a factorization of the sparse matrix $A$ . These methods produce an exact solution to the system $Ax = b$ if all computations are performed exactly and are generally considered more robust than inexact and iterative methods. As we shall see, there are fundamental limits on the speed of certain sparse direct methods, which make iterative methods very appealing for some problems.

Note from the outset that our presentation will be on illustrating the big ideas rather than presenting the careful step-by-step details needed to actually code a sparse direct method yourself. An excellent reference for the latter is Tim Davis’ wonderful book Direct Methods for Sparse Linear Systems.

Let us begin by reviewing how $LU$ factorization works for general matrices. Suppose that the $(1,1)$ entry of $A$ is nonzero. Then, $LU$ factorization proceeds by subtracting scaled multiples of the first row from the other rows to zero out the first column. If one keeps track of these scaling, then one can write this process as a matrix factorization, which we may demonstrate pictorially as

(3) $\begin{equation*} \underbrace{\begin{bmatrix} * & * & \cdots & * \\ * & * & \cdots & * \\ \vdots & \vdots & \ddots & \vdots \\ * & * & \cdots & *\end{bmatrix}}_{=A} = \begin{bmatrix} 1 & & &\\ * & 1 & & \\ \vdots & & \ddots & \\ * & & & 1\end{bmatrix} \begin{bmatrix} * & * & \cdots & * \\ & * & \cdots & * \\ & \vdots & \ddots & \vdots \\ & * & \cdots & *\end{bmatrix}. \end{equation*}$

Here, $*$ ‘s denote nonzero entries and blanks denote zero entries. We then repeat the process on the $(N-1)\times (N-1)$ submatrix in the bottom right (the so-called Schur complement). Continuing in this way, we eventually end up with a complete $LU$ factorization

(4) $\begin{equation*} \underbrace{\begin{bmatrix} * & * & \cdots & * \\ * & * & \cdots & * \\ \vdots & \vdots & \ddots & \vdots \\ * & * & \cdots & *\end{bmatrix}}_{=A} = \underbrace{\begin{bmatrix} 1 & & &\\ * & 1 & & \\ \vdots & \vdots & \ddots & \\ * &* &\cdots & 1\end{bmatrix}}_{=L} \underbrace{\begin{bmatrix} * & * & \cdots & * \\ & * & \cdots & * \\ & & \ddots & \vdots \\ & & & *\end{bmatrix}}_{=U}. \end{equation*}$

In the case that $A$ is symmetric positive definite (SPD), one has that $U = DL^\top$ for $D$ a diagonal matrix consisting of the entries on $U$ . This factorization $A = LDL^\top$ is a Cholesky factorization of $A$ .³ For general non-SPD matrices, one needs to incorporate partial pivoting for Gaussian elimination to produce accurate results.⁴

Let’s try the same procedure for a sparse matrix. Consider a sparse matrix with the following sparsity pattern:

(5) $\begin{equation*} A = \begin{bmatrix} * & * & * & & * \\ * & * & & * & \\ * & & * & & \\ & * & & * & \\ * & & & & * \end{bmatrix}. \end{equation*}$

When we eliminate the $(1,1)$ entry, we get the following factorization:

(6) $\begin{equation*} \underbrace{\begin{bmatrix} * & * & * & & * \\ * & * & & * & \\ * & & * & & \\ & * & & * & \\ * & & & & * \end{bmatrix}}_{=A} = \begin{bmatrix} 1 & & & & \\ * & 1& & & \\ * & & 1 & & \\ & & & 1 & \\ * & & & & 1 \end{bmatrix} \begin{bmatrix} * & * & * & & * \\ & * & \bullet & * & \bullet \\ & \bullet & * & & \bullet \\ & * & & * & \\ & \bullet & \bullet & & * \end{bmatrix} \end{equation*}$

Note that the Schur complement has new additional nonzero entries (marked with a $\bullet$ ) not in the original sparse matrix $A$ . The Schur complement of $A$ is denser than $A$ was; there are new fill-in entries. The worst-case scenario for fill-in is the arrowhead matrix:

(7) $\begin{equation*} \underbrace{\begin{bmatrix} * & * & * & \cdots & * \\ * & * & & & \\ * & & * & & \\ \vdots & & & \ddots & \\ * & & & & * \end{bmatrix}}_{=A} = \begin{bmatrix} 1 & & & & \\ * & 1& & & \\ * & & 1 & & \\ \vdots & & & \ddots & \\ * & & & & 1 \end{bmatrix} \begin{bmatrix} * & * & * & * & * \\ & * & \bullet & \cdots & \bullet \\ & \bullet & * & \cdots & \bullet \\ & \vdots & \vdots & \ddots & \vdots \\ & \bullet & \bullet & \cdots & * \end{bmatrix} \end{equation*}$

After one step of Gaussian elimination, we went from a matrix with $\mathcal{O}(N)$ nonzeros to a fully dense Schur complement! However, the arrowhead matrix also demonstrates a promising strategy. Simply construct a permutation matrix which reorders the first entry to be the last⁵ and then perform Gaussian elimination on the symmetrically permuted matrix $PAP^\top$ instead. In fact, the entire $LU$ factorization can be computed without fill-in:

(8) $\begin{equation*} \underbrace{\begin{bmatrix} * & & & & * \\ & * & & & * \\ & & * & & * \\ & & & \ddots & \vdots \\ * & * & * & \cdots & * \end{bmatrix}}_{=PAP^\top} = \begin{bmatrix} 1 & & & & \\ & 1& & & \\ & & 1 & & \\ & & & \ddots & \\ * & * & * & \cdots & 1 \end{bmatrix} \begin{bmatrix} * & & & & * \\ & * & & & * \\ & & * & & * \\ & & & \ddots & \vdots \\ & & & & * \end{bmatrix}. \end{equation*}$

This example shows the tremendous importance of reordering of the rows and columns when computing a sparse $LU$ factorization.

The Best Reordering

As mentioned above, when computing an $LU$ factorization of a dense matrix, one generally has to reorder the rows (and/or columns) of the matrix to compute the solution accurately. Thus, when computing the $LU$ factorization of a sparse matrix, one has to balance the need to reorder for accuracy and to reorder to reduce fill-in. For these reasons, for the remainder of this post, we shall focus on computing Cholesky factorizations of SPD sparse matrices, where reordering for accuracy is not necessary.⁶ Since we want the matrix to remain SPD, we must restrict ourselves to symmetric reordering strategies where $A$ is reordered to $PAP^\top$ where $P$ is a permutation matrix.

Our question is deceptively simple: what reordering produces the least fill-in? In matrix language, what permutation $P$ minimizes $\operatorname{nnz}(L)$ where $LDL^\top = PAP^\top$ is the Cholesky factorization of $PAP^\top$ ?

Note that, assuming no entries in the Gaussian elimination process exactly cancel, then the Cholesky factorization depends only on the sparsity pattern of $A$ (the locations of the zeros and nonzeros) and not on the actual numeric values of $A$ ‘s entries. This sparsity structure is naturally represented by a graph $\mathcal{G}(A)$ whose nodes are the indices $\{1,\ldots,N\}$ with an edge between $i \ne j$ if, and only if, $a_{ij} \ne 0$ .

Now let’s see what happens when we do Gaussian elimination from a graph point-of-view. When we eliminate the $(1,1)$ entry from matrix, this results in all nodes of the graph adjacent to $1$ becoming connected to each other.⁷

This shows why the arrowhead example is so bad. By eliminating the a vertex connected to every node in the graph, the eliminated graph becomes a complete graph.

Reordering the matrix corresponds to choosing in what order the vertices of the graph are eliminated. Choosing the elimination order is then a puzzle game; eliminate all the vertices of the graph in the order that produces the fewest fill-in edges (shown red).⁸

Finding the best elimination ordering for a sparse matrix (graph) is a good news/bad news situation. For the good news, many graphs possess a perfect elimination ordering, in which no fill-in is produced at all. There is a simple algorithm to determine whether a graph (sparse matrix) possesses a perfect elimination ordering and if so, what it is.⁹ Some important classes of graphs can be eliminated perfectly (for instance, trees). More generally, the class of all graphs which can be eliminated perfectly is precisely the set of chordal graphs, which are well-studied in graph theory.

Now for the bad news. The problem of finding the best elimination ordering (with the least fill-in) for a graph is NP-Hard. This means, assuming the widely conjectured result that ${\rm P} \ne {\rm NP}$ , that finding the best elimination ordering would be a hard computational problem than the worst-case $\mathcal{O}(N^3)$ complexity for doing Gaussian elimination in any ordering! One should not be too pessimistic about this result, however, since (assuming ${\rm P} \ne {\rm NP}$ ) all it says is that there exists no polynomial time algorithm guaranteed to produce the absolutely best possible elimination ordering when presented with any graph (sparse matrix). If one is willing to give up on any one of the bolded statements, further progress may be possible. For instance, there exists several good heuristics, which find reasonably good elimination orderings for graphs (sparse matrices) in linear $\mathcal{O}(\operatorname{nnz}(A))$ (or nearly linear) time.

Can Sparse Matrices be Eliminated in Linear Time?

Let us think about the best reordering question in a different way. So far, we have asked the question “Can we find the best ordering for a sparse matrix?” But another question is equally important: “How efficiently can we solve a sparse matrix, even with the best possible ordering?”

One might optimistically hope that every sparse matrix possesses an elimination ordering such that its Cholesky factorization can be computed in linear time (in the number of nonzeros), meaning that the amount of time needed to solve $Ax = b$ is proportional to the amount of data needed to store the sparse matrix $A$ .

When one tests a proposition like this, one should consider the extreme cases. If the matrix $A$ is dense, then it requires $\mathcal{O}(N^3)$ operations to do Gaussian elimination,¹⁰ but $A$ only has $\operatorname{nnz}(A) = N^2$ nonzero entries. Thus, our proposition cannot hold in unmodified form.

An even more concerning counterexample is given by a matrix $A$ whose graph $\mathcal{G}(A)$ is a $\mathcal{O}(\sqrt{N}) \times \mathcal{O}(\sqrt{N})$ 2D grid graph.

Sparse matrices with this sparsity pattern (or related ones) appear all the time in discretized partial differential equations in two dimensions. Moreover, they are truly sparse, only having $\operatorname{nnz}(A) = \mathcal{O}(N)$ nonzero entries. Unforunately, no linear time elimination ordering exists. We have the following theorem:

Theorem: For any elimination ordering for a sparse matrix $A$ with $\mathcal{G}(A)$ being a $\sqrt{N}\times \sqrt{N}$ 2D grid graph, in any elimination ordering, the Cholesky factorization $PAP^\top = LDL^\top$ requires $\Omega(N^{3/2})$ operations and satisfies $\operatorname{nnz}(L)= \Omega(N\log N)$ .¹¹

The proof is contained in Theorem 10 and 11 (and the ensuing paragraph) of classic paper by Lipton, Rose, and Tarjan. Natural generalizations to $d$ -dimensional grid graphs give bounds of $\Omega(N^{3(d-1)/d})$ time and $\operatorname{nnz}(L)=\Omega(N^{2(d-1)/d})$ for $d > 2$ . In particular, for 2D finite difference and finite element discretizations, sparse Cholesky factorization takes $\Omega(N^{3/2})$ operations and produces a Cholesky factor with $\operatorname{nnz}(L)= \Omega(N\log N)$ in the best possible ordering. In 3D, sparse Cholesky factorization takes $\Omega(N^{2})$ operations and produces a Cholesky factor with $\operatorname{nnz}(L)= \Omega(N^{4/3})$ in the best possible ordering.

Fortunately, at least these complexity bounds are attainable: there is an ordering which produces a sparse Cholesky factorization with $PAP^\top = LDL^\top$ requiring $\Theta(N^{3/2})$ operations and with $\operatorname{nnz}(L)= \Theta(N\log N)$ nonzero entries in the Cholesky factor.¹² One such asymptotically optimal ordering is the nested dissection ordering, one of the heuristics alluded to in the previous section. The nested dissection ordering proceeds as follows:

Find a separator $S$ consisting of a small number of vertices in the graph $\mathcal{G}(A)$ such that when $S$ is removed from $\mathcal{G}(A)$ , $\mathcal{G}(A)$ is broken into a small number of edge-disjoint and roughly evenly sized pieces $\mathcal{G}_1,\ldots,\mathcal{G}_k$ .¹³
Recursively use nested dissection to eliminate each component $\mathcal{G}_1,\ldots, \mathcal{G}_k$ individually.
Eliminate $S$ in any order.

For example, for the 2D grid graph, if we choose $S$ to be a cross through the center of the 2D grid graph, we have a separator of size $|S| = \Theta(\sqrt{N})$ dividing the graph into $4$ roughly $\sqrt{N}/2\times \sqrt{N}/2$ pieces.

Let us give a brief analysis of this nested dissection ordering. First, consider the sparsity of the Cholesky factor $\operatorname{nnz}(L)$ . Let $S(N)$ denote the number of nonzeros in $L$ for an elimination of the $\sqrt{N}\times \sqrt{N}$ 2D grid graph using the nested dissection ordering. Then step 2 of nested dissection requires us to recursively eliminate four $\sqrt{N/4} \times \sqrt{N/4}$ 2D grid graphs. After doing this, for step 3, all of the vertices of the separator might be connected to each other, so the separator graph will potentially have as many as $\mathcal{O}(|S|^2) = \mathcal{O}(N)$ edges, which result in nonzero entries in $L$ . Thus, combining the fill-in from both steps, we get

(9) $\begin{equation*} S(N) = 4S\left(\frac{N}{4}\right) + \mathcal{O}(N). \end{equation*}$

Solving this recurrence using the master theorem for recurrences gives $\operatorname{nnz}(L) = S(N) = \mathcal{O}(N\log N)$ . If one instead wants the time $T(N)$ required to compute the Cholesky factorization, note that for step 3, in the worst case, all of the vertices of the separator might be connected to each other, leading to a $\sqrt{N}\times \sqrt{N}$ dense matrix. Since a $\sqrt{N}\times \sqrt{N}$ matrix requires $\mathcal{O}((\sqrt{N})^3) = \mathcal{O}(N^{3/2})$ , we get the recurrence

(10) $\begin{equation*} T(N) = 4T\left(\frac{N}{4}\right) + \mathcal{O}(N^{3/2}), \end{equation*}$

which solves to $T(N) = \mathcal{O}(N^{3/2})$ .

Conclusions

As we’ve seen, sparse direct methods (as exemplified here by sparse Cholesky) possess fundamental scalability challenges for solving large problems. For the important class of 2D and 3D discretized partial differential equations, the time to solve $Ax = b$ scales like $\Theta(N^{3/2})$ and $\Theta(N^2)$ , respectively. For truly large-scale problems, these limitations may be prohibitive for using such methods.

This really is the beginning of the story, not the end for sparse matrices however. The scalability challenges for classical sparse direct methods has spawned many exciting different approaches, each of which combats the scalability challenges of sparse direct methods for a different class of sparse matrices in a different way.

Upshot: Sparse matrices occur everywhere in applied mathematics, and many operations on them can be done very fast. However, the speed of computing an $LU$ factorization of a sparse matrix depends significantly on the arrangement of its nonzero entries. Many sparse matrices can be factored quickly, but some require significant time to factor in any reordering.

4 thoughts on “Big Ideas in Applied Math: Sparse Matrices”

Pingback: Big Ideas in Applied Math: The Fast Fourier Transform – Ethan Epperly
Pingback: Big Ideas in Applied Math: Concentration Inequalities – Ethan Epperly
Pingback: Does Sketching Work? – Ethan N. Epperly
Pingback: Which Sketch Should I Use? – Ethan N. Epperly