是否有适用于大数的快速卷积？

Question

我正在对许多相当长的数组进行卷积。尽管它们在 float64 的范围内，但它们中的值也很大。不幸的是，scipy 的 convolve 在使用 fft 方法时给出了错误的答案。这是 MWE：

from scipy.signal import convolve
import numpy as np

A = np.array([0] + [1e100]*10000)
convolve(A, A)

这给出了错误的输出：

array([ 9.15304445e+187, -7.04080342e+187,  1.00000000e+200, ...,
        3.00000000e+200,  2.00000000e+200,  1.00000000e+200])

虽然速度很快：

%timeit convolve(A, A)
458 µs ± 5.09 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

您可以通过以下方式获得正确答案：

convolve(A, A, method="direct")
array([0.e+000, 0.e+000, 1.e+200, ..., 3.e+200, 2.e+200, 1.e+200])

但是这要慢得多：

%timeit convolve(A, A, method="direct")
23.4 ms ± 511 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

有什么办法可以快速得到正确答案吗？我很高兴使用任何免费提供的库。

---

我在不太极端的输入大小下遇到了同样的问题：

A = np.array([0] + [1e10]*10000)
convolve(A, A)
array([1.96826156e+08, 1.35742176e+08, 1.00000000e+20, ...,
       3.00000000e+20, 2.00000000e+20, 1.00000000e+20])

赏金

赏金是针对不太极端的例子

 A = np.array([0] + [1e10]*10000)

Answer 1

没有其他评论进来，我不妨将我的评论变成一个答案，尽管不是一个非常令人满意的答案：

解释 FFT 方法结果

这种情况只是浮点计算精度有限的情况，加上其数值范围很大，导致在面对大绝对值时会产生很大的绝对误差。对于更详细的细节和解释，我们当然可以转向 Goldberg 的“每个计算机科学家应该了解浮点运算”。

简化一点：双精度及其 53 位尾数为我们提供了大约 16 位十进制数字 (

log10(2**53)

)。这意味着即使在最好的情况下，对 1e0 左右的值进行计算也会产生 1e-16 左右的数值误差，并且这些误差当然可能会累积。我们可以通过缩小卷积来看到这一点

A = np.array([0] + [1e0]*10000)
convolve(A, A)
array([-1.10378972e-12, -2.39154439e-12,  1.00000000e+00, ...,
        3.00000000e+00,  2.00000000e+00,  1.00000000e+00])

1e-12 的大致错误是我们可以预料到的，对浮点算术有基本了解的人不会对这种数字错误的计算眨眼。

现在我们将结果缩放至 1e200。会发生什么？我们仍然有相同的有限尾数。因此，误差也随之变化。这正是我们所看到的。

A = np.array([0] + [1e100]*10000)
convolve(A, A)
array([1.97142496e+188, 3.52040171e+187, 1.00000000e+200, ...,
      3.00000000e+200, 2.00000000e+200, 1.00000000e+200])

注意结果和误差之间仍然存在大约 1e12 的因子。相对误差（或信噪比，如果您愿意的话）是相同的。

解释直接法结果

为什么直接卷积不受影响？仅仅是因为这是一个有点人工的数据集。当您开始认真处理这些数据集或使用更多可变数据时，直接卷积的表面优势就消失了。我们可以通过以单精度卷积随机值并将结果与双精度卷积进行比较来证明这一点。

A = np.random.random(10000).astype('f4')
fourier = convolve(A, A)
direct = convolve(A, A, method='direct')
A = A.astype('f8')
reference = convolve(A, A, method='direct')
abs(direct - reference).sum()
23.32324818390154
abs(fourier - reference).sum()
3.2571286909226727

单精度的fft方法比单精度的直接方法更接近双精度的直接方法。这表明 fft 方法通常会更准确，这是有道理的，因为它涉及较少的算术运算，因此添加舍入误差的机会也较少。

缓解问题

有什么我们可以做的吗？并不真地。建议使用长双精度。然而，80 位 x86 格式只提供 64 位尾数，即 20 位十进制数字。所以我们的错误仍然在 1e184 左右。 IEEE-754 中最大的格式是具有 237 位尾数的 binary256。这给出了 72 位数字，将我们的错误变成 1e132。绝对数量仍然很大。

需要 664 位尾数才能准确表示 1e200 到 1e0 范围内的值。那时我们已经深入到任意精度数字的领域，问题是我们是否能够找到快速且准确的解决方案。所以，不。

您应该做的是查看相对值。 1e188 很大，但与 1e200 相比却很小。按绝对最大值缩放结果，不要回头。

旁注

说到缩放结果，您的计算非常接近指数的范围限制。由于卷积和许多其他操作可以简单地按常数缩放，因此我建议您将数字降低到更合理的值范围。问题是您面临中间结果可能溢出指数的风险。

为了实现精确缩放，永远不会改变尾数中的位（非规范化数字除外）并且可以精确反转，我建议使用 2 的幂。

def scale_down(arr):
    arr = np.asarray(arr)
    scale = abs(arr).max()
    _, power = math.frexp(scale)
    scale = math.ldexp(1., power)
    arr = arr * (1. / scale)
    return arr, scale

A, scale = scale_down(A)
convolve(A, A) * scale**2

# or for two:
A, ascale = scale_down(A)
B, bscale = scale_down(B)
convolve(A, B) * (ascale * bscale)

一个实际示例，其中对于当前值来说这是必要的，即采用卷积结果的向量范数。

A = np.array([0] + [1e200]*10000)
np.linalg.norm(A)
inf
A, scale = scale_down(A)
np.linalg.norm(A) * scale
1e202