对numpy数组进行快速迭代以获得平方残差

Question

我喜欢使用许多已知信号形状的最小二乘匹配数据(一个numpy数组的浮点数).我的代码有效,但对于我打算做的很多次运行来说太慢了:

import numpy
import time

samples = 50000
width_signal = 100
data = numpy.random.normal(0, 1, samples)
signal = numpy.random.normal(0, 1, width_signal)  # Placeholder

t0 = time.clock()
for i in range(samples - width_signal):
    data_chunk = data[i:i + width_signal]
    residuals = data_chunk - signal
    squared_residuals = residuals**2
    summed_residuals = numpy.sum(squared_residuals)
t1 = time.clock()
print('Time elapsed (sec)', t1-t0)

编辑:纠正了一个错误:第一个方形残差,然后将它们相加.

在我的机器上运行大约需要0.2秒.由于我有许多数据集和信号形状,这太慢了.我的具体问题不允许使用典型的MCMC方法,因为信号形状太不同了.它必须是蛮力.

典型的数据量为50,000浮点数据和100个信号浮点数.这些可以变化几个因素.

我的测试表明:

• 数据的总和占numpy.sum(residuals)90%的时间.我尝试过Python sum(residuals),对于小型阵列(〜<50个元素)来说速度更快,对于更大的阵列来说速度更慢.我应该插入一个if条件吗？
• 我尝试过numpy.roll()而不是直接获取数据,而且.roll()速度较慢.

问题:

• 加速是否有合理的改进？
• 是否有更快的方法来汇总数组？我不知道C,但如果它快得多,我可以试试.
• GPU能帮忙吗？我有很多事要做.如果是这样,我在哪里可以找到执行此操作的代码段？

Answer 1

基于提出的各种方法Compute mean squared, absolute deviation and custom similarity measure - Python/NumPy,我们期待在这里解决我们的案例.

方法#1

我们可以np.lib.stride_tricks.as_strided基于杠杆scikit-image's view_as_windows来获得滑动窗口,从而在这里有我们的第一个解决方案,就像这样 -

from skimage.util import view_as_windows

d = view_as_windows(data,(width_signal))-signal # diffs
out = np.einsum('ij,ij->i',d,d)

上使用的更多信息的as_strided基础view_as_windows.

方法#2

再次根据答案文章中的矩阵乘法技巧,我们可以更好地表现,就像这样 -

def MSD_strided(data, signal):
    w = view_as_windows(data,(width_signal))
    return (w**2).sum(1) + (signal**2).sum(0) - 2*w.dot(signal)

方法#3

我们将通过引入均匀过滤和卷积来改进方法#2 -

from scipy.ndimage.filters import uniform_filter 

def MSD_uniffilt_conv(data, signal):
    hW = width_signal//2
    l = len(data)-len(signal)+1
    parte1 = uniform_filter(data**2,width_signal)[hW:hW+l]*width_signal
    parte3 = np.convolve(data, signal[::-1],'valid')    
    return parte1 + (signal**2).sum(0) - 2*parte3

标杆

发布样本的时间 -

In [117]: %%timeit
     ...: for i in range(samples - width_signal + 1):
     ...:     data_chunk = data[i:i + width_signal]
     ...:     residuals = data_chunk - signal
     ...:     squared_residuals = residuals**2
     ...:     summed_residuals = numpy.sum(squared_residuals)
1 loop, best of 3: 239 ms per loop

In [118]: %%timeit
     ...: d = view_as_windows(data,(width_signal))-signal
     ...: np.einsum('ij,ij->i',d,d)
100 loops, best of 3: 11.1 ms per loop

In [209]: %timeit MSD_strided(data, signal)
10 loops, best of 3: 18.4 ms per loop

In [210]: %timeit MSD_uniffilt_conv(data, signal)
1000 loops, best of 3: 1.71 ms per loop

~140x 加快第三个!