C++最小二乘拟合 (高阶最小二乘拟合)（附C++代码）

配置Eigen矩阵运算库

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

拟合原理

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

$a\ast X^2+b\ast X+c=y$

如要要拟合的4个点为：

(3, 6)

(5, 8)

(6, 12)

(7, 4)

则有：

即有：

$AW=B$

$A^{T}AW=A^{T}B$

$(A^{T}A{)}^{-1}{A}^{T}AW=(A^{T}A{)}^{-1}{A}^{T}B$

$W=(A^{T}A{)}^{-1}{A}^{T}B$

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

准备数据

由于Python的Matplotlib库可以很好地可视化数据，所以选用Python来生成数据文件——points.txt。

代码如下，程序中设置的系数为a = 2.2, b = 1.4, c = -1.3：

import numpy as np
import matplotlib.pyplot as plt


# 设置点数
num = 100

# 设置噪声系数
rate = 0

# 生成数据
x = np.linspace(-10, 10, num)
y = []
a = 2.2
b = 1.4
c = -1.3

for i in x:
    noise = np.random.random()
    temp = i + noise * rate
    y.append(a * pow(temp, 2) + b * temp + c)

# 绘制曲线
plt.plot(x, y)
plt.show()

# 写入文件
with open("./points.txt", 'w') as f:
    for i in range(num):
        f.write(str(x[i]))
        f.write("\t")
        f.write(str(y[i]))
        f.write("\n")

曲线为：

程序会在当前目录下生成points.txt文件。

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;
}

结果

程序输出为

2.2
1.4
-1.3

与设置的系数相吻合，但这是建立在噪声为0的情况下的结果，想要体验不同噪声对拟合结果的朋友可以在给出的Python程序中对rate参数进行修改。