在 numpy 中创建 tilted/successively 移位矩阵的最快方法

Fastest way to create a tilted/successively shifted matrix in numpy

我需要创建一个如下所示的 (w,N)-矩阵:

w//2............N-1,N-1
.             \     N-1
.              \    N-1
.               \   N-1
1...............N-1,N-1
0...................N-1
00..................N-2
. \                 N-3
.  \                .
.   \               .
000000..............N-w//2

这是一个 (w,N) 矩阵,w 为奇数。中间行是从 0 到 N 的范围。对于中间行上方的每个行索引,该行像 scipy.ndimage.shift(mode='nearest') 一样向左移动,对于中间行下方的每一行,它向右移动同样的方法。

N 通常在 10^4 左右,w 通常在 10 到 10^2 之间。

我想出了两种方法:

from scipy.ndimage import shift
middle = np.arange(0, N)
final = np.vstack(
    [shift(middle, i, mode='nearest') for i in range(-w//2, 0)] + 
    [middle] + 
    [shift(middle, i, mode='nearest') for i in range(1, w//2)] ) 

需要 0.035 秒才能达到 运行。

np.vstack([
        np.maximum(
            0,
            np.minimum(
                N-1,
                np.arange(-step, N-step)
            )
        )
        for step in range(-w//2, w//2)
    ])

需要 0.021 秒才能达到 运行。

这些数字是 N=10^3 和 w=21。

我真的很想尽可能降低这些数字,最好降低到 1 毫秒左右。

我尝试了多处理,但这并没有真正帮助,开销太大而无法从并发中获得一些东西。我也知道我可以将此结果存储在某个地方,但这需要此函数的调用者进行重大更改,因此稍后会完成。

是否有任何数学关系可以表示这样的tilt/shift操作?我想不出一个,但如果有的话,numpy 可能会利用它来击败我的结果。

是的,有什么让我的代码更快的想法吗?

初始化一个具有适当形状和水平值的数组,从 0N(含)

w, N = 11, 10
arr = np.empty(shape= [w, N], dtype= int)
 
arr[:] = np.arange(N)
arr

>>>   [[0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.]]

从每一行中减去一个适当的值

arr += np.arange(w).reshape([-1, 1])[::-1] - (1+w//2)
arr

>>>   [[ 5.,  6.,  7.,  8.,  9., 10., 11., 12., 13., 14.],
       [ 4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12., 13.],
       [ 3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12.],
       [ 2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11.],
       [ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.],
       [ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9.],
       [-1.,  0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.],
       [-2., -1.,  0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.],
       [-3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.,  6.],
       [-4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.],
       [-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.]]

其中值交叉限制值重新分配给他们的限制值

arr[arr<0] = 0
arr[arr>N-1] = N-1

arr

>>>   [[5., 6., 7., 8., 9., 9., 9., 9., 9., 9.],
       [4., 5., 6., 7., 8., 9., 9., 9., 9., 9.],
       [3., 4., 5., 6., 7., 8., 9., 9., 9., 9.],
       [2., 3., 4., 5., 6., 7., 8., 9., 9., 9.],
       [1., 2., 3., 4., 5., 6., 7., 8., 9., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 0., 1., 2., 3., 4., 5., 6., 7., 8.],
       [0., 0., 0., 1., 2., 3., 4., 5., 6., 7.],
       [0., 0., 0., 0., 1., 2., 3., 4., 5., 6.],
       [0., 0., 0., 0., 0., 1., 2., 3., 4., 5.],
       [0., 0., 0., 0., 0., 0., 1., 2., 3., 4.]]

Edit

尝试为脚本计时

import timeit

script = '''
w, N   = 21, 10**3
arr    = np.empty(shape= [w, N], dtype= int)
arr[:] = np.arange(N)
arr   += np.arange(w).reshape([-1, 1])[::-1] - (1+w//2)

arr[arr<0] = 0
arr[arr>N-1] = N-1
'''

time = timeit.timeit(script, number= 100000, setup= 'import numpy as np') / 100000

time

>>> 0.00019059010320999733    # 0.19 ms