Sparse Matrices

One of the most important enabling developments in scientific computing is sparse matrix technology. This technology includes the data structures to represent the matrices, the techniques for manipulating them, the algorithms used, and the efficient mapping of the data structures and algorithms to high performance and parallel machines.

Sparse matrices occur in many applications including solving partial differential equations (PDEs), text-document matrices used for latent semantic indexing (LSI), linear and nonlinear optimization, and manipulating network and graph models. The type of matrices used in this course will be similar to those from PDE applications, using various discretization methods: finite differences, finite volumes, finite elements, or spectral elements.

A sparse matrix is one with most of its entries equal to 0.0 (zero), or (paraphrasing Iain Duff) one with enough zeros to take advantage of. The opposite of sparse is dense, meaning that all of the entries of the matrix are explicitly stored, even if some of them are exactly equal to zero. Almost always 1D vectors are stored as dense 1D arrays, even when the 2D matrices involved are sparse.

The notation used from now on is:

For some applications a small number of entries equal to zero are actually stored, so a better definition of nnz is "the number of interesting entries" or "the number of entries explicitly stored". However, the notation nnz is now standard in almost all of the literature and codes involving sparse matrices, so it also will be used here.

Categories of Linear System Solvers

Solution methods for linear systems of equations can be placed in three general categories:

  1. direct methods factor the coefficient matrix into 2 or more matrices of simpler type, which can be handled separately. E.g., Gaussian elimination factors A = LU, where L and U are triangular. Direct methods are accurate and well-understood in their numerical behaviour, but for PDEs they can involve "fill-in" and hence memory usage. Also, direct methods typically do not give any part of the solution until the very end of the algorithm. With LU factorization, >O(n3) (more precisely, >2*n3/3)) operations must be completed to get the LU factors before the O(n2) triangular solves can be done. Only in the last phase, when performing the upper triangular solves with U that any part of the solution vector is generated.
  2. iterative methods generate a sequence of approximations x0, x1, x2, ... that converge to a solution. E.g., SOR iteration, conjugate gradient method. Iterative methods give an approximate solution at every iteration, unlike direct methods. Iterative methods typically use much less memory by not requiring the fill-in generated by factorization methods. However, iterative methods can fail to converge and for PDEs of modern research interest, can behave erratically.
  3. hybrid methods mix direct and iterative. Examples: Gaussian elimination followed by "iterative refinement", or an iterative method "preconditioned" by a direct method to generate an approximate factorization.
All robust practical methods are hybrid in some way or another. This is a rule applicable across most of computer science, not just scientific computing! Whenever you see multiple methods with complementary strengths and weaknesses, start looking for a hybrid approach. As a general rule, this idea is a great Ph.D. thesis generator.