Approximate minimum degree permutation algorithm in Mathematica

While the MinCut (and also the MinimumBandwidth command in the GraphUtilities package) may be useful in many situations, this is a very general problem and might also benefit from more general solutions. For example, a while ago I had some data about the relationship between 120 items in the form of a distance matrix, basically a measurement of the distance between each object. With the data labelled arbitrarily, the matrix looks like this:

Black means zero distance (which you see along the diagonal because each item is zero distance from itself) while white represents the maximum possible distance between the items. Obviously, the matrix is completely scrambled. Laboriously using prior knowledge of the relationships between the various items, we were able to manually reorder the distance matrix into the form

which nicely shows the clustering along the diagonal. Though we were quite happy with this, it nagged at me: shouldn't there be some way to automatically reorder the data and hence rearrange the matrix? I tried two commands from Matlab (symamd -- a variant of amd, and symrcm which implements the minimum bandwidth algorithm of Cuthill-McKee). I tried both MinCut and MinimumBandwidth ordering in the graph utilities in Mathematica. None came close to what we could do manually -- that is -- the reordered matrices, while more organized than the first matrix above, were still quite random looking.

At this point, someone suggested that what we really needed was a clustering algorithm because what we were trying for was to cluster the elements into small groups. Accordingly, I tried

clusters = FindClusters[scrambled -> Range[Length[scrambled]], 10, 
           DistanceFunction -> (Norm[#1 - #2, 1] &)];
ind = Flatten[clusters];
ArrayPlot[scrambled[[ind, ind]]]

where scrambled is the unsorted (original) distance matrix. This was by far the best automated results we were able to obtain. The beauty of it is that by changing the DistanceFunction and the number of requested clusters, we could get many different plausible clusterings.

Finally I got a feedback from Wolfram support on the AMD algorithm. It turned out that there is (almost as usual) an undocumented implementation of the AMD algorithm within Mathematica. The algorithm is exactly identical to the MATLAB implementation, thus exactly what I was looking for. By calling

SparseArray`ApproximateMinimumDegree[m_Matrix]

one gets the order index array for the preconditioned matrix through

G=G[[order,order]]

The algorithm results in a moderate speedup of the CholeskyDecomposition for very large (rank > 10000) in the range of a few percent. It seems that the Decomposition is in general already somehow optimized. Biggest drawback of the CholeskyDecomposition routine is that it does not return SparseArrays yet. However support told me they are working on an sparse implementation.

There is an alternative solution to implement AMD into Mathematica, by using the MinCut function from the GraphUtilities Package. This function is not just working on Graphs but is also usable for Matrices.

The algorithm is reordering the Matrix into blocks and effectively reducing the off diagonal elements. It can be involved rather easily. If one starts with a matrix G to be preconditioned the following commands are involved.

Needs["GraphUtilities`"];
order = MinCut[G, Quotient[Length[G], 2]];
G = G[[Flatten[order], Flatten[order]]];

Here Length[G]/2 block diagonals are defined resulting in an almost diagonalized matrix. This matrix is much better suited for CholeskyDecomposition.