26
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – Bill Barth Sep 6 '12 at 19:55
  • 2
    $\begingroup$ 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. $\endgroup$ – Daniel Velkov Sep 6 '12 at 20:53
  • $\begingroup$ @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? $\endgroup$ – none Sep 6 '12 at 21:25
  • $\begingroup$ 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? $\endgroup$ – Daniel Velkov Sep 6 '12 at 21:36
  • 4
    $\begingroup$ What are you going to do with the permuted matrix? It may be better to simply permute the vector when applying the operator. $\endgroup$ – Jed Brown Sep 8 '12 at 15:22
9
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ @none is hack 2. what you call 'manually swapping rows and columns'? $\endgroup$ – Stefano M Sep 6 '12 at 23:22
  • 1
    $\begingroup$ 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. $\endgroup$ – Erik P. Sep 7 '12 at 3:36
  • $\begingroup$ @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. $\endgroup$ – Stefano M Sep 7 '12 at 7:45
  • 2
    $\begingroup$ This is a really good canidate for a cython function. Shoudln't be more than 10 lines . . . want me to give it a crack? $\endgroup$ – meawoppl Sep 9 '12 at 0:29
  • $\begingroup$ Lol. I started to Cython this, then found the right answer in a function that I use all the time. Doh. See my posted answer. $\endgroup$ – meawoppl Jan 14 '14 at 23:57
5
$\begingroup$

Warning: The below example works properly, but using the full set of parameters suggested at the post end exposes a bug, or at least an "undocumented feature" in the numpy.take() function. See comments below for details. Bug report filed.

You can do this in-place with numpy's take() function, but it requires a bit of hoop jumping.

Here is an example of doing a random permutation of an identity matrix's rows:

import numpy as np
i = np.identity(10)
rr = range(10)
np.random.shuffle(rr)
np.take(i, rr, axis=0)
array([[ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.]])

To do it in-place, all you need to do is specify the "out" parameter to be the same as the input array AND you have to set the mode="clip" or mode="wrap". If you don't set the mode it will make a copy to restore the array state on a Python exception (see here).

On a final note, take seems to be a array method, so instead of

np.take(i, rr, axis=0)

you could call

i.take(rr, axis=0)

if that is more to your taste. So in total you call should look something like the following:

#Inplace Rearrange
arr = makeMyBixMatrix()
pVec0, pVec1 = calcMyPermutationVectors()
arr.take(pVec0, axis=0, out=arr, mode="clip")
arr.take(pVec1, axis=1, out=arr, mode="clip")

To permute both rows and columns I think you either have to run it twice, or pull some ugly shenanigans with numpy.unravel_index that hurts my head to think about.

$\endgroup$
  • $\begingroup$ As said, in place algorithms are hard. Your solution does not work with numpy 1.6.2. and 1.7.1 (duplicate rows/columns). Had no time to check if 1.8.x fixes this problem $\endgroup$ – Stefano M Jan 20 '14 at 16:16
  • $\begingroup$ Hmmm. Can you post test code somewhere? In my head, I feel as though there needs to be a sort operation on the indices that happens first before the plucking. I will investigate more this PM. $\endgroup$ – meawoppl Jan 21 '14 at 23:13
  • 1
    $\begingroup$ When I run this code I get 1.6.2, test take, not overwriting: True, test not-in-place take: True, test in-place take: False, rr [3, 7, 8, 1, 4, 5, 9, 0, 2, 6], arr [30 70 80 70 40 50 90 30 80 90], ref [30 70 80 10 40 50 90 0 20 60]. So np.take at least for numpy 1.6.2 is not aware of doing a in-place permutation and messes things up. $\endgroup$ – Stefano M Jan 23 '14 at 21:51
  • $\begingroup$ Yeouch. Well demonstrated. This probably qualifies as a bug IMHO. At the very least the docs should say input and output can not be the same array, probably check to see and except if it is. $\endgroup$ – meawoppl Feb 12 '14 at 7:14
  • $\begingroup$ Agreed on the bug: maybe you should add a note to your post to warn readers that your solution can produce wrong results. $\endgroup$ – Stefano M Feb 12 '14 at 7:57
2
$\begingroup$

If you have a sparse matrix stored in COO format, the following might be helpful

    A.row = perm[A.row];
    A.col = perm[A.col];

assuming that A contains the COO matrix, and perm is a numpy.array containing the permutation. This will only have $m$ memory overhead, where $m$ is the number of non-zero elements of the matrix.

$\endgroup$
  • $\begingroup$ but what's the memory overhead for storing a full dense matrix as a sparse C00 matrix in te first place? $\endgroup$ – Federico Poloni Jan 14 '14 at 2:54
  • $\begingroup$ As the number of elements is equal in both a sparse and a (full) dense representation, the memory difference is merely a constant (2 ints and 1 float in sparse representation per element as to a single float in the dense representation). But the memory overhead of this method will be $n^2$ in a dense case, so one could probably then better stick to regular numpy.ndarrays. $\endgroup$ – Vincent Traag Jan 28 '14 at 9:06
1
$\begingroup$

I don't have enough reputation to comment, but I think the following SO question might be helpful: https://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.

$\endgroup$
  • $\begingroup$ The OP is asking for a permutation, and this is not possible with basic slicing. $\endgroup$ – Stefano M Apr 13 '13 at 19:40
  • $\begingroup$ 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. $\endgroup$ – hadsed Apr 14 '13 at 13:40
-1
$\begingroup$

What about

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

$\endgroup$
  • 1
    $\begingroup$ This constructs a temporary, which is exactly what he wants to avoid. $\endgroup$ – Michael Grant Apr 1 '13 at 5:35

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.