在 NUMPY 中映射 2D 延迟向量和 1D 向量

Question

我正在尝试从一维向量生成一个 2D 向量，其中元素每行沿行移动一个增量。

我希望我的输入看起来像这样：

input:
t = [t1, t2, t3, t4, t5]
out = 
[t5,  0,  0,  0,  0]
[t4, t5,  0,  0,  0]
[t3, t4, t5,  0,  0]
[t2, t3, t4, t5,  0]
[t1, t2, t3, t4, t5]
[ 0, t1, t2, t3, t4]
[ 0,  0, t1, t2, t3]
[ 0,  0,  0, t1, t2]
[ 0,  0,  0,  0, t1]

我不知道在不使用 for 循环的情况下执行此操作的方法，并且计算效率对于 im 使用它的任务很重要。有没有办法在没有 for 循环的情况下做到这一点？

这是我使用 for 循环的代码：

import numpy as np
t = np.linspace(-3, 3, 7)
z = np.zeros((2*len(t) - 1, len(t)))
diag = np.arange(len(t))
for index, val in enumerate(np.flip(t, 0)):
z[diag + index, diag] = val
print(z)

Answer 1

你在这里要求的东西被称为Toeplitz矩阵，它是：

从左到右的每个降序对角线都是恒定的矩阵

一个区别是你想要矩阵的下三角形。

你碰巧很幸运，你可以使用scipy.linalg.toeplitz来构造你的矩阵，然后np.tril访问下三角形。

import numpy as np
from scipy.linalg import toeplitz
v = np.array([1, 2, 3, 4, 5])
t = np.pad(v[::-1], (0, 4), mode='constant')

求解矩阵并访问下三角形：

np.tril(toeplitz(t, v))

还有我们的输出！

array([[5, 0, 0, 0, 0],
[4, 5, 0, 0, 0],
[3, 4, 5, 0, 0],
[2, 3, 4, 5, 0],
[1, 2, 3, 4, 5],
[0, 1, 2, 3, 4],
[0, 0, 1, 2, 3],
[0, 0, 0, 1, 2],
[0, 0, 0, 0, 1]])

要推广此方法，只需从v的形状计算t所需的填充：

v = # any one dimension array
t = np.pad(v[::-1], (0, v.shape[0]-1), mode='constant')

Answer 2

不知道非循环方法，但您可以使用roll和column_stack加快速度。

v = np.array([1, 2, 3, 4, 5])
t = np.pad(v[::-1], (0, 4), mode='constant')
np.column_stack([np.roll(t, i) for i in range(len(v))]) 
array([[5, 0, 0, 0, 0],
[4, 5, 0, 0, 0],
[3, 4, 5, 0, 0],
[2, 3, 4, 5, 0],
[1, 2, 3, 4, 5],
[0, 1, 2, 3, 4],
[0, 0, 1, 2, 3],
[0, 0, 0, 1, 2],
[0, 0, 0, 0, 1]])

Answer 3

这是一种利用strides-tricks的矢量化方法，它只是将视图使用到零填充数组中，并且视图的内存效率很高，因此性能很高 -

def map2D(a):
n = len(a)
p = np.zeros(n-1,dtype=a.dtype)
a_ext = np.r_[p,a,p]
s0 = a_ext.strides[0]
strided = np.lib.stride_tricks.as_strided
return strided(a_ext[-n:], (len(a_ext)-n+1,n), (-s0,s0), writeable=False)

示例运行 -

In [81]: a = np.array([2,5,6,7,9])
In [82]: map2D(a)
Out[82]: 
array([[9, 0, 0, 0, 0],
[7, 9, 0, 0, 0],
[6, 7, 9, 0, 0],
[5, 6, 7, 9, 0],
[2, 5, 6, 7, 9],
[0, 2, 5, 6, 7],
[0, 0, 2, 5, 6],
[0, 0, 0, 2, 5],
[0, 0, 0, 0, 2]])

如果需要输出具有自己的内存空间，请使用.copy()

5k元素数组的时序 -

In [83]: a = np.random.randint(0,9,(5000))
# From this post's soln
In [84]: %timeit map2D(a)
10000 loops, best of 3: 26.3 µs per loop
# If you need output with its own memory space
In [97]: %timeit map2D(a).copy()
10 loops, best of 3: 43.6 ms per loop
# @user3483203's soln
In [87]: %%timeit
...: t = np.pad(a[::-1], (0, len(a)-1), mode='constant')
...: out = np.tril(toeplitz(t, a))
1 loop, best of 3: 264 ms per loop
# @coldspeed's soln
In [89]: %%timeit
...: t = np.pad(a[::-1], (0, len(a)-1), mode='constant')
...: out = np.column_stack([np.roll(t, i) for i in range(len(a))])
1 loop, best of 3: 336 ms per loop