Efficient Way to Permutate a Symmetric Square Matrix in Numpy
<p>What's the best way to do the following in Numpy when dealing with <strong>symmetric square matrices</strong> (<code>NxN</code>) where <code>N &gt; 20000</code>?    </p>

<pre><code>&gt;&gt;&gt; a = np.arange(9).reshape([3,3])
&gt;&gt;&gt; a = np.maximum(a, a.T)
&gt;&gt;&gt; a
array([[0, 3, 6],
       [3, 4, 7],
       [6, 7, 8]])
&gt;&gt;&gt; perm = np.random.permutation(3)
&gt;&gt;&gt; perm
array([1, 0, 2])
&gt;&gt;&gt; shuffled_arr = a[perm, :][:, perm]
&gt;&gt;&gt; shuffled_arr
array([[4, 3, 7],
       [3, 0, 6],
       [7, 6, 8]])
</code></pre>

<p>This takes about 6-7 secs when N is about 19K. While the same opertation in Matlab takes less than a second:</p>

<pre><code>perm = randperm(N);
shuffled_arr = arr(perm, perm);
</code></pre>

<pre><code>In [703]: N=10000
In [704]: a=np.arange(N*N).reshape(N,N);a=np.maximum(a, a.T)
In [705]: perm=np.random.permutation(N)
</code></pre>

<p>One indexing step is quite a bit faster:</p>

<pre><code>In [706]: timeit a[perm[:,None],perm]   # same as `np.ix_...`
1 loop, best of 3: 1.88 s per loop

In [707]: timeit a[perm,:][:,perm]
1 loop, best of 3: 8.88 s per loop

In [708]: timeit np.take(np.take(a,perm,0),perm,1)
1 loop, best of 3: 1.41 s per loop
</code></pre>

<p><code>a[perm,perm[:,None]]</code> is in the 8s category.</p>

python matlab numpy