Numpy submatrix operations

Question

I'm looking for an efficient way to perform submatrix operations over a larger matrix without resorting to for loops.

I'm currently doing the operation (for a 3x3 window):

newMatrix = numpy.zeros([numRows, numCols])
for i in range(1, numRows-1):
    for j in range(1, numCols-1):
        sub = matrix[i-1:i+2, j-1:j+2]
        newMatrix[i][j] = ... #do things with sub matrix

This is considerably slower than normal operations with numpy matrices. Is there anything numpy has to offer to solve this, or is that hoping for too much?

Edit: Specific example

xWeight = numpy.array([[-1./8, 0, 1./8], [-2./8, 0, 2./8], [-1./8, 0, 1./8]])
yWeight = numpy.array([[1./8, 2./8, 1./8], [0, 0, 0], [-1./8, -2./8, -1./8]])

Inside loop:

        dz_dx = numpy.sum(xWeight * sub)
        dz_dy = numpy.sum(yWeight * sub)

Answer 1

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.

Answer 2

I found a solution in numpy.lib.stride_tricks

from numpy.lib.stride_tricks import as_strided

In the method:

    expansion = stride.as_strided(matrix, shape = (numRows-2, numCols-2, 3, 3), strides = matrix.strides * 2)
    xWeight = numpy.array([[-1./8, 0, 1./8], [-2./8, 0, 2./8], [-1./8, 0, 1./8]])
    yWeight = numpy.array([[1./8, 2./8, 1./8], [0, 0, 0], [-1./8, -2./8, -1./8]])

    dx = xWeight * expansion
    dy = yWeight * expansion

    dx = numpy.sum(numpy.sum(dx, axis=3), axis=2)
    dy = numpy.sum(numpy.sum(dy, axis=3), axis=2)

There may well be a better solution, but this is sufficiently simple and general purpose for what I was after. This went through a 1600x1200 matrix in 3.41 seconds, vs 188.47 seconds using for loops.

(Feel free to offer said better solution, if you have it)

Answer 3

Use scipy instead for image processing operations:

http://scipy-lectures.github.io/advanced/image_processing/

Answer 4

It seems you can use np.ix_, see this example from the documentation:

a = np.arange(10).reshape(2, 5)
#array([[0, 1, 2, 3, 4],
#       [5, 6, 7, 8, 9]])

ixgrid = np.ix_([0,1], [2,4])

a[ixgrid]
#array([[2, 4],
#       [7, 9]])