为什么我的矢量化 Numpy 代码比非矢量化代码花费的时间更长

Why is my vectorized Numpy code taking longer than the non-vectorized code

所以我正在使用大量数据计算泊松分布。我有一个形状数组 (2666667,19) - “尖峰”,和一个形状数组 (19,100) - “placefields”。我曾经有一个循环遍历 2666667 维度的 for 循环,大约需要 60 秒才能完成。然后,我了解到如果我对循环进行矢量化,它会变得更快,所以我尝试这样做。矢量化形式工作并输出相同的结果,但是,现在需要 120 秒:/

这是原始循环(60 秒):

def compute_probability(spikes,placefields):
    nTimeBins = len(spikes[0])
    probability = np.empty((nTimeBins, 99)) #empty probability matrix
    for i in range(nTimeBins):
        nspikes = np.tile(spikes[:,i],(99))
        nspikes = np.swapaxes(nspikes,0,1)
        maxL = stats.poisson.pmf(nspikes,placefields)
        maxL = maxL.prod(axis=0)
        probability[i,:] = maxL
    return probability

这里是矢量化形式 (120s)

def compute_probability(spikes,placefields):  

    placefields = np.reshape(placefields,(19,99,1))
    #prepared placefields

    nspikes = np.tile(spikes, (99,1,1))
    nspikes = np.swapaxes(nspikes,0,1)
    #prepared nspikes

    probability = stats.poisson.pmf(nspikes,placefields)
    probability = np.swapaxes(probability.prod(axis=0),0,1)
    return probability

为什么这么慢。我认为可能是矢量化形式创建的平铺数组非常庞大,它们占用了大量内存。我怎样才能让它走得更快? 下载 samplespikes 和 sampleplacefields(如评论所建议)- https://mega.nz/file/lpRF1IKI#YHq1HtkZ9EzYvaUdlrMtBwMg-0KEwmhFMYswxpaozXc

编辑: 问题是虽然它是矢量化的,但这个巨大的数组占用了太多的内存。我已将计算分成块,现在效果更好:

placefields = np.reshape(placefields,(len(placefields),99,1))
nspikes = np.swapaxes(np.tile(spikes, (xybins,1,1)),0,1)
probability = np.empty((len(spikes[0]), xybins))

chunks = len(spikes[0])//20
n = int(len(spikes[0])/chunks)
for i in range(0,len(nspikes[0][0]),n):
    nspikes_chunk = nspikes[:,:,i:i+n]
    probability_chunk = stats.poisson.pmf(nspikes_chunk,placefields)
    probability_chunk = np.swapaxes(probability_chunk.prod(axis=0),0,1)
    if len(probability_chunk)<(len(spikes)//chunks):
        probability[i:] = probability_chunk
    else:
        probability[i:i+len(probability_chunk)] = probability_chunk

这可能是由于 memory/cache 影响

第一个代码适用于适合 CPU 缓存的小数组。这不是很好,因为每个 Numpy 函数调用都需要一些时间。第二个代码解决了这个问题。但是,它 allocate/fill 巨大的数组 在几个 GiB 的内存中。在 CPU 缓存中工作比在主内存 (RAM) 中工作要快得多。当工作数组仅使用一次时尤其如此(因为昂贵的 OS page-faults),这在您的代码中似乎就是这种情况。如果您没有足够的内存,OS 将 read/write 临时数据存储在 SSD/HDD 与 RAM 和 CPU 缓存相比非常慢的存储设备中。

最好的解决方案可能是在块上工作,以便操作既矢量化(减少 Numpy 函数调用的开销)又适合 CPU 缓存(减少 RAM 的成本 reads/writes).请注意,如今主流 PC 处理器上的最后一级缓存的大小通常只有几 MiB。

要点是矢量化并不总能使事情变得更快。为了获得更好的性能,应该关心被操纵数据块的大小,以便它们适合 CPU 缓存。

PS:注意,如果不太在意精度,可以用simple-precision(np.float32)代替double-precision(np.float64)来加速稍微计算一下。