title: Compute commutation matrix in numpy/scipy efficiently
body: <p>I am trying to compute the <a href="https://en.wikipedia.org/wiki/Commutation_matrix" rel="nofollow noreferrer">commutation matrix</a> in python for a large dataset. I wrote the following code but found it performs terribly (and runs into memory errors for examples of around 500 by 500). In my code <code>a</code> and <code>b</code> are equivilent to the <code>m</code> and <code>n</code> notation in the linked wikipedia page. Can anyone provide a quicker and more memory efficient alternative to my current attempt?</p>
  
  <pre><code>def vec(matrix):
  
      #Return vectorised matrix
      return(matrix.transpose().reshape(matrix.shape[0]*matrix.shape[1],1))
  
  
  def commutation(a, b):
  
      # Example matrix with unique elements
      m = np.arange(a*b).reshape(a,b)
  
      # Vec(m) 
      vecm = vec(m)
      vecm = vecm.reshape(vecm.shape[0])
  
      # Get row inds
      rowInds = np.arange(a*b)
  
      # Get column inds
      colInds = np.argsort(vecm)
      colInds = colInds.reshape(colInds.shape[0])
  
      # Work out mapping between them.
      K = scipy.sparse.csr_matrix((np.ones(a*b),(rowInds,colInds)))
  
      return(K)
  </code></pre>