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)

What exactly do you do in "do things with sub matrix"? If you do not specify what are you doing, we cannot help you vectorize your code. — Viktor Kerkez, Aug 15 '13 at 21:22
Jim: No, but I imagine there's cross over. Brionius: No, that the old matrix the submatrix is extracted from. Viktor: I intentionally left it vague, because I'm looking for general guidance on this type of problem, of which there are several, not looking for a specific fix. But for example see edits. — James Jenkinson, Aug 15 '13 at 22:05
Where do the `dz` elements go? Stackoverflow is really center on discrete examples and specific fixes. — Daniel, Aug 15 '13 at 22:21
They're used to solve various further operations that can easily be accomplished with simple numpy operations. This method is what's slowing everything down. — James Jenkinson, Aug 15 '13 at 22:27

score 3 · Answer 1 · answered Aug 16 '13 at 15:05

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)

score 3 · Accepted Answer · answered Aug 16 '13 at 15:25

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.

This looks like a neater solution for this specific problem. — James Jenkinson, Aug 17 '13 at 09:00

score 2 · Answer 3 · answered Aug 15 '13 at 21:21

2

Use scipy instead for image processing operations:

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

answered Aug 15 '13 at 21:21

Jiminion

5,080
1
31
54

score 1 · Answer 4 · answered Aug 16 '13 at 04:34

1

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]])

answered Aug 16 '13 at 04:34

Saullo G. P. Castro

56,802
26
179
234

Numpy submatrix operations

4 Answers4