多项式的卷积 returns 意外值

Question

我正在使用 R 3.6.2（平台 = x86_64-w64-mingw32）

在以下两个多项式系数向量的卷积中，我希望第一个条目正好是 1.0，但 convolve 函数不同：

g <- c(1, -49, 1155, -17441, 189700, -1583071, 10545901, -57608692, 
263063351, -1018546561, 3380085631, -9693547553, 24176423345, 
-52691112850)
u <- c(1, -6, 11, -6)
convolve(g, rev(u), type = 'o')
# output
 [1]  1.000172e+00 -5.500020e+01  1.460000e+03 -2.491600e+04
 [5]  3.073450e+05 -2.920052e+06  2.223567e+07 -1.394361e+08
 [9]  7.342188e+08 -3.293898e+09  1.273071e+10 -4.275645e+10
[13]  1.256299e+11 -3.246592e+11  6.402486e+11 -7.246608e+11
[17]  3.161467e+11

请注意，结果中的第一个条目是 1.000172，而不是 1.0。

在 Python 3.7.4 中执行相同的卷积提供了预期的答案：

import numpy as np
g = [1, -49, 1155, -17441, 189700, -1583071, 10545901, -57608692, 263063351, -1018546561, 3380085631, -9693547553, 24176423345, -52691112850]
u = [1, -6, 11, -6]
np.convolve(g,u)
array([            1,           -55,          1460,        -24916,
              307345,      -2920052,      22235673,    -139436079,
           734218840,   -3293897685,   12730714010,  -42756453616,
        125629929970, -324659189789,  640248619213, -724660781420,
        316146677100], dtype=int64)

另外，当我使用 Rcpp vignette 中的 convolveCpp 示例时，我得到的结果与上面 Python 中的相同。

convolve 或基础 fft 是否存在舍入或精度问题？

Answer 1

R convolve function 确实在其计算中使用了 FFT。如果我在 MATLAB 中使用 FFT 复制您的实验，我也会得到一个不准确的结果：

format long

g = [1, -49, 1155, -17441, 189700, -1583071, 10545901, -57608692, ...
263063351, -1018546561, 3380085631, -9693547553, 24176423345, -52691112850];
u = [1, -6, 11, -6];
r = ifft(fft([g,zeros(1,numel(u)-1)]) .* fft(up[u,zeros(1,numel(g)-1)]);

这里的第一个值 r(1) 是 0.999971277573529。该值比 R 中 convolve 的结果更接近 1。如果 R 使用较少的 FFT 实现，您看到的差异很可能完全是由于 FFT 中的数值不精确.

请注意，如果我将输入转换为单精度浮点数，r(1) 将变为 4.6261e+04，这意味着此特定问题确实需要高精度以避免灾难性错误。

Python 的 np.convolve 与 MATLAB 的 conv 非常相似，在其计算中不使用 FFT，因此能够产生精确的结果。

多项式的卷积 returns 意外值

Convolution of polynomials returns unexpected value

floating-point

r

fft

rounding

convolution