Numpy methods are going to beat python loops almost always, so I am going to skip your 1.

As for 2, in this particular case the following works:

a = np.array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11],
              [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]])
a = a.reshape(2, 3, 4)
>>> a
array([[[ 0,  1,  2,  3],
        [ 4,  5,  6,  7],
        [ 8,  9, 10, 11]],

       [[12, 13, 14, 15],
        [16, 17, 18, 19],
        [20, 21, 22, 23]]])
>>> np.mean(a, axis=1)
array([[  4.,   5.,   6.,   7.],
       [ 16.,  17.,  18.,  19.]])

The trick is in the reshape. For a general case where you want blocks of n columns, the following is an option

a = a.reshape((a.shape[0], -1, n))

Your concerns in 3 are mostly unwarranted. reshape returns a view of the original array, not a copy, so the conversion to 3D only requires altering the shape and strides attributes of the array, without having to copy any of the actual data.

EDIT To be sure that reshaping does not copy the array, but returns a view, do the reshape as

a.shape = a = a.reshape((a.shape[0], -1, n))

The example in the docs goes along the lines of:

>>> a = np.arange(12).reshape(3,4)
>>> b = a.T
>>> b.shape = (12,)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: incompatible shape for a non-contiguous array

And in general there are only problems if you have been doing transpose, rollaxis, swapaxes or the like on your array.

To compute some operation slice by slice in numpy, it is very often useful to reshape your array and use extra axes.

To explain the process we'll use here: you can reshape your array, take the mean, reshape it again and take the mean again. Here I assume blocksize is 2

t = np.array([[0, 1, 2, 3, 4, 5], [0, 1, 2, 3, 4, 5],[0, 1, 2, 3, 4, 5],[0, 1, 2, 3, 4, 5],[0, 1, 2, 3, 4, 5],[0, 1, 2, 3, 4, 5],])
t = t.reshape([6, 3, 2])
t = np.mean(t, axis=2)
t = t.reshape([3, 2, 3])
np.mean(t, axis=1)

outputs

array([[ 0.5,  2.5,  4.5],
       [ 0.5,  2.5,  4.5],
       [ 0.5,  2.5,  4.5]])

Now that it's clear how this works, you can do it in one pass only:

t = t.reshape([3, 2, 3, 2])
np.mean(t, axis=(1, 3))

works too (and should be quicker since means are computed only once - I guess). I'll let you substitute height/blocksize, width/blocksize and blocksize accordingly.

See @askewcan nice remark on how to generalize this to any dimension.

mats_row=np.array([[A1,A2,A3,...,A_n]]) #Create array of matrices with shape (1,M,N,N)
mats_column=np.transpose(mats_row,(1,0,2,3)) #Make a copy with shape (M,1,N,N)
block=.5*(mats_row+mats_column) #Add, broadcasting to a (M,M,N,N) array

This works by utilizing Numpy's broadcasting capabilities. The basic idea is similar to if you added a matrix where each row was to a matrix where each column was , but broadcasting allows you to do this without ever explicitly forming these intermediates matrices.

One approach would be to reshape splitting each of those two axes into two more and find mean along the latter of those two, giving us the average values. We will keep the dimensions with keepdims=True, so as to facilitate replicating along the reduced axes later on with np.repeat for the final output.

Thus, one implementation would be -

b = a.reshape(2,2,2,2).mean((1,3), keepdims=1)
out = np.repeat(np.repeat(b,(2),axis=(1)),(2), axis=3).reshape(4,4)

Sample run -

In [17]: a
Out[17]: 
array([[ 1,  2,  3,  4],
       [ 5,  6,  7,  8],
       [ 9, 10, 11, 12],
       [13, 14, 15, 16]])

In [18]: b = a.reshape(2,2,2,2).mean((1,3), keepdims=1)

In [19]: np.repeat(np.repeat(b,(2),axis=(1)),(2), axis=3).reshape(4,4)
Out[19]: 
array([[  3.5,   3.5,   5.5,   5.5],
       [  3.5,   3.5,   5.5,   5.5],
       [ 11.5,  11.5,  13.5,  13.5],
       [ 11.5,  11.5,  13.5,  13.5]])

Generically put, the solution would look like this -

m,n = a.shape
p,q = (2,2) # Block size
b = a.reshape(m//p,p,n//q,q).mean((1,3), keepdims=1)
out = np.repeat(np.repeat(b,(p),axis=(1)),(q), axis=3).reshape(a.shape)

Performance boost

We could replace that replicating part with an initialization based section, like so -

out = np.empty((m//p,p,n//q,q),dtype=float)
out[:] = b
out.shape = a.shape

If your array arr has a length divisible by 3:

np.mean(arr.reshape(-1, 3), axis=1)

Reshaping to a higher dimensional array and then performing some form of reduce operation on one of the additional dimensions is a staple of numpy programming.

You can do:

# sample data
a=np.arange(24).reshape((6,4))

rows, cols = a.shape    
a.reshape(rows//2, 2, cols//2, 2).mean(axis=(1,-1))

Output:

array([[ 2.5,  4.5],
       [10.5, 12.5],
       [18.5, 20.5]])

If you are open to other packages, Pandas as a convenient groupby function:

out = (pd.Series(a.ravel(), 
                  index = pd.MultiIndex.from_product((pairs,pairs)))
         .groupby(level=(0,1)).mean()
         .unstack().to_numpy()
      )

Output:

array([[5.25      , 5.        , 3.5       ],
       [5.33333333, 3.11111111, 3.33333333],
       [6.5       , 3.33333333, 4.        ]])