如何在numpy中高效地实现x[i][j] = y[i+j]?

Question

如何在numpy中高效地实现x[i][j] = y[i+j]?

4

假设x是一个形状为(A,B)的矩阵，y是一个大小为A+B-1的数组。

for i in range(A):
    for j in range(B):
        x[i][j] = y[i+j]

如何使用numpy中的函数高效实现等价代码？

- codebomb

看起来 y.reshape 是你想要的，这基本上只是将一个一维数组转换为二维数组。 - Andrew

通常情况下，在使用numpy数组进行元素访问时，你应该使用x[i,j]而不是x[i][j]，以获得更快的访问速度。 - JoshAdel

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Divakar · Accepted Answer

方法一 使用Scipy的Hankel函数 -

from scipy.linalg import hankel

x = hankel(y[:A],y[A-1:]

方法二 使用NumPy广播 -

（本段主要介绍使用 NumPy 广播的方法）

x = y[np.arange(A)[:,None] + np.arange(B)]

方法三 使用 NumPy步幅技术 -

n = y.strides[0]
x = np.lib.stride_tricks.as_strided(y, shape=(A,B), strides=(n,n))

运行时测试 -

In [93]: def original_app(y,A,B):
    ...:     x = np.zeros((A,B))
    ...:     for i in range(A):
    ...:         for j in range(B):
    ...:             x[i][j] = y[i+j]
    ...:     return x
    ...: 
    ...: def strided_method(y,A,B):
    ...:     n = y.strides[0]
    ...:     return np.lib.stride_tricks.as_strided(y, shape=(A,B), strides=(n,n))
    ...: 

In [94]: # Inputs
    ...: A,B = 100,100
    ...: y = np.random.rand(A+B-1)
    ...: 

In [95]: np.allclose(original_app(y,A,B),hankel(y[:A],y[A-1:]))
Out[95]: True

In [96]: np.allclose(original_app(y,A,B),y[np.arange(A)[:,None] + np.arange(B)])
Out[96]: True

In [97]: np.allclose(original_app(y,A,B),strided_method(y,A,B))
Out[97]: True

In [98]: %timeit original_app(y,A,B)
100 loops, best of 3: 5.29 ms per loop

In [99]: %timeit hankel(y[:A],y[A-1:])
10000 loops, best of 3: 114 µs per loop

In [100]: %timeit y[np.arange(A)[:,None] + np.arange(B)]
10000 loops, best of 3: 60.5 µs per loop

In [101]: %timeit strided_method(y,A,B)
10000 loops, best of 3: 22.4 µs per loop

基于strides的其他方法 -

似乎在一些基于图像处理的模块中已经使用了strides技术： extract_patches 和 view_as_windows。因此，我们还有两种方法 -

from skimage.util.shape import view_as_windows
from sklearn.feature_extraction.image import extract_patches

x = extract_patches(y,(B))
x = view_as_windows(y,(B))

In [151]: np.allclose(original_app(y,A,B),extract_patches(y,(B)))
Out[151]: True

In [152]: np.allclose(original_app(y,A,B),view_as_windows(y,(B)))
Out[152]: True

In [153]: %timeit extract_patches(y,(B))
10000 loops, best of 3: 62.4 µs per loop

In [154]: %timeit view_as_windows(y,(B))
10000 loops, best of 3: 108 µs per loop