在numpy数组上快速迭代以求残差平方 - python

我喜欢最小二乘匹配具有许多已知信号形状的数据(浮点数的数组)。我的代码可以运行,但是对于我计划执行的许多运行来说太慢了:

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可以帮助吗?我有很多事情要做。如果是这样,我在哪里可以找到一个代码片段来做到这一点?
  • 参考方案

    基于 Compute mean squared, absolute deviation and custom similarity measure - Python/NumPy 中提出的各种方法,我们希望在这里解决问题。

    方法#1

    我们可以利用基于 np.lib.stride_tricks.as_stridedscikit-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)
    

    More info on use of as_strided based 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 加快了第三个速度!

    Python:如何停止多线程的numpy? - python

    我知道这似乎是一个荒谬的问题,但是我必须在与部门中其他人共享的计算服务器上定期运行作业,当我开始10个作业时,我真的希望它只占用10个核心而不是更多;我不在乎每次运行一个内核所需的时间是否更长:我只是不想让它侵犯其他人的领土,这将需要我放弃工作等等。我只想拥有10个核心,仅此而已。更具体地说,我在基于Python 2.7.3和numpy 1.6.1的Redh…

    numpy.savetxt“元组索引超出范围”? - python

    我试图在文本文件中写几行,这是我使用的代码:import numpy as np # Generate some test data data = np.arange(0.0,1000.0,50.0) with file('test.txt', 'w') as outfile: outfile.write('…

    Python / Scipy过滤器离散化 - python

    我目前正在尝试从Matlab转向Python,并在多个方面取得了成功。但是,我经常使用的Matlab信号处理工具箱中的一个函数是impinvar函数,用于从其模拟版本计算数字滤波器。在Scipy.signal中,我仅发现bilinear函数可以执行类似的操作。但是,与Matlab bilinear function相比,它不需要可选参数来对频率进行一些预变形…

    Python GPU资源利用 - python

    我有一个Python脚本在某些深度学习模型上运行推理。有什么办法可以找出GPU资源的利用率水平?例如,使用着色器,float16乘法器等。我似乎在网上找不到太多有关这些GPU资源的文档。谢谢! 参考方案 您可以尝试在像Renderdoc这样的GPU分析器中运行pyxthon应用程序。它将分析您的跑步情况。您将能够获得有关已使用资源,已用缓冲区,不同渲染状态上…

    Python:图像处理可产生皱纹纸效果 - python

    也许很难描述我的问题。我正在寻找Python中的算法,以在带有某些文本的白色图像上创建皱纹纸效果。我的第一个尝试是在带有文字的图像上添加一些真实的皱纹纸图像(具有透明度)。看起来不错,但副作用是文本没有真正起皱。所以我正在寻找更好的解决方案,有什么想法吗?谢谢 参考方案 除了使用透明性之外,假设您有两张相同尺寸的图像,一张在皱纹纸上明亮,一张在白色背景上有深…