优化"差函数"的计算

Question

我的代码调用了许多"差异函数"来计算" Yin算法 "(基频提取器).

差异函数(论文中的等式6)定义为:

这是我对差异函数的实现:

def differenceFunction(x, W, tau_max):
   df = [0] * tau_max
   for tau in range(1, tau_max):
      for j in range(0, W - tau):
          tmp = long(x[j] - x[j + tau])
          df[tau] += tmp * tmp
   return df

例如:

x = np.random.randint(0, high=32000, size=2048, dtype='int16')
W = 2048
tau_max = 106
differenceFunction(x, W, tau_max)

有没有办法优化这个双循环计算(仅使用python,最好没有其他库而不是numpy)？

编辑:更改代码以避免索引错误(j循环,@艾略特答案)

EDIT2:更改代码使用x [0](j循环,@ hynekcer评论)

Answer 1

首先,您应该考虑数组的边界.您最初编写的代码将得到一个IndexError.您可以通过矢量化内循环来获得显着的加速

import numpy as np

# original version
def differenceFunction_2loop(x, W, tau_max):
   df = np.zeros(tau_max, np.long)
   for tau in range(1, tau_max):
      for j in range(0, W - tau): # -tau eliminates the IndexError
          tmp = np.long(x[j] -x[j + tau])
          df[tau] += np.square(tmp)
   return df

# vectorized inner loop
def differenceFunction_1loop(x, W, tau_max):
    df = np.zeros(tau_max, np.long)
    for tau in range(1, tau_max):
        tmp = (x[:-tau]) - (x[tau:]).astype(np.long)
        df[tau] = np.dot(tmp, tmp)
    return df

x = np.random.randint(0, high=32000, size=2048, dtype='int16')
W = 2048
tau_max = 106
twoloop = differenceFunction_2loop(x, W, tau_max)
oneloop = differenceFunction_1loop(x, W, tau_max)

# confirm that the result comes out the same. 
print(np.all(twoloop == oneloop))
# True

现在进行一些基准测试.在ipython我得到以下

In [103]: %timeit twoloop = differenceFunction_2loop(x, W, tau_max)
1 loop, best of 3: 2.35 s per loop

In [104]: %timeit oneloop = differenceFunction_1loop(x, W, tau_max)
100 loops, best of 3: 8.23 ms per loop

所以,加速大约300倍.

Answer 2

编辑:将速度提高到220μs - 请参见最后的编辑 - 直接版本

可以通过自相关函数或类似地通过卷积来容易地评估所需的计算.Wiener-Khinchin定理允许用两个快速傅里叶变换(FFT)计算自相关,时间复杂度为O(n log n).我使用卷积加速而功能fftconvolve从SciPy的包.一个优点是,它很容易解释为什么它的工作原理.一切都是矢量化的,没有Python解释器级别的循环.

from scipy.signal import fftconvolve

def difference_by_convol(x, W, tau_max):
    x = np.array(x, np.float64)
    w = x.size
    x_cumsum = np.concatenate((np.array([0.]), (x * x).cumsum()))
    conv = fftconvolve(x, x[::-1])
    df = x_cumsum[w:0:-1] + x_cumsum[w] - x_cumsum[:w] - 2 * conv[w - 1:]
    return df[:tau_max + 1]

• 与Elliot答案中的differenceFunction_1loop功能相比:使用FFT更快:430μs与原始1170μs相比.它开始变得更快.数值精度很高.与精确整数结果相比,最高相对误差小于1E-14.(因此可以很容易地舍入到完全长整数解.)tau_max >= 40
• 该参数tau_max对于算法并不重要.它最终只限制输出.索引0处的零元素被添加到输出中,因为索引应该在Python中以0开头.
• 该参数W在Python中并不重要.尺寸最好是内省的.
• 最初将数据转换为np.float64以防止重复转换.它快了一半.任何小于np.int64的类型都会因溢出而无法接受.
• 所需的差函数是双能减去自相关函数.这可以通过卷积来评估:correlate(x, x) = convolve(x, reversed(x).
• "从Scipy v0.19开始,正常情况会convolve自动选择此方法或基于估算速度更快的直接方法." 因为卷积计算更也就是说启发式不足以这种情况下tau比tau_max,它必须比直接法更快FFT来抵消.
• 它也可以通过Numpy ftp模块在没有Scipy的情况下通过重写答案来计算在Matlab中使用FFT计算自相关到Python(下面的结尾).我认为上面的解决方案可以更容易理解.

证明:(对于Pythonistas :-)

最初的天真实现可以写成:

df = [sum((x[j] - x[j + t]) ** 2   for j in range(w - t))  for t in range(tau_max + 1)]

哪里tau_max < w.

按规则推导 (a - b)**2 == a**2 + b**2 - 2 * a * b

df = [  sum(x[j] ** 2 for j in range(w - t))
      + sum(x[j] ** 2 for j in range(t, w))
      - 2 * sum(x[j] * x[j + t] for j in range(w - t))
      for t in range(tau_max + 1)]

在帮助下替换前两个元素x_cumsum = [sum(x[j] ** 2 for j in range(i)) for i in range(w + 1)]可以在线性时间内轻松计算.sum(x[j] * x[j + t] for j in range(w - t))用conv = convolvefft(x, reversed(x), mode='full')具有大小输出的卷积代替len(x) + len(x) - 1.

df = [x_cumsum[w - t] + x_cumsum[w] - x_cumsum[t]
      - 2 * convolve(x, x[::-1])[w - 1 + t]
      for t in range(tau_max + 1)]

通过矢量表达式进行优化:

df = x_cumsum[w:0:-1] + x_cumsum[w] - x_cumsum[:w] - 2 * conv[w - 1:]

每个步骤也可以通过测试数据进行数字测试和比较

编辑: Numpy FFT直接实施解决方案.

def difference_fft(x, W, tau_max):
    x = np.array(x, np.float64)
    w = x.size
    tau_max = min(tau_max, w)
    x_cumsum = np.concatenate((np.array([0.]), (x * x).cumsum()))
    size = w + tau_max 
    p2 = (size // 32).bit_length()
    nice_numbers = (16, 18, 20, 24, 25, 27, 30, 32)
    size_pad = min(x * 2 ** p2 for x in nice_numbers if x * 2 ** p2 >= size)
    fc = np.fft.rfft(x, size_pad)
    conv = np.fft.irfft(fc * fc.conjugate())[:tau_max]
    return x_cumsum[w:w - tau_max:-1] + x_cumsum[w] - x_cumsum[:tau_max] - 2 * conv

它比我之前的解决方案快两倍以上,因为卷积的长度被限制在最接近的"漂亮"数字之后W + tau_max,其中小的素数因子被评估为满2 * W.也没有必要像使用`fftconvolve(x,reversed(x))那样将相同的数据转换两次.

In [211]: %timeit differenceFunction_1loop(x, W, tau_max)
1.1 ms ± 4.51 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [212]: %timeit difference_by_convol(x, W, tau_max)
431 µs ± 5.69 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [213]: %timeit difference_fft(x, W, tau_max)
218 µs ± 685 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

对于tau_max> = 20,最新的解决方案比Eliot的difference_by_convol更快.由于类似的开销成本比率,该比率并不太依赖于数据大小.