Tell me more ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.
up vote 12 down vote favorite

I want to modify a dense square transition matrix in-place by changing the order of several of its rows and columns, using python's numpy library. Mathematically this corresponds to pre-multiplying the matrix by the permutation matrix P and post-multiplying it by P^-1 = P^T, but this is not a computationally reasonable solution.

Right now I am manually swapping rows and columns, but I would have expected numpy to have a nice function f(M, v) where M has n rows and columns, and v has n entries, so that f(M, v) updates M according to the index permutation v. Maybe I am just failing at searching the internet.

Something like this might be possible with numpy's "advanced indexing" but my understanding is that such a solution would not be in-place. Also for some simple situations it may be sufficient to just separately track an index permutation, but this is not convenient in my case.

Added:
Sometimes when people talk about permutations, they only mean the sampling of random permutations, for example as part of a procedure to obtain p-values in statistics. Or they mean counting or enumerating all possible permutations. I'm not talking about these things.

Added:
The matrix is small enough to fit into desktop RAM but big enough that I do not want to copy it thoughtlessly. Actually I would like to use matrices as large as possible, but I don't want to deal with the inconvenience of not being able to hold them in RAM, and I do O(N^3) LAPACK operations on the matrix which would also limit the practical matrix size. I currently copy matrices this large unnecessarily, but I would hope this could be easily avoided for permutation.

share|improve this question
It would be good if you could update the question to give the size of your matrices. "Gigantic" does not mean the same thing to all people. – Bill Barth Sep 6 '12 at 19:55
2  
You are right that the advanced (or so called fancy) indexing creates a copy. But if you accept to live with that fact then your code is just M[v] to permute the rows. – Daniel Velkov Sep 6 '12 at 20:53
@daniel: And it would be M[v, :][:, v] to do the whole permutation? Would this be the best way to get the permutation using fancy indexing? And would it use 3x the matrix memory, including the size of the original matrix, the row+column permuted matrix, and the temporary row permuted matrix? – none Sep 6 '12 at 21:25
That's correct, you would have your original matrix and 2 copies. Btw why do you need to permute both rows and columns at the same time? – Daniel Velkov Sep 6 '12 at 21:36
2  
What are you going to do with the permuted matrix? It may be better to simply permute the vector when applying the operator. – Jed Brown Sep 8 '12 at 15:22

3 Answers

up vote 6 down vote

According to the docs, there is no in-place permutation method in numpy, something like ndarray.sort.

So your options are (assuming that M is a $N\times N$ matrix and p the permutation vector)

  1. implementing your own algorithm in C as an extension module (but in-place algorithms are hard, at least for me!)

  2. $N$ memory overhead

    for i in range(N):
    M[:,i] = M[p,i]
for i in range(N):
    M[i,:] = M[i,p]

  3. $N^2$ memory overhead

    M[:,:] = M[p,:]
M[:,:] = M[:,p]

Hope that these suboptimal hacks are useful.

share|improve this answer
@none is hack 2. what you call 'manually swapping rows and columns'? – Stefano M Sep 6 '12 at 23:22
1  
I would combine options 1 and 2: write C code that uses a buffer of order N to write each permuted column to, then writes it back to where it came from; then do the same for rows. As @Stefano writes, this takes only $O(N)$ extra memory, which you are already spending to store the permutation $p$ in the first place. – Erik P. Sep 7 '12 at 3:36
@ErikP. for a C implementation $O(N)$ extra memory is reasonable and for sure your scatter write to temp and copy back approach is sound. The interesting question however is if there are more efficient algorithms, given $O(N)$ extra memory. The answer is hard I think, since we should take into account processor architecture, memory access patterns, cache hits, ... This said I would follow your advice and go for a simple and easy to implement algorithm. – Stefano M Sep 7 '12 at 7:45
This is a really good canidate for a cython function. Shoudln't be more than 10 lines . . . want me to give it a crack? – meawoppl Sep 9 '12 at 0:29
up vote 2 down vote

I don't have enough reputation to comment, but I think the following SO question might be helpful: http://stackoverflow.com/questions/4370745/view-onto-a-numpy-array

The basic points are that you can use basic slicing and that will create a view on to the array without copying, but if you do advanced slicing/indexing then it will create a copy.

share|improve this answer
The OP is asking for a permutation, and this is not possible with basic slicing. – Stefano M Apr 13 at 19:40
You are correct of course. I thought it would be useful for the OP to understand what was happening with slicing (in case they didn't know) since they were concerned about when copies would be happening. If he used something from your answer, I think that would be good to know since you use them inside your loops. – hadsed Apr 14 at 13:40
up vote -1 down vote

What about

my_array[:,[0, 1]] = my_array[:,[1, 0]]

share|improve this answer
This constructs a temporary, which is exactly what he wants to avoid. – Michael C. Grant Apr 1 at 5:35

Your Answer

 
discard

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.