将 a[i+j] 数组映射到 A[i,j] 矩阵

Question

我有一个数组a，并想创建一个新矩阵A，用于A[i,j] = a[i+j]。示例：

import numpy as np

a = np.random.rand(3)
A = np.zeros((2,2));

for i in range(2):
    for j in range(2):
        A[i,j] = a[i+j]

有没有办法在没有 for 循环的情况下做到这一点？ （使用 numpy）

Answer 1

使用stride_tricks.as_strided

这将是stride_tricks 的完美用例：

from np.lib.stride_tricks import as_strided

将步幅设置为 (8, 8)（即 (1, 1) 就插槽而言）。这样，我们基本上将结果数组A 映射为i, j -> k = i + j。更详细的描述是：我们将每个i, j 对映射到一个自然数k，由步幅定义为k = i*s_i + j*s_j，其中s_i 和s_j 是步幅设置为1s，ie k = i + j。这样就得到了想要的结果：A[i, j] = a[k] = a[i + j]。

>>> a
array([0.53954179, 0.51789927, 0.33982179])

>>> A = as_strided(a, shape=(2,2), strides=(8,8))
array([[0.53954179, 0.51789927],
       [0.51789927, 0.33982179]])

其他注意事项

更通用的解决方案是从a 的元数据中获取形状和步幅。

• A 的形状由(len(a)//2+1,)*2 给出。

• 正如@Daniel F 所指出的，内存插槽大小并不总是等于8，这确实取决于您的数组的dtype。最好从a 的步幅定义strides，而不是：a.strides*2。

这归结为：

>>> A = as_strided(a, shape=(len(a)//2+1,)*2, strides=a.strides*2)

---

使用网格索引

或者，您可以构建一个坐标网格（您可以使用itertools.product 这样做）然后将适当的值从a 复制到A：

i, j = np.array(list(product(range(2), range(2)))).T

然后初始化A和copy：

>>> A = np.zeros((2,2))
>>> A[i, j] = a[i + j]

然而，与as_strided 方法相比，这将使内存使用量翻倍。

Answer 2

我认为以下内容比使用stride_tricks 更具可读性。 tril_indices,triu_indices 为您提供下/上三角形的索引。所以你可以直接给他们写信。

A = np.zeros((2,2))

A[np.tril_indices(2)] = a
A[np.triu_indices(2)] = a

如果你担心我可能效率低下，因为我们写了两次对角线。我们可以写上三角形，将上面的非对角值复制到下面的非对角值：

A = np.zeros((2,2))
A[np.tril_indices(2)] = a
A[np.triu_indices(2,1)] = B[np.tril_indices(2,-1)]

Answer 3

In [40]: a=np.random.randint(1,100, (9,))
In [41]: a
Out[41]: array([ 6, 35, 69, 60, 63, 51, 72, 57, 22])

可广播的索引效果很好：

In [42]: i,j=np.arange(3)[:,None], np.arange(3)
In [43]: i,j
Out[43]: 
(array([[0],
        [1],
        [2]]),
 array([0, 1, 2]))
In [44]: i+j
Out[44]: 
array([[0, 1, 2],
       [1, 2, 3],
       [2, 3, 4]])

将其应用于一维数组：

In [45]: a[i+j].reshape(3,3)
Out[45]: 
array([[ 6, 35, 69],
       [35, 69, 60],
       [69, 60, 63]])

它们也可以用于分配

A[i,j] = a[i+j]

In [46]: _[i,j]
Out[46]: 
array([[ 6, 35, 69],
       [35, 69, 60],
       [69, 60, 63]])

as_strided 工作正常，但更难理解，并且容易出错。最近的版本鼓励我们使用sliding_window_view。

In [47]: np.lib.stride_tricks.as_strided(a, shape=(3,3), strides=(8,8))
Out[47]: 
array([[ 6, 35, 69],
       [35, 69, 60],
       [69, 60, 63]])

只要跨步数组只是“读取”，用作视图就更快，但后续操作可能需要复制。

其实对于这个小案例来说，视图生成并没有更快：

In [48]: timeit np.lib.stride_tricks.as_strided(a, shape=(3,3), strides=(8,8))
12.4 µs ± 43.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [49]: %%timeit
    ...: i,j=np.arange(3)[:,None], np.arange(3)
    ...: a[i+j].reshape(3,3)
    ...: 
    ...: 
8.9 µs ± 189 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

这些也产生所需的i,j：

np.ix_(range(3),range(3))
np.array(list(itertools.product(range(3), range(3)))).T