On the reordering of sparse matrices

Question

I have been reading on different techniques used to reorder sparse matrices to achieve better performance, the most popular being the Cuthill-McKee or Reverse Cuthill-McKee algorithm. Most of those techniques focus on reducing the bandwidth of the matrix, which is defined as the furthest distance of a non-zero entry from the main diagonal.

If I understand correctly, the reordering algorithm optimizations are purely computational: the goal is to come up with a storage of the matrix more adequate for caching effects. Is this true or does the reordering affect the convergence/accuracy of the iterative algorithm?

Also, why do people mostly focus on the matrix bandwidth and not for example at the average bandwidth of each matrix row? If most of the points are located next to the diagonal except for 1, the matrix bandwidth will still be high but the performance should be nearly optimal (assuming matrix size n >> 1).

Answer 1

When using an iterative method, you will typically use a preconditioner that speeds up convergence. A good example is the incomplete LU factorization (ILU).

When you take the LU factorization of a sparse matrix, the L and U factors might lose some of its sparsity, the extra entries are called fill in. The ILU will ignore some of this fill in to form a approximate factorization.

When you reorder the matrix, the main goal is to reduce the amount of fill in (banded matrices have almost no fill in, however, just one entry can really cause a lot of fill in if far from the diagonal). This means the approximation of the ILU will be more correct and you need less storage.

Edit: this is (to my knowledge) the most common use of reordering when solving linear systems. It can of course also serve other purposes. If the memory access pattern is better, multiplications will typically also be better.

Answer 2

For direct factorization, you would ideally want to minimize the total fill-in. However, this is an NP-Hard combinatorial optimization problem that is intractable to solve for matrices of interesting size. Reducing the bandwidth does reduce an upper bound on the fill-in, which is a surrogate for what you really want.

In iterative methods, incomplete LU factorization (ILU) is often used as a preconditioning strategy- this is, an approximate LU factorization $A \approx LU$ is computed, and then $M^{-1}=U^{-1}L^{-1}$ (actually solving systems of equations involving $U$ and $L$) is used to precondition the iterative method.

In the simplest version of ILU, called ILU(0), no fill-in is allowed in the approximate factorization- the sparsity pattern of $A$ is used. In the ILU(0) scheme there isn't any need to reorder for sparsity.

There are also more complicated ILU schemes. For example, ILU(k) produces an approximate LU factorization with the sparsity pattern of $A^{k}$. Another common scheme is to drop fill-in entries below some size tolerance. For these more complicated variations on ILU, it can help to reorder A to reduce fill-in.

Answer 3

You have found Cuthill McKee, but there is also the "minimum degree" method as well as a bunch of others. Here are some considerations.

1. You state that these method improve caching. Well, this stuff was invented long before there were caches. The actual motivation for Cuthill McKee is reduction of the bandwidth.
2. You ask about iterative methods: these schemes are originally intended for full LU factorization, which is a direct method. Analyzing the effect of reordering on ILU is much harder.
3. Minimum degree and "multiple minimum degree" try to find independent points which benefit vectorization and parallelism.
4. You have missed "spectral" reordering methods, which use the "Fiedler vector" or "Perron vector" to split the variables into two sets with minimal connection. This is done for parallelism, especially if you do it recursively.
5. Finding sets taht can be handled independently is also an old idea. For Cartesian domains this was done by Alan George in 1971 or so, later there were papers by Tarjan that did it for general graphs. The motivation here is initially fill-in reduction, but later this was of course applied to parallelism: domain decomposition.
6. Since you're asking about bandwidth and paralleism: orderings for parallelism often have a large bandwidth, more or less equal to the whole matrix. However, this does not matter since fill-in estimates based on bandwidth are only upper bounds.
7. Other graph-based domain decomposition / graph splitting methods do not have method names attached, but can be found in software packages such as Metis, Zoltan, Chaco.
8. In most cases the importance of parallelism is such that one takes a slight increase in number of iterations for granted, but there are papers by Benzi and Szyld that show that re-orderings can sometimes reduce the number of iterations.
9. You can also order the matrix by sorting the rows by length. This can benefit the matrix-vector product substantially. D'Azevedo showed this for the Cray X1 vector computer, and this idea was also applied by many authors to GPU implementation. Intel even tried to patent these 50 year old ideas.