I'm not sure if I am addressing your question on vectorizing directly, but consider the following code, which does (maybe?) what your asking. import scipy from numpy import reshape,ones, zeros, arange, array A=reshape(arange(100),[10,10]) nr,nc=A.shape B=zeros(A.shape) #initialize array #calculate averages filt=ones([3,3])/9.0 #the weighting matrix B[1:-1,1:-1]=(scipy.signal.convolve2d(A,filt, mode='same'))[1:-1,1:-1] #initialize variance matrix for interior elements B_var =zeros(array(A.shape)-2) for i in arange(-1,2): for j in arange(-1,2): B_var+=(A[(1+i):(nr-1+i),(1+j):(nc-1+j)]-B[1:-1,1:-1])**2 B_var/=9.0 #variance of 8x8 interior elements This may or may not be the best or fastest way to get your answer. In particular, it may be just as fast to use the same loop structure for B as for B_var since convolve2D is not known as a speed demon. On the whole, though, I usually like to use scipy/numpy functions when available. If you need to do this calculation a lot, it would be pretty straightforward to parallelize with c or cuda. You'd just need to take care with minimizing memory accesses (by doing things like calculating average and variance after doing one fetch of the 9 elements).