本文介绍使用C++结合Eigen库实现多项式曲线拟合的方法。内容涵盖最小二乘法原理、范德蒙矩阵构建、正规方程与SVD求解策略对比、数据预处理（归一化、异常值剔除）及模型评估（RMSE、R²）。通过模块化设计与命令行交互，提供从理论到工程落地的完整解决方案，强调数值稳定性与过拟合防范。
在实际业务场景中，常面临数据规律不明确的问题。需要构建简洁的数学模型去逼近复杂的现实。多项式拟合是这个过程最基础也最关键的工具之一。
为什么要用多项式？因为它'看起来非线性'，但'算起来却是线性的'。
假设我们要拟合的数据满足这样一个关系：
$$ f(x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_n x^n $$
这个函数整体上看当然是非线性的。可注意！它对参数 $\mathbf{a} = [a_0, a_1, …, a_n]^T$ 的依赖却是完全线性的。
构造所谓的 。
$$ \mathbf{X} = \begin{bmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^n \ 1 & x_2 & x_2^2 & \cdots & x_2^n \ \vdots & \vdots & \vdots & \ddots & \vdots \ 1 & x_m & x_m^2 & \cdots & x_m^n \ \end{bmatrix} $$
每一行对应一个数据点，每一列代表一个幂次项。
std::vector<std::vector<double>> buildVandermonde(const std::vector<double>& x, int degree) {
int m = x.size();
int n = degree + 1;
std::vector<std::vector<double>> X(m, std::vector<double>(n));
for (int i = 0; i < m; ++i) {
for (int j = 0; j < n; ++j) {
X[i][j] = std::pow(x[i], j);
}
}
return X;
}
当 $x_i$ 比较大时，高阶项可能导致矩阵元素数量级差异巨大，称为'病态矩阵'。解决办法是归一化先行。
graph TD
A[原始数据点 xi, yi] --> B{是否需要预处理？}
B -->|是 | C[归一化/去噪]
B -->|否 | D[直接构造矩阵]
C --> D
D --> E[初始化空矩阵 X]
E --> F[遍历每个 xi]
F --> G[计算 1, xi, xi², ..., xin]
G --> H[填入矩阵第i行]
H --> I{是否所有点处理完毕？}
I -->|否 | F
I -->|是 | J[输出范德蒙矩阵 X]
现实中更常见的情况是：阶数越高，训练误差越低，但预测能力反而下降。
| 模型类型 | 偏差 | 方差 | 表现 |
|---|---|---|---|
| 低阶（如线性） | 高 | 低 | 欠拟合，错过趋势 |
| 中阶（如3~6次） | 适中 | 适中 | 最佳平衡区 |
| 高阶（>8次） | 低 | 高 | 过拟合，记住噪声 |
可以使用信息准则自动判断：
$$ \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;
}
常用方法有两个：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;
}
}
到了这一步，你已经有了干净的数据、合理的阶数、稳定的矩阵。接下来，就是见证奇迹的时刻：求解线性系统。
中小规模（$n<50$）还能应付，但面对真实工业数据时，效率和稳定性都不够看。
当矩阵接近奇异、多重共线性严重、或者根本秩亏时，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
最后，让我们把这些模块组装成一个可用的工具。
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 精度。
而这，也正是所有数据科学工作的缩影。
下次当你面对一堆混乱的数据时，不妨问问自己：我是真的发现了规律，还是只是记住了噪声？
毕竟，最优模型通常并非最复杂者，而是最能经得起质疑的那个。
微信公众号「极客日志」，在微信中扫描左侧二维码关注。展示文案：极客日志 zeeklog
使用加密算法（如AES、TripleDES、Rabbit或RC4）加密和解密文本明文。 在线工具，加密/解密文本在线工具，online
生成新的随机RSA私钥和公钥pem证书。 在线工具，RSA密钥对生成器在线工具，online
基于 Mermaid.js 实时预览流程图、时序图等图表，支持源码编辑与即时渲染。 在线工具，Mermaid 预览与可视化编辑在线工具，online
将字符串编码和解码为其 Base64 格式表示形式即可。 在线工具，Base64 字符串编码/解码在线工具，online
将字符串、文件或图像转换为其 Base64 表示形式。 在线工具，Base64 文件转换器在线工具，online
将 Markdown（GFM）转为 HTML 片段，浏览器内 marked 解析；与 HTML 转 Markdown 互为补充。 在线工具，Markdown 转 HTML在线工具，online