Numpy 子矩阵运算答案

【问题标题】：Numpy submatrix operationsNumpy 子矩阵运算
【发布时间】：2013-08-18 04:36:48
【问题描述】：

我正在寻找一种有效的方法来对更大的矩阵执行子矩阵运算，而无需使用 for 循环。

我目前正在进行操作（对于 3x3 窗口）：

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

这比使用 numpy 矩阵的正常操作要慢得多。有什么 numpy 可以解决这个问题，还是希望太高？

编辑：具体例子

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

内循环：

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

【问题讨论】：

您是否尝试使用 numpy 进行图像处理？
第4行的matrix应该是newMatrix吗？
“用子矩阵做事”到底是做什么的？如果您不指定您在做什么，我们将无法帮助您矢量化您的代码。
吉姆：不，但我想有交叉。 Brionius：不，是从旧矩阵中提取子矩阵。 Viktor：我故意含糊其辞，因为我正在寻找有关此类问题的一般指导，其中有几个，而不是寻找特定的解决方法。但例如查看编辑。
dz 元素去哪了？ Stackoverflow 真正以离散示例和特定修复为中心。

标签： python numpy submatrix

【解决方案1】：

在我看来，您正在尝试做一个简单的卷积？

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

没有详细说明这个标志，但希望这能让你走上正轨。

【讨论】：

这看起来像是针对这个特定问题的更简洁的解决方案。

【解决方案2】：

我在 numpy.lib.stride_tricks 中找到了解决方案

from numpy.lib.stride_tricks import as_strided

在方法中：

    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)

可能有更好的解决方案，但这对于我所追求的而言足够简单和通用。这在 3.41 秒内完成了 1600x1200 矩阵，而使用 for 循环则为 188.47 秒。

（如果您有更好的解决方案，请随时提供）

【讨论】：

【解决方案3】：

看来你可以用np.ix_，看这个例子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]])

【讨论】：

【解决方案4】：

使用 scipy 代替图像处理操作：

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

【讨论】：