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$ – Brian Borchers 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$ – Johan Falkenjack 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
  • 1
    $\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$ – Brian Borchers 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.

| cite | improve this answer | |
  • 1
    $\begingroup$ Thanks! Today this seems totally obvious but last night my brain apparently took a leave of absence :) $\endgroup$ – Johan Falkenjack Aug 25 '16 at 11:09
  • 1
    $\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$ – w1th0utnam3 Mar 27 '17 at 15:09

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.