R 中的拉格朗日插值与 poly.calc 和自定义实现之间的差异

Question

我正在尝试将拉格朗日插值多项式拟合到一组点：

x0 <- c(576020169, 576020229, 576020296, 576020366)
y0 <- c(4816.391, 4908.896, 5011.254, 5116.875)

我正在尝试两种不同的方法。第一个是通过polynom的poly.calc:

library(polynom)
lagrangeInterpolation1 <- as.function(poly.calc(x0, y0))

第二个是通过我在这个答案上找到的一些代码，它使用了函数工厂：

lagrange <- function(x0, y0) {
    f <- function (x) {
        sum(y0 * sapply(seq_along(x0), function(j) {
            prod(x - x0[-j])/prod(x0[j] - x0[-j])
        }))
    }
    return(Vectorize(f, "x"))
}
lagrangeInterpolation2 <- lagrange(x0, y0)

应用两者会产生截然不同的结果，通过自定义代码获得的结果是正确的。例如，对于 x=576020225 的值（位于插值区间内）：

lagrangeInterpolation1(576020225) # 2.486671e+18
lagrangeInterpolation2(576020225) # 4902.472 , reasonable and correct

知道这里可能发生什么吗？

Answer 1

我大胆猜测失败的多项式拟合者头脑有点简单，并且由于 X 值中的 DC 偏移较大而建立了一个朴素的多项式拟合矩阵，其条件数很糟糕。任何旨在准确的合理拟合算法都会预先将 X 数据范围重新调整到 -1,1（或至少减去 Xn 的平均值）。

从所有 X 值中减去 576020300，然后再次尝试拟合可能会解决问题。对于 Excel 来说确实如此（尽管我承认我原以为 Excel 图表实际上是正确的！）。

计算具有大加法常数 A 的 (A+x)^N 可能是一场灾难。

顺便说一句，您的一组点看起来几乎完全是线性的，这不是拉格朗日插值的特别好的功能测试。

FWIW 我使用 Excel 得到线性方程

y = 1.51927651798024x - 8.7512909703529400E+08
R² = 9.9997709635592700E-01

但是如果我尝试在那里拟合二次方程，我会得到一个完全无意义的“答案”，它具有看似合理但毫无意义的 R² 值 - 我有点惊讶，因为我认为它比这更好。我其实已经预料到它会起作用！

Excel 使用具有大偏移量的原始 X 值进行二次拟合失败：

y = -1.2529309982056700E-04x2 + 1.4434425517145300E+05x - 4.1573047570104200E+13
R² = 9.9999994440119700E-01

Excel 成功进行二次拟合并修改了 X' = X-567020300 :

y = -1.2529309982056700E-04x2 + 1.5172783075004900E+00x + 5.0172911254711800E+03
R² = 9.9999994440119700E-01

在我看来，多项式库需要一些关注！