最小二乘法拟合直线和曲线的C++代码

　　用最小二乘法拟合直线或者多项式曲线的C++代码如下（Windows平台）。ReadFile函数用于读入待拟合的离散数据，二维$(x,y)$坐标。LeastSquareFit函数用于拟合，它的输入是待拟合的离散数据向量和多项式曲线的阶数（order），即要拟合需要事先知道曲线的阶数，拟合得到的参数保存在ret向量中。

#include 
#include 
#include 
#include 

#include "Eigen-3.3/Eigen/Dense"

void ReadFile(const std::string& filename, 
    std::vector* x_vec,
    std::vector* y_vec) {
    std::ifstream in(filename, std::ios::in);
    std::string str;
    while (getline(in, str)) {
        std::stringstream input(str);
        double x, y;
        input >> x;
        x_vec->push_back(x);
        input >> y;
        y_vec->push_back(y);
        std::cout << std::fixed << std::setprecision(9);
        std::cout << " x = " << x << ", y = " << y << std::endl;
    }
    std::cout.unsetf(std::ios::fixed);
}

bool LeastSquareFit(std::vector& x, std::vector& y,
    size_t orders, Eigen::VectorXd& ret) {
    int ret_size = orders + 1;
    ret.resize(ret_size);

    for (int i = 0; i < ret_size; i++) {
        ret[i] = 0;
    }
    // input check
    if (x.size() < 2 || y.size() < 2 || x.size() != y.size() || orders < 1) {
        std::cout << "Error: data size less than 2 or x y size not equal" << std::endl;
        return false;
    }

    // map sample data from STL vector to eigen vector
    Eigen::Map sample_x(x.data(), x.size());
    Eigen::Map sample_y(y.data(), y.size());

    // Vandermonde matrix of x-axis coordinate vector of sample data
    Eigen::MatrixXd vandermonde_matrix(x.size(), orders + 1);
    Eigen::VectorXd vandermonde_column = sample_x;  // Vandermonde column

    // construct Vandermonde matrix column by column
    for (size_t i = 0; i < orders + 1; ++i) {
        if (0 == i) {
            vandermonde_matrix.col(0) = Eigen::VectorXd::Constant(x.size(), 1, 1);
            continue;
        }
        if (1 == i) {
            vandermonde_matrix.col(1) = vandermonde_column;
            continue;
        }

        vandermonde_column = vandermonde_column.array() * sample_x.array();
        vandermonde_matrix.col(i) = vandermonde_column;
    }

    ret = vandermonde_matrix.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(sample_y);
    return true;
}

int main()
{
    std::vector sample_x;
    std::vector sample_y;
    Eigen::VectorXd result;
    const std::string file_path = "C:/"; 
    const std::string filename = file_path + "data.txt";
    ReadFile(filename, &sample_x, &sample_y);
    if (LeastSquareFit(sample_x, sample_y, 1, result)) {
        std::cout << "result = " << result.transpose() << std::endl;
    }
    return 0;
}

　　我们用一些数据来验证一下，首先是简单的直线方程：$y=x+2$。
　　拟合结果如下，黑点为待拟合的离散数据（共51个点），红色直线为生成数据的直线$y=x+2$，绿线为拟合得到的直线。二者基本吻合，证明C++程序正确。

　　下面试一下曲线：$y=0.4x^3-2.3x^2+1.5x-2$。
　　拟合结果如下，黑点为待拟合的离散数据，红色曲线为生成数据的方程曲线，绿色曲线为根据拟合得到的参数画出的，结果同样比较接近，证明最小二乘法同样适用于曲线的情况。

　　该代码只依赖Eigen，所以可以在Linux下运行，只需要修改找到Eigen的头文件即可，例如在Ubuntu下Eigen一般安装在/usr/include/eigen3路径下，所以引用改为#include "eigen3/Eigen/Dense"即可。