Tutorial: Bandwidth reduction - The CutHill-McKee Algorithm

1. The problem

I've recently come across an interesting NP-complete problem :D. It's the bandwidth reduction problem, also known in the field of sparse matrix applications as the bandwidth minimization problem (or BMP in short):

For a given simetric sparse matrix, M(nxn), the problem is to reduce its bandwidth B by permuting rows and columns such as to move all the nonzero elements of M in a band as close as possible to the diagonal.

In other words, the problem consists of transforming through successive row and column permutations this (8x8 input matrix):

into, possibly, this (8x8 output matrix):

Now here's an interesting thing regarding this problem. If you consider the given sparse matrix (the input matrix) and the resulting matrix (the output matrix) as adjacency matrices for some given graphs, this is what you end up with:

.

The graph associated with the input matrix.

The graph associated with the output matrix.

As you can easily see, the two graph structures are identical, the only thing that is different is the node (vertex) labeling. In other words the bandwidth reduction problem can also be viewed as a graph labeling problem:

Find the node labeling that minimizes the bandwidth B of the adjacency matrix of the graph G(n) , where we can formally define: B=max|Li-Lj|, i,j=1..n and Li is the label of node i, Lj is the label of node j and nodes i and j are adjacent or neighbours (i.e., there is an edge in graph G(n) connecting nodes i and j).

2. The solution

Back in the 60's and 70's many graph algorithms were proposed for "solving" the bandwidth reduction problem. I wrote "solving" beacause the problem is after all NP-complete and none of the algorithms described claimed to find the exact solution(s) regardless of the input matrix, but they were more or less succesfull in finding a relatively good solution in a resonable amount of time (the so called heuristic approach to NP-complete problem solving).
.
Initially the most widely used of these heuristics was the Reverse Cuthill-McKee algorithm (RCM), a modification by Alan George of the original algorithm developed by CutHill and McKee in 1969. In 1976 the GPS algorithm (named after its developers Gibbs, Poole and Stockmeyer) was proposed. On average GPS and RCM find solutions of relative equal quality, but GPS is up to ten times as fast.

One last thing I should remind you before I continue with the description of the RCM algorithm: the degree of a node in a graph is the number of nodes adjacent to it. So here is the RCM algorithm for a given graph G(n).

Step 0: Prepare an empty queue Q and an empty result array R.

Step 1: Select the node in G(n) with the lowest degree (ties are broken arbitrarily) that hasn't previously been inserted in the result array. Let us name it P (for Parent).

Step 2: Add P in the first free position of R.

Step 3: Add to the queue all the nodes adjacent with P in the increasing order of their degree.

Step 4.1: Extract the first node from the queue and examine it. Let us name it C (for Child).

Step 4.2: If C hasn't previously been inserted in R, add it in the first free position and add to Q all the neighbours of C that are not in R in the increasing order of their degree.

Step 5: If Q is not empty repeat from Step 4.1 .

Step 6: If there are unexplored nodes (the graph is not connected) repeat from Step 1 .

Step 7: Reverse the order of the elements in R. Element R[i] is swaped with element R[n+1-i].

The result array will be interpreted like this: R[L] = i means that the new label of node i (the one that had the initial label of i) will be L.

As you can see steps 1-6 describe nothing but a classic Breadth First Search (BFS) with the minor modification that the nodes are explored in the increasing order of their degree. Step 7 is not mandatory, it is the modification introduced by George to the initial algorithm (it has the purpose of further reducing the profile of a matrix, but that is another story). If you don't reverse the result array you end up with the original Cuthill-McKee algorithm.

Here is how the RCM algorithm should perform on the 8x8 input matrix presented at the start of this post:
- At first we must choose a starting node. Nodes initially labeled with 1, 4 and 7 are equally good choices as all of them have the degree of 1.
- Let us start with node 1. We'll add 1 in R and its neighbours in Q in the increasing order of their degree: R=[1]; Q=[5];
- We extract the first element from Q. As 5 isn't in R we add it. We also add all the neighbours of 5, that are not in R, in Q in the increasing order of their degree: R[1,5]; Q=[3];
- We extract the first element from Q. 3 is not in R so we add 3 to R and all the neighbours of 3, that are not in R, to Q in the increasing order of their degree: R=[1,5,3]; Q=[2];
- We extract the first element from Q. 2 is not in R so we add 2 to R and all the neighbours of 2, that are not in R, to Q in the increasing order of their degree: R=[1,5,3,2]; Q=[6,8];

- We extract the first element from Q. 6 is not in R so we add 6 to R and all the neighbours of 6, that are not in R, to Q in the increasing order of their degree: R=[1,5,3,2,6]; Q=[8,8];

- We extract the first element from Q. 8 is not in R so we add 8 to R and all the neighbours of 8, that are not in R, to Q in the increasing order of their degree: R=[1,5,3,2,6,8]; Q=[8];
- We extract the first element from Q. 8 is in R so we do nothing: R=[1,5,3,2,6,8]; Q=[];
- Q is empty so we check the size of R. If the size of R equals 8 then we are done, otherwise there are still unexplored nodes (the graph is not connected). As the size of R is 6 we must restart from Step 1.
- We must again choose a starting node. Nodes labeled 4 and 7 are equally good start solutions as they both have a degree of 1.
- We choose node 4. We'll add 4 in R and its neighbours in Q in the increasing order of their degree: R=[1,5,3,2,6,8,4]; Q=[7];

- We extract the first element from Q. 7 is not in R so we add 7 to R and all the neighbours of 7, that are not in R, to Q in the increasing order of their degree: R=[1,5,3,2,6,8,4,7]; Q=[];

- Q is empty so we check the size of R. The size of R is 8 so all the nodes have been relabeled.
- We reverse R and obtain the solution R=[7,4,8,6,2,3,5,1] with the following interpretation:

If you look closely at the two graph images at the start of the post you can see how the node which had the initial label of 7 will be relabeled with 1, the node which had the initial label of 4 will be relabeled with 2, the node which had the initial label of 8 will be relabeled with 3, and so on.

.
The only thing left now is to build the adjacency matrix of the relabeled graph. You should obtain the 8x8 output matrix presented at the start of this tutorial.  

P.S. For a Delphi implementation (source code included) of the Cuthill-McKee and RCM algorithms, as well as a comparison of these methods with a more modern aproach based on a genetic algorithm and simulated annealing please see this post.

P.P.S. Help required / Collaboration opportunity

I have recently (February 2013) been contacted by Dr. Lubomir ONDRIS (lubomir.ondris(at)gmail.com) who is a scientist from Slovakia that has dedicated a lot of effort over the years into studying the bandwidth reduction problem and improving the available algorithms that attempt to solve it.

He created a personal hybrid algorithm that is quite competitive on the Everstine and other mesh collections (i.e., common benchmarks for testing the performance of bandwidth reduction routines) and he would need access to non-artificial, fairly big (>500 nodes) and complicated (inter-connected) meshes (matrices) in order to further test and improve his algorithm. Typically, such matrices are encountered when solving very large equation systems that arise in finite element analysis.

If you can help Dr. Ondris in any way with regard to his data requirements, please contact him directly via the e-mail address provided above. Here is a short description, in Dr. Ondris' own words, of what he requires:

"Again, what I need are parts of standard input data for finite element programs: Finite elements with their node numbers and nodes with their degrees of freedoms. And, of course, the bandwidth and profile achieved by the renumbering routine used."

Hope that you can help & thank you very much!