r - 多项式卷积返回意外值
问题描述
我正在使用 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 vignetteconvolveCpp
中的示例时,我得到的结果与上面的 Python 相同。
是否存在舍入或精度问题convolve
或基础fft
?
解决方案
Rconvolve
函数确实在其计算中使用了 FFT。如果我使用 FFT 在 MATLAB 中复制你的实验,我也会得到一个不精确的结果:
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
。此值比convolve
R 中的结果更接近 1。如果 R 使用较少的 FFT 实现,您看到的差异很可能完全是由于 FFT 中的数值不精确造成的。
请注意,如果我将输入转换为单精度浮点数,r(1)
则变为4.6261e+04
,这意味着这个特定问题确实需要高精度以避免灾难性错误。
Pythonnp.convolve
与 MATLAB 非常相似conv
,在其计算中不使用 FFT,因此能够产生精确的结果。