简介：最小二乘法是一种广泛应用于数据建模与曲线拟合的数学方法，通过使数据点到拟合曲线的误差平方和最小，实现最优逼近。本文介绍在C++环境下实现多项式曲线拟合的完整流程，涵盖数据准备、模型定义、矩阵构建、线性系统求解及结果评估等关键步骤。结合PolynomialFit.h与PolynomialFit.cpp文件，利用Eigen等数值库高效求解，支持可视化分析与拟合优化，适用于线性回归、信号处理和工程计算等实际场景。
1. 最小二乘法的基本原理与数学推导
最小二乘法的核心思想与数学建模
最小二乘法通过最小化观测值与模型预测值之间的 误差平方和 ，寻找最优参数估计。设观测向量为 $ \mathbf{y} \in \mathbb{R}^n $，设计矩阵为 $ \mathbf{A} \in \mathbb{R}^{n \times m} $，待估参数为 $ \mathbf{c} \in \mathbb{R}^m $，则拟合残差为 $ \mathbf{r} = \mathbf{y} - \mathbf{A}\mathbf{c} $，目标函数为：
J(\mathbf{c}) = |\mathbf{r}|^2 = (\mathbf{y} - \mathbf{A}\mathbf{c})^T(\mathbf{y} - \mathbf{A}\mathbf{c})
对该二次函数求极小值，对 $ \mathbf{c} $ 求梯度并令其为零：
\nabla_{\mathbf{c}} J = -2\mathbf{A}^T\mathbf{y} + 2\mathbf{A}^T\mathbf{A}\mathbf{c} = 0
得到 法方程 （Normal Equation）：
\mathbf{A}^T\mathbf{A}\mathbf{c} = \mathbf{A}^T\mathbf{y}
当 $ \mathbf{A}^T\mathbf{A} $ 可逆时，唯一解为：
\mathbf{c} = (\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T\mathbf{y}
该解在设计矩阵满秩时存在且唯一，体现了最小二乘估计的 无偏性 与 最小方差性 （高斯-马尔可夫定理），为后续多项式拟合提供理论基础。
2. 多项式模型构建与数据预处理
在最小二乘法的工程应用中，构建合理的数学模型是实现有效拟合的前提。多项式模型因其形式简洁、可微性强以及对非线性趋势具有良好的逼近能力，成为曲线拟合中最常用的函数类之一。然而，一个成功的拟合过程不仅依赖于模型的选择，更取决于前期的数据质量控制与数值稳定性保障。本章将系统阐述如何从实际观测数据出发，构建适配的多项式模型，并通过一系列预处理技术提升后续计算的准确性与鲁棒性。
2.1 多项式模型的定义与阶数选择
多项式模型作为最基础且广泛应用的函数逼近工具，在科学计算和工程建模中占据核心地位。其本质是利用一组幂函数的线性组合来逼近未知的真实关系 $ y = f(x) $。该模型的形式灵活，可通过调节阶数适应不同复杂度的趋势特征。但阶数过高易导致过拟合，阶数过低则无法捕捉关键变化趋势。因此，合理选择多项式阶数是平衡拟合精度与泛化性能的关键环节。
2.1.1 多项式函数的形式化表达
设输入变量为 $ x \in \mathbb{R} $，输出响应为 $ y \in \mathbb{R} $，则一个 $ n $ 阶多项式模型可表示为：
y(x) = c_0 + c_1 x + c_2 x^2 + \cdots + c_n x^n = \sum_{i=0}^{n} c_i x^i
其中，$ c_0, c_1, …, c_n \in \mathbb{R} $ 是待估计的系数向量 $ \mathbf{c} \in \mathbb{R}^{n+1} $，$ n $ 为多项式的最高阶次。此表达式构成了一组基函数 $ {1, x, x^2, …, x^n} $ 的线性展开，属于广义线性模型范畴。
这种形式的优势在于其解析性质优良——无限可微、连续、易于积分与求导，适合用于物理系统的动态建模或信号趋势提取。例如，在传感器校准中常使用二次多项式修正非线性偏移；在轨迹规划中采用五次多项式保证加速度连续。
值得注意的是，尽管多项式本身是非线性的函数（关于 $ x $），但在参数 $ \mathbf{c} $ 上是线性的，这使得我们可以将其嵌入线性最小二乘框架进行高效求解。
2.1.2 不同阶数对拟合能力的影响机制
多项式阶数直接影响模型的表达能力。低阶模型（如一次或二次）只能描述简单趋势，而高阶模型理论上可以逼近任意光滑函数（根据魏尔斯特拉斯逼近定理）。但实践中需权衡“拟合能力”与“数值稳定性”。
下表展示了不同阶数多项式的主要特性对比：
| 阶数 | 函数形式 | 拟合能力 | 典型应用场景 |
|---|---|---|---|
| 1 | $ c_0 + c_1x $ | 线性趋势 | 趋势分析、斜率估计 |
| 2 | $ c_0 + c_1x + c_2x^2 $ | 单峰/单谷曲线 | 抛物运动、热膨胀补偿 |
| 3 | 三次多项式 | S形曲线、拐点 | 控制系统响应拟合 |
| 4-5 | 四至五次 | 多个极值点 | 轨迹插值、平滑过渡 |
| ≥6 | 高阶多项式 | 极强局部拟合能力 | 易引发震荡、过拟合 |
随着阶数增加，模型自由度上升，能够更好地匹配训练数据中的细节波动。然而，这也带来了严重的副作用： 龙格现象（Runge’s phenomenon） ——即在区间端点附近出现剧烈振荡，尤其当节点等距分布时更为显著。
// 示例：生成不同阶数多项式拟合效果（伪代码）
#include <vector>
#include <cmath>
std::vector<double> evaluate_polynomial(const std::vector<double>& coeffs,
const std::vector<double>& x_vals) {
std::vector<double> y_vals;
for (double x : x_vals) {
double y = 0.0;
for (size_t i = 0; i < coeffs.size(); ++i) {
y += coeffs[i] * pow(x, i); // 计算 c_i * x^i
}
y_vals.push_back(y);
}
return y_vals;
}
代码逻辑逐行解读：
- 第5行：定义函数
evaluate_polynomial，接收系数向量coeffs和输入点集x_vals。 - 第7–13行：遍历每个输入点 $ x $，计算对应的多项式输出值。
- 第10–11行：内层循环累加每一项 $ c_i x^i $，完成多项式求值。
- 参数说明：
-
coeffs: 长度为 $ n+1 $ 的向量，对应 $ [c_0, c_1, …, c_n] $ -
x_vals: 待预测的输入点集合 - 返回值：对应每个 $ x $ 的预测输出 $ y $
该函数可用于可视化不同阶数下的拟合曲线形态，辅助判断是否出现异常波动。
2.1.3 模型复杂度与泛化性能的权衡分析
选择多项式阶数本质上是在 偏差-方差权衡（Bias-Variance Tradeoff） 中寻找最优平衡点。低阶模型偏差大、方差小；高阶模型偏差小、方差大。
我们可以通过构建 验证误差 vs. 阶数 曲线来指导选择：
graph LR
A[多项式阶数 n=1] --> B[训练误差较高<br>验证误差较高]
C[n=2] --> D[训练误差下降<br>验证误差最低]
E[n=4] --> F[训练误差很低<br>验证误差回升]
G[n=6+] --> H[训练误差接近零<br>验证误差急剧上升]
style C stroke:#2ecc71,stroke-width:2px
style E stroke:#e74c3c,stroke-width:2px
上图表明：当阶数适中时，验证误差达到最小值，此时模型具备最佳泛化能力。超过该阈值后，虽然训练误差持续下降，但模型开始记忆噪声而非学习规律，导致测试表现恶化。
为此，推荐采用以下策略确定最优阶数：
1. 使用交叉验证（Cross-validation）评估各阶模型的平均预测误差；
2. 引入信息准则如 AIC 或 BIC 进行自动选择；
3. 结合领域知识设定上限（如机械系统通常不超过四阶）；
最终目标不是使训练误差最小，而是让模型在新数据上保持稳定预测能力。
2.2 数据准备与预处理技术
高质量的输入数据是获得可靠拟合结果的基础。原始采集数据往往包含缺失值、异常点、量纲差异等问题，若不加以处理，可能导致设计矩阵病态、求解失败或结果失真。因此，必须实施系统的数据清洗与变换流程，确保数据满足最小二乘法的基本假设：独立、同分布、零均值误差。
2.2.1 原始数据的采集与清洗策略
真实世界的数据常来自传感器、日志文件或数据库记录，存在多种质量问题。典型问题包括：
- 缺失值（NaN 或空字段）
- 重复采样点
- 时间戳错乱或非单调
- 单位不统一（如温度混用摄氏与华氏）
清洗步骤应遵循如下流程：
flowchart TD
Start[开始] --> Load[加载原始数据]
Load --> CheckMissing[检查缺失值]
CheckMissing -- 存在 --> Impute[插值填补或删除]
CheckMissing -- 无 --> RemoveDup[去除重复点]
RemoveDup --> SortData[按x排序]
SortData --> ValidateRange[检查x/y范围合理性]
ValidateRange --> Output[输出清洗后数据]
对于缺失值，常见处理方式有：
- 删除含有缺失的样本（适用于少量缺失）
- 使用线性插值、样条插值或均值填充（适用于时间序列）
示例代码（C++风格）展示如何剔除无效数据：
#include <vector>
#include <algorithm>
#include <cmath>
struct DataPoint {
double x, y;
};
std::vector<DataPoint> clean_data(std::vector<DataPoint> raw_data) {
std::vector<DataPoint> cleaned;
for (const auto& pt : raw_data) {
if (!std::isnan(pt.x) && !std::isnan(pt.y)) { // 排除NaN
if (std::isfinite(pt.x) && std::isfinite(pt.y)) { // 排除非有限值
cleaned.push_back(pt);
}
}
}
// 按x升序排列
std::sort(cleaned.begin(), cleaned.end(),
[](const DataPoint& a, const DataPoint& b) {
return a.x < b.x;
});
return cleaned;
}
代码逻辑逐行解读：
- 第9–14行：遍历原始数据，仅保留 x 和 y 均为有效数值的点。
- std::isnan() 判断是否为 NaN， std::isfinite() 排除无穷大。
- 第17–21行：使用 lambda 表达式按 x 排序，便于后续构造范德蒙德矩阵。
- 输出为结构化数据列表，便于后续矩阵构造。
此步骤虽简单，却极大提升了后续计算的稳定性。
2.2.2 异常值检测与剔除方法（如3σ准则）
异常值（outliers）会严重扭曲最小二乘解，因其以平方误差为目标函数，对大残差极为敏感。一种经典且高效的检测方法是 3σ准则（Three-Sigma Rule） ，基于正态分布假设：约99.7%的数据落在均值±3倍标准差范围内。
算法步骤如下：
1. 对初始拟合残差 $ r_i = y_i - \hat{y}_i $ 计算均值 $ \mu_r $ 和标准差 $ \sigma_r $
2. 标记满足 $ |r_i - \mu_r| > 3\sigma_r $ 的点为异常值
3. 可选迭代：重新拟合并再次检测，称为“迭代3σ剔除”
#include <vector>
#include <numeric>
#include <cmath>
std::vector<int> detect_outliers_3sigma(const std::vector<double>& residuals,
double threshold_factor = 3.0) {
// 计算残差均值
double mean_res = std::accumulate(residuals.begin(), residuals.end(), 0.0) / residuals.size();
// 计算标准差
double var = 0.0;
for (double r : residuals) {
var += (r - mean_res) * (r - mean_res);
}
double std_dev = std::sqrt(var / residuals.size());
// 找出超出 threshold_factor * σ 的索引
std::vector<int> outlier_indices;
for (size_t i = 0; i < residuals.size(); ++i) {
if (std::abs(residuals[i] - mean_res) > threshold_factor * std_dev) {
outlier_indices.push_back(i);
}
}
return outlier_indices;
}
参数说明：
- residuals : 输入残差向量，长度为 $ m $
- threshold_factor : 默认为3，可根据需求调整（如2.5用于宽松检测）
- 返回值：异常点在原数据中的索引列表
该方法计算开销小，适合实时系统。但对于非高斯噪声或多模态分布可能误判，建议结合箱线图（IQR）等鲁棒方法联合使用。
2.2.3 数据归一化与标准化在数值稳定性中的作用
高阶多项式涉及 $ x^k $（$ k $ 较大）运算，若原始 $ x $ 取值范围较宽（如 $ x \in [0, 1000] $），则 $ x^n $ 会导致数量级爆炸，造成设计矩阵元素差异巨大，进而引起 病态条件数（ill-conditioning） ，严重影响法方程求解精度。
解决方案是对输入变量进行 归一化（Normalization）或标准化（Standardization） ：
| 方法 | 公式 | 特性 |
|---|---|---|
| 最小-最大归一化 | $ x’ = \frac{x - x_{\min}}{x_{\max} - x_{\min}} $ | 映射到 [0,1] 区间 |
| Z-score标准化 | $ x’ = \frac{x - \mu_x}{\sigma_x} $ | 均值为0，方差为1 |
推荐使用 Z-score标准化 ，因其对异常值相对稳健，且更适合后续正则化扩展。
示例代码实现标准化：
std::vector<double> standardize(const std::vector<double>& data) {
double mean = std::accumulate(data.begin(), data.end(), 0.0) / data.size();
double variance = 0.0;
for (double val : data) {
variance += (val - mean) * (val - mean);
}
variance /= data.size();
double stddev = std::sqrt(variance);
std::vector<double> normalized;
for (double val : data) {
normalized.push_back((val - mean) / stddev);
}
return normalized;
}
优势分析：
- 经过标准化后，$ x $ 分布集中于 $[-2, 2]$ 附近，避免 $ x^6 $ 等项过大；
- 提高矩阵 $ \mathbf{A}^T\mathbf{A} $ 的条件数，减少舍入误差累积；
- 加速梯度类优化算法收敛（若后续引入迭代求解）；
注意：变换后的模型参数需反变换回原始尺度才能解释物理意义。
2.3 设计矩阵与观测向量的构造实践
最小二乘法的核心是将非线性拟合问题转化为线性代数问题。这一转化依赖于正确构造设计矩阵 $ \mathbf{A} $ 和观测向量 $ \mathbf{y} $。尤其在多项式拟合中，设计矩阵即为著名的 范德蒙德矩阵（Vandermonde Matrix） ，其结构直接决定求解难度。
2.3.1 从样本点生成范德蒙德矩阵的过程详解
设有 $ m $ 个采样点 $ (x_1, y_1), (x_2, y_2), …, (x_m, y_m) $，拟合 $ n $ 阶多项式，则设计矩阵 $ \mathbf{A} \in \mathbb{R}^{m \times (n+1)} $ 定义为：
\mathbf{A} =
\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}
每一行对应一个样本点的幂次展开，列数等于模型参数个数。
C++ 实现如下：
#include <vector>
#include <Eigen/Dense>
Eigen::MatrixXd build_vandermonde(const std::vector<double>& x_vals, int degree) {
int m = x_vals.size(); // 样本数
int n = degree + 1; // 列数 = 阶数+1
Eigen::MatrixXd A(m, n);
for (int i = 0; i < m; ++i) {
A(i, 0) = 1.0; // 常数项
for (int j = 1; j < n; ++j) {
A(i, j) = A(i, j-1) * x_vals[i]; // 递推计算 x^j
}
}
return A;
}
逻辑分析：
- 使用 Eigen::MatrixXd 存储矩阵，便于后续线性代数操作；
- 内层循环采用递推方式计算幂次（$ x^j = x^{j-1} \cdot x $），减少重复调用 pow() 函数；
- 时间复杂度 $ O(mn) $，空间复杂度 $ O(mn) $；
该矩阵一旦构造完成，即可参与法方程构建：$ (\mathbf{A}^T\mathbf{A})\mathbf{c} = \mathbf{A}^T\mathbf{y} $
2.3.2 观测向量Y的组织方式与维度匹配原则
观测向量 $ \mathbf{y} \in \mathbb{R}^m $ 是所有输出测量值组成的列向量：
\mathbf{y} =
\begin{bmatrix}
y_1 \ y_2 \ \vdots \ y_m
\end{bmatrix}
必须确保其长度与设计矩阵行数一致。常见错误包括：
- 数据清洗后未同步更新 $ \mathbf{y} $
- 归一化仅作用于 $ x $ 而忽略 $ y $ 的缩放（影响残差尺度）
使用 Eigen 构造示例：
Eigen::VectorXd build_observation_vector(const std::vector<double>& y_vals) {
return Eigen::Map<const Eigen::VectorXd>(y_vals.data(), y_vals.size());
}
维度匹配原则总结：
| 对象 | 维度 | 说明 |
|---|---|---|
| $ \mathbf{A} $ | $ m \times (n+1) $ | $ m $: 样本数，$ n+1 $: 参数数 |
| $ \mathbf{y} $ | $ m \times 1 $ | 必须与 $ \mathbf{A} $ 行数相同 |
| $ \mathbf{c} $ | $ (n+1) \times 1 $ | 输出系数向量 |
任何维度不匹配都将导致矩阵乘法失败或求解异常。
2.3.3 高阶多项式带来的病态矩阵问题及其初步应对
随着多项式阶数升高，范德蒙德矩阵趋于 高度病态（ill-conditioned） ，表现为 $ \text{cond}(\mathbf{A}) \gg 1 $，轻微扰动即可引起解的巨大变化。
例如，当 $ x \in [0,10] $，$ n=8 $ 时，$ \mathbf{A} $ 中元素跨度可达 $ 1 $ 到 $ 10^8 $，导致浮点精度损失严重。
缓解措施包括：
1. 输入变量标准化 （已在前节说明）
2. 使用正交多项式基 （如勒让德多项式）替代幂基
3. 引入正则化项 （见第四章）
当前阶段，优先执行标准化是最简便有效的手段。
综上，本节完整展示了从原始数据到矩阵构造的全流程，强调了每一个环节对最终结果的影响。唯有严谨对待每一步预处理与建模决策，方能获得稳健可靠的拟合结果。
3. 法方程求解与数值计算实现
在最小二乘法的完整流程中，构建出设计矩阵 $ A $ 与观测向量 $ \mathbf{y} $ 后，核心任务转化为求解一个关键线性系统—— 法方程（Normal Equation） 。该方程的形式为：
(A^T A)\mathbf{c} = A^T \mathbf{y}
其中 $ \mathbf{c} $ 是待估计的多项式系数向量。虽然形式简洁，但其数值求解过程涉及矩阵代数、条件稳定性以及算法效率等多个层面的技术挑战。本章将从法方程的数学结构出发，深入分析其代数特性与潜在的数值问题，并系统比较主流线性求解器的性能差异，最终通过现代C++高性能数值库（如Eigen）实现端到端的拟合计算流程。
3.1 法方程的代数结构与求解路径
法方程是连接模型构造与参数估计的关键桥梁。它并非直接来自原始数据关系，而是通过对误差平方和函数进行梯度优化推导而来。理解其代数本质有助于我们判断解的存在性、唯一性以及数值稳定性。
3.1.1 推导 (A^T * A) * C = A^T * Y 的完整过程
考虑一组观测数据点 $ (x_i, y_i), i=1,\dots,n $，假设我们使用一个 $ m $ 阶多项式进行拟合：
f(x) = c_0 + c_1 x + c_2 x^2 + \cdots + c_m x^m
对于每个样本点 $ x_i $，预测值可表示为：
\hat{y} i = \sum {j=0}^{m} c_j x_i^j
定义误差项为 $ e_i = y_i - \hat{y}_i $，则总误差平方和为：
S(\mathbf{c}) = \sum_{i=1}^{n} e_i^2 = | \mathbf{y} - A\mathbf{c} |^2
其中 $ A \in \mathbb{R}^{n \times (m+1)} $ 是范德蒙德矩阵，每一行为 $ [1, x_i, x_i^2, \dots, x_i^m] $，$ \mathbf{y} \in \mathbb{R}^n $ 为观测值向量，$ \mathbf{c} \in \mathbb{R}^{m+1} $ 为未知系数向量。
为了使 $ S(\mathbf{c}) $ 最小化，对 $ \mathbf{c} $ 求偏导并令其为零：
\nabla_{\mathbf{c}} S(\mathbf{c}) = -2A^T(\mathbf{y} - A\mathbf{c}) = 0
整理得：
A^T A \mathbf{c} = A^T \mathbf{y}
这就是著名的 法方程 。当 $ A^T A $ 可逆时，存在唯一解：
\mathbf{c} = (A^T A)^{-1} A^T \mathbf{y}
这一推导揭示了最小二乘解的闭式表达，也为后续的数值实现提供了理论依据。
注意 ：此闭式解仅在 $ A^T A $ 正定时成立，且实际中往往不推荐直接求逆，原因将在下文详述。
3.1.2 对称正定矩阵的特性及其在求解中的优势
矩阵 $ G = A^T A $ 被称为 Gram矩阵 或 信息矩阵 ，具有如下重要性质：
| 性质 | 描述 |
|---|---|
| 对称性 | $ G^T = (A^T A)^T = A^T A = G $ |
| 半正定性 | 对任意非零向量 $ \mathbf{v} $，有 $ \mathbf{v}^T G \mathbf{v} = |A\mathbf{v}|^2 \geq 0 $ |
| 正定性条件 | 当且仅当 $ A $ 列满秩（即rank($A$) = $m+1$），$ G $ 正定，从而可逆 |
由于 $ G $ 是对称矩阵，可以利用专门针对对称结构的求解器（如Cholesky分解），显著提升计算效率和数值精度。
例如，若 $ G $ 正定，则存在唯一的下三角矩阵 $ L $，使得：
G = LL^T
然后通过前向代入和后向代入即可高效求解：
// 示例伪代码：Cholesky分解求解 Ax=b
Matrix G = A.transpose() * A;
Vector b = A.transpose() * y;
// 分解 G = LL^T
LLT<MatrixXd> llt(G);
Vector c = llt.solve(b); // 内部自动执行前/后代入
逻辑分析 ：
A.transpose() * A构造 Gram 矩阵；LLT<MatrixXd>是 Eigen 库提供的 Cholesky 分解类；.solve(b)自动判断是否已分解，并调用相应的三角求解过程；- 参数说明：输入必须是对称正定矩阵，否则分解可能失败或产生错误结果。
这种结构化处理相比通用高斯消元能减少约一半的浮点运算量，是工程实践中常用的加速手段。
mermaid 流程图：法方程求解决策路径
graph TD
A[开始求解法方程] --> B{A是否列满秩?}
B -- 是 --> C[G = A^T A 是否正定?]
C -- 是 --> D[使用Cholesky分解]
C -- 否 --> E[尝试LDL^T分解]
B -- 否 --> F[矩阵奇异, 需正则化或降阶]
D --> G[求解 Lz = b; L^T x = z]
E --> G
G --> H[输出系数向量 c]
该流程图展示了面对不同矩阵属性时应采取的策略分支，强调了“先验判断”在数值计算中的重要性。
3.1.3 数值不稳定性的来源与条件数分析
尽管法方程形式优美，但在高阶多项式拟合中极易出现严重的 数值病态 问题。其根本原因在于 $ A^T A $ 的 条件数（Condition Number） 过大。
条件数定义为：
\kappa(G) = |G| \cdot |G^{-1}|
通常采用谱条件数：$ \kappa_2(G) = \frac{\sigma_{\max}}{\sigma_{\min}} $，其中 $ \sigma $ 为奇异值。
当 $ \kappa(G) \gg 1 $ 时，微小的数据扰动会导致解的巨大偏差。例如，在等距节点上构造高阶范德蒙德矩阵时，随着阶数增加，条件数呈指数级增长：
| 多项式阶数 $ m $ | 条件数 $ \kappa_2(A^T A) $ （近似） |
|---|---|
| 5 | $ 10^3 $ |
| 10 | $ 10^7 $ |
| 15 | $ 10^{12} $ |
| 20 | $ >10^{16} $（双精度失效） |
这意味着即使使用双精度浮点数（有效位约16位），当阶数超过15时，有效数字几乎全部丢失。
解决方案方向 包括：
- 数据归一化（如将 $ x \in [-1,1] $）
- 使用正交多项式基代替单项式基（Legendre、Chebyshev）
- 改用QR或SVD分解避免显式构造 $ A^T A $
这些方法将在后续章节展开，此处仅指出： 盲目使用法方程闭式解在高维场景下是危险的 。
3.2 线性系统的经典求解算法比较
虽然法方程提供了一个明确的目标方程，但如何高效、稳定地求解 $ (A^T A)\mathbf{c} = A^T \mathbf{y} $，仍是数值线性代数的核心课题。不同的分解方法在稳定性、速度和内存占用方面各有优劣。
3.2.1 高斯消元法：原理与编程实现步骤
高斯消元法是最基础的线性方程组求解方法，通过行变换将增广矩阵化为上三角形式，再回代求解。
设我们要解 $ G\mathbf{c} = \mathbf{b} $，其中 $ G = A^T A $，$ \mathbf{b} = A^T \mathbf{y} $。
算法步骤 如下：
- 构造增广矩阵 $ [G|\mathbf{b}] $
- 对每列 $ k $ 执行主元选择（部分主元）
- 将第 $ k $ 行以下所有行减去第 $ k $ 行的倍数，消去下方元素
- 得到上三角矩阵后，从最后一行开始回代求解
以下是简化的C++实现片段（忽略边界检查）：
void gaussianElimination(MatrixXf& G, VectorXf& b) {
int n = G.rows();
for (int k = 0; k < n - 1; ++k) {
// 部分主元选取
int maxRow = k;
for (int i = k + 1; i < n; ++i)
if (fabs(G(i,k)) > fabs(G(maxRow,k)))
maxRow = i;
G.row(k).swap(G.row(maxRow));
std::swap(b(k), b(maxRow));
// 消元
for (int i = k + 1; i < n; ++i) {
float factor = G(i,k) / G(k,k);
G.block(i,k+1,1,n-k-1) -= factor * G.block(k,k+1,1,n-k-1);
b(i) -= factor * b(k);
G(i,k) = 0; // 显式置零
}
}
// 回代
for (int i = n - 1; i >= 0; --i) {
for (int j = i + 1; j < n; ++j)
b(i) -= G(i,j) * b(j);
b(i) /= G(i,i);
}
}
逐行解读与参数说明 ：
MatrixXf和VectorXf：Eigen 中动态大小的单精度矩阵/向量类型；- 外层循环控制主元列 $ k $，内层寻找最大绝对值行以提高稳定性；
swap()操作交换整行，防止除以接近零的主元；block(i,k+1,1,n-k-1)提取子块进行向量化操作，避免逐个赋值；- 回代阶段从底向上逐步解出每个 $ c_i $；
- 时间复杂度为 $ O(n^3) $，空间复杂度 $ O(n^2) $。
缺点 ：直接作用于 $ A^T A $ 会加剧舍入误差；无法利用对称性；不适合稀疏结构。
3.2.2 LU分解：带部分主元的分解策略
LU分解将矩阵分解为一个下三角矩阵 $ L $ 和上三角矩阵 $ U $，即 $ G = PLU $，其中 $ P $ 为置换矩阵（用于主元重排）。
相较于高斯消元，LU分解可分离“分解”与“求解”阶段，适合多次右端项求解。
在Eigen中调用方式如下：
PartialPivLU<MatrixXd> lu(G);
VectorXd c = lu.solve(b);
其内部流程如下：
graph LR
A[原始矩阵 G] --> B[PA = LU 分解]
B --> C{是否成功?}
C -- 是 --> D[Ly = Pb]
D --> E[Ux = y]
E --> F[返回解 x]
C -- 否 --> G[报错: 奇异矩阵]
优点分析 ：
- 支持部分主元，增强稳定性；
- 分解后可用于多个右端项快速求解；
- 相比直接求逆更安全；
局限性 ：
- 仍需构造 $ A^T A $，继承其病态性；
- 不适用于非方阵（而最小二乘本身常为超定系统）；
因此，LU更适合中小型、良态的问题。
3.2.3 QR分解：基于Householder变换的稳定求解方案
QR分解是解决最小二乘问题最稳健的方法之一。它将设计矩阵 $ A \in \mathbb{R}^{n \times p} $ 分解为：
A = QR
其中 $ Q \in \mathbb{R}^{n \times n} $ 正交（$ Q^T Q = I $），$ R \in \mathbb{R}^{n \times p} $ 上梯形（前 $ p $ 行为上三角，其余为零）。
代入最小二乘目标：
| \mathbf{y} - A\mathbf{c} |^2 = | \mathbf{y} - QR\mathbf{c} |^2 = | Q^T \mathbf{y} - R\mathbf{c} |^2
令 $ Q^T \mathbf{y} = \begin{bmatrix} d \ e \end{bmatrix} $，则最优解满足：
R_1 \mathbf{c} = d
其中 $ R_1 $ 为 $ R $ 的前 $ p $ 行组成的上三角矩阵。
该方法的优势在于：
- 无需构造 $ A^T A $，从根本上规避病态放大；
- Householder反射保证高度数值稳定性；
- 适用于秩亏情况（结合截断SVD）
Eigen 实现示例：
HouseholderQR<MatrixXd> qr(A);
VectorXd c = qr.solve(y);
参数说明与逻辑分析 ：
HouseholderQR使用 Householder 反射构造正交矩阵；.solve(y)自动执行 $ Q^Ty $ 并求解上三角系统；- 即使 $ A $ 不满秩，也能返回最小范数解（若启用
ColPivHouseholderQR则更强健）；- 计算成本略高于LU，但稳定性远胜于法方程。
综上， 在高阶多项式拟合中，强烈建议优先采用QR分解而非法方程闭式解 。
3.3 基于C++数值库的高效实现
现代科学计算依赖于高度优化的底层线性代数库。C++生态中， Eigen 是最广泛使用的开源矩阵库，支持表达式模板、SIMD加速和BLAS/LAPACK接口绑定，非常适合嵌入式与高性能场景。
3.3.1 使用Eigen库进行矩阵运算的接口调用
以下是一个完整的最小二乘拟合实现示例，涵盖从数据输入到系数求解全过程：
#include <Eigen/Dense>
#include <vector>
#include <iostream>
using namespace Eigen;
using namespace std;
VectorXd polynomialFit(const vector<double>& x,
const vector<double>& y,
int degree) {
int n = x.size();
MatrixXd A(n, degree + 1);
// 构造范德蒙德矩阵
for (int i = 0; i < n; ++i) {
A(i, 0) = 1.0;
for (int j = 1; j <= degree; ++j) {
A(i, j) = A(i, j-1) * x[i];
}
}
VectorXd Y(y.data(), n);
// 使用QR分解求解
return A.householderQr().solve(Y);
}
int main() {
vector<double> x = {0, 1, 2, 3, 4};
vector<double> y = {1, 3, 7, 13, 21}; // 接近 y = x^2 + x + 1
int deg = 2;
VectorXd coeffs = polynomialFit(x, y, deg);
cout << "Fitted coefficients:\n" << coeffs.transpose() << endl;
return 0;
}
代码逐行解析 ：
VectorXd：动态长度双精度向量；MatrixXd：动态大小双精度矩阵；- 范德蒙德矩阵逐行构造，利用递推 $ x^j = x^{j-1} \cdot x $ 减少幂运算开销；
householderQr()返回 QR 分解对象，.solve()完成求解；- 输出应接近
[1, 1, 1]，对应 $ 1 + x + x^2 $；性能提示 ：
- 若数据量大，可用
Map<VectorXd>避免拷贝；- 开启编译器优化
-O3 -march=native启用向量化；- 使用
ColPivHouseholderQR提高秩亏鲁棒性。
3.3.2 BLAS/LAPACK底层函数在大规模计算中的优化支持
Eigen 在后台可链接 Intel MKL 或 OpenBLAS，从而调用高度优化的 Fortran 实现。例如，矩阵乘法 A*B 会被映射到底层 dgemm_ 函数，充分利用CPU缓存层级和多核并行。
典型链接配置（CMake）：
find_package(OpenBLAS REQUIRED)
target_link_libraries(myapp Eigen3::Eigen ${OPENBLAS_LIBRARIES})
此时， A.transpose() * A 等操作将获得数量级的速度提升，尤其在 $ n > 1000 $ 时效果显著。
| 操作规模 $ n \times m $ | Eigen原生（ms） | OpenBLAS加速（ms） | 加速比 |
|---|---|---|---|
| 1000×100 | 85 | 12 | 7.1× |
| 2000×200 | 680 | 89 | 7.6× |
这表明，在工业级应用中，合理集成BLAS/LAPACK是必不可少的性能保障。
3.3.3 实际代码片段展示：从矩阵构建到系数求解的全过程
综合以上技术，封装一个健壮的拟合函数：
#include <Eigen/QR>
#include <stdexcept>
struct PolynomialFitter {
static VectorXd fit(const Ref<const VectorXd>& x,
const Ref<const VectorXd>& y,
int degree) {
if (x.size() != y.size())
throw invalid_argument("x and y size mismatch");
if (degree < 0 || degree >= x.size())
throw invalid_argument("invalid degree");
int n = x.size();
MatrixXd A(n, degree + 1);
A.col(0).setOnes();
for (int j = 1; j <= degree; ++j)
A.col(j) = A.col(j-1).array() * x.array();
HouseholderQR<MatrixXd> qr(A);
if (qr.rank() < degree + 1)
throw runtime_error("Design matrix rank deficient");
return qr.solve(y);
}
};
扩展能力说明 ：
- 使用
Ref<>允许传入数组、子向量等多种格式；- 添加异常处理防止非法输入；
- 可进一步加入正则化项形成岭回归；
- 支持模板泛型以兼容 float/double 类型。
该模块已具备投入生产的潜力，体现了现代C++科学编程的最佳实践。
总结性观察 ：
本章揭示了从法方程到实际代码之间的鸿沟——理论上简单的 $ (A^TA)^{-1}A^Ty $ 在实践中充满陷阱。唯有深入理解矩阵结构、条件数影响及各类分解算法特性，才能构建出既准确又高效的拟合系统。下一章将进一步探讨如何评估这些解的质量，并引入正则化等高级技巧应对过拟合问题。
4. 拟合结果评估与模型优化策略
在构建多项式模型并完成最小二乘求解后，如何科学地评估拟合质量、识别潜在问题，并进一步提升模型的泛化能力，是工程实践中至关重要的环节。一个看似“完美”贴合训练数据的高阶多项式，可能实际上已陷入过拟合陷阱，失去对未知数据的预测能力。因此，必须引入系统化的评估指标和优化机制，从多个维度审视模型表现。本章将深入探讨拟合质量的量化方法，揭示过拟合现象的本质特征及其诊断手段，并重点介绍正则化技术在抑制模型复杂度、增强稳定性方面的关键作用。
4.1 拟合质量的量化评估指标
为了客观衡量多项式拟合的效果，仅凭肉眼观察曲线是否“贴近”原始数据点远远不够。我们需要借助统计学中的量化指标，从误差大小、解释能力以及泛化性能等角度全面评价模型。常用的评估指标包括决定系数 $ R^2 $、均方根误差（RMSE），以及基于数据划分的交叉验证方法。这些工具不仅帮助我们判断当前模型的好坏，还为不同阶数或结构的模型之间提供可比较的标准。
4.1.1 决定系数 $ R^2 $ 的意义与局限性
决定系数 $ R^2 $ 是衡量模型对目标变量变异程度解释比例的核心指标，其定义如下：
R^2 = 1 - \frac{SS_{\text{res}}}{SS_{\text{tot}}}
其中：
- $ SS_{\text{res}} = \sum_{i=1}^{n}(y_i - \hat{y} i)^2 $：残差平方和（Residual Sum of Squares）
- $ SS {\text{tot}} = \sum_{i=1}^{n}(y_i - \bar{y})^2 $：总平方和（Total Sum of Squares）
- $ y_i $：第 $ i $ 个观测值
- $ \hat{y}_i $：对应的拟合值
- $ \bar{y} $：所有观测值的均值
$ R^2 $ 的取值范围通常在 [0, 1] 之间，越接近 1 表示模型解释了越多的数据变异性，拟合效果越好。当模型完全拟合所有数据点时，$ SS_{\text{res}} = 0 $，此时 $ R^2 = 1 $；若模型预测恒等于均值，则 $ R^2 = 0 $。
然而，$ R^2 $ 存在一个显著缺陷：它会随着模型复杂度（如多项式阶数）的增加而单调上升，即使新增项并无实际意义。例如，在极端情况下使用 $ n-1 $ 阶多项式拟合 $ n $ 个点，必然得到 $ R^2 = 1 $，但这并不代表模型具有良好的外推能力。
下面通过 Python 实现 $ R^2 $ 的计算过程：
import numpy as np
def compute_r_squared(y_true, y_pred):
ss_res = np.sum((y_true - y_pred) ** 2)
ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
r2 = 1 - (ss_res / ss_tot)
return r2
# 示例数据
y_true = np.array([1.2, 2.3, 3.5, 4.1, 5.0])
y_pred = np.array([1.1, 2.4, 3.4, 4.2, 4.9])
r2_score = compute_r_squared(y_true, y_pred)
print(f"R² Score: {r2_score:.4f}")
代码逻辑逐行解析：
1. ss_res = np.sum((y_true - y_pred) ** 2) ：计算残差平方和，反映模型预测偏差。
2. ss_tot = np.sum((y_true - np.mean(y_true)) ** 2) ：计算总平方和，代表数据本身的波动总量。
3. r2 = 1 - (ss_res / ss_tot) ：根据公式计算 $ R^2 $，体现模型相对于基准（均值）的改进程度。
4. 返回浮点数值，便于后续比较分析。
该函数可用于任意回归任务的结果评估，但在高维或多参数场景下应结合调整后的 $ R^2 $ 使用以避免误导。
参数说明表：
| 参数 | 类型 | 含义 |
|---|---|---|
y_true | np.ndarray | 真实观测值向量 |
y_pred | np.ndarray | 模型预测输出向量 |
| 返回值 | float | 决定系数 $ R^2 $，越大表示解释力越强 |
⚠️ 注意事项：当模型欠拟合严重时，$ R^2 $ 可能为负值，表明模型比简单均值更差。
4.1.2 均方根误差 RMSE 的物理含义与计算方式
均方根误差（Root Mean Square Error, RMSE）是衡量预测值与真实值之间平均差异的重要指标，尤其适用于关注误差幅度的应用场景。其数学表达为：
\text{RMSE} = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}
相比 $ R^2 $，RMSE 具有明确的物理单位（与原数据一致），能直观反映误差的“典型规模”。例如，在温度拟合中，RMSE 为 0.5°C 表示平均偏差约为半度，便于工程师理解。
RMSE 对异常值敏感，因其平方操作放大了大误差的影响，这既是优点也是缺点——既能快速暴露模型在某些区域的失效，也可能被个别噪声点扭曲整体评估。
以下是 RMSE 的实现代码：
def compute_rmse(y_true, y_pred):
mse = np.mean((y_true - y_pred) ** 2)
rmse = np.sqrt(mse)
return rmse
rmse_value = compute_rmse(y_true, y_pred)
print(f"RMSE: {rmse_value:.4f}")
逐行分析：
1. (y_true - y_pred) ** 2 ：逐点计算误差平方。
2. np.mean(...) ：求平均值得到 MSE（均方误差）。
3. np.sqrt(mse) ：开方获得 RMSE，还原至原始量纲。
相较于 $ R^2 $，RMSE 更适合用于损失函数设计和模型调优过程中作为优化目标。
RMSE 与其他误差指标对比表：
| 指标 | 公式 | 特点 |
|---|---|---|
| MAE | $\frac{1}{n}\sum|y_i-\hat{y}_i|$ | 对异常值不敏感，解释性强 |
| MSE | $\frac{1}{n}\sum(y_i-\hat{y}_i)^2$ | 易于微分，常用于梯度下降 |
| RMSE | $\sqrt{\text{MSE}}$ | 单位一致，突出大误差影响 |
| MAPE | $\frac{1}{n}\sum\left | \frac{y_i-\hat{y}_i}{y_i}\right |
选择何种指标取决于具体应用需求：若强调稳健性，可用 MAE；若需惩罚大偏差，则优先 RMSE。
4.1.3 训练集与测试集划分下的交叉验证方法
仅在训练数据上评估模型会导致乐观偏误，无法反映真实泛化能力。为此，必须采用数据分割策略进行交叉验证（Cross Validation）。最基础的方法是留出法（Hold-out Method），即将数据随机分为训练集和测试集。
from sklearn.model_selection import train_test_split
X = np.linspace(0, 10, 100)
y = 2 * X + np.sin(X) + np.random.normal(0, 0.5, size=X.shape)
X_train, X_test, y_train, y_test = train_test_split(
X.reshape(-1, 1), y, test_size=0.3, random_state=42
)
# 构建范德蒙德矩阵（以3阶为例）
def vander_matrix(x, degree):
return np.vander(x.flatten(), degree + 1, increasing=True)
V_train = vander_matrix(X_train, 3)
V_test = vander_matrix(X_test, 3)
# 最小二乘求解
coeffs = np.linalg.solve(V_train.T @ V_train, V_train.T @ y_train)
y_train_pred = V_train @ coeffs
y_test_pred = V_test @ coeffs
# 分别计算训练集与测试集的 RMSE
train_rmse = compute_rmse(y_train, y_train_pred)
test_rmse = compute_rmse(y_test, y_test_pred)
print(f"Train RMSE: {train_rmse:.4f}, Test RMSE: {test_rmse:.4f}")
流程图（Mermaid 格式）展示数据划分与评估流程：
graph TD
A[原始数据集] --> B{数据预处理}
B --> C[划分训练集/测试集]
C --> D[构建训练设计矩阵]
D --> E[求解最小二乘系数]
E --> F[生成训练预测值]
F --> G[计算训练RMSE/R²]
C --> H[构建测试设计矩阵]
H --> I[应用相同模型预测]
I --> J[计算测试RMSE/R²]
J --> K[比较性能差异]
K --> L{是否存在过拟合?}
L -->|测试误差显著更大| M[降低模型复杂度或引入正则化]
L -->|接近| N[模型合理]
代码逻辑说明：
- train_test_split 将数据按 70%/30% 划分，确保独立验证。
- 使用相同的多项式基函数构造训练与测试的设计矩阵。
- 系数仅由训练集求得，测试集仅用于评估，模拟真实预测场景。
- 若测试 RMSE 明显高于训练 RMSE，则提示过拟合风险。
此方法虽简单有效，但对于小样本数据可能存在偶然性。进阶做法可采用 k 折交叉验证（k-Fold CV），提升评估稳定性。
4.2 过拟合现象的识别与诊断
尽管高阶多项式能够精确通过每一个训练点，但这种“完美拟合”往往是以牺牲平滑性和泛化能力为代价的。过拟合表现为模型过度适应训练数据中的噪声或局部波动，导致在新数据上的表现急剧下降。识别和诊断过拟合是模型优化的前提，需结合可视化工具与定量分析手段共同完成。
4.2.1 高阶多项式导致的震荡行为观察
当多项式阶数过高时，尤其是在非均匀采样或边界区域，会出现剧烈振荡现象，即所谓的“龙格现象”（Runge’s Phenomenon）。即使在内部插值点附近表现良好，边缘部分仍可能出现极大偏差。
考虑如下实验：用不同阶数多项式拟合函数 $ f(x) = \frac{1}{1+25x^2} $ 在区间 [-1, 1] 上的等距采样点。
import matplotlib.pyplot as plt
def runge_function(x):
return 1 / (1 + 25 * x**2)
x_plot = np.linspace(-1, 1, 200)
y_plot = runge_function(x_plot)
plt.figure(figsize=(10, 6))
plt.plot(x_plot, y_plot, 'k-', label="True Function", linewidth=2)
for deg in [5, 10, 15]:
x_sample = np.linspace(-1, 1, deg+1)
y_sample = runge_function(x_sample)
V = np.vander(x_sample, deg+1, increasing=True)
coeffs = np.linalg.solve(V.T @ V, V.T @ y_sample)
y_fit = np.polyval(coeffs[::-1], x_plot)
plt.plot(x_plot, y_fit, '--', label=f'Poly Degree {deg}')
plt.scatter(x_sample, y_sample, c='red', s=20, zorder=5)
plt.title("Runge's Phenomenon in High-Degree Polynomial Fitting")
plt.xlabel("x"); plt.ylabel("y")
plt.legend(); plt.grid(True)
plt.show()
图像分析：
- 低阶（如5阶）拟合较为平稳，虽不能完全捕捉峰谷，但趋势正确。
- 随着阶数增至10、15，曲线开始出现剧烈摆动，尤其在两端超出真实函数范围。
- 所有采样点均被准确穿过，但中间区域产生虚假波动。
这说明单纯追求训练误差最小化是危险的，必须引入约束机制控制模型复杂度。
4.2.2 残差分布图与 Q-Q 图的可视化判据
残差 $ e_i = y_i - \hat{y}_i $ 提供了关于模型未解释部分的信息。理想情况下，残差应呈现均值为零、方差恒定、独立同分布的特性。绘制残差图有助于发现系统性偏差。
residuals = y_train - y_train_pred
plt.subplot(1, 2, 1)
plt.scatter(y_train_pred, residuals)
plt.axhline(0, color='red', linestyle='--')
plt.xlabel("Fitted Values")
plt.ylabel("Residuals")
plt.title("Residual vs Fitted Plot")
plt.grid(True)
from scipy import stats
stats.probplot(residuals, dist="norm", plot=plt.subplot(1, 2, 2))
plt.title("Q-Q Plot of Residuals")
plt.show()
解读：
- 左图若残差呈明显趋势（如U形或漏斗状），说明模型遗漏重要变量或存在异方差。
- 右图 Q-Q 图中点偏离对角线，表明残差不服从正态分布，影响统计推断有效性。
上述图形构成诊断过拟合的基础工具链。
4.2.3 模型复杂度与误差趋势的关系曲线分析
绘制训练误差与测试误差随模型阶数变化的趋势图，是识别最优复杂度的关键手段。
degrees = range(1, 16)
train_errors, test_errors = [], []
for d in degrees:
V_train_d = vander_matrix(X_train, d)
V_test_d = vander_matrix(X_test, d)
# 正规方程求解
try:
coeffs_d = np.linalg.solve(V_train_d.T @ V_train_d, V_train_d.T @ y_train)
y_train_p = V_train_d @ coeffs_d
y_test_p = V_test_d @ coeffs_d
train_errors.append(compute_rmse(y_train, y_train_p))
test_errors.append(compute_rmse(y_test, y_test_p))
except np.linalg.LinAlgError:
train_errors.append(np.inf)
test_errors.append(np.inf)
plt.figure(figsize=(9, 5))
plt.plot(degrees, train_errors, 'b-o', label='Training RMSE')
plt.plot(degrees, test_errors, 'r-s', label='Testing RMSE')
plt.xlabel("Polynomial Degree")
plt.ylabel("RMSE")
plt.title("Bias-Variance Tradeoff Curve")
plt.legend(); plt.grid(True)
plt.xticks(degrees); plt.show()
表格：各阶数下误差对比（示例）
| Degree | Train RMSE | Test RMSE |
|---|---|---|
| 1 | 0.8912 | 0.9015 |
| 3 | 0.4521 | 0.4678 |
| 5 | 0.3210 | 0.3402 |
| 7 | 0.2103 | 0.2891 |
| 9 | 0.1021 | 0.3785 |
| 12 | 0.0032 | 0.6210 |
可见，测试误差先降后升，形成典型的 U 形曲线，最低点对应最佳模型阶数。
4.3 正则化技术在最小二乘中的扩展应用
面对病态设计矩阵和过拟合风险，传统最小二乘难以胜任。正则化通过在目标函数中加入惩罚项，限制参数幅度，从而提升模型稳定性和泛化能力。岭回归（L2）和 Lasso 回归（L1）是最具代表性的两类方法。
4.3.1 岭回归（L2正则化）的数学表达与作用机理
岭回归的目标函数为：
\min_{\mathbf{c}} | \mathbf{A}\mathbf{c} - \mathbf{y} |^2 + \lambda | \mathbf{c} |^2
其解为：
\mathbf{c}_{\text{ridge}} = (\mathbf{A}^T\mathbf{A} + \lambda \mathbf{I})^{-1} \mathbf{A}^T \mathbf{y}
其中 $ \lambda > 0 $ 为正则化参数，控制惩罚强度。加入单位矩阵使得 $ \mathbf{A}^T\mathbf{A} + \lambda \mathbf{I} $ 始终可逆，解决了秩亏或病态问题。
def ridge_regression(A, y, lam):
n_features = A.shape[1]
ATA = A.T @ A
regularized_ATA = ATA + lam * np.eye(n_features)
coeff = np.linalg.solve(regularized_ATA, A.T @ y)
return coeff
# 应用示例
lam = 1e-4
V_high = vander_matrix(X_train, 10)
c_ridge = ridge_regression(V_high, y_train, lam)
y_ridge_pred = vander_matrix(X_test, 10) @ c_ridge
rmse_ridge = compute_rmse(y_test, y_ridge_pred)
优势：
- 改善矩阵条件数，防止数值不稳定。
- 缩小系数幅值，减少过拟合。
- 所有变量保留，适合相关性强的情形。
4.3.2 Lasso回归（L1正则化）对稀疏解的促进效果
Lasso 使用 L1 惩罚：
\min_{\mathbf{c}} | \mathbf{A}\mathbf{c} - \mathbf{y} |^2 + \lambda | \mathbf{c} |_1
其特点是能将部分系数压缩至零，实现自动特征选择。适用于高维稀疏场景。
from sklearn.linear_model import Lasso
lasso = Lasso(alpha=0.01, max_iter=10000)
lasso.fit(V_high, y_train)
y_lasso_pred = lasso.predict(vander_matrix(X_test, 10))
对比岭回归与Lasso特点：
| 特性 | 岭回归（L2） | Lasso（L1） |
|---|---|---|
| 解的性质 | 连续收缩 | 可产生稀疏解 |
| 变量选择 | 不进行 | 自动剔除无关变量 |
| 数值求解 | 解析解存在 | 需迭代算法（如坐标下降） |
| 适用场景 | 多重共线性 | 特征筛选 |
4.3.3 正则化参数 λ 的选择方法：岭迹法与交叉验证
选择合适的 $ \lambda $ 至关重要。常用方法包括：
- 岭迹法（Ridge Trace） ：绘制不同 $ \lambda $ 下系数变化路径，观察趋于稳定的拐点。
- 交叉验证 ：使用
sklearn.model_selection.GridSearchCV或RidgeCV自动搜索最优值。
from sklearn.linear_model import RidgeCV
alphas = np.logspace(-6, -1, 30)
ridge_cv = RidgeCV(alphas=alphas, cv=5)
ridge_cv.fit(V_high, y_train)
best_lambda = ridge_cv.alpha_
print(f"Optimal λ: {best_lambda}")
最终选定的 $ \lambda $ 应使验证误差最小且系数稳定。
综上所述，合理的评估体系与正则化策略相结合，构成了现代曲线拟合中不可或缺的技术支柱。
5. 工程实践中的曲线拟合应用与系统封装
5.1 PolynomialFit类的设计与接口抽象
在实际工程项目中，将最小二乘多项式拟合封装为一个可复用、高内聚的C++类是提升代码可维护性和扩展性的关键。 PolynomialFit 类的设计应遵循面向对象原则，兼顾数值稳定性、异常处理和接口简洁性。
5.1.1 类成员变量规划：设计矩阵、系数向量、阶数控制
class PolynomialFit {
private:
int degree; // 多项式阶数
bool fitted; // 拟合状态标志
Eigen::VectorXd coefficients; // 拟合系数向量 [c0, c1, ..., cn]
Eigen::MatrixXd A; // 设计矩阵（范德蒙德矩阵）
std::vector<double> x_data; // 原始输入数据x
std::vector<double> y_data; // 原始输出数据y
double r_squared; // 决定系数R²缓存
double rmse; // 均方根误差缓存
};
其中 degree 控制模型复杂度； coefficients 存储求解得到的多项式系数； A 是通过输入 x_data 构造的范德蒙德矩阵，形式如下：
A =
\begin{bmatrix}
1 & x_0 & x_0^2 & \cdots & x_0^n \
1 & x_1 & x_1^2 & \cdots & x_1^n \
\vdots & \vdots & \vdots & \ddots & \vdots \
1 & x_{m-1} & x_{m-1}^2 & \cdots & x_{m-1}^n \
\end{bmatrix}
该结构便于后续调用 QR 分解或正则化方法进行稳健求解。
5.1.2 核心接口定义：fit()、predict()、getR2()等功能模块
类对外暴露的标准接口包括：
| 接口函数 | 功能说明 |
|---|---|
bool fit(const std::vector<double>& x, const std::vector<double>& y, int n) | 执行拟合并返回是否成功 |
double predict(double x) | 给定输入x，返回多项式预测值 |
std::vector<double> predict(const std::vector<double>& x) | 批量预测 |
double getR2() const | 获取决定系数 R² |
double getRMSE() const | 获取均方根误差 |
std::vector<double> getCoefficients() const | 返回拟合系数 |
示例 fit() 函数实现逻辑：
bool PolynomialFit::fit(const std::vector<double>& x, const std::vector<double>& y, int n) {
if (x.size() != y.size() || x.empty()) return false;
if (n < 0 || n >= x.size()) return false;
degree = n;
x_data = x; y_data = y;
int m = x.size();
// 构造范德蒙德矩阵 A ∈ R^(m × (n+1))
A.resize(m, n + 1);
for (int i = 0; i < m; ++i) {
A(i, 0) = 1.0;
for (int j = 1; j <= n; ++j) {
A(i, j) = A(i, j-1) * x[i];
}
}
// 转换观测向量 Y
Eigen::VectorXd Y(m);
for (int i = 0; i < m; ++i) Y(i) = y[i];
// 使用QR分解求解最小二乘问题
try {
coefficients = A.colPivHouseholderQr().solve(Y);
fitted = true;
// 计算评估指标
computeMetrics(Y);
} catch (...) {
fitted = false;
}
return fitted;
}
使用 Eigen 的 colPivHouseholderQr() 可有效应对秩亏或近奇异情况，提高鲁棒性。
5.1.3 异常处理机制：秩亏、奇异矩阵等错误状态捕获
在构造高阶多项式时，若样本点分布集中或存在重复值，可能导致设计矩阵列相关甚至秩亏。为此需加入以下检测机制：
- 在
fit()中检查A.rank()是否小于min(m, n+1) - 判断条件数
A.conditionNumber()是否超过阈值（如 1e8） - 对
solve()抛出的异常进行捕获并设置fitted = false
可通过如下方式增强诊断信息：
if (A.colPivHouseholderQr().rank() < A.cols()) {
std::cerr << "Warning: Design matrix is rank-deficient. Consider reducing polynomial degree." << std::endl;
return false;
}
此类机制确保系统在非理想数据下仍能安全运行，并提供调试线索。
5.2 拟合结果的可视化分析
5.2.1 使用Matplotlib-CPP或GNUPLOT接口绘制原始数据与拟合曲线
借助 matplotlib-cpp 封装库，可在C++中实现高质量绘图。首先准备数据：
#include "matplotlibcpp.h"
namespace plt = matplotlibcpp;
// 生成密集采样点用于绘制平滑曲线
std::vector<double> x_plot(200);
std::vector<double> y_plot(200);
double dx = (xmax - xmin)/199;
for(int i=0; i<200; ++i) x_plot[i] = xmin + i*dx;
auto y_pred = poly.predict(x_plot);
plt::figure();
plt::plot(x_data, y_data, "bo", {{"label", "Noisy Data"}});
plt::plot(x_plot, y_pred, "r-", {{"linewidth", "2"}, {"label", "Fitted Curve"}});
plt::xlabel("x"); plt::ylabel("y");
plt::legend(); plt::grid(true);
plt::title("Polynomial Fit Result (degree = " + std::to_string(degree) + ")");
plt::show();
5.2.2 多阶拟合对比图与残差图的联合呈现
采用子图布局展示不同阶数下的拟合效果及残差分布：
plt::subplot(2, 1, 1);
plt::plot(x_data, y_data, "ko");
plt::plot(x_plot, y_pred, "-", {{"label", "deg="+std::to_string(d)}});
plt::legend(); plt::grid(true);
plt::subplot(2, 1, 2);
auto residuals = compute_residuals();
plt::stem(x_data, residuals, {{"markerfmt", "r."}});
plt::axhline(0, {{"color", "gray"}, {"linestyle", "--"}});
plt::ylabel("Residual");
plt::grid(true);
残差图有助于识别系统偏差或异方差性。
5.2.3 动态更新图形界面支持实时调试
结合 ImGui 或 Qt 可构建交互式拟合工具，允许用户动态调整阶数并即时刷新图像，适用于传感器校准等场景。例如：
graph TD
A[用户拖动滑块改变阶数] --> B{触发fit()}
B --> C[重新计算系数]
C --> D[更新预测曲线]
D --> E[刷新图表显示]
E --> F[显示R²/RMSE变化趋势]
此流程支持快速原型验证与参数调优。
5.3 最小二乘拟合在实际场景中的典型应用
5.3.1 工程传感器数据的趋势提取与校准
在温度、压力传感器中，原始读数常包含漂移趋势。设某热电偶采集数据含缓慢上升基线：
| 时间(s) | 原始读数(V) | 真实值(V) |
|---|---|---|
| 0 | 0.102 | 0.100 |
| 10 | 0.208 | 0.200 |
| 20 | 0.325 | 0.300 |
| 30 | 0.440 | 0.400 |
| 40 | 0.570 | 0.500 |
| 50 | 0.695 | 0.600 |
| 60 | 0.830 | 0.700 |
| 70 | 0.970 | 0.800 |
| 80 | 1.120 | 0.900 |
| 90 | 1.270 | 1.000 |
| 100 | 1.430 | 1.100 |
使用 PolynomialFit 对“原始读数 - 真实值”误差序列建模，拟合二阶多项式后反向补偿，显著降低系统误差。
5.3.2 通信信号中噪声背景下的基线恢复
在EEG或雷达信号处理中，低频基线漂移会掩盖有效特征。对每段窗口数据局部拟合一次或二次多项式，再从原信号中减去拟合趋势项，即可实现基线校正。
5.3.3 时间序列预测中的局部多项式拟合策略
对于非平稳时间序列（如股价波动），全局拟合易失效。采用滑动窗口 + 局部最小二乘拟合（如LOESS思想），仅对未来短时段外推，结合置信区间评估风险，可作为高频交易辅助模型。
简介：最小二乘法是一种广泛应用于数据建模与曲线拟合的数学方法，通过使数据点到拟合曲线的误差平方和最小，实现最优逼近。本文介绍在C++环境下实现多项式曲线拟合的完整流程，涵盖数据准备、模型定义、矩阵构建、线性系统求解及结果评估等关键步骤。结合PolynomialFit.h与PolynomialFit.cpp文件，利用Eigen等数值库高效求解，支持可视化分析与拟合优化，适用于线性回归、信号处理和工程计算等实际场景。
扫一扫
882
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
