I have applied SciPy's implementation of the Cuthill-McKee algorithm to a $48 \times 48$ sparse non-symmetric matrix in Compressed Sparse Row (CSR) format and the output is an array of length $48$ which the documentation calls a "permutation array".

How do I use this "permutation" array?

Obviously, the array is in some way representing how my matrix is to be reordered, but a permutation array of length $48$ is only enough to reorder one axis. What am I missing?

  • $\begingroup$ Is your matrix symmetric and positive definite? In that case you would normally apply the permutation symmetrically to both rows and columns of the matrix. $\endgroup$ Aug 24 '16 at 20:50
  • $\begingroup$ Positive definite, but not symmetric, but apparently scipy adds the matrix and its transpose to make it symmetric before applying RCM. $\endgroup$ Aug 25 '16 at 11:10
  • $\begingroup$ @BrianBorchers, I'm not aware that there is any requirement on positive definiteness to use RCM. Is there a reason why you say this? I thought it just needs to be symmetric in structure for the algorithm to "make sense". $\endgroup$
    – Nick C.
    Aug 25 '16 at 14:55
  • 2
    $\begingroup$ If you're going to use RCM with a sparse Cholesky factorization then the matrix needs to be positive definite. That's the most common use-case. $\endgroup$ Aug 25 '16 at 16:02
  • $\begingroup$ Thanks, I see. My only experience with RCM is in the context of reordering the nodes of a finite element mesh, and then you are really just thinking about how those nodes are "connected" on a graph and the values (in the context of the final assembled FE system) don't matter, from what I recall. $\endgroup$
    – Nick C.
    Aug 26 '16 at 3:17

The reverse Cuthill-McKee algorithm produces a reordering that applies to both the rows and columns. This is because it works by considering matrices as graphs of (undirected) connected nodes. According to the function's documentation in SciPy, the output array is the permuted row/column indices, so you can simply do the following

perm = reverse_cuthill_mckee(arr_orig)
arr_perm = arr_orig[perm, perm]

and numpy's fancy indexing functionality will take care of it. Of course, this involves making a copy of your array, and I assume that's acceptable.

  • 3
    $\begingroup$ For larger sparse matrices it is much faster to apply the permutation directly to the sparse matrix index array: mat.indices = perm.take(mat.indices) than convert the matrix to csc: mat = mat.tocsc() and apply the permutation again: mat.indices = perm.take(mat.indices). If you need the matrix in csr you can convert it back again. Using the advanced indexing is much slower than this - even when the matrix is converted twice between the formats. Tested with 25k x 25k matrices. $\endgroup$ Mar 27 '17 at 15:09

I'm not sure if Nick's answer is right or maybe it was working in the previous version but it does not really in this one.

First of all, reverse_cuthill_mckee does not take a NumPy array as an input but it takes a csr_matrix.

I think it should be:

arr = [
    [0, 1 , 2, 0],
    [0, 0, 0, 1],
    [2, 0, 0, 3],
    [0, 0, 0, 0]
graph = csr_matrix(arr)
perm = reverse_cuthill_mckee(graph)
reord_arr = csr_matrix.toarray(graph[perm])

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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