優化「差分函數」的計算

Question

我的代碼調用許多「差分函數」來計算「Yin algorithm」（基頻提取器）。

的不同功能（在本文EQ 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) 

有沒有一種方法來優化這個雙循環補償utation（只有python，最好沒有其他庫，而不是numpy）？

編輯：改變的代碼，以避免指數誤差（j循環，@Elliot回答）

EDIT2：改變的代碼，以用x [0]（j循環，@hynekcer評語）

Answer 1

編輯：改進的速度，以220微秒 - 看到在端編輯 - 直接版本

所需計算可以通過Autocorrelation function或類似地通過卷積能夠容易地評價。維納 - 欽欽定理允許用兩個快速傅立葉變換（FFT）計算自相關，其時間複雜度爲O（n  log  n）。 我使用來自Scipy包的加密卷積函數fftconvolve。一個好處是，這裏很容易解釋爲什麼它的工作原理。一切都是矢量化的，在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] 

• 與differenceFunction_1loop功能Elliot's answer相比：它是FFT更快：430  微秒相比原來的1170微秒。它的啓動速度大約爲tau_max >= 40。數值精度很高。與精確整數結果相比，最高相對誤差小於1E-14。 （因此它可以很容易地四捨五入到精確的長整數解決方案。）
• 參數tau_max對算法不重要。它最終只限制輸出。索引0處的零元素被添加到輸出中，因爲索引應該在Python中以0開頭。
• 參數W在Python中不重要。尺寸更好地被內省。
• 數據最初轉換爲np.float64以防止重複轉換。它快了一半。任何小於np.int64的類型都會因爲溢出而無法接受。
• 所需的差分函數是雙能量減自相關函數。這可以通過卷積來評估：correlate(x, x) = convolve(x, reversed(x)。
• 「自從Scipy v0.19正常convolve自動選擇這種方法或基於估計速度更快的直接方法。」這種啓發式方法不適用於這種情況，因爲卷積比tau_max估計得多得多，並且它必須比直接方法快得多。
• 它也可以通過Nippy的ftp模塊來計算，不需要Scipy，通過將Python的答案Calculate autocorrelation using FFT in matlab重寫（最後在下面）。我認爲上面的解決方案可以更容易理解。

---

證明：（用於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)] 

優化由向量的表達式：通過numpy的FFT實施溶液直接：

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

每一步也可以測試並比較由測試數據數值

---

EDIT。

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) 

最新解決方案比艾略特的difference_by_convol爲tau_max> = 20更快這比不依賴由於間接費用的類似比率多的數據大小。

Answer 2

我不知道如何找到替代你的嵌套循環問題，但對於算術函數，你可以使用numpy庫。它比手動操作更快。

import numpy as np 
tmp = np.subtract(long(x[j] ,x[j + tau]) 

Answer 3

我會做這樣的事情：

>>> x = np.random.randint(0, high=32000, size=2048, dtype='int16') 
>>> tau_max = 106 
>>> res = np.square((x[tau_max:] - x[:-tau_max])) 

但是我相信這是不是這樣做的最快方式。

Answer 4

首先，你應該考慮數組的邊界。你原來編寫的代碼將得到一個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 5

在相對優化的算法可以用numba.jit優化解釋：

import timeit 

import numpy as np 
from numba import jit 


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 = int(x[j] - x[j + tau]) 
      df[tau] += tmp * tmp 
    return df 


@jit 
def differenceFunction2(x, W, tau_max): 
    df = np.ndarray(shape=(tau_max,)) 
    for tau in range(1, tau_max): 
     for j in range(0, W - tau): 
      tmp = int(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) 


print('old', 
     timeit.timeit('differenceFunction(x, W, tau_max)', 'from __main__ import differenceFunction, x, W, tau_max', 
        number=20)/20) 
print('new', 
     timeit.timeit('differenceFunction2(x, W, tau_max)', 'from __main__ import differenceFunction2, x, W, tau_max', 
        number=20)/20) 

結果：

old 0.18265145074453273 
new 0.016223197058214667 

您可以更好的結果結合算法的優化和numba.jit。

Answer 6

下面是使用列表理解另一種方法。它大約不到原始功能所用時間的十分之一，但不會超過Elliot's answer。無論如何，只要把它放在那裏。

import numpy as np 
import time 

# 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 

# with list comprehension 
def differenceFunction_1loop_listcomp(x, W, tau_max): 
    df = [sum(((x[:-tau]) - (x[tau:]).astype(np.long))**2) for tau in range(1, tau_max)] 
    return [0] + df[:] 

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

s = time.clock() 
twoloop = differenceFunction_2loop(x, W, tau_max) 
print(time.clock() - s) 

s = time.clock() 
oneloop = differenceFunction_1loop(x, W, tau_max) 
print(time.clock() - s) 

s = time.clock() 
listcomprehension = differenceFunction_1loop_listcomp(x, W, tau_max) 
print(time.clock() - s) 

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

性能測試結果（大約）：

differenceFunction_2loop() = 0.47s 
differenceFunction_1loop() = 0.003s 
differenceFunction_1loop_listcomp() = 0.033s