加快矩阵计算（循环子数组）[numpy]答案

【问题标题】：Speed up matrix calculations(looping over subarrays) [numpy]加快矩阵计算（循环子数组）[numpy]
【发布时间】：2017-03-14 20:47:55
【问题描述】：

我正在尝试为我的学生作业编写一个算法，它运行良好。但是，计算需要很长时间，尤其是对于大数组。这部分代码正在减慢所有程序。

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

【问题讨论】：

这并不完整 - 将所有变量 (F,A,X) 放在一起，以便人们可以使用某些东西。如果对数组进行迭代，通常最好转换为 python 列表，因为它非常慢 - 使用向量操作更快。
@kabanus 我不能在程序工作期间生成它们。
我建议打印一次，然后将开头的结果粘贴到这里。
@Divakar 我刚刚添加了信息：logBN -- 矩阵和 s2 -- 标量。

标签： python algorithm performance numpy matrix

【解决方案1】：

在通过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] 标量的除法效率不高，因此将其替换为乘以它们的倒数。因此，在进入循环之前，作为预处理存储它们的倒数。

【讨论】：

太棒了！谢谢！我有改变算法的想法，但是您使用 einsum 的解决方案的运行速度快了 8 倍。虽然我不明白为什么它可以更快地运行普通的 numpy 函数。
@Acapello 相信我，我必须研究代数函数，因此帖子顶部的评论。是一个艰苦的会议！ :) 这更快，因为我们执行的操作要少得多。

【解决方案2】：

只是试图理解你的内部循环

    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 在里面。您正在对 1st 2 dim 上的值求和，因此 term1 和 term2 的形状均为 X.shape[2]，您将其保存在 q 中。

那个遮罩/窗口是 20x20。

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

【讨论】：