How to assign square subarrays in a large matrix without loops in numpy

Question

How can I vectorize this loop that fills the two square submatrices of the larger matrix (also preserves the greater matrix symmetry) using numpy arrays:

for x in range(n): assert m[x].shape == (n,) M[i:i+n,j+x] = m[x] M[j+x,i:i+n] = m[x] 

This is tempting, but not consistent with the above loop (see example of disagreement below):

 assert m.shape == (n,n) M[i:i+n,j:j+n] = m M[j:j+n,i:i+n] = m 

Here is a small example (failure for n> 1):

 from numpy import arange,empty,NAN from numpy.testing import assert_almost_equal for n in (1,2,3,4): # make the submatrix m = (10 * arange(1, 1 + n * n)).reshape(n, n) N = n # example below, submatrix is the whole thing # M1 using loops, M2 "vectorized" M1 = empty((N, N)) M2 = empty((N, N)) M1.fill(NAN) M2.fill(NAN) i,j = 0,0 # not really used when (n == N) # this results in symmetric matrix for x in range(n): assert m[x].shape == (n,) M1[i:i+n,j+x] = m[x] M1[j+x,i:i+n] = m[x] # this does not work as expected M2[i:i+n,j:j+n] = m M2[j:j+n,i:i+n] = m assert_almost_equal(M1,M1.transpose(),err_msg="M not symmetric?") print "M1\n",M1,"\nM2",M2 assert_almost_equal(M1,M2,err_msg="M1 (loop) disagrees with M2 (vectorized)") 

As a result, we get:

 M1 = [10 30 30 40] # symmetric M2 = [10 20 30 40] # ie m 

Answer 1

Your test is incorrect: for i, j = 0,0, your M2 [] = assignments just overwrite the same matrix block.

The reason you get the symmetric matrix when using M1 is because you assign the values ​​of M1 in one cycle.

if you divide the loop into two:

 for x in range(n): M1[i:i+n,j+x] = m[x] for x in range(n): M1[j+x,i:i+n] = m[x] 

M1 will obviously be the same as M2.

So, to summarize, the following code works (equivalent to your calculation of M2), but! it will only work if there is no overlap between the blocks above and below the diagonal. if you need to decide what to do there

 xs=np.arange(4).reshape(2,2) ys=np.zeros((7,7)) ys[i:i+n,j:j+n]=xs ys[j:j+n,i:i+n]=xs.T print ys >> array([[ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 1., 0., 0.], [ 0., 0., 0., 2., 3., 0., 0.], [ 0., 0., 2., 0., 0., 0., 0.], [ 0., 1., 3., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.]])