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

Question

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

我想知道是否有实现以下一个快速和有效的方式：

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

我想排在Z偏移添加：在

import numpy as np 

period = [1, 2, 3] 

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

z = np.outer(x,y) 

print(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 ]] 

最终，我会是厕所国王在垂直/列轴（轴= 1）上求和。我需要一个最终的阵列类似如下：

[ 0.0 0.2 2.4 24.6 46.0 60.0] 

Answer 1

您也可以先计算指标，并分配一次：

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 2

循环在第一层面的工作：

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

Answer 3

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

是一个衣衫褴褛的名单。我们可以用viectorized阵列魔法来构建它，至少不是用普通的东西。

为了解决这个问题，我们需要或者变平或垫这种结构

[[ 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)可能会避免这种情况，在不扩展它的情况下对其进行编译迭代。处理非常大的数组时，这可能会有所不同。

结果相反，但这很容易修复。我甚至可以在施工期间修复它。