Numpy：矩阵数组移位/按索引插入

Question

我有一个对象，我已将一个大型 for 循环方法转换为一系列矢量化 numpy 数组（大约快 50 倍）。现在我正在尝试添加一个新方法，我需要处理一个 numpy 矩阵，然后根据矩阵中的数组索引“移位”子数组内容（即插入值）。我知道我可以使用 for 循环来完成此操作，但我试图通过使用矢量数学来实现加速收益来避免这种情况。

我想知道是否有一种快速有效的方法来完成以下任务：

import numpy as np

period = [1, 2, 3]

x = [1, 10, 100]
y = [.2, .4, .6]

z = np.outer(x,y)

print(z)

结果：

[[  0.2   0.4   0.6]
 [  2.    4.    6. ]
 [ 20.   40.   60. ]]

我想移动z中的行以添加基于句点的零数作为z中的行索引，基本上如下：

[[   0.0   0.2   0.4    0.6 ]
 [   0.0   0.0   2.0    4.0    6.0 ]
 [   0.0   0.0   0.0   20.0   40.0   60.0 ]]

最终，我希望在垂直/列轴 (axis=1) 上求和。我需要一个如下所示的最终数组：

[   0.0   0.2   2.4   24.6   46.0   60.0]

Answer 1

[[   0.0   0.2   0.4    0.6 ]
 [   0.0   0.0   2.0    4.0    6.0 ]
 [   0.0   0.0   0.0   20.0   40.0   60.0 ]]

是一个参差不齐的列表。我们可以用矢量化数组魔法来构建它，至少不能用普通的东西。

为了解决这个问题，我们需要展平或填充这个结构

[[   0.0   0.2   0.4    0.6    0.0    0.0]
 [   0.0   0.0   2.0    4.0    6.0    0.0 ]
 [   0.0   0.0   0.0   20.0   40.0   60.0 ]]

或

[   0.0   0.2   0.4    0.6   0.0   0.0   2.0    4.0    6.0   0.0   0.0   0.0   20.0   40.0   60.0 ] 

sum.reduceat 让我们对平面数组的块求和，但你想要一个跳过求和。我想我可以探索展平转置。

我的第一个想法是填充数组看起来像一个对角线，[.2,2,20] 放在对角线上，[.4,4,40] 放在下一个偏移量上，依此类推。我知道sparse 可以从一个矩阵和一组偏移量构建一个矩阵，但我认为numpy 中没有这样的功能。它们都一次使用一个偏移量。

但它看起来也像stride_tricks 可以产生的那种偏移量。

让我们探索一下：

In [458]: as_strided =np.lib.index_tricks.as_strided

In [459]: Z=np.pad(z,[[0,0],[3,3]],mode='constant')
In [460]: Z
Out[460]: 
array([[  0. ,   0. ,   0. ,   0.2,   0.4,   0.6,   0. ,   0. ,   0. ],
       [  0. ,   0. ,   0. ,   2. ,   4. ,   6. ,   0. ,   0. ,   0. ],
       [  0. ,   0. ,   0. ,  20. ,  40. ,  60. ,   0. ,   0. ,   0. ]])

In [461]: Z.strides
Out[461]: (72, 8)       # prod an offset with (72+8, 8)
In [462]: as_strided(Z,shape=(3,6),strides=(80,8))
Out[462]: 
array([[  0. ,   0. ,   0. ,   0.2,   0.4,   0.6],
       [  0. ,   0. ,   2. ,   4. ,   6. ,   0. ],
       [  0. ,  20. ,  40. ,  60. ,   0. ,   0. ]])

这是我们想要的那种转变，但方向错了；所以让我们翻转Z：

In [463]: Z1=Z[::-1,:].copy()
In [464]: as_strided(Z1,shape=(3,6),strides=(80,8))
Out[464]: 
array([[  0. ,   0. ,   0. ,  20. ,  40. ,  60. ],
       [  0. ,   0. ,   2. ,   4. ,   6. ,   0. ],
       [  0. ,   0.2,   0.4,   0.6,   0. ,   0. ]])
In [465]: as_strided(Z1,shape=(3,6),strides=(80,8)).sum(0)
Out[465]: array([  0. ,   0.2,   2.4,  24.6,  46. ,  60. ])

可以将概括留给读者。

是否有任何速度优势尚不清楚。可能不是这个小案例，也许是一个非常大的案例。

---

这会清理填充和跨步

In [497]: Z=np.pad(z,[[0,0],[1,4]],mode='constant')
In [498]: Z.strides
Out[498]: (64, 8)
In [499]: as_strided(Z,shape=(3,6),strides=(64-8,8))
Out[499]: 
array([[  0. ,   0.2,   0.4,   0.6,   0. ,   0. ],
       [  0. ,   0. ,   2. ,   4. ,   6. ,   0. ],
       [  0. ,   0. ,   0. ,  20. ,  40. ,  60. ]])

---

这忽略了z 的构造方式。如果外积是问题的核心，我可能会尝试在 1d y 上跨步，并使用 x 进行加权求和。

In [553]: x=np.array([1,10,100]); y=np.array([.2,.4,.6])
In [554]: z=np.concatenate(([0,0],y[::-1],[0,0,0]))
In [555]: z
Out[555]: array([ 0. ,  0. ,  0.6,  0.4,  0.2,  0. ,  0. ,  0. ])
In [556]: Z=as_strided(z,shape=(3,6), strides=(8,8))
In [557]: Z
Out[557]: 
array([[ 0. ,  0. ,  0.6,  0.4,  0.2,  0. ],
       [ 0. ,  0.6,  0.4,  0.2,  0. ,  0. ],
       [ 0.6,  0.4,  0.2,  0. ,  0. ,  0. ]])
In [558]: np.dot(x,Z)
Out[558]: array([ 60. ,  46. ,  24.6,   2.4,   0.2,   0. ])

在这个构造中Z 是z 的一个视图，因此比前面的Z 小。但我确定dot 在将其发送到已编译代码时会复制一份。 np.einsum('i,ij',x,Z) 可能会避免这种情况，在不扩展视图的情况下对其进行编译迭代。这在处理非常大的数组时可能会有所不同。

结果是相反的，但这很容易解决。我什至可以在施工期间修复它。

Answer 2

您也可以先计算索引并立即分配：

a = np.array(
    [[0.2 ,  0.4 ,  0.6],
     [2.,    4.,    6. ],
     [20.,   40.,   60. ]])

s0, s1 = a.shape
rows = np.repeat(np.arange(s0), s1).reshape(a.shape)
cols = (np.add.outer(np.arange(0, s0), np.arange(s1)) + 1)
res = np.zeros((s0, s0 + s1))
res[rows, cols] = a
np.sum(res,axis=0)

>>> np.sum(res,axis=0)
array([  0. ,   0.2,   2.4,  24.6,  46. ,  60. ])

Answer 3

循环第一个维度有效：

a = np.array(
    [[0.2 ,  0.4 ,  0.6],
     [2.,    4.,    6. ],
     [20.,   40.,   60. ]])

​
s0, s1 = a.shape
res = np.zeros((s0, s0 + s1))
for i in range(s1):
    res[i, i + 1: i + s0 + 1] = a[i] 

>>> np.sum(res,axis=0)
array([  0. ,   0.2,   2.4,  24.6,  46. ,  60. ])