多项式最小二乘插值法-C++【可直接复制粘贴/欢迎评论点赞】

目录

背景

优点

缺点

 代码

---

背景

多项式最小二乘插值法是一种数学方法，用于在给定的数据点集合上找到最佳拟合的多项式。与简单的多项式插值（如拉格朗日插值或牛顿插值）不同，最小二乘法不强制多项式精确通过每一个数据点，而是寻求一种多项式，使得所有数据点到多项式上对应点的垂直距离的平方和最小。这种方法在数据点受到噪声影响或数据点分布不均匀时尤其有用。

优点

1. 灵活性：最小二乘法不强制多项式精确通过每一个数据点，而是寻求一种多项式，使得所有数据点到多项式上对应点的垂直距离的平方和最小。这种灵活性使得它适用于各种类型的数据分布，包括噪声较大的数据。

2. 计算效率高：最小二乘法的计算过程相对简单，特别是对于线性最小二乘问题，存在高效的数值解法（如高斯消元法、QR分解等）。这使得多项式最小二乘插值法在处理大规模数据集时具有较高的计算效率。

3. 易于实现：在多种编程环境中，如Python、MATLAB等，都有现成的库和函数可以实现多项式最小二乘插值法，如Python的NumPy库中的polyfit函数。这大大降低了实现的难度和复杂度。

4. 广泛的应用领域：多项式最小二乘插值法被广泛用于回归分析、信号处理、图像处理、物理建模等多个领域。通过拟合数据中的趋势和模式，它可以帮助我们理解和预测变量的行为，并为决策提供有力的支持。

5. 可以处理非线性关系：虽然多项式本身是线性的（在多项式系数确定后），但通过选择合适的多项式次数，它可以用来拟合非线性关系的数据。

缺点

1. 过拟合和欠拟合问题：如果多项式的阶数选择不当，可能会出现过拟合或欠拟合的问题。过拟合是指多项式过于复杂，以至于它捕捉到了数据中的噪声而不是真正的趋势；欠拟合则是指多项式过于简单，无法充分描述数据的真实分布。

2. 对异常值敏感：最小二乘法对异常值（即远离其他数据点的数据点）较为敏感。异常值可能会显著影响多项式的拟合结果，导致拟合曲线偏离真实的数据分布。

3. 可能无法找到全局最优解：虽然多项式最小二乘插值法可以找到使得残差平方和最小的多项式，但这并不一定是全局最优解。在某些情况下，可能存在多个局部最优解，而算法可能陷入其中一个局部最优解。

4. 需要选择合适的多项式次数：选择合适的多项式次数是一个挑战。次数过低可能导致欠拟合，次数过高则可能导致过拟合。这通常需要根据具体的数据分布和问题的需求来进行权衡。

5. 可能不适用于所有类型的数据：虽然多项式最小二乘插值法具有广泛的应用领域，但它并不适用于所有类型的数据。例如，对于具有复杂非线性关系或周期性变化的数据，可能需要采用其他更复杂的模型来进行拟合。

 代码

 SolverEqGauss.h：

#include<iostream>
#include<iomanip>
#include<math.h>
using namespace std;
void SolverEqGauss(double** A, double* b, int n, double* x)
{
    int m, i, j, k;
    double** a;
    double* temp;
    //double temp[15];
    double d;
    a = new double* [n];
    temp = new double[n];
    for (i = 0; i < n; i++)
    {
        a[i] = new double[n + 1];
    }
    m = n;
    //赋值
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            a[i][j] = A[i][j];
        }
        a[i][n] = b[i];
    }
    
    for (k = 0; k < n - 1; k++) //找列主元最大值
    {
        double max = 0;
        int hang = 0, num = 0;
        for (i = k; i < n; i++)
        {
            if (fabs(a[i][k]) > max)
            {
                max = fabs(a[i][k]);
                hang = i;
            }
        }
        
        if (k != hang) //换行
        {
            for (i = 0; i < m + 1; i++)
            {
                temp[i] = a[k][i];
                a[k][i] = a[hang][i];
                a[hang][i] = temp[i];
            }
        }
        for (i = k + 1; i < m; i++) //消元
        {
            d = a[i][k] / a[k][k];
            for (j = 0; j < n + 1; j++)
            {
                a[i][j] = a[i][j] - d * a[k][j];
            }
        }
    }
    //memset(temp, 0, 15 * sizeof(double)); //将temp清0，准备存放解向量
    for (i = 0; i < n; i++)
    {
        temp[i] = 0;
    }
    for (i = m - 1; i >= 0; i--) //求解向量
    {
        d = 0;
        for (k = 0; k < n; k++)
        {
            d = d + temp[k] * a[i][k];
        }
        temp[i] = (a[i][n] - d) / a[i][i];
    }
    //cout << "此方程组的解向量转置为：("; //输出解向量
    for (i = 0; i < m; i++)
    {
        x[i] = temp[i];
        //cout << " " << fixed << setprecision(4) << temp[i];//5位小数
    }
    //cout << endl;
    //delete[]temp;
    for (int i = 0; i < n; i++)
    {
        delete a[i];
    }
    delete[]a;
}

 main.c：

#include<iostream>
#include "SolverEqGauss.h"
#include<iomanip>
#include<math.h>
using namespace std;
double Polyfit(double arrX[],double arrY[], int order, int n, double x, double Coeffi[]){
//其中arrX和arrY表示待插值数据点列表，order表示多项式阶次，n表示数据点个数，x表示待求点，Coeffi表示多项式系数列表。
int i,j,k;
double** A = new double* [order + 1];
    double* b = new double[order + 1];
    for (int i = 0; i < order + 1; i++)
    {
        A[i] = new double[order + 1];
    }
//求系数矩阵A
	for (i = 0; i <= order; i++)
	{
		for (j = 0; j <= order; j++)
		{

            for (int k = 0; k < n; k++)
		{
            A[i][j]+=pow(arrX[k], i + j);//不要A[i][j]=pow(arrX[k], i + j);会覆盖，变成无穷大
        }
		}
	}
//求方程右边的常数向量b
	for (i = 0; i <= order; i++)
	{
		b[i] = 0;
		for (int k = 0; k < n; k++)
		{

            b[i]+=arrY[k]*pow(arrX[k], i );


		}
	}
SolverEqGauss(A, b, order + 1, Coeffi);//调用高斯列主元消去法解法方程组(线性方程组)

	double res = 0;  //用拟合后的多项式计算校验点数据 
	for (int i = 0; i<order + 1; i++)
	{


            res += Coeffi[i] * pow(x, i); //Coeffi[i] 是多项式的第 i 个系数（按照降幂排列），pow(x, i) 是 x 的 i 次方，这两者的乘积就是多项式中 x 的 i 次项的数值

	}
//释放内存
    for (int i = 0; i < order + 1; i++)
    {
        delete A[i];
    }
    delete[]A;
    delete[]b;
    return res;

}


int main()
{
    int i, j; int order=0;double xs , ys = 0;
    cout << fixed << setprecision(6);
    int count ;//数据个数
    cin >> count;
    double x[200] = { };
    double y[200] = { };
    for (i = 0; i < count; i++) {
        cin >> x[i];
    }
    for (i = 0; i < count; i++) {
        cin >> y[i];
    }
    cin >> order;//需拟合次数
    cin >> xs;//输入待求点值
    double* CoeArray = new double[order + 1];
    ys = Polyfit(x, y, order, count, xs, CoeArray);
    for (i = 0; i <= order; i++)
    {
        cout << CoeArray[i] <<" ";
    }
    cout << ys;
    delete[]CoeArray;

return 0;
}