基于C++的多项式曲线拟合代码实现与实战

可以使用信息准则自动判断：

$$ \text{AIC} = 2k - 2\ln(L),\quad \text{BIC} = k\ln(m) - 2\ln(L) $$

其中 $k = n+1$ 是参数个数，$L$ 是似然函数最大值。

最小二乘法：优雅背后的数学逻辑

它的思想极其朴素：让预测值和真实值之间的总体偏差尽可能小。

目标函数成了：

$$ S(\mathbf{a}) = \sum_{i=1}^m \left( y_i - \sum_{k=0}^n a_k x_i^k \right)^2 = | \mathbf{Y} - \mathbf{X}\mathbf{a} |^2 $$

正规方程推导：微积分遇上线性代数

为了最小化 $S(\mathbf{a})$，我们对每个 $a_j$ 求偏导并令其为零：

$$ \frac{\partial S}{\partial a_j} = -2 \sum_{i=1}^m \left( y_i - \sum_{k=0}^n a_k x_i^k \right) x_i^j = 0 $$

整理一下得到：

$$ \mathbf{X}^T \mathbf{X} \mathbf{a} = \mathbf{X}^T \mathbf{Y} $$

最终解为：

$$ \mathbf{a} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{Y} $$

中间那个 $(\mathbf{X}^T \mathbf{X})$ 很容易变得'又扁又长'，条件数爆炸。建议能不用正规方程就不用。

#include <Eigen/Dense>
Eigen::VectorXd solve_normal_equation(const Eigen::MatrixXd& X, const Eigen::VectorXd& Y) {
    Eigen::MatrixXd XtX = X.transpose() * X;
    Eigen::VectorXd XtY = X.transpose() * Y;
    return XtX.ldlt().solve(XtY); // 使用 LDLT 分解提升稳定性
}

凸优化视角：为什么你能相信这个解？

最小二乘的目标函数是一个关于 $\mathbf{a}$ 的二次函数，其 Hessian 矩阵为 $2\mathbf{X}^T\mathbf{X}$。只要 $\mathbf{X}$ 列满秩，$\mathbf{X}^T\mathbf{X}$ 就是正定的，意味着目标函数严格凸。

这意味着存在唯一的全局最小值。

flowchart TD
P[最小二乘问题] --> Q{X是否列满秩？}
Q -->|是 | R[XtX正定 → 凸函数]
Q -->|否 | S[存在无穷多解或无解]
R --> T[梯度为零 ⇒ 唯一最优解]
S --> U[需引入正则化或SVD]

数据处理实战：让理论真正跑起来

再完美的理论，碰上垃圾数据也会崩盘。

文件读取：别让格式毁了你的努力

我们需要一个健壮的文件解析器：

bool read_csv(const std::string& filename, std::vector<double>& x_vec, std::vector<double>& y_vec) {
    std::ifstream file(filename);
    if (!file.is_open()) {
        throw std::runtime_error("Cannot open file: " + filename);
    }
    std::string line;
    int line_num = 0;
    while (std::getline(file, line)) {
        ++line_num;
        std::istringstream ss(line);
        std::string cell_x, cell_y;
        if (!std::getline(ss, cell_x, ',')) continue;
        if (!std::getline(ss, cell_y, ',')) {
            throw std::invalid_argument("Incomplete data at line " + std::to_string(line_num));
        }
        try {
            double x = std::stod(cell_x);
            double y = std::stod(cell_y);
            x_vec.push_back(x);
            y_vec.push_back(y);
        } catch (...) {
            throw std::invalid_argument("Invalid number format at line " + std::to_string(line_num));
        }
    }
    if (x_vec.empty()) {
        throw std::length_error("No valid data found in file.");
    }
    return true;
}

异常值剔除：别让 outliers 毁了整锅汤

常用方法有两个：Z-score 和 IQR。推荐使用 IQR：

void remove_outliers_iqr(std::vector<double>& x_data, std::vector<double>& y_data) {
    if (y_data.size() < 4) return;
    std::vector<double> y_sorted = y_data;
    std::sort(y_sorted.begin(), y_sorted.end());
    double Q1 = percentile(y_sorted, 0.25);
    double Q3 = percentile(y_sorted, 0.75);
    double IQR = Q3 - Q1;
    double lower_bound = Q1 - 1.5 * IQR;
    double upper_bound = Q3 + 1.5 * IQR;
    std::vector<double> new_x, new_y;
    for (size_t i = 0; i < y_data.size(); ++i) {
        if (y_data[i] >= lower_bound && y_data[i] <= upper_bound) {
            new_x.push_back(x_data[i]);
            new_y.push_back(y_data[i]);
        }
    }
    x_data = std::move(new_x);
    y_data = std::move(new_y);
}

归一化：拯救病态矩阵的最后一道防线

标准做法是 Z-score 标准化：

$$ x' = \frac{x - \mu_x}{\sigma_x} $$

C++ 实现如下：

void normalize_data(std::vector<double>& x_data) {
    double sum = 0.0;
    for (double x : x_data) sum += x;
    double mean = sum / x_data.size();
    double var_sum = 0.0;
    for (double x : x_data) var_sum += (x - mean) * (x - mean);
    double stddev = std::sqrt(var_sum / x_data.size());
    if (stddev == 0.0) stddev = 1.0;
    for (double& x : x_data) {
        x = (x - mean) / stddev;
    }
}

高效求解策略：不只是快，更是稳

到了这一步，你已经有了干净的数据、合理的阶数、稳定的矩阵。接下来，就是见证奇迹的时刻：求解线性系统。

高斯消元：教学经典 vs 实战局限

中小规模（$n<50$）还能应付，但面对真实工业数据时，效率和稳定性都不够看。

SVD：终极武器，专治各种不服

当矩阵接近奇异、多重共线性严重、或者根本秩亏时，SVD 是唯一靠谱的选择。

其核心思想是分解：

$$ \mathbf{V} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^T $$

则最小二乘解为：

$$ \mathbf{a} = \mathbf{V} \boldsymbol{\Sigma}^{-1} \mathbf{U}^T \mathbf{y} $$

Eigen 中只需一行：

Eigen::JacobiSVD<Eigen::MatrixXd> svd(V, Eigen::ComputeThinU | Eigen::ComputeThinV);
return svd.solve(y);

虽然耗时稍长，但它能在其他方法失败的地方成功。

flowchart TD
Start[开始求解] --> IsStable{矩阵是否良态？}
IsStable -- 是 | UseQR[使用 QR 分解]
IsStable -- 否 | UseSVD[使用 SVD 分解]
UseQR --> Result1[快速获得近似解]
UseSVD --> Result2[获得最稳健解]
Result1 --> End
Result2 --> End

完整系统集成：从命令行到可视化

最后，让我们把这些模块组装成一个可用的工具。

模块化设计：CurveFitter 类登场

class CurveFitter {
private:
    std::vector<double> x, y;
    Eigen::VectorXd coefficients;
    int degree;
public:
    CurveFitter(const std::vector<double>& x_data, const std::vector<double>& y_data, int deg)
        : x(x_data), y(y_data), degree(deg) {}
    void build_vandermonde_matrix(Eigen::MatrixXd& V);
    bool solve_normal_equations();
    double evaluate(double x_val) const;
    double compute_rmse() const;
    double compute_r_squared() const;
    void save_results(const std::string& output_path) const;
};

RAII + Eigen 的组合拳，既保证资源安全，又获得极致性能。

命令行交互：让用户掌控节奏

while ((opt = getopt(argc, argv, "i:d:o:h")) != -1) {
    switch (opt) {
        case 'i': input_file = optarg; break;
        case 'd': degree = atoi(optarg); break;
        case 'o': output_file = optarg; break;
        case 'h': cout << "Usage: " << argv[0] << " -i data.txt -d 3 -o result.csv" << endl; return 0;
    }
}

支持 -i, -d, -o 参数，无需修改代码即可批量运行不同配置。

自动绘图：一键生成图表

输出 CSV 文件供 Gnuplot 使用：

set title "Polynomial Fit (degree $1)"
set xlabel "x"
set ylabel "y"
set grid
set terminal png size 800,600
set output 'fit_plot.png'
plot '$2' using 1:2 with points pt 7 ps 0.8 title "Data",\
     '$2' using 1:3 with lines lw 2 title "Fitted Curve"
EOF

C++ 中调用：

std::system("gnuplot plot.gp");

从此告别手动画图，真正实现自动化报告生成。

模型诊断：防止自我欺骗的关键一步

最后提醒一句：不要盲目相信你的拟合结果。

系统应具备基本的自检能力：

if (degree >= 8 && rmse_on_test_region > 2 * global_rmse) {
    std::cout << "[Warning] Possible overfitting detected. "
              << "Consider reducing polynomial degree." << std::endl;
}
double r2 = fitter.compute_r_squared();
if (r2 < 0.8) {
    std::cout << "[Suggestion] Low R² (" << r2 << "). Model may underfit. Try higher degree or check data noise." << std::endl;
}

还可以预留接口引入正则化：

Eigen::MatrixXd regularized_A = V.transpose() * V + lambda * Eigen::MatrixXd::Identity(degree+1, degree+1);
coefficients = regularized_A.ldlt().solve(V.transpose() * y_vector);

岭回归虽简单，却能在关键时刻提供保障。

写在最后：拟合不仅是技术，更是思维方式

回过头看，曲线拟合这件事，本质上是一场在有限信息中寻找秩序的努力。

我们用多项式作为语言，用最小二乘作为推理规则，用正则化作为防错机制。每一步都在平衡：表达力 vs 稳定性，拟合优度 vs 泛化能力，速度 vs 精度。

而这，也正是所有数据科学工作的缩影。

下次当你面对一堆混乱的数据时，不妨问问自己：我是真的发现了规律，还是只是记住了噪声？

毕竟，最优模型通常并非最复杂者，而是最能经得起质疑的那个。