具有相等数量的数据点和系数的多项式拟合

问题描述

我目前正在尝试使用jupyter进行多项式拟合。给定xs中具有对应ys的数据点，下面的函数将返回次数为m的最小二乘多项式。

doc

假设我有以下六个数据点，并且拟合了5级的多项式：

from numpy import *
from matplotlib import pyplot as plt

def findas(m,xs,ys):
    A = array([[0]*(m+1)]*(m+1))
    b = array([0]*(m+1))
    for k in range(m+1):
        b[k] = sum(ys*xs**k)
        for i in range(m+1):
            A[k,i] = sum(xs**(k+i))
    coefs = linalg.solve(A,b)
    print(coefs)
    def fit(x):
        return sum(coefs*(x**array(range(len(coefs)))))
    return fit

根据我的理解，结果曲线应该精确地通过每个数据点（实际上，拉格朗日多项式应该是结果）。

xs = array([1,2,3,4,5,6])
ys = array([-5.21659,2.53152,2.05687,14.1135,20.9673,33.5652])
ft = findas(5,ys)

样本输出：

但是，事实并非如此。曲线距离很远！这里发生了什么？这与舍入误差有关吗？预先感谢！

解决方法

似乎存在截断错误！代码块

A = array([[0]*(m+1)]*(m+1))
b = array([0]*(m+1))
for k in range(m+1):
...

应为：

A = array([[0.]*(m+1)]*(m+1))
b = array([0.]*(m+1))
for k in range(m+1):
...

即我们必须将零指定为浮点数。

此外，舍入误差会在矩阵求逆过程中放大。当我们要求逆的矩阵的特征值在数量级上明显不同时，情况尤其如此。

,

您的代码显示正确；您重新发现了尝试用有限精度算术求近似奇异矩阵的问题。矩阵A看起来像这样

[[       6       21       91      441     2275    12201]
 [      21       91      441     2275    12201    67171]
 [      91      441     2275    12201    67171   376761]
 [     441     2275    12201    67171   376761  2142595]
 [    2275    12201    67171   376761  2142595 12313161]
 [   12201    67171   376761  2142595 12313161 71340451]]

请注意最大值和最小值之间的差异有多大。其特征值由

给出

[7.35326515e+07 1.98781177e+04 5.75757797e+01 1.74921341e+00
 5.89932892e-02 1.37532926e-04]

请注意，最大与最小的比率约为10 ^ 11。因此，矩阵在理论上不是奇异的，而在数值计算中则是奇异的。它的求反会导致非常大的舍入误差，就像除以很小的数字会导致最终结果的精度大大下降一样。

更多详细信息here和相关链接