Algorithms on Sparse Matrices

In general, the "best" representation of a sparse matrix depends on what operations are to be done with it. For iterative methods, the dominant operation is ...

Sparse matrix-vector multiplication

In the discussion of each data structure shown, sparse matrix-vector products are assumed to be the prototypical operation, because that is the most common operation in iterative methods for solving linear systems of equations. The operation is multiplying a sparse matrix times a dense vector:

w = A * d

where w and d are n×1 vectors. Recall that the linear system being solved is denoted as

A*x = b

which is why different vectors (d and w) are used here. Looking ahead, generally the estimate of the solution vector x is updated on each iteration by some multiple of d, and the residual vector r = b - A*x is updated by some multiple of w. Because the matrix vector product is needed on each iteration, it accounts for the majority of computational time for iterative solvers.

You should be able to easily write the code in one of C, C++, or Fortran that does this multiplication when A is a dense matrix stored in a 2D array. If not, stop reading this until you can do so.

Transposed sparse matrix-vector multiplication

The transposed matrix-vector product
w = AT * d
also occurs in some iterative solvers, but less frequently than the (untransposed) matrix-vector product. For most of the data structures described below, the w = A * d operation is used to illustrate manipulating a sparse matrix, but you should try also writing down the transposed matrix-vector product to make sure you understand the data structures.

Direct solvers

No single optimal data structure exists for direct solvers. E.g., For LU factorization, "banded storage" can be the most efficient choice in terms of computational rates for regular mesh finite difference matrices from PDE's. For more irregular problems with complicated geometries, skyline storage is best for do-it-yourself direct solvers (e.g., LU factorization), if pivoting is unnecessary. All fill-in will be within the "skyline" or "envelope" of nonzeros, but pivoting destroys that feature.

For direct methods (LU factorization, QR factorization) even higher computational rates are possible with frontal and supernodal methods, but those beyond the scope of this course. Those methods allow using BLAS level 2 and sometimes BLAS level 3 methods, with concommittent efficiencies. Implementations and descriptions of those can be found in the packages SuperLU (from Xiaoye Li, Lawrence Berkeley Lab) and UMFPACK (from Tim Davis at the University of Florida). The most general and widely used data structure is compressed sparse column (CSC), and it is the one that many libraries for sparse matrix problems use. With the increase in memory sizes available, most linear solver packages now allow the user to input a sparse matrix in COO format and then convert it internally to a more efficient data structure, often CSR. Before, this was not possible - two copies of the sparse matrix could not be held in memory, so any conversion required complicated in-place methods.

Next page: Fill-in and LU Factorization