Python实战：用最小二乘法搞定曲线拟合（附Eigen库对比代码）

Python实战：用最小二乘法搞定曲线拟合（附Eigen库对比代码）

最近在整理一个传感器数据校准的项目，发现无论数据看起来多“干净”，直接用原始读数做预测总是不太准。和团队里的算法工程师聊了聊，他一句话点醒了我：“你这就是个典型的曲线拟合问题，试试最小二乘法，从原理到代码，半小时就能上手。” 我一开始还将信将疑，毕竟“最小二乘”听起来像是教科书里那些复杂的数学推导。但真正动手用Python实现了一遍，并且对比了C++里常用的Eigen库写法后，才发现它的核心思想异常直观，而且对于处理实验数据、校准模型、甚至机器学习中的线性回归，都是绕不开的基石。这篇文章，我就从一个实践者的角度，带你重新认识最小二乘法，不止于公式，更聚焦于如何用代码把它用起来，解决真实世界中的数据拟合难题。

无论是刚开始接触数据科学的学生，还是需要快速验证算法原型的工程师，理解并实现最小二乘法都能让你在面对杂乱数据时，多一份笃定。我们会从最基础的原理切入，用Python的NumPy一步步构建求解器，然后引入更健壮的加权与迭代加权版本以应对异常值，最后，我会展示如何用C++的Eigen库实现同样的功能，并对比两种语言在性能和易用性上的差异。更重要的是，我们会深入探讨一个实际应用中无法回避的问题：如何选择合适的多项式阶数？盲目追求高阶拟合带来的“过拟合”陷阱，远比想象中更常见。

1. 最小二乘法的核心：从几何直觉到矩阵求解

很多人第一次接触最小二乘法，都会被那一串求和符号和偏导公式吓退。其实，我们可以换个更直观的方式来理解。想象你在二维平面上有一系列散点，你想找一条直线，让所有这些点到这条直线的垂直距离的平方和最小。这个“距离平方和最小”，就是“最小二乘”中“二乘”（Least Squares）的含义。它比直接用距离之和更好，因为平方操作能放大较大误差的影响，使得求解过程更稳定，并且数学上可导，便于计算。

1.1 代数形式与矩阵形式的抉择

早期的教材多从代数形式推导。对于拟合直线y = a*x + b，我们需要求解参数a和b。通过建立误差平方和函数，并分别对a和b求偏导令其为零，可以得到一个二元一次方程组。这个方法清晰，但一旦模型变得复杂，比如拟合一个三次多项式y = a*x³ + b*x² + c*x + d，代数推导就会变得异常繁琐。

注意：代数推导在低维（1-2个参数）时有助于理解原理，但绝不适用于实际编程。它的计算复杂度高，且容易在推导中出错。

现代计算几乎无一例外地采用矩阵形式。它将拟合问题优雅地转化为一个线性方程组求解问题。对于有n个数据点，拟合m阶多项式的情况，我们可以构造如下矩阵方程：

A * x = B

这里：

• A被称为设计矩阵。它的每一行对应一个数据点，每一列对应多项式的一项。例如，对于二次拟合y = p0 + p1*x + p2*x²，第i行就是[1, x_i, x_i²]。
• x是我们要求解的参数向量，即[p0, p1, p2]^T。
• B是观测值向量，即所有y_i组成的列向量。

我们的目标是找到参数向量x，使得A*x尽可能接近B。根据线性代数知识，当方程数（数据点）多于未知数（参数）时，这是一个“超定方程组”，通常没有精确解。最小二乘给出的就是那个最优的近似解，其标准解（在(A^T A)可逆的情况下）为：

x = (A^T * A)^-1 * A^T * B

这个公式是核心中的核心。A^T * A是一个方阵，求逆后即得到参数。下面我们用Python的NumPy库来直观感受一下实现是多么简洁。

1.2 Python实现：用NumPy构建最小二乘求解器

NumPy提供了强大的线性代数工具，让我们可以几乎“直译”上面的矩阵公式。假设我们有一组数据点x_data和y_data，想要拟合一个三次多项式。

import numpy as np import matplotlib.pyplot as plt def ordinary_least_squares(x_data, y_data, degree): """ 普通最小二乘法 (OLS) 多项式拟合 参数: x_data: 自变量数据，一维数组 y_data: 因变量数据，一维数组 degree: 多项式阶数 返回: coeff: 多项式系数，从低次到高次 [c0, c1, ..., c_degree] """ # 1. 构建设计矩阵 A # 使用范德蒙德矩阵 (Vandermonde matrix) 是最简单的方式 A = np.vander(x_data, degree + 1, increasing=True) # 列顺序为 x^0, x^1, ..., x^degree # 2. 直接利用公式求解 x = (A^T A)^-1 A^T y # np.linalg.lstsq 是更数值稳定的专业函数，但这里为了演示公式，我们手动计算 ATA = np.dot(A.T, A) ATb = np.dot(A.T, y_data) # 使用求解线性方程组的方式，比直接求逆更稳定 coeff = np.linalg.solve(ATA, ATb) return coeff # 生成带噪声的示例数据 np.random.seed(42) x_raw = np.linspace(0, 10, 30) y_true = 2.5 * x_raw + 1.8 # 真实关系：直线 y_noise = y_true + np.random.randn(len(x_raw)) * 4 # 加入高斯噪声 # 进行线性拟合 (degree=1) coeff_linear = ordinary_least_squares(x_raw, y_noise, degree=1) print(f"拟合的直线参数 (截距, 斜率): {coeff_linear}") # 进行三次多项式拟合 coeff_cubic = ordinary_least_squares(x_raw, y_noise, degree=3) print(f"拟合的三次多项式参数: {coeff_cubic}")

运行这段代码，你会立刻得到拟合的参数。np.vander函数巧妙地为我们构建了设计矩阵。为了验证拟合效果，我们可以计算预测值并绘图：

# 生成平滑曲线用于绘制拟合结果 x_plot = np.linspace(x_raw.min(), x_raw.max(), 300) A_plot = np.vander(x_plot, 2, increasing=True) # 用于直线 y_linear_fit = np.dot(A_plot, coeff_linear) A_plot_cubic = np.vander(x_plot, 4, increasing=True) # 用于三次曲线 y_cubic_fit = np.dot(A_plot_cubic, coeff_cubic) plt.figure(figsize=(10, 6)) plt.scatter(x_raw, y_noise, alpha=0.6, label='带噪声的数据点') plt.plot(x_plot, y_linear_fit, 'r-', linewidth=2, label=f'线性拟合: y={coeff_linear[1]:.2f}x+{coeff_linear[0]:.2f}') plt.plot(x_plot, y_cubic_fit, 'g--', linewidth=2, label='三次多项式拟合') plt.xlabel('X') plt.ylabel('Y') plt.legend() plt.grid(True, linestyle='--', alpha=0.5) plt.title('普通最小二乘法 (OLS) 拟合对比') plt.show()

通过图像，你可以清晰地看到直线和曲线对同一组数据的拟合差异。这就是最小二乘法的基本应用。然而，ordinary_least_squares函数中的手动求解(A^T A)^-1 A^T b在数值计算上可能存在隐患。

提示：当A^T A的条件数很大时（即矩阵接近奇异），直接求逆或求解会带来巨大的数值误差。在实际项目中，强烈建议使用np.linalg.lstsq(A, y_data, rcond=None)这个经过高度优化的函数，它使用更稳定的数值算法（如SVD）来求解最小二乘问题。

2. 应对现实挑战：加权与迭代重加权最小二乘法

普通最小二乘法有一个很强的假设：所有数据点的误差是独立同分布的。这意味着每个点对最终拟合结果的“话语权”是相同的。但在现实中，这个假设经常被打破。例如，有些数据来自高精度传感器，有些则来自噪声很大的旧设备；或者在拟合时，我们希望近期数据比历史数据权重更高。

2.1 加权最小二乘法 (WLS)：赋予数据不同权重

加权最小二乘法的思想非常直接：为每个数据点分配一个权重。越可靠的数据点，权重越大，它在拟合过程中的影响力也就越大。这在金融时间序列分析、传感器融合等领域非常有用。

其矩阵公式只是在OLS的基础上引入了一个对角权重矩阵 W，其中对角线元素W[i, i]就是第i个数据点的权重。优化目标变为最小化加权误差平方和。解的形式变为：

x = (A^T * W^T * W * A)^-1 * A^T * W^T * W * B

由于W是对角阵，W^T * W就是一个对角线元素为w_i²的矩阵。实现起来只需要在构建A和B时，将每一行都乘以对应的权重w_i即可。

def weighted_least_squares(x_data, y_data, weights, degree): """ 加权最小二乘法 (WLS) 参数: weights: 权重向量，长度与数据点相同。权重越大，该点影响力越大。 """ if len(weights) != len(x_data): raise ValueError("权重向量长度必须与数据点长度一致") # 1. 构建设计矩阵 A A = np.vander(x_data, degree + 1, increasing=True) # 2. 构建权重矩阵 W (这里用对角阵的平方根 W_sqrt，方便计算) W_sqrt = np.diag(np.sqrt(weights)) # 因为公式中是 W^T W，我们直接用 sqrt(W) 乘到A和B上 # 3. 对A和B进行加权 A_weighted = np.dot(W_sqrt, A) B_weighted = np.dot(W_sqrt, y_data) # 4. 求解加权后的最小二乘问题 # 使用 np.linalg.lstsq 更稳定 coeff, residuals, rank, s = np.linalg.lstsq(A_weighted, B_weighted, rcond=None) return coeff # 示例：假设中间部分的数据更可靠 weights = np.ones_like(x_raw) weights[10:20] = 5.0 # 将第10到20个数据点的权重设为5倍 coeff_wls = weighted_least_squares(x_raw, y_noise, weights, degree=1) print(f"加权最小二乘拟合参数: {coeff_wls}")

通过调整weights数组，你可以控制拟合曲线更“偏向”哪些区域的数据。

2.2 迭代重加权最小二乘法 (IRLS)：智能降噪的利器

WLS要求我们事先知道权重，但这通常很难。迭代重加权最小二乘法 (IRLS) 是一个更强大的工具，它能根据拟合的残差（误差）动态调整权重。其核心思想是：

1. 先用OLS进行一次拟合，得到初始参数和残差。
2. 根据残差大小计算新的权重：残差大的点，赋予较小的权重（因为它可能是离群点）。
3. 用新的权重进行WLS拟合，得到新的参数和残差。
4. 重复步骤2-3，直到参数收敛或达到迭代次数。

常用的权重函数是残差绝对值的p-2次方（p是范数阶数，常用1或2）。当p=2时，权重恒为1，退化为OLS；当p=1时，权重为1 / |残差|，能有效抑制离群点的影响，因此IRLS常用于稳健回归。

def iterative_reweighted_least_squares(x_data, y_data, degree, p=1.2, max_iter=20, tol=1e-6): """ 迭代重加权最小二乘法 (IRLS) 参数: p: 范数参数，通常 1 <= p < 2。p=2等价于OLS，p接近1对离群点更鲁棒。 max_iter: 最大迭代次数 tol: 参数收敛容差 """ n_samples = len(x_data) # 初始权重全为1 weights = np.ones(n_samples) coeff = ordinary_least_squares(x_data, y_data, degree) # 初始解 A = np.vander(x_data, degree + 1, increasing=True) for i in range(max_iter): coeff_old = coeff.copy() # 计算当前残差 residual = y_data - np.dot(A, coeff) # 防止除零，加一个小常数 epsilon = 1e-8 # 根据残差更新权重，这里使用Huber-like权重函数 # 当p=1时，权重 = 1 / max(|residual|, epsilon) abs_res = np.abs(residual) weights = 1.0 / np.maximum(abs_res, epsilon) ** (2 - p) # 常见的权重更新公式 # 进行本轮加权最小二乘拟合 W_sqrt = np.diag(np.sqrt(weights)) A_weighted = np.dot(W_sqrt, A) B_weighted = np.dot(W_sqrt, y_data) coeff, _, _, _ = np.linalg.lstsq(A_weighted, B_weighted, rcond=None) # 检查收敛 if np.linalg.norm(coeff - coeff_old) < tol: print(f'IRLS 在第 {i+1} 次迭代后收敛。') break else: print(f'IRLS 在 {max_iter} 次迭代后未完全收敛。') return coeff # 人为添加一个严重的离群点 x_with_outlier = np.append(x_raw, 8.0) y_with_outlier = np.append(y_noise, 60.0) # 在x=8, y=60处加入一个异常高的点 # 对比OLS和IRLS coeff_ols_outlier = ordinary_least_squares(x_with_outlier, y_with_outlier, degree=1) coeff_irls_outlier = iterative_reweighted_least_squares(x_with_outlier, y_with_outlier, degree=1, p=1.2) print(f"存在离群点时，OLS拟合结果: {coeff_ols_outlier}") print(f"存在离群点时，IRLS拟合结果: {coeff_irls_outlier}")

运行这段代码并绘图，你会直观地看到，OLS拟合的直线会被那个离群点明显“拉偏”，而IRLS拟合的直线则几乎不受影响，更贴近大多数数据点形成的趋势。这就是稳健回归的价值所在。

3. C++ Eigen库实现：性能与工程化的考量

Python在原型验证和数据分析上得天独厚，但在对性能有苛刻要求的嵌入式系统、高频交易或大型数值计算中，C++仍是首选。Eigen是一个C++模板库，提供了媲美Matlab的线性代数API，语法简洁，运行效率极高。下面我们看看如何用Eigen实现同样的最小二乘求解。

首先，你需要确保开发环境已包含Eigen。它通常只是一个头文件库，下载后包含路径即可。

3.1 普通最小二乘法的Eigen实现

Eigen的语法非常直观，几乎可以看作是NumPy的C++版本。

#include <iostream> #include <Eigen/Dense> // 核心矩阵运算模块 using namespace Eigen; VectorXd ordinaryLeastSquaresEigen(const VectorXd& x_data, const VectorXd& y_data, int degree) { int n = x_data.size(); // 1. 构建设计矩阵 A (n x (degree+1)) MatrixXd A(n, degree + 1); for (int i = 0; i < n; ++i) { for (int j = 0; j <= degree; ++j) { A(i, j) = std::pow(x_data(i), j); // 第j列是 x^j } } // 2. 使用Eigen的BDCSVD分解求解最小二乘问题，这是最稳健的方法之一 // BDCSVD是分治SVD算法，适合大型矩阵 VectorXd coeff = A.bdcSvd(ComputeThinU | ComputeThinV).solve(y_data); return coeff; } int main() { // 示例：生成数据 int n = 30; VectorXd x_data = VectorXd::LinSpaced(n, 0, 10); VectorXd y_true = 2.5 * x_data.array() + 1.8; VectorXd noise = VectorXd::Random(n) * 4; // Eigen的Random生成[-1,1]均匀分布，这里模拟噪声 VectorXd y_data = y_true + noise; int degree = 1; VectorXd coeff = ordinaryLeastSquaresEigen(x_data, y_data, degree); std::cout << "Eigen OLS 拟合参数 (截距， 斜率):\n" << coeff.transpose() << std::endl; // 计算R²分数 VectorXd y_pred = A * coeff; // 这里A需要重新计算或传递 double sse = (y_data - y_pred).squaredNorm(); // 误差平方和 double sst = (y_data.array() - y_data.mean()).square().sum(); // 总平方和 double r_squared = 1.0 - sse / sst; std::cout << "R² 分数: " << r_squared << std::endl; return 0; }

Eigen的BDCSVD分解求解器 (.bdcSvd().solve()) 是求解最小二乘问题的推荐方法，它在数值稳定性上优于直接计算(A^T A)^-1 A^T B。对于加权最小二乘，思路与Python一致：构建加权后的矩阵A_weighted和向量B_weighted，然后调用同一个求解器。

3.2 Python与C++ Eigen的性能与生态对比

将两种实现方式放在一起对比，能帮助我们做出更合适的技术选型。

我的经验是，在项目初期或进行探索性数据分析时，Python是无可争议的首选。你可以用几行代码快速验证想法，并用Matplotlib看到即时结果。而当算法定型，需要集成到对延迟和资源有严格限制的生产系统（比如自动驾驶的感知模块或工业控制器的实时校准）时，用C++和Eigen进行重写和优化就变得必要。很多时候，团队会采用“Python研发，C++部署”的混合模式。

4. 拟合的陷阱：过拟合与模型选择实战

掌握了工具，下一个关键问题是：我的多项式到底该选几阶？用最小二乘法拟合一个高阶多项式，总能让误差平方和变得更小，甚至对训练数据达到“完美”拟合。但这通常意味着灾难——过拟合。模型不仅学到了数据背后的规律，也记住了数据中的噪声，导致在新数据上表现极差。

4.1 可视化过拟合：从欠拟合到完美“插值”

让我们用一个例子来感受。假设真实的数据生成模型是一个带有噪声的二次曲线y = 1 + 2x - 0.5x² + noise。我们用不同阶数的多项式去拟合它。

def generate_sample_data(): np.random.seed(123) x = np.linspace(0, 5, 25) # 25个训练点 y_true = 1 + 2*x - 0.5*x**2 y_noisy = y_true + np.random.randn(len(x)) * 1.5 # 加入噪声 return x, y_noisy, y_true def fit_and_plot(x_train, y_train, max_degree=8): """拟合不同阶数多项式并绘图""" x_plot = np.linspace(x_train.min(), x_train.max(), 300) plt.figure(figsize=(14, 10)) train_errors = [] # 生成一个平滑区间外的测试点，用于观察外推能力 x_test_extrap = np.array([-1, 6]) # 训练数据区间是[0,5]，这里测试-1和6 y_test_extrap_true = 1 + 2*x_test_extrap - 0.5*x_test_extrap**2 for idx, degree in enumerate(range(1, max_degree+1)): coeff = np.polyfit(x_train, y_train, degree) # NumPy内置的多项式拟合函数 p = np.poly1d(coeff) # 生成多项式函数 y_plot = p(x_plot) y_pred_train = p(x_train) # 计算训练集上的均方误差 (MSE) mse_train = np.mean((y_pred_train - y_train)**2) train_errors.append(mse_train) # 计算外推点的误差 y_pred_extrap = p(x_test_extrap) mse_extrap = np.mean((y_pred_extrap - y_test_extrap_true)**2) plt.subplot(3, 3, idx+1) plt.scatter(x_train, y_train, s=30, alpha=0.7, label='训练数据') plt.plot(x_plot, y_plot, 'r-', linewidth=2, label=f'Deg {degree}') plt.scatter(x_test_extrap, y_test_extrap_true, c='green', s=100, marker='*', label='外推真值') plt.scatter(x_test_extrap, y_pred_extrap, c='darkred', s=100, marker='X', label='外推预测') plt.title(f'Degree {degree}\nTrain MSE: {mse_train:.3f}\nExtrap MSE: {mse_extrap:.3f}') plt.legend(loc='best', fontsize='x-small') plt.grid(True, linestyle='--', alpha=0.3) plt.ylim(y_train.min()-5, y_train.max()+5) plt.tight_layout() plt.show() # 绘制训练误差随阶数变化的曲线 plt.figure(figsize=(8, 4)) plt.plot(range(1, max_degree+1), train_errors, 'bo-', linewidth=2, markersize=8) plt.xlabel('多项式阶数') plt.ylabel('训练集均方误差 (MSE)') plt.title('训练误差 vs. 模型复杂度 (过拟合的典型信号)') plt.grid(True) plt.axvline(x=2, color='r', linestyle='--', alpha=0.5, label='真实阶数 (2)') plt.legend() plt.show() x_train, y_train, y_true = generate_sample_data() fit_and_plot(x_train, y_train, max_degree=8)

运行这段代码，你会看到一系列子图。在阶数较低（1阶，直线）时，模型过于简单，无法捕捉数据的弯曲趋势，这是欠拟合。在阶数为2或3时，曲线平滑地穿过了数据点聚集的区域，与真实二次曲线形状吻合。从第4阶开始，曲线开始出现不必要的波动，试图穿过每一个数据点，包括噪声点。到了第7、8阶，曲线在训练数据区间内剧烈震荡，虽然训练误差（MSE）已经非常小，但在区间外的预测点（绿色星号和红色X号）上，预测值与真实值相差甚远，外推误差激增，这就是典型的过拟合。

4.2 如何选择最佳阶数？交叉验证与信息准则

避免过拟合不能只靠“看感觉”。有两个系统性的方法：

1. 交叉验证 (Cross-Validation):将数据分成训练集和验证集。只用训练集拟合不同阶数的模型，然后在未见过的验证集上评估误差。选择在验证集上误差最小的模型。更稳健的方法是使用K折交叉验证。

from sklearn.model_selection import KFold from sklearn.metrics import mean_squared_error import warnings warnings.filterwarnings('ignore') def select_degree_by_cv(x_data, y_data, max_degree=10, n_splits=5): """使用5折交叉验证选择最佳多项式阶数""" kf = KFold(n_splits=n_splits, shuffle=True, random_state=42) degree_candidates = range(1, max_degree+1) cv_scores = {deg: [] for deg in degree_candidates} for train_idx, val_idx in kf.split(x_data): x_train_fold, x_val_fold = x_data[train_idx], x_data[val_idx] y_train_fold, y_val_fold = y_data[train_idx], y_data[val_idx] for deg in degree_candidates: coeff = np.polyfit(x_train_fold, y_train_fold, deg) p = np.poly1d(coeff) y_val_pred = p(x_val_fold) mse = mean_squared_error(y_val_fold, y_val_pred) cv_scores[deg].append(mse) # 计算每个阶数的平均MSE avg_scores = {deg: np.mean(scores) for deg, scores in cv_scores.items()} best_degree = min(avg_scores, key=avg_scores.get) return best_degree, avg_scores best_deg, scores = select_degree_by_cv(x_train, y_train, max_degree=8) print(f"交叉验证建议的最佳阶数: {best_deg}") print("各阶数对应的平均验证MSE:") for deg, score in scores.items(): print(f" 阶数 {deg}: {score:.4f}")

2. 信息准则 (Information Criterion):如赤池信息量准则 (AIC)或贝叶斯信息准则 (BIC)。它们在衡量模型拟合优度的同时，加入了对于参数数量的惩罚，从而倾向于选择更简洁的模型。statsmodels等库在线性回归结果中会直接给出AIC/BIC值。

import statsmodels.api as sm def calculate_aic_bic(x_data, y_data, max_degree=8): """计算不同阶数多项式回归的AIC和BIC""" results = [] for deg in range(1, max_degree+1): A = np.vander(x_data, deg + 1, increasing=True) model = sm.OLS(y_data, A).fit() results.append((deg, model.aic, model.bic)) return results ic_results = calculate_aic_bic(x_train, y_train, max_degree=8) print("阶数 | AIC | BIC") print("-" * 40) for deg, aic, bic in ic_results: print(f"{deg:4d} | {aic:12.2f} | {bic:12.2f}")

通常，AIC和BIC值最小的模型就是权衡了拟合能力和复杂度的“最佳”模型。在实际项目中，我通常会结合交叉验证和信息准则，再辅以业务逻辑的考量（例如，在物理模型中，阶数可能有理论依据），来做出最终决定。记住，没有免费的午餐，更复杂的模型并不总是更好的模型。最小二乘法给了你强大的拟合能力，但如何明智地使用这种能力，避免落入过拟合的陷阱，才是区分新手和专家的关键。