首页 > 解决方案 > 矩阵元素乘法与移位列

问题描述

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

逐元素乘法定义如下: 在此处输入图像描述

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

在此处输入图像描述

标签: pythonnumpymatrixlinear-algebramatrix-multiplication

解决方案


我们可以A与它自己的切片之一连接,然后得到那些滑动窗口。要获得这些窗口,我们可以利用np.lib.stride_tricks.as_strided基于scikit-image's view_as_windows. B然后,将这些窗口与最终输出相乘。有关使用based 的更多信息as_stridedview_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]]])

我们还可以利用multi-cores模块numexpr进行最后一步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

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

如果您真的从基于跨步视图的方法中挤出了最好的,请使用基于原始视图的方法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更好,对于大于那些的东西,基于跨步的应该是要走的路。


推荐阅读