用于高动态范围的快速准确的卷积算法(如 FFT)?

信息处理 fft 卷积 算法 傅立叶 快速卷积
2021-12-29 14:52:12

由于评估单位根周围的所有内容,基于 FFT 的卷积似乎受到浮点分辨率有限的影响,正如您在此 Python 代码1014

from scipy.signal import convolve, fftconvolve
a = [1.0, 1E-15]
b = [1.0, 1E-15]
convolve(a, b)     # [  1.00000000e+00,   2.00000000e-15,   1.00000000e-30]
fftconvolve(a, b)  # [  1.00000000e+00,   2.11022302e-15,   1.10223025e-16]

有没有不受此问题困扰的快速卷积算法?
还是直接(二次时间)卷积是获得准确解决方案的唯一方法?

(这么小的数字是否足够重要而不是砍掉我的观点。)

4个回答

免责声明:我知道这个话题比较老,但是如果有人正在寻找“快速准确的卷积高动态范围”或类似的东西,这是少数几个不错的结果中的第一个。我想分享我对这个主题的见解,以便将来可能对某人有所帮助。如果我在回答中使用了错误的术语,我深表歉意,但是我在这个主题上发现的所有内容都相当模糊,即使在这个线程中也会导致混淆。无论如何,我希望读者能理解。

直接卷积对每个点的机器精度大多是准确的,即对于结果的每个点的双精度,相对误差通常大致或接近 1.e-16。每个点有 16 个正确数字。不过,舍入误差对于非典型的大卷积可能很重要,严格来说,应该小心取消并使用Kahan summation和足够高精度的数据类型,但实际上误差几乎总是最优的。

除了舍入误差之外,FFT 卷积的误差是“全局相对”误差,这意味着每个点的误差取决于机器精度和结果的峰值。例如,如果结果的峰值为2.e9,则每个点的绝对误差为因此,如果结果中的值应该非常小,比如说21091016=2107109,该点的相对误差可能很大。如果您在结果尾部需要较小的相对误差,FFT 卷积基本上是无用的,例如,您的数据有某种指数衰减并且尾部需要准确的值。有趣的是,如果 FFT 卷积不受该误差的限制,与直接卷积相比,它的舍入误差要小得多,因为你显然做的加法/乘法更少。这实际上就是为什么人们经常声称 FFT 卷积更准确的原因,并且在某种意义上他们几乎是正确的,因此他们可以非常坚定。

不幸的是,没有简单的通用解决方案来获得快速准确的卷积,但根据您的问题,可能有一个......我发现了两个:

如果您有可以通过尾部多项式很好地逼近的平滑内核,那么带有切比雪夫插值的黑盒快速多极子方法可能对您来说很有趣。如果您的内核“不错”,这实际上可以完美运行:您将获得线性(!)计算复杂度和机器精度精度。如果这适合您的问题,您应该使用它。然而实现起来并不容易。

对于某些特定的内核(我认为是凸函数,通常来自概率密度),您可以使用“指数移位”在结果尾部的某些部分获得最佳误差。有一篇PHD 论文和一个带有 python 实现的 github系统地使用了它,作者将其称为准确 FFT 卷积然而,在大多数情况下,这并不是很有用,因为它要么回归到直接卷积,要么无论如何你都可以使用 FFT 卷积。尽管代码是自动执行的,这当然很好。

- - - - - - - - - - 编辑: - - - - - - - - - -

我看了一点Karatsuba算法(实际上我做了一个小实现),在我看来,它通常具有类似 FFT 卷积的错误行为,即相对于结果的峰值会出现错误。由于算法的分而治之的性质,结果尾部的一些值实际上有更好的误差,但我没有看到一种简单的系统方法来判断哪些值或无论如何如何使用这种观察。太糟糕了,起初我认为 Karatsuba 可能在直接卷积和 FFT 卷积之间有用。但我没有看到 Karatsuba 应该优于常见的两种卷积算法的常见用例。

再加上我上面提到的指数偏移:在很多情况下,您可以使用它来改善卷积的结果,但它又不是通用的解决方法。实际上,我将它与 FFT 卷积一起使用以获得非常好的结果(在所有输入的一般情况下:在最坏的情况下,与正常 FFT 卷积相同的误差,在每个点上与机器精度的最好相对误差)。但同样,这仅对特定的内核和数据非常有效,但对我来说,内核和数据或在某种程度上呈指数衰减。

一个候选者是Karatsuba 算法,它在时间内运行。它不是基于转换的。Music-DSP Source Code Archive 中还有一些带有解释的代码,看起来像是对类似算法的独立发现。O(Nlog23)O(N1.5849625)

使用您问题中的数字测试Karatsuba 算法的Python 实现(由sudo pip install karatsuba

import numpy as np
from karatsuba import *
k = make_plan(range(2), range(2))
l = [np.float64(1), np.float64(1E-15)]
np.set_printoptions(formatter={'float': lambda x: format(x, '.17E')})
print "Karatsuba:"
print(k(l, l)[0:3])
print "Direct:"
print(np.convolve(l, l)[0:3])

打印:

Karatsuba:
[1.0, 1.9984014443252818e-15, 1.0000000000000001e-30]
Direct:
[1.00000000000000000E+00 2.00000000000000016E-15 1.00000000000000008E-30]

与其放弃快速卷积算法,为什么不使用具有更高动态范围的 FFT?

这个问题的答案显示了如何使用具有提升多精度的 Eigen FFT 库。

我相信 Cordic 算法的精度可以随心所欲地扩展,如果是这样,请使用整数 DFT 和适合您问题的字长。

直接卷积也是如此,使用非常长的整数。