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
– Viktor Kerkez, Commented Aug 15, 2013 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
– James Jenkinson, Commented Aug 15, 2013 at 22:05
Where do the dz elements go? Stackoverflow is really center on discrete examples and specific fixes. — Daniel
– Daniel, Commented Aug 15, 2013 at 22:21

Viktor Kerkez · Accepted Answer · 2013-08-16 15:25:52Z

3

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.

answered Aug 16, 2013 at 15:25

Viktor Kerkez

46.8k13 gold badges109 silver badges88 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

James Jenkinson Over a year ago

This looks like a neater solution for this specific problem.

James Jenkinson · Accepted Answer · 2013-08-16 15:05:55Z

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)

Jiminion · Accepted Answer · 2013-08-15 21:21:01Z

2

Use scipy instead for image processing operations:

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

answered Aug 15, 2013 at 21:21

Jiminion

5,1701 gold badge35 silver badges56 bronze badges

Comments

Saullo G. P. Castro · Accepted Answer · 2013-08-16 04:34:30Z

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, 2013 at 4:34

Saullo G. P. Castro

59.4k28 gold badges191 silver badges244 bronze badges

Collectives™ on Stack Overflow

Numpy submatrix operations

4 Answers 4

1 Comment

Comments

Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

4 Answers 4

1 Comment

Comments

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related