加快矩阵计算（在子阵列上循环）[numpy]

我正在尝试为我的学生作业编写一个算法，它运行良好。但是，计算需要很长时间，尤其是对于大数组。这部分代码会减慢所有程序。加快矩阵计算（在子阵列上循环）[numpy]

Shapes: X.shape = mask.shape = logBN.shape = (500,500,1000), 
     F.shape = (20,20), 
     A.shape = (481,481), 
     s2 -- scalar.

我应该如何改变这种代码，使其更快？

h = F.shape[0] 
w = F.shape[1] 
q = np.zeros((A.shape[0], A.shape[1], X.shape[2])) 
for i in range(A.shape[0]): 
    for j in range(A.shape[1]): 
     mask[:,:,:] = 0 
     mask[i:i + h,j:j + w,:] = 1 
     q[i,j,:] = ((logBN*(1 - mask)).sum(axis=(0,1)) + 
        (np.log(norm._pdf((X[i:i + h,j:j + w,:]-F[:,:,np.newaxis])/s2)/s2)).sum(axis=(0,1))

来源

2016-11-01 Acapello

这是不完整的 - 将所有变量（F，A，X）放在一起，以便人们可以使用某些东西。如果在数组上迭代，通常最好转换为python列表，因为它非常慢 - 使用向量操作最快。 – kabanus

@kabanus我不能在程序工作期间生成它们。 – Acapello

我建议打印一次，并在开始时粘贴结果。 – kabanus

只是试图让你的内环感

mask[:,:,:] = 0 
    mask[i:i + h,j:j + w,:] = 1 
    q[i,j,:] = ((logBN*(1 - mask)).sum(axis=(0,1)) + 
       (np.log(norm._pdf((X[i:i + h,j:j + w,:]-F[:,:,np.newaxis])/s2)/s2)).sum(axis=(0,1))

看起来像

idx = (slice(i,i+h), slice(j,j_w), slice(None)) 
mask = np.zeros(X.shape) 
mask(idx) = 1 
mask = 1 - mask 
# alt mask=np.ones(X.shape);mask[idx]=0 
term1 = (logBN*mask).sum(axis=(0,1)) 
term2 = np.log(norm._pdf((X[idx] - F[...,None])/s2)/s2).sum(axis=(0,1)) 
q[i,j,:] = term1 + term2

所以idx和mask在A定义一个子阵列。您正在阵列外部使用logBN;里面还有term。您正在对第1个2个dim进行求和，因此term1和term2的形状为X.shape[2]，您可以将其保存在q中。

该面具/窗口是20x20。

作为第一次切割，我试图一次计算所有i,j的term2。这看起来像一个典型的滑动窗口问题。我也尝试将term1表示为减法 - 整个logBN减去此窗口。

来源

2016-11-01 17:42:48 hpaulj

后通过log，exp，power代数运算重杂耍，这一切来到这个 -

# Params 
m,n = F.shape[:2] 
k1 = 1.0/(s2*np.sqrt(2*np.pi)) 
k2 = -0.5/s2**2 
k3 = np.log(k1)*m*n 

out = np.zeros((A.shape[0], A.shape[1], X.shape[2])) 
for i in range(A.shape[0]): 
    for j in range(A.shape[1]): 
     mask[:] = 1 
     mask[i:i + h,j:j + w,:] = 0 
     XF = (X[i:i + h,j:j + w,:]-F[:,:,np.newaxis])   
     p1 = np.einsum('ijk,ijk->k',logBN,mask) 
     p2 = k2*np.einsum('ijk,ijk->k',XF,XF) 
     out[i,j,:] = p1 + p2 
out += k3

使用几件事情是 -

1] norm._pdf基本上是：norm.pdf(x) = exp(-x**2/2)/sqrt(2*pi)。所以，我们可以内嵌实现并在脚本级别对其进行优化。

2]由标量划分将不会有效，所以这些被它们的倒数乘法代替。所以，在进入循环之前，作为一个预处理存储他们的倒数。

来源

2016-11-01 18:25:42 Divakar

这就是炉排！谢谢！我有想法改变一个算法，但是你用'einsum'的解决方案工作速度快8倍。虽然我不明白为什么它运行速度更快，普通的numpy功能。 – Acapello

@Acapello相信我，我不得不通过代数函数，因此在帖子顶部的评论。是一个艰苦的会议！ :)这更快，因为我们正在做更小的操作。 – Divakar

加快矩阵计算（在子阵列上循环）[numpy]

回答

相关问题