python - 带移位列的矩阵元素明智乘法

Question

假设我有两个数组，A 和 B。

元素乘法定义如下:

我想以类似卷积的方式进行逐元素乘法，即将每一列向右移动一步，例如，第 1 列现在是第 2 列，第 3 列现在是第 1 列。这应该产生一个( 2 x 3 x 3 )数组(所有 3 种可能性的 2x3 矩阵)

Answer 1

我们可以将 A 与它自己的切片之一连接起来，然后得到那些滑动窗口。要获得这些窗口，我们可以利用 np.lib.stride_tricks.as_strided基于 scikit-image's view_as_windows .然后，将这些窗口与 B 相乘以获得最终输出。 More info on use of as_strided based view_as_windows .

因此，我们将有一个像这样的矢量化解决方案 -

In [70]: from skimage.util.shape import view_as_windows

In [71]: A1 = np.concatenate((A,A[:,:-1]),axis=1)

In [74]: view_as_windows(A1,A.shape)[0]*B
Out[74]: 
array([[[1, 0, 3],
        [0, 0, 6]],

       [[2, 0, 1],
        [0, 0, 4]],

       [[3, 0, 2],
        [0, 0, 5]]])

我们还可以利用多核 numexpr module对于 broadcasted-multiplication 的最后一步，这在较大的数组上应该更好。因此，对于示例案例，它将是 -

In [53]: import numexpr as ne

In [54]: w = view_as_windows(A1,A.shape)[0]

In [55]: ne.evaluate('w*B')
Out[55]: 
array([[[1, 0, 3],
        [0, 0, 6]],

       [[2, 0, 1],
        [0, 0, 4]],

       [[3, 0, 2],
        [0, 0, 5]]])

比较所提出的两种方法的大型阵列计时 -

In [56]: A = np.random.rand(500,500)
    ...: B = np.random.rand(500,500)

In [57]: A1 = np.concatenate((A,A[:,:-1]),axis=1)
    ...: w = view_as_windows(A1,A.shape)[0]

In [58]: %timeit w*B
    ...: %timeit ne.evaluate('w*B')
1 loop, best of 3: 422 ms per loop
1 loop, best of 3: 228 ms per loop

---

挤出最好的基于跨步的方法

如果您真的从基于跨 View 的方法中挤出最好的，请使用原始的 np.lib.stride_tricks.as_strided基于一个避免功能开销的 view_as_windows -

def vaw_with_as_strided(A,B):
    A1 = np.concatenate((A,A[:,:-1]),axis=1)
    s0,s1 = A1.strides
    S = (A.shape[1],)+A.shape
    w = np.lib.stride_tricks.as_strided(A1,shape=S,strides=(s1,s0,s1))
    return w*B

与 @Paul Panzer's array-assignment 比较基于一个，交叉似乎在 19x19 形状的数组 -

In [33]: n = 18
    ...: A = np.random.rand(n,n)
    ...: B = np.random.rand(n,n)

In [34]: %timeit vaw_with_as_strided(A,B)
    ...: %timeit pp(A,B)
10000 loops, best of 3: 22.4 µs per loop
10000 loops, best of 3: 21.4 µs per loop

In [35]: n = 19
    ...: A = np.random.rand(n,n)
    ...: B = np.random.rand(n,n)

In [36]: %timeit vaw_with_as_strided(A,B)
    ...: %timeit pp(A,B)
10000 loops, best of 3: 24.5 µs per loop
10000 loops, best of 3: 24.5 µs per loop

因此，对于任何小于 19x19 的东西，array-assignment 似乎更好，而对于比那些更大的东西，strided-based 应该是可行的方法。