The nice thing about array slicing in numpy is you don't need the for loops that you are using. Also the reason that it is only putting the center element is because you only put a single element there (c1[c1mid,c1mid] is a single number) here is what you could do:

    z[7:12,7:12] = c1
    z[7:12,27:32] = c2
    z[26:33,6:14] = c3
    z[25:34,25:33] = c4

There are several ways to get submatrix in numpy:

In [35]: ri = [0,2]
    ...: ci = [2,3]
    ...: a[np.reshape(ri, (-1, 1)), ci]
Out[35]: 
array([[ 2,  3],
       [10, 11]])

In [36]: a[np.ix_(ri, ci)]
Out[36]: 
array([[ 2,  3],
       [10, 11]])

In [37]: s=a[np.ix_(ri, ci)]

In [38]: np.may_share_memory(a, s)
Out[38]: False

note that the submatrix you get is a new copy, not a view of the original mat.

It seems to me that you're trying to do a simple convolution?

def do(m):
    rows, cols = m.shape
    newMatrix = np.zeros_like(m)
    for i in range(1, rows-1):
        for j in range(1, cols-1):
            sub = matrix[i-1:i+2, j-1:j+2]
            newMatrix[i][j] = numpy.sum(xWeight * sub)
    return newMatrix[1:-1, 1:-1]
>>> res1 = do(matrix)
>>> res2 = scipy.signal.convolve2d(matrix, xWeight)[2:-2,2:-2]
>>> np.allclose(np.abs(res1), np.abs(res2))
True

Didn't went into details about the sign, but that should hopefully put you on the right track.

to answer your question s is only a 1d array ... (even if you did actually transpose it ... which you did not)

>>> u,s,v = linalg.svd(A)
>>> s
array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])
>>>

for selecting a submatrix I think this does what you want ... there may be a better way

>>> rows = range(10,15)
>>> cols = range(5,8)
>>> A[rows][:,cols]
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])

or probably better

>>> A[15:32, 2:7]
array([[ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])

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 M1 values in a single loop.

if you would split 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]

The M1 will be obviously the same as M2.

So Summarizing, the following code works (equivalent to your M2 calculation) but ! it will only work if there is no overlap between blocks above and below the diagonal. if there is you have 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.]])

Here is how you can do it:

>>> A[3:5, 3:5] = B
>>> A
array([[ 1. ,  1. ,  1. ,  1. ,  1. ],
       [ 1. ,  1. ,  1. ,  1. ,  1. ],
       [ 1. ,  1. ,  1. ,  1. ,  1. ],
       [ 1. ,  1. ,  1. ,  0.1,  0.2],
       [ 1. ,  1. ,  1. ,  0.3,  0.4]])

Is n potentially large, so the result is a large sparse matrix with nonzero values concentrated along the diagonal? Sparse matrices are designed with this kind of matrix in mind (from FD and FE PDE problems). I did this a lot in MATLAB, and some with the scipy sparse module.

That module has a block definition mode that might work, but what I'm more familiar with is the coo to csr route.

In the coo format, nonzero elements are defined by 3 vectors, i, j, and data. You can collect all the values for A, B, etc in these arrays (applying the appropriate offset for the values in B etc), without worrying about overlaps. Then when that format is converted to csr (for matrix calculations) the overlapping values are summed - which is exactly what you want.

I think the sparse documentation has some simple examples of this. Conceptually the simplest thing to do is iterate over the n submatrices, and collect the values in those 3 arrays. But I also worked out a more complex system whereby it can be done as one big array operation, or by iterating over a smaller dimension. For example each submatrix has 16 values. In a realistic case 16 will be much smaller than n.

I'd have play around with code to give a more concrete example.

==========================

Here's a simple example with 3 blocks - functional, but not the most efficient

Define 3 blocks:

In [620]: A=np.ones((4,4),int)    
In [621]: B=np.ones((4,4),int)*2
In [622]: C=np.ones((4,4),int)*3

lists to collect values in; could be arrays, but it is easy, and relatively efficient to append or extend lists:

In [623]: i, j, dat = [], [], []

In [629]: def foo(A,n):
   # turn A into a sparse, and add it's points to the arrays
   # with an offset of 'n'
   ac = sparse.coo_matrix(A)
   i.extend(ac.row+n)
   j.extend(ac.col+n)
   dat.extend(ac.data)


In [630]: foo(A,0)

In [631]: i
Out[631]: [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3]    
In [632]: j
Out[632]: [0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3]

In [633]: foo(B,1)
In [634]: foo(C,2)  # do this in a loop in the real world

In [636]: M = sparse.csr_matrix((dat,(i,j)))

In [637]: M
Out[637]: 
<6x6 sparse matrix of type '<class 'numpy.int32'>'
    with 30 stored elements in Compressed Sparse Row format>

In [638]: M.A
Out[638]: 
array([[1, 1, 1, 1, 0, 0],
       [1, 3, 3, 3, 2, 0],
       [1, 3, 6, 6, 5, 3],
       [1, 3, 6, 6, 5, 3],
       [0, 2, 5, 5, 5, 3],
       [0, 0, 3, 3, 3, 3]], dtype=int32)

If I've done this right, overlapping values of A,B,C are summed.

More generally:

In [21]: def foo1(mats):
      i,j,dat = [],[],[]
      for n,mat in enumerate(mats):
          A = sparse.coo_matrix(mat)
          i.extend(A.row+n)
          j.extend(A.col+n)
          dat.extend(A.data)
      M = sparse.csr_matrix((dat,(i,j)))
      return M
   ....:   

In [22]: foo1((A,B,C,B,A)).A
Out[22]: 
array([[1, 1, 1, 1, 0, 0, 0, 0],
       [1, 3, 3, 3, 2, 0, 0, 0],
       [1, 3, 6, 6, 5, 3, 0, 0],
       [1, 3, 6, 8, 7, 5, 2, 0],
       [0, 2, 5, 7, 8, 6, 3, 1],
       [0, 0, 3, 5, 6, 6, 3, 1],
       [0, 0, 0, 2, 3, 3, 3, 1],
       [0, 0, 0, 0, 1, 1, 1, 1]], dtype=int32)

Coming up with a way of doing this more efficiently may depend on how the individual submatrices are generated. If they are created iteratively, you might as well collect the i,j,data values iteratively as well.

==========================

Since the submatrices are dense, we can get the appropriate i,j,data values directly, without going through a coo intermediary. And without looping if the A,B,C are collected into one larger array.

If I modify foo1 to return a coo matrix, I see the i,j,data lists (as arrays) as given, without summation of duplicates. In the example with 5 matrices, I get 80 element arrays, which can be reshaped as

In [110]: f.col.reshape(-1,16)
Out[110]: 
array([[0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3],
       [1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4],
       [2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5],
       [3, 4, 5, 6, 3, 4, 5, 6, 3, 4, 5, 6, 3, 4, 5, 6],
       [4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7]], dtype=int32)

In [111]: f.row.reshape(-1,16)
Out[111]: 
array([[0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
       [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4],
       [2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5],
       [3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6],
       [4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7]], dtype=int32)

In [112]: f.data.reshape(-1,16)
Out[112]: 
array([[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
       [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3],
       [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])

I should be able generate those without a loop, especially the row and col.

In [143]: mats=[A,B,C,B,A]

the coordinates for the elements of an array

In [144]: I,J=[i.ravel() for i in np.mgrid[range(A.shape[0]),range(A.shape[1])]] 

replicate them with offset via broadcasting

In [145]: x=np.arange(len(mats))[:,None]
In [146]: I=I+x    
In [147]: J=J+x

Collect the data into one large array:

In [148]: D=np.concatenate(mats,axis=0)

In [149]: f=sparse.csr_matrix((D.ravel(),(I.ravel(),J.ravel())))

or as a compact function

def foo3(mats):
    A = mats[0]
    n,m = A.shape
    I,J = np.mgrid[range(n), range(m)]
    x = np.arange(len(mats))[:,None]
    I = I.ravel()+x
    J = J.ravel()+x
    D=np.concatenate(mats,axis=0)
    f=sparse.csr_matrix((D.ravel(),(I.ravel(),J.ravel())))
    return f

In this modest example the 2nd version is 2x faster; the first scales linearly with the length of the list; the 2nd is almost independent of its length.

In [158]: timeit foo1(mats)
1000 loops, best of 3: 1.3 ms per loop

In [160]: timeit foo3(mats)
1000 loops, best of 3: 653 µs per loop

If you want to add B to A with the upper left-hand corner of B going to index (r, c) in A, you can do it using the index and the shape attribute of B:

A[r:r+B.shape[0], c:c+B.shape[1]] += B

If you want to just set the elements (overwrite instead of adding), replace += with =. In your particular example:

>>> A = np.zeros((5, 6), dtype=int)
>>> B = np.r_[np.arange(2, 10), 3].reshape(3, 3)

>>> r, c = 1, 2

>>> A[r:r+B.shape[0], c:c+B.shape[1]] += B
>>> A
array([[0, 0, 0, 0, 0, 0],
       [0, 0, 2, 3, 4, 0],
       [0, 0, 5, 6, 7, 0],
       [0, 0, 8, 9, 3, 0],
       [0, 0, 0, 0, 0, 0]])

The indexing operation produces a view into A since it is simple indexing, meaning that the data is not copied, which makes the operation fairly efficient for large arrays.

It probably is more useful to work with the transpose of W rather than W itself, both for human-readability and to facilitate writing the code. This means that the entries that affect each S_i are grouped together in one of the inner parentheses of W, i.e. in a row of W rather than a column as you have it now.

Then, S_i = np.array[S[j,:] for j in np.shape(S)[0] if W_T[i,j] == 1], where W_T is the transpose of W. If you need/want to stick with W as is, you need to reverse the indices i and j.

As for the outer loop, you could try to nest this in another similar comprehension without an if statement--however this might be awkward since you aren't actually building one output matrix (the S_i can easily be different dimensions, unless you're somehow guaranteed to have the same number of 1s in every column of W). This in fact raises the question of what you want--a list of these arrays S_i? Otherwise if they are separate variables as you have it written, there's no good way to refer to them in a generalizable way as they don't have indices.