logo标签列表

将数据拟合到三次多项式

c++mathstatistics
13

我正在编写一个C++程序，在其中有独立和依赖数据的向量，希望将它们拟合成一个三次函数。然而，我遇到了生成能够拟合我的数据的多项式的问题。

问题的一部分是我不能使用各种数值包，例如GSL（长话短说）；甚至可能对我的情况来说过于复杂。我不需要一个非常通用的最小二乘拟合解决方案。我特别想把我的数据拟合成一个三次函数。我可以访问Sony的向量库，该库支持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的矩阵？我想将上述算法改写为适用于拟合三次多项式的简化情况。我不是在寻找代码解决方案，但如果有人能指点我正确的方向，我将不胜感激。

即使您无法使用 GSL 或者 uBLAS，您仍然可以查看它们的源代码，以了解如何实现。 - Lightness Races in Orbit
1
“查看源代码”并没有真正解决我的问题，完全错过了重点。这些库已经建立在强大的矩阵库之上，并使用我试图避免的相同技术。我正在寻找将我的问题分解为4x4矩阵或更简单的东西。 - lhumongous
2
我需要一个支持任意维度矩阵求逆的矩阵库。这是一个很大的要求，因为非方阵矩阵没有定义逆矩阵。 - John
3个回答
12

这里是我正在使用的页面，虽然该页面本身并没有直接回答你的问题。我的回答概括如下：

如果您无法直接使用 Nx4 矩阵，则手动进行这些矩阵计算，直到将问题简化为仅涉及 4x4 或更小的矩阵。在本回答中，我将概述如何手动执行您需要的特定矩阵计算。

--

假设你有一堆数据点 (x1,y1)...(xn,yn)，并且你正在寻找最适合这些点的三次方程式 y = ax^3 + bx^2 + cx + d。
然后按照上面的链接，你会写下这个方程：

enter image description here

我将用A，x和B表示这些矩阵。然后，按照上面的链接，您需要乘以A的转置，这将给您4x4矩阵AT*A，您可以对其进行反演。在等式中，以下是计划：

A * x = B .................... [我们开始的内容]

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

x = (AT * A)-1 * AT * B ... [乘以AT * A的逆]

您说您很满意反转4x4矩阵，因此，如果我们可以编写一种方法来获取这些矩阵，而不实际使用矩阵对象，那么我们应该没问题。

因此，这里有一种方法，尽管它可能有点像为您的口味制作自己的矩阵库。:)

  • 为4x4矩阵的16个元素编写显式方程。第个元素（从（0,0）开始）由以下公式给出：x1i * x1j + x2i * x2j + ... + xNi * xNj。

  • 使用您的矩阵库反转该4x4矩阵。即 (AT * A)-1。

  • 现在我们只需要 AT * B，这是一个4x1矩阵。它的第个元素由以下公式给出：x1i * y1 + x2i * y2 + ... + xNi * yN。

  • 将手动创建的4x4矩阵(AT * A)-1与手动创建的4x1矩阵AT * B相乘，以获得最小二乘系数的4x1矩阵。

祝好运！

9
是的，我们可以将问题限制为使用“4x4矩阵”进行计算。即使对于M个数据点的三次最小二乘拟合，也只需要解四个未知数的线性方程组。假设所有的x坐标都不同，则系数矩阵是可逆的，因此原则上可以通过求解系数矩阵的逆来解决该系统。我们假设M大于4，这通常是最小二乘拟合的情况。
这里有一个Maple的写作Fitting a cubic to data，它几乎完全隐藏了正在解决的细节。一阶最小条件（关于系数的一阶导数作为平方误差的参数）给出了四个线性方程，通常称为normal equations。
您可以在代码中“组装”这四个方程，然后应用矩阵求逆或更复杂的解决策略。显然，您需要以某种形式存储数据点。一种可能性是两个线性数组，一个用于x坐标，一个用于y坐标，长度均为M，即数据点的数量。

NB: 我将讨论此矩阵组装，以1为基的数组下标。多项式系数实际上是0为基的数组下标更简洁，但是将其重写成使用0为基的下标的C或任何其他语言留给读者作为练习。

正常方程组最容易用矩阵形式表示，通过引用一个Mx4数组A，其条目为x坐标数据的幂：

A(i,j) = 第i个数据点的x坐标的(j-1)次方

让A'表示A的转置，这样A'A就是一个4x4的矩阵。

如果我们让d是包含数据点y坐标（按给定顺序）的长度为M的列，则正常方程组就是这样的：

A'Au = A'd

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

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

你的反对意见似乎集中在存储和/或操作Mx4数组A或其转置方面。因此，我的答案将侧重于如何组装矩阵A'A和列A'd，而无需一次性明确存储所有A。换句话说，我们将隐式地进行指示的矩阵-矩阵和矩阵-向量乘法，以获得一个4x4系统，您可以解决：
Cu = f
如果您考虑条目C(i,j)是A'的第i行与A的第j列的乘积，再加上A'的第i行实际上只是A的第i列的转置，那么应该清楚：
C(i,j) = 所有数据点上的SUM x^(i+j-2)
这当然是使用基于0的下标简化表述的地方之一！
累加矩阵C的条目可能是有意义的，这些条目仅取决于i+j的值，即所谓的Hankel矩阵，并将它们放入长度为7的线性数组中，例如：
W(k) = 所有数据点上的SUM x^k

其中k = 0,..,6。4x4矩阵C具有“条纹”结构，这意味着只出现这七个值。循环遍历数据点的x坐标列表，可以在W的适当条目中累积每个数据点的每个幂的适当贡献。

类似的策略可用于组装列f = A' d，即循环遍历数据点并累积以下四个求和：

f(k) = 求和（x^k)*y所有数据点

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

注意事项：这满足仅使用4x4矩阵的目标。但是，通常尝试避免显式形成正规方程系数矩阵，因为这些矩阵在数值分析中经常被称为病态矩阵。特别是x坐标紧密间隔的情况在尝试通过反转系数矩阵来解决系统时可能会导致困难。

一种更为复杂的解决这些正规方程的方法是使用正规方程的共轭梯度法，可以通过计算逐个条目的矩阵-向量乘积A u和A' v的代码来实现（使用上述关于A的条目的内容）。
由于其自然的迭代方法，特别是当可以使用略微额外精度计算所需的点积时，共轭梯度法的准确性通常是令人满意的。
哈哈...你比我快了46秒。 :) - Chris Cunningham
@Chris Cunningham：嘿，赞一个漂亮的布局！实际上早先有一个回答，现在已经被删除了，只是简单地提到了拉格朗日插值。我忽略了它，继续前进... - hardmath
这个实现并不太难，而且似乎给了我想要的系数。但是数值精度是一个问题。我的目标硬件只支持单精度浮点数也没有帮助。下一步我会研究共轭梯度法。 - lhumongous
3

出于稳定性原因，您不应该进行完全矩阵求逆。而是应该进行LU分解和前向-后向替换。除此之外，其他解决方案都是正确的。

你能为你所说的内容添加一个参考文献吗？我想阅读一下。 - Chris Cunningham
“Numerical Recipes” 这本书是首选，任何一本关于数值方法的好书也可以。 - duffymo
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接