具有相同数量的数据点和系数的多项式拟合答案

【问题标题】：Polynomial fitting with equal number of data points and coefficients具有相同数量的数据点和系数的多项式拟合
【发布时间】：2021-02-19 05:42:51
【问题描述】：

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

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

假设我有以下六个数据点并拟合一个 5 次多项式：

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,xs,ys)

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

xdense = arange(1,6.1,0.1)
ydense = [ft(x) for x in xdense]   

plt.plot(xdense,ydense)
plt.plot(xs,ys,'rx')
plt.show()

样本输出：

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

【问题讨论】：

@Mr.T 现在好点了吗？抱歉，不要用这么多。
@Mr.T 真的很抱歉！我只是注意到粘贴代码时所有标识都消失了！是的，我正在使用 numpy，并且我已经进行了上面的编辑。
这很可能是 chi^2 中的局部最小值...拟合是局部优化而不是全局优化
@mikuszefski 这是有道理的。假设样本点后面的“真实”函数为 f(x) = 3x+2，假设样本点数为 n，回归多项式的次数为 N。如果 n=N>=1，这样的局部-但某些 N 仍会出现非全局最小值？
等等...我有点糊涂了。这是完全线性的，因此最坏的情况可能比问题条件严重的点更高。

标签： python regression curve-fitting polynomials data-fitting

【解决方案1】：

似乎有截断错误！代码块

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):
...

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

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

【讨论】：

【解决方案2】：

您的代码显示正确；您重新发现了尝试使用有限精度算术反转近奇异矩阵的问题。矩阵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 和相关链接

【讨论】：

非常感谢您的评论！我才知道，当我把 A = array([[0.]*(m+1)]*(m+1)), b = array([0.]*(m+1))... 加上 0. 而不仅仅是 0 时，它似乎工作正常！
哦，这很有趣。我所说的要点仍然正确，但这里的额外转折是计算是在整数而不是精度低得多的浮点数中执行的，所以我描述的这些问题在前面已经体现出来了。它实际上想知道为什么你只用 6 分就遇到这些问题。您可能希望将其写成您自己的答案，以便其他人更容易找到它
嗯。 Runge 对我来说是新的，但是按照链接中的链接，它指出它类似于 Gibb 的现象，这与舍入误差无关，而是一个基本的固有属性。寻找invert nearly singular matrices的可能性我基本上回到了我上面提到的SVD。