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

Question

让x是一个形状为（A，B）的矩阵，并且是大小为A + B-1的阵列。

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

如何在numpy中使用函数有效实现等效代码？

Answer 1

方法＃1使用Scipy's hankel -

from scipy.linalg import hankel 

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

方法＃2使用NumPy broadcasting -

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

方法3使用NumPy strides technique -

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