C++最小二乘拟合 (高阶最小二乘拟合)

原文链接：
https://blog.csdn.net/weixin_44344462/article/details/88850409

https://blog.csdn.net/xsz591541060/article/details/107532611

配置Eigen矩阵运算库
后续计算需要利用矩阵运算来求解拟合系数，用到Eigen矩阵运算库，配置方法可自行搜索或MacOs可以参考Mac配置Eigen库进行配置。

拟合原理
以二次曲线拟合为例，拟合曲线应满足以下方程：

即有：

由上面最后一个等式利用矩阵的转置与求逆，则可以得出其拟合曲线的系数W矩阵。

C++最小二乘拟合代码
以下代码理论上可以完成N阶的曲线拟合，需要对参数N进行设置，并提供正确的数据点。

#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Geometry>
#include <Eigen/Eigenvalues>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <cmath>

using namespace std;


int main() {
    // txt点数据文件路径
    const string fileName = "./points.txt";
    
    // 设置是几次拟合
    const int N = 2;
    
    // 创建两个vector
    vector<float> x, y;
    
    // 读取文件
    ifstream f(fileName);
    
    if (!f) {
        cout << "数据文件打开失败" << endl;
        exit(EXIT_FAILURE);
    }
    
    float tempx, tempy;
    
    while (f >> tempx >> tempy) {
        x.push_back(tempx);
        y.push_back(tempy);
    }
    
    if (x.size() != y.size()) {
        cout << "数据文件内容有误" << endl;
        exit(EXIT_FAILURE);
    }
    
    // 创建A矩阵
    Eigen::MatrixXd A(x.size(), N + 1);
    
    for (unsigned int i = 0; i < x.size(); ++i) {  // 遍历所有点
        
        for (int n = N, dex = 0; n >= 1; --n, ++dex) {  // 遍历N到1阶
            A(i, dex) = pow(x[i], n);
        }
        
        A(i, N) = 1;  //
    }
    
    // 创建B矩阵
    Eigen::MatrixXd B(y.size(), 1);
    
    for (unsigned int i = 0; i < y.size(); ++i) {
        B(i, 0) = y[i];
    }
  	
    // 创建矩阵W
    Eigen::MatrixXd W;
    W = (A.transpose() * A).inverse() * A.transpose() * B;
    
  	// 打印W结果
    cout << W << endl;
}

与设置的系数相吻合，但这是建立在噪声为0的情况下的结果。

#include <QCoreApplication>
#include <QtCore>
#include <QObject>
#include <QVector>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Geometry>
#include <Eigen/Eigenvalues>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <cmath>

using namespace Eigen;
using namespace std;


bool polyfit(const QVector<double> &x, const QVector<double> &y, quint8 dim, QVector<double> &z)
{
    Q_ASSERT(x.size() == y.size());
    quint8 N = 3;
    // 设置是几次拟合
    if(dim < 0 && dim > 20)
    {
        return false;
    }
    N = dim;
        // 创建A矩阵
        Eigen::MatrixXd A(x.size(), N + 1);

        for (unsigned int i = 0; i < x.size(); ++i) {  // 遍历所有点

            for (int n = N, dex = 0; n >= 1; --n, ++dex) {  // 遍历N到1阶
                A(i, dex) = pow(x[i], n);
            }

            A(i, N) = 1;  //
        }

        // 创建B矩阵
        Eigen::MatrixXd B(y.size(), 1);

        for (unsigned int i = 0; i < y.size(); ++i) {
            B(i, 0) = y[i];
        }

        // 创建矩阵W
        Eigen::MatrixXd W;
        W = (A.transpose() * A).inverse() * A.transpose() * B;



        // 打印W结果
        cout << W << endl;

        cout << "rows:"<<W.rows() << " cols:" << W.cols() << endl;

        double t;
        z.clear();
        for(int i = 0; i< W.rows()*W.cols(); i++)
        {
            t = W(i,0);
            z.push_back(t);
        }

    return true;
}

int main(int argc, char *argv[])
{
    QCoreApplication a(argc, argv);

    // 创建两个vector
    QVector<double> x, y, z;
    x.push_back(1);
    x.push_back(2);
    x.push_back(3);
    x.push_back(4);
    x.push_back(5);
    x.push_back(6);

    y.push_back(10);
    y.push_back(13);
    y.push_back(15);
    y.push_back(21);
    y.push_back(22);
    y.push_back(28);

    bool ret = polyfit(x, y, 5, z);

    qDebug()<<"QVector z: "<<z;


    return a.exec();
}