编辑:改进的速度,以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更快这比不依赖由于间接费用的类似比率多的数据大小。
你想更快的方法吗?或者使用矩阵寻找矢量化解决方案。 – Dark
我要找使用python – PatriceG
一种方式最快的方法,如果你想保留的for循环,并获得atmost速度使用'numba'。你是否使用其他库? – Dark