将数据拟合到3次多项式

Question

我正在编写一个C++程序,其中我有独立和依赖数据的向量,我想要适合三次函数.但是,我在生成可以适合我的数据的多项式时遇到了麻烦.

部分问题是我不能使用各种数字包,例如GSL(长篇故事); 对我的情况来说甚至可能有点矫枉过正.对于最小二乘拟合,我不需要非常通用的解决方案.我特别希望将我的数据拟合为三次函数.我可以访问索尼的矢量库,它支持4x4矩阵,可以计算它们的反转等.

在Scilab中进行原型设计时,我使用了以下函数:

function p = polyfit(x, y, n)
    m = length(x);
    aa = zeros(m, n+1)
    aa(:,1) = ones(m,1)
    for k = 2:n+1
        aa(:,k) = x.^(k-1)
    end
    p = aa\y
endfunction

不幸的是,这不能很好地映射到我目前的环境.以上示例需要支持M x N + 1维的矩阵.就我而言,那是M x 4,其中M取决于我有多少样本数据.还有左分裂的问题.我需要一个支持任意维矩阵的逆矩阵库.

是否存在最小二乘算法,我可以避免计算aa\y,或者至少将其限制为4x4矩阵？我想我正在尝试将上述算法重写为一个更简单的情况,适用于拟合三次多项式.我不是在寻找代码解决方案,但如果有人能指出我正确的方向,我会很感激.

Answer 1

这是我正在处理的页面,尽管该页面本身并未直接解决您的问题.我的答案摘要是:

如果你不能直接使用Nx4矩阵,那么"手动"进行那些矩阵计算,直到问题变成只有4x4或更小矩阵的东西.在这个答案中,我将概述如何"手动"进行所需的特定矩阵计算.

-

假设您有一堆数据点,(x1,y1)...(xn,yn)并且您正在寻找y = ax^3 + bx^2 + cx + d最适合这些点的三次方程式.

然后按照上面的链接,你写下这个等式:

我会写A,x并B为那些矩阵.然后按照上面的链接,你想要乘以转置A,它将为你提供可以反转的4x4矩阵A^T*A.在方程式中,以下是计划:

A*x = B .................... [我们开始的时候]

(A ^T*A)*x = A ^T*B ..... [乘以A ^T ]

x =(A ^T*A)^-1*A ^T*B ... [乘以A ^T*A 的倒数]

你说你对反转4x4矩阵感到满意,所以如果我们能编写一种方法来获得这些矩阵而不实际使用矩阵对象,我们应该没问题.

所以,这是一种方法,虽然它可能有点像为你的口味制作自己的矩阵库.:)

• 为4x4矩阵的16个条目中的每个条目写一个显式方程.该(i,j)条目(我以(0,0)开始)由x ₁ ⁱ*x ₁ ^j + x ₂ ⁱ*x ₂ ^j + ... + x _N ⁱ*x _N ^j给出.

• 使用矩阵库反转4x4矩阵.那是(A^T*A) - ¹.

• 现在我们所需要的只是A^T*B,这是一个4x1矩阵.它的第i个条目由x ₁ ⁱ*y ₁ + x ₂ ⁱ*y ₂ + ... + x _N ⁱ*y _N给出.

• 增加我们的手工创建的4x4矩阵(A^牛逼*A)^-1由我们手工创建4X1矩阵A^牛逼*B获得最小二乘系数的4X1矩阵的立方.

祝好运!

Answer 2

是的,我们可以将问题限制为使用"4x4矩阵"进行计算.即使对于M个数据点,立方体的最小二乘拟合仅需要四个未知数中的四个线性方程的解.假设所有x坐标都是不同的,则系数矩阵是可逆的,因此原则上系统可以通过反转系数矩阵来求解.我们假设M大于4,通常是最小二乘拟合的情况.

这是Maple的一篇文章,它将一个立方体拟合到数据中,几乎完全隐藏了所解决的细节.一阶最小标准(关于系数的一阶导数作为平方和误差的参数)得到四个线性方程,通常称为正规方程.

您可以在代码中"组合"这四个方程,然后应用矩阵逆或更复杂的解决方案策略.显然,您需要以某种形式存储数据点.一种可能性是两个线性阵列,一个用于x坐标,一个用于y坐标,长度M都是数据点的数量.

注意:我将根据基于1的数组下标来讨论这个矩阵汇编.多项式系数实际上是一个应用程序,其中基于0的数组下标使事情变得更清晰和简单,但是用C或任何其他有利于基于0的下标的语言重写它仍然是读者的练习.

通过参考Mx4阵列A,正则方程的线性系统最容易以矩阵形式表示,其中Mx4阵列的条目是x坐标数据的幂:

A(i,j)=第i个数据点的x坐标上升到幂j-1

设A'表示A的转置,因此A'A是4×4矩阵.

如果我们让d是一个长度为M的列,包含数据点的y坐标(按给定的顺序),那么正规方程组就是这样的:

A'A u = A'd

其中u = [p0,p1,p2,p3]'是具有最小二乘拟合的三次多项式的系数列:

P(x)= p0 + p1*x + p2*x ^ 2 + p3*x ^ 3

你的异议似乎集中在存储和/或操纵Mx4阵列A或其转置的困难.因此,我的答案将集中在如何组装矩阵A'A和A列而不一次明确地存储所有A.换句话说,我们将隐式地执行指示的矩阵 - 矩阵和矩阵 - 向量乘法,以获得可以解决的4x4系统:

C u = f

如果你认为条目C(i,j)是A的第i行与A的第j列的乘积,再加上A'的第i行实际上只是A的第i列的转置这一事实,应该很清楚:

所有数据点上的C(i,j)= SUM x ^(i + j-2)

这肯定是通过使用基于0的下标简化博览会的地方!

累积矩阵C的条目可能是有意义的,矩阵C仅取决于i + j的值,即所谓的Hankel矩阵,在长度为7的线性阵列中,使得:

所有数据点上的W(k)= SUM x ^ k

其中k = 0,..,6.4x4矩阵C具有"条纹"结构,这意味着只显示这七个值.循环遍历数据点的x坐标列表,可以在W的相应条目中累积每个数据点的每个幂的适当贡献.

可以使用类似的策略来组合列f = A'd,即循环数据点并累积以下四个总和:

所有数据点上的f(k)= SUM(x ^ k)*y

其中k = 0,1,2,3.[当然在上面的总和中,值x,y是公共数据点的坐标.]

注意事项: 这满足了仅使用4x4矩阵的目标.然而,人们通常试图避免显式形成正规方程的系数矩阵,因为这些矩阵通常是数值分析中被称为病态的.特别是当x坐标紧密间隔的情况下,当试图通过反转系数矩阵来解决系统时会导致困难.

求解这些正规方程的一种更复杂的方法是正规方程的共轭梯度法,可以用代码计算矩阵向量乘积A u和A'v一次一个条目(使用我们上面所说的)关于A)的条目.

共轭梯度法的准确性通常是令人满意的,因为它具有自然的迭代方法,尤其是 当人们可以用一点额外的精度计算出所需的点积.