Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

In MATLAB there is a nice implementation of the so called AMD (approximate minimum degree permutation) algorithm named amd (see Online MATLAB Documentation). There is an alternative algorithm called COLAMD (see Online MATLAB Documentation). This algorithm is especially interesting when one wants to model large resistor networks e.g. calculate path resistances or reduce a network of resistors (for the theory see, e.g., Efficient Methods for Large Resistor Networks). There the AMD step is a preconditioning step for a large matrix before a Cholesky decomposition (available in Mathematica through CholeskyDecomposition) is done.

The original fortran coded algorithms can be found on the webpage from the autors.

Is anybody aware of an implementation of the AMD or the COLAMD algorithm in Mathematica?

share|improve this question
If you do have access to MATLAB but you prefer to do your work in Mathematica because you're much more fluent in it (like me), then consider using MATLink to call these MATLAB functions directly. Disclosure: I'm one of the MATLink authors. –  Szabolcs Dec 11 '13 at 21:28
Would it make sense to instead use MinCut to permute rows/columns? –  Daniel Lichtblau Dec 11 '13 at 21:56
@Szabolcs: It crossed my mind to use MATLink :-) , that would be my last resort if there is no Mathematica implementation available. Since I'm planning to do a broad statistical analysis on the networks (thousands of different resistor networks topographies) I hesitate to involve MATLink if not absolutely necessary. –  Rainer Dec 12 '13 at 5:30
There is an example in the documentation for MinCut that applies to permuting rows/columns of a matrix, for the purpose of decreasing fill-in. Which is to say, it is a precondition-by-permutation example. As such, it could be an alternative to AMD. –  Daniel Lichtblau Dec 12 '13 at 15:51
Also I wonder if just reordering based on vertex degree would suffice as a preconditioner? For a sparse matrix that should be fairly efficient, though I don't know how effective the resulting preconditioning would be in practice. –  Daniel Lichtblau Dec 12 '13 at 15:53
show 3 more comments

3 Answers

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.

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.

share|improve this answer
add comment

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:

enter image description here

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

enter image description here

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]]]

enter image description here

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.

share|improve this answer
add comment
up vote 2 down vote accepted

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


one gets the order index array for the preconditioned matrix through


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.

share|improve this answer
add comment

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.