On the reordering of sparse matrices

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).

• Do you realize the Cuthill-Mckee algorithm was published back in 1969? In numerical analysis, 50-year-old methods are mostly "ancient history", and that is certainly true for a simple bandwidth minimization algorithm. Even if you want to minimize bandwidth (which is irrelevant for modern equation solvers), there are better methods to do it. Sep 13 '20 at 18:08
• In my finite element application, I have a 3D mesh in which each node has an integer index that points to coordinates. The indexing of these nodes is totally arbitrary, except it affects the sparsity pattern of the system. After feeding a connectivity matrix into METIS and reordering the nodes (I have zero knowledge of the theory behind what it does), my sparse matrix-vector products (intel MKL) were about 30% faster, with the same numerical results and number of iterations for my conjugate gradient solver. In short, the reordering did not affect the calculations, except that they were faster. Sep 14 '20 at 3:18
• @alephzero in HPC applications, the RCM method is still widely used and has been shown in a recent study to be more efficient than other existing methods for some applications on massively parallel computers Sep 14 '20 at 8:06
• @CharlieS Are you referring to METIS' fill reducing matrix ordering or mesh partitioning function? Not sure if you're talking about a serial or parallel solver Sep 14 '20 at 8:14
• @solalito I used METIS_NodeND, described as "ﬁll reducing orderings of sparse matrices using the multilevel nested dissection algorithm". I use the parallel Intel MKL library for sparse matrix-vector products (mkl_sparse_?_mv), so you could say the bottleneck of my conjugate gradient solver is parallelized. Sep 14 '20 at 11:50

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.

• Thanks for the answer. Are there special considerations to be taken when solving the system in parallel? In this case, I'm guessing the minimum theoretical bandwidth ts the number of elements per core since the neighboring elements are on another core. Sep 14 '20 at 8:12
• See my answer for a bunch of considerations about parallelism. You need to take a much more global view. And bandwidth has very little to do with parallelis. In fact, orderings that are better for parallelism typically have a very large bandwidth. Sep 15 '20 at 0:25

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.

• The OP asked about convergence, and then re-ordering of ILU(0) can definitely make a big difference. Sep 28 '20 at 14:58

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.