多项式的卷积返回了意料之外的值。

Question

多项式的卷积返回了意料之外的值。

5

我正在使用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文档的convolveCpp示例时，我得到了与上面Python代码中相同的结果。

convolve或其底层的fft是否存在舍入或精度问题？

- user2474226

这是“fft”近似的问题。在Python中，from scipy import signal signal.fftconvolve(g, u)也会给出一个近似的答案，它也使用了“fft”。 - Bas

OP的例子是有效的，但系数的大小掩盖了问题。只需使用convolve(c(1, 2), c(3, 4), type="open")就足以显示convolve的结果与多项式卷积不同-- convolve返回c(4, 11, 6)，而多项式卷积（例如由pracma::conv实现）返回c(3, 10, 8)。 - Robert Dodier

@RobertDodier 序列 [1,2] 和 [3,4] 的卷积应该返回 [3,10,8]。我不知道为什么你的 convolve 返回 [4,11,6]，如果这是 R 语言，它支持 OP 的问题。也许 R 的 convolve 函数有问题。 - Cris Luengo

1

我在混淆问题，我认为OP的问题确实与数值精度有关。仔细重新阅读问题陈述后，我发现它说“convolve(g, rev(u), type="o")”（即反转“u”），结果产生了多项式卷积，而“convolve(g, u, type="o")”则是其他东西（可能更像相关而不是卷积）。这似乎是一个选择不当的默认设置，但OP似乎已经意识到了，所以我的评论没有帮助。 - Robert Dodier

1

我现在明白了，重新阅读? convolve，它说：“请注意，两个序列‘x’和‘y’的卷积的通常定义是由‘convolve(x, rev(y), type = "o")’给出的。”再一次，我没有足够的启示去阅读文档...我仍然认为这是一个非常糟糕的默认行为。 - Robert Dodier

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Cris Luengo · Accepted Answer

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，因此能够产生精确的结果。