矩阵元素乘法与移位列

Question

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

逐元素乘法定义如下：

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

Answer 1

针对列数运行 for 循环并使用 np.roll() 围绕轴 =1，移动列并进行矩阵乘法。

请参阅此参考中接受的answer。 希望这会有所帮助。

Answer 2

我们可以将A 与它自己的切片之一连接起来，然后得到那些滑动窗口。要获得这些窗口，我们可以利用基于scikit-image's view_as_windows 的np.lib.stride_tricks.as_strided。然后，将这些窗口与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]]])

我们还可以利用multi-cores 和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

---

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

如果你真的从基于跨步视图的方法中挤出最好的东西，请使用基于 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 似乎更好，对于大于这些的东西，基于跨步的应该是要走的路。

Answer 3

我实际上可以用 2 列从其两侧填充数组（以获得 2x5 数组） 并以'b'作为内核运行conv2，我认为它更有效

Answer 4

请注意view_as_windows/as_strided。尽管这些函数很简洁，但知道它们具有相当明显的恒定开销是很有用的。以下是 @Divakar 的基于 view_as_windows 的解决方案 (vaw) 与我的基于复制重塑的方法之间的比较。

如您所见，vaw 在中小型操作数上的速度不是很快，仅在数组大小 30x30 以上时才开始发光。

代码：

from simple_benchmark import BenchmarkBuilder, MultiArgument
import numpy as np
from skimage.util.shape import view_as_windows

B = BenchmarkBuilder()

@B.add_function()
def vaw(A,B):
    A1 = np.concatenate((A,A[:,:-1]),axis=1)
    w = view_as_windows(A1,A.shape)[0]
    return w*B

@B.add_function()
def pp(A,B):
    m,n = A.shape
    aux = np.empty((n,m,2*n),A.dtype)
    AA = np.concatenate([A,A],1)
    aux.reshape(-1)[:-n].reshape(n,-1)[...] = AA.reshape(-1)[:-1]
    return aux[...,:n]*B

@B.add_arguments('array size')
def argument_provider():
    for exp in range(4, 16):
        dim_size = int(1.4**exp)
        a = np.random.rand(dim_size,dim_size)
        b = np.random.rand(dim_size,dim_size)
        yield dim_size, MultiArgument([a,b])

r = B.run()
r.plot()

import pylab
pylab.savefig('vaw.png')