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

【问题标题】：Use numpy to create matrix based on conditional indices of an array使用 numpy 根据数组的条件索引创建矩阵
【发布时间】：2020-07-08 04:50:47
【问题描述】：

我想基于 n + m 长度数组的元素创建一个 n by m 矩阵。

这里有一个简单的 double 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]

我发现理解比双循环更快

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

有没有更快的方法？

另一个好处是合并模运算符。我实际上只需要i+j % 2 == 0 所在的索引。在双 for 循环中，取模方法似乎要快一点，但这对于通过 numpy 生成此矩阵可能不方便或权宜之计。

最好不要这样做，因为矩阵乘法会发生在之后，必要的元素无论如何都会乘以零。我提到模数只是在这种情况下，这会导致更快的解决方案。

对于这个 MWE

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

注意：

我要求一个 numpy 解决方案，假设 numpy 会最快，但任何纯 python（包括 numpy/scipy）解决方案都可以，只要它比纯 python double for 循环更快

动机：

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

【问题讨论】：

arr[:, None]+ arr
这似乎不起作用。我不明白它怎么知道要做成什么形状？最终结果不会是 n x m 矩阵吗？这不会构成 (n+m) x (n+m) 矩阵吗？

标签： python python-3.x numpy matrix

【解决方案1】：

您可以对array 使用高级索引。为了提高效率，您可以将模板数组中已经存在的奇数位置归零。

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)

这是一个“压缩”视图，如果您需要修改条目，您必须进行复制（它仍然会更快）。

【讨论】：

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

【解决方案2】：

您可以创建一个hankel 矩阵：

>>> 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 处的元素设置为零，您可以这样做 (original post)：

>>> 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]])

【讨论】：

IIRC 那些特殊的矩阵函数虽然方便，但并不是很快。
你保罗。如果你继续给我反馈我的答案，我会称你为大师！ :) 我认为它应该很快，但当然这取决于 scipy.linalg.toeplitz 是否由实习生完成。您对这些功能的缓慢有一些参考吗？
参见，例如stackoverflow.com/a/43735362/7207392 看看那里的 cmets，可能是 scipy 已被修补，并且功能不再缓慢。再备注：可以直接使用hankel，不用翻转toeplitz
即使你最后复制as_strided至少比我花哨的索引方法快。
好的，我通读了整个评论线程，@Divakar 本人确实修补了特殊矩阵（并且 PR 接受了 scipy）。所以你应该表现得很好。无法 +1，因为已经这样做了。

【解决方案3】：

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

定义一个函数：

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

将其用作np.fromfunction的参数：
```
 arr = np.fromfunction(tVal, (n, m))
```

对于您的目标形状 (6 * 4)，结果是：

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 实际上不是为每个数组单独调用元素。相反，它只被称为 once，具有 2 个 arrays（r 和 c）形状作为目标数组，为每个单元格填充各自的参数。所以这个函数对这些数组进行操作（而不是每个数组的单个值单元格索引）。

这就是为什么这个函数必须包含 where，而不是 if 用于 r 和 c 特定单元格的值。

关于变量名的说明：matrix is a class in Numpy （ndarray 的一个子类型），所以使用变量不是一个好习惯同名。使用其他名称，就像我在示例中所做的那样。

【讨论】：

这适用于原始数组生成为 np.arange 的情况，但这不是真实情况，我只是使用它以便很容易判断最终矩阵应该是什么。通常，数组中的值变化很大。这可以轻松应用于该场景吗？
您可以使用 fromfuncttion，只要您能够编写在 2 个数组上运行的函数，每个元素的值（索引）填充。在更复杂的情况下，您必须遍历数组元素。阅读nditer。它是一个在 Numpy 数组上操作的 iterator，但也允许访问当前元素的索引和 write 访问元素本身。跨度>

【解决方案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)

【讨论】：

【解决方案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，我会尝试查看整体模式。这就是小例子特别有价值的地方。

【讨论】：