python - 矩阵元素乘法与移位列
问题描述
假设我有两个数组,A 和 B。
我想以类似卷积的方式进行逐元素乘法,即将每一列向右移动一步,例如,第 1 列现在是第 2 列,第 3 列现在是第 1 列。这应该会产生 ( 2 by 3 by 3) 数组(所有 3 种可能性的 2x3 矩阵)
解决方案
我们可以A
与它自己的切片之一连接,然后得到那些滑动窗口。要获得这些窗口,我们可以利用np.lib.stride_tricks.as_strided
基于scikit-image's view_as_windows
. B
然后,将这些窗口与最终输出相乘。有关使用based 的更多信息as_strided
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
进行最后一步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
更好,对于大于那些的东西,基于跨步的应该是要走的路。
推荐阅读
- memcached - 有人可以解释一下 memcached 中的平板是什么吗
- sql-server - 我正在使用终端使用 VS Code 从现有数据库 .Net Core 3.0 生成模型,但它不起作用
- r - R sink() 消息并输出到同一个文件 - 完整性检查
- angular - 将清晰度数据网格的行滚动到顶部并以编程方式扩展它的更好方法是什么?
- c - fscanf 从头开始读取
- java - 排序和过滤
- ios - 具有典型菜单栏样式导航的应用程序是如何编码的?
- django - Django 显示选择值
- swift - 如何在不重新启动页面的情况下更新导航控制器上一页的信息?
- python-3.x - 为什么我在 Jupyter 中使用 TensorFlow 时遇到问题?