使用numpy根据数组的条件索引创建矩阵

Question

我希望创建一个基于长度为n + m的数组元素来构建大小为n by m矩阵。

这里可以使用简单的双重for循环来完成，但我希望有一种更快捷的解决方案。该矩阵相对较小。

n = 4
m = 6
s = n + m

array = np.arange(s)  # note: arange is only for example. real array varies.
matrix = np.zeros((n,m))

for i in range(n):
    for j in range (m):
        matrix[i,j] = array[i+j]

我发现理解速度比双重for循环快。

matrix3 = [[array[i+j] for i in range(m)] for j in range(n)]

有更快的方法吗？

另一个额外的奖励是使用模运算符。我实际上只需要满足 i+j % 2 == 0 的索引。在双重循环中，模运算方法似乎要快一些，但这可能对于通过numpy生成此矩阵来说不方便或不明智。

如果不这样做也可以，因为矩阵乘法将在之后发生，必要的元素无论如何都将被乘以零。我提到模除仅是出于考虑到这可能导致更快的解决方案。

针对这个MWE（最小化工作示例）

for i in range(n):
    for j in range (m):
        if (i + j) % 2 == 0:
            matrix[i,j] = array[i+j]

注意:

我希望使用numpy解决方案，并且认为这将是最快的，但只要比纯Python双重循环更快，任何纯Python（包括numpy/scipy）解决方案都可以接受

动机:

我正在尝试从双重循环中消除所有对数组的依赖项，以便我可以使用广播而不是双重循环。这是剩下的唯一一个数组

- Charlie Crown

---

arr[:, None]+ arr - hpaulj

这似乎不起作用。我不知道它如何知道要制作什么形状？最终结果不会是一个n乘m的矩阵，对吗？这样做不会制作一个（n+m）乘（n+m）的矩阵吗？ - Charlie Crown

Answer 1

您可以创建一个汉克尔矩阵：

>>> from scipy.linalg import hankel
>>> matrix = hankel(array[0:n], array[n:s])
>>> matrix
array([[0, 1, 2, 3, 4, 6],
       [1, 2, 3, 4, 6, 7],
       [2, 3, 4, 6, 7, 8],
       [3, 4, 6, 7, 8, 9]])

如果您绝对想要将满足 (i+j)%2==1 的元素设置为零，可以这样做（原始帖子）：

>>> matrix[::2, 1::2] = 0
>>> matrix[1::2, ::2] = 0
>>> matrix
array([[0, 0, 2, 0, 4, 0],
       [0, 2, 0, 4, 0, 7],
       [2, 0, 4, 0, 7, 0],
       [0, 4, 0, 7, 0, 9]])

您还可以将array的每个其他值设置为零，然后构建的矩阵将在所需位置处有零：

>>> array[1::2]=0
>>> hankel(array[0:n], array[n:s])
array([[0, 0, 2, 0, 4, 6],
       [0, 2, 0, 4, 6, 0],
       [2, 0, 4, 6, 0, 8],
       [0, 4, 6, 0, 8, 0]])

- bousof

---

如果我没记错的话，那些特殊的矩阵函数虽然方便，但并不是很快。 - Paul Panzer

嗨，保罗。如果你继续对我的答案提供反馈，我会称呼你为大师！ :) 我认为它应该很快，但当然这取决于scipy.linalg.toeplitz是否由实习生完成。你有关于这些函数速度缓慢的参考资料吗？ - bousof

例如，查看 https://dev59.com/TqHia4cB1Zd3GeqPSmJ3#43735362 在那里浏览评论，可能是scipy已经被修补过了，这些函数不再缓慢。另外一条备注：您可以直接使用“hankel”而不是翻转“toeplitz”。 - Paul Panzer

即使您最后进行复制，as_strided 的速度仍然比我花哨的索引方法要快。 - Paul Panzer

好的，我已经阅读了整个评论线程，并且特殊矩阵确实由@Divakar本人修补（并已被合并到scipy中）。所以就性能而言，你应该没问题了。不能再点赞因为之前已经点过了。 - Paul Panzer

@PaulPanzer 没问题。很高兴知道他们提高了性能。对于科学代码的可读性，使用专用矩阵名称创建矩阵甚至更加重要，我个人认为。 - bousof

Answer 2

您可以使用高级索引进入数组。为了提高效率，您可以在模板数组中零化奇数位置。

np.where(np.arange(m+n)&1,0,array)[sum(np.ogrid[:n,:m])]
# array([[0, 0, 2, 0, 4, 0],
#        [0, 2, 0, 4, 0, 6],
#        [2, 0, 4, 0, 6, 0],
#        [0, 4, 0, 6, 0, 8]])

或者（更快）

template = np.where(np.arange(m+n)&1,0,array)
np.lib.stride_tricks.as_strided(template,(n,m),2*template.strides)

这是一种“压缩”视图，如果您需要修改条目，则必须制作副本（这仍然会更快）。

- Paul Panzer

---

@bousof <strided_view>.copy()将创建一个连续布局的数组版本。 - Paul Panzer

Answer 3

创建表格的更简单方法是：

1. Define a function:

    def tVal(r, c):
        sm = r + c
        return np.where(sm % 2 == 0, sm, 0)
   

2. Use it as an argument of np.fromfunction:

    arr = np.fromfunction(tVal, (n, m))
   

针对您的目标形状（6 * 4），结果为：

（保留HTML标记）。

array([[0., 0., 2., 0., 4., 0.],
       [0., 2., 0., 4., 0., 6.],
       [2., 0., 4., 0., 6., 0.],
       [0., 4., 0., 6., 0., 8.]])

请注意，tVal 并不是为每个数组元素单独调用的。它只被调用一次，使用两个数组（r 和 c）作为目标数组，填充各自的参数以处理这些数组（而不是每个单元格索引的单个值）。

因此，此函数必须包含where，而不是特定单元格的r和c值的if。

关于变量名称的备注：在Numpy中，matrix是一个ndarray的子类型，因此使用相同名称的变量是不好的做法。最好使用其他名称，就像我在示例中所做的那样。

- Valdi_Bo

---

这对于原始数组由np.arange生成的情况可以工作，但这不是真实情况，我只是使用这个来说明最终矩阵应该是什么。一般来说，数组中的值变化很大。这个方法是否可以轻松地应用于这种情况？ - Charlie Crown

只要您能编写操作两个填充有值（索引）的数组中每个元素的函数，就可以使用fromfuncttion。在更复杂的情况下，您必须迭代数组元素。请阅读关于nditer的相关内容。它是一个在Numpy数组上操作的迭代器，但也允许访问当前元素的索引和对元素本身进行写入访问。 - Valdi_Bo

Answer 4

我会直接在numpy级别上处理：

matrix = np.arange(n * m).reshape(n,m)
matrix = matrix // m + matrix % m             # matrix // m is i and matrix % m is j

对于n, m = 4, 6，它按预期给出：

array([[0, 1, 2, 3, 4, 5],
       [1, 2, 3, 4, 5, 6],
       [2, 3, 4, 5, 6, 7],
       [3, 4, 5, 6, 7, 8]], dtype=int32)

- Serge Ballesta

Answer 5

你的第一个例子：

In [30]: arr=np.arange(24)                                                              
In [31]: [[arr[i+j] for i in range(6)] for j in range(4)]                               
Out[31]: 
[[0, 1, 2, 3, 4, 5],
 [1, 2, 3, 4, 5, 6],
 [2, 3, 4, 5, 6, 7],
 [3, 4, 5, 6, 7, 8]]

利用“广播”的方式：

In [32]: np.arange(4)[:,None]+np.arange(6)                                              
Out[32]: 
array([[0, 1, 2, 3, 4, 5],
       [1, 2, 3, 4, 5, 6],
       [2, 3, 4, 5, 6, 7],
       [3, 4, 5, 6, 7, 8]])

外部的i循环被替换为一个(n,1)数组；内部的j循环被替换为(m,)数组；两者结合的结果是一个(n,m)数组。

您更为详细的情况：

In [35]: arr = np.arange(24) 
    ...: res = np.zeros((4,6),int) 
    ...: for i in range(4): 
    ...:     for j in range(6): 
    ...:         if (i+j)%2 ==0: 
    ...:             res[i,j] = arr[i+j] 
    ...:                                                                                
In [36]: res                                                                            
Out[36]: 
array([[0, 0, 2, 0, 4, 0],
       [0, 2, 0, 4, 0, 6],
       [2, 0, 4, 0, 6, 0],
       [0, 4, 0, 6, 0, 8]])

这是原始版本，仅设置了偶数值。

In [37]: Out[32]                                                                        
Out[37]: 
array([[0, 1, 2, 3, 4, 5],
       [1, 2, 3, 4, 5, 6],
       [2, 3, 4, 5, 6, 7],
       [3, 4, 5, 6, 7, 8]])

找出赔率：

In [38]: Out[32]%2                                                                      
Out[38]: 
array([[0, 1, 0, 1, 0, 1],
       [1, 0, 1, 0, 1, 0],
       [0, 1, 0, 1, 0, 1],
       [1, 0, 1, 0, 1, 0]])

乘法：

In [39]: Out[32]*(Out[32]%2==0)                                                         
Out[39]: 
array([[0, 0, 2, 0, 4, 0],
       [0, 2, 0, 4, 0, 6],
       [2, 0, 4, 0, 6, 0],
       [0, 4, 0, 6, 0, 8]])

一般来说，为了最大程度地利用 numpy，我尝试看到整体模式。这就是小例子特别有价值的地方。

- hpaulj