从最小二乘法到非线性拟合：数据建模中的核心算法与实践

1. 从“差不多”到“刚刚好”：为什么我们需要拟合算法

做数学建模，尤其是处理数据的时候，我们经常会遇到一个场景：手里有一堆实验或者观测得到的数据点，它们看起来似乎遵循某种规律，比如像一条直线，或者一个抛物线。我们的直觉是，如果能找到一个数学公式，能完美地穿过所有这些点，那不就万事大吉了吗？但现实往往很骨感。数据点通常不会那么“听话”，它们会因为测量误差、环境干扰等各种原因，散落在我们理想中的曲线周围。

这时候，如果你强行找一个复杂的公式（比如一个高次多项式）去穿过每一个点，这就是“插值”。插值出来的曲线会完美经过所有数据点，但往往会在点与点之间产生剧烈的、不合理的波动，这种现象被称为“龙格现象”。想象一下，你根据几个月的温度数据，硬是画出一条穿过每个日均温点的复杂曲线，用它来预测明天的温度，结果可能极其离谱。这就像是为一件量身定做的衣服，为了完全贴合身体每一个细微的起伏，把衣服做得满是褶皱，反而失去了整体合身的美感和预测性。

所以，我们需要的不是一件“完全贴身”但可能奇形怪状的衣服，而是一件“整体合身”的成衣。这就是拟合算法的核心思想：我们不要求曲线必须经过每一个数据点，而是要求它从整体趋势上，与所有数据点最为“接近”。这个“接近”的程度，需要一个量化的标准来衡量，最常用、最经典的就是最小二乘法。它的目标很简单：找到一条曲线，使得所有数据点到这条曲线的垂直距离的平方和最小。为什么是平方和？因为距离有正有负，直接相加会相互抵消，无法真实反映总体偏差；取绝对值在数学上不好处理（不可导）；而取平方既能消除正负影响，又保持了函数的良好性质（光滑可导），便于我们通过求导等数学工具找到那个最优解。

简单来说，拟合就是在承认数据有“噪声”的前提下，去寻找背后那个更简洁、更稳定、更能反映普遍规律的数学模型。从“清风数学建模”的语境来看，掌握拟合算法，意味着你拿到了从杂乱数据中提炼科学规律的钥匙，这是解决预测、趋势分析、参数估计等大量建模问题的基本功。

2. 最小二乘法的“灵魂”：目标函数与求解原理

我们先把问题说具体。假设我们有一组数据点(x_i, y_i), i=1,2,...,n。我们认为y和x之间存在某种函数关系y = f(x)，其中f(x)的形式是已知的，但里面有一些待定的参数。比如，我们认为它是直线关系，那么f(x) = kx + b，其中k和b就是待求参数；如果认为是二次关系，那么f(x) = ax^2 + bx + c，参数就是a, b, c。

拟合的目标是找到一组参数值，使得函数f(x)在所有这些x_i点上的计算值f(x_i)，与真实观测值y_i的总体差距最小。这个“总体差距”就是我们的目标函数，也叫损失函数。最小二乘法定义这个差距为残差平方和RSS(Residual Sum of Squares)：

RSS = Σ [y_i - f(x_i)]^2，其中求和Σ是从i=1到n。

这里的y_i - f(x_i)就是第i个点的残差，即垂直距离。我们的任务就是调整f(x)中的参数，让RSS这个值达到最小。这是一个典型的无约束优化问题。

那么，如何找到这组使RSS最小的参数呢？这里就体现出最小二乘法的巧妙之处。对于参数是线性组合的函数形式（即线性拟合），我们可以通过严格的数学推导得到解析解。

以最简单也最常用的一元线性拟合y = kx + b为例。我们的目标函数是：RSS(k, b) = Σ (y_i - k*x_i - b)^2

这是一个关于k和b的二元二次函数。根据微积分，要使RSS最小，它对k和b的偏导数应该同时为零。这就得到了所谓的正规方程组：

对b求偏导：∂RSS/∂b = -2 Σ (y_i - k*x_i - b) = 0->Σ y_i = k Σ x_i + n*b对k求偏导：∂RSS/∂k = -2 Σ [x_i*(y_i - k*x_i - b)] = 0->Σ (x_i*y_i) = k Σ (x_i^2) + b Σ x_i

这是一个关于k和b的二元一次方程组，直接求解即可得到：

k = [n*Σ(x_i*y_i) - Σx_i * Σy_i] / [n*Σ(x_i^2) - (Σx_i)^2]b = [Σy_i * Σ(x_i^2) - Σx_i * Σ(x_i*y_i)] / [n*Σ(x_i^2) - (Σx_i)^2]

这就是最小二乘线性拟合的解析解，也是许多教科书和基础教程里给出的公式。它的优势是计算确定、快速。对于多项式拟合（如y = a0 + a1*x + a2*x^2 + ...），虽然参数更多，但因其对参数仍是线性的，同样可以通过构建并求解一个更大的正规方程组来获得解析解。

注意：这里说的“线性”指的是参数是线性的，而不是x是线性的。y = a*sin(x) + b*exp(x)对参数a, b来说也是线性的，同样可以用上述方法求解。判断标准是：目标函数RSS对每个待求参数的偏导数方程是否仍然是这些参数的线性方程组。

3. 当公式失效时：非线性拟合与数值优化算法

现实世界的数据关系往往没那么“规矩”。很多模型的参数是非线性的。例如：

• 指数衰减/增长模型：y = a * exp(b*x)或y = a * exp(-b*x)
• 幂律模型：y = a * x^b
• 饱和增长模型（如 Logistic 模型）：y = L / (1 + exp(-k*(x-x0)))
• 自定义的复杂物理、化学、生物模型。

对于这类模型，RSS对参数的偏导数方程不再是线性方程组，我们无法像上一节那样直接推导出一个简洁的解析解公式。这时，我们就需要借助数值优化算法来寻找使RSS最小的那组参数。这个过程可以形象地理解为：在一个由参数构成的多维“山地”中，RSS是海拔高度，我们的目标是找到最低的那个山谷。数值算法就是我们的“登山向导”（其实是“下山向导”）。

常用的数值优化算法有很多，在拟合中常见的有：

3.1 梯度下降法及其变种这是最直观的优化思想。既然我们想找最低点，那就沿着当前点最陡的下坡方向走一步。这个“最陡的下坡方向”就是目标函数RSS在当前参数值处的负梯度方向。梯度是一个向量，指向函数值增长最快的方向，那么负梯度自然就是下降最快的方向。

• 步骤：
  1. 随机初始化一组参数值。
  2. 计算当前参数下的RSS及其梯度。
  3. 沿负梯度方向更新参数：新参数 = 旧参数 - 学习率 * 梯度。
  4. 重复步骤2-3，直到RSS的变化小于某个阈值，或达到最大迭代次数。
• 关键点：“学习率”是一个超参数，步子太小收敛慢，步子太大可能跨过最低点甚至发散。现代深度学习中常用的 Adam、RMSprop 等优化器，都是梯度下降的改进版，能自适应调整每个参数的学习率。

3.2 高斯-牛顿法与列文伯格-马夸尔特法这两种方法是专门为最小二乘问题设计的更高效的算法。

• 高斯-牛顿法：它是对牛顿法的近似。牛顿法需要计算目标函数的二阶导数（海森矩阵），计算量大。高斯-牛顿法利用最小二乘问题的特殊结构，用雅可比矩阵（一阶偏导数矩阵）的乘积来近似海森矩阵，大大减少了计算量。它通常比梯度下降收敛更快。
• 列文伯格-马夸尔特法：可以看作是梯度下降和高斯-牛顿法的融合与改进。它在参数空间中进行一种“信任域”搜索：当当前近似很好时，它更像高斯-牛顿法，快速收敛；当近似不好时，它更像梯度下降，采取更保守的步骤。LM算法非常鲁棒，是许多科学计算库（如 SciPy, MATLAB）中非线性最小二乘拟合的默认或首选算法。

实操心得：在实际建模中，我们很少需要自己从头实现这些优化算法。像 Python 的SciPy.optimize.curve_fit函数，或者 MATLAB 的lsqcurvefit函数，内部已经集成了强大的 LM 算法。我们的工作重点是：1. 正确定义模型函数；2. 为算法提供合理的参数初始猜测值。初始值选得好，算法收敛快且容易找到全局最优；初始值差得太远，算法可能陷入局部最优或直接失败。提供初始值时，可以基于对物理背景的理解或通过绘制数据图进行粗略估计。

4. 拟合好坏的“裁判”：评价指标与过拟合陷阱

找到拟合曲线后，我们立刻会面临两个问题：1. 这条曲线拟合得到底“好不好”？2. 是不是用的模型越复杂越好？

4.1 常用的评价指标

• 误差平方和 / 均方误差：SSE = Σ(y_i - ŷ_i)^2，这就是我们最小化的目标RSS本身。MSE = SSE / n。它们衡量的是总体误差的绝对大小，数值越小越好，但其数值依赖于y本身的量纲，不同数据集之间难以直接比较。
• 确定系数：R^2。这是最常用、最直观的指标。它表示模型能够解释的数据波动的比例。计算公式为R^2 = 1 - SSE/SST，其中SST = Σ(y_i - ȳ)^2是总平方和，ȳ是y的平均值。R^2越接近1，说明模型对数据的解释能力越强。在线性拟合中，R^2就是相关系数的平方。
• 调整后的确定系数：Adjusted R^2。当模型参数（变量）增多时，R^2会天然地增大，即使新增的变量没有实际解释力。为了惩罚不必要的复杂度，引入了调整R^2：Adjusted R^2 = 1 - [(1-R^2)*(n-1)/(n-p-1)]，其中p是参数个数。在比较不同复杂度的模型时，调整R^2比普通R^2更公平。
• 均方根误差：RMSE = sqrt(MSE)。它的量纲和原始数据y一致，可以直观理解为“平均每个点的预测误差大概是多少个单位”，在实际业务解释中非常有用。

4.2 过拟合：模型“学傻了”过拟合是拟合（乃至整个机器学习）中最核心的陷阱。它指的是模型在训练数据上表现极好（R^2很高，误差很小），但在新的、未见过的测试数据上表现很差。

• 为什么会产生？模型过于复杂，强大到不仅学到了数据背后的普遍规律，还把数据中的随机噪声、特定采样偏差也当作规律“记”了下来。就像学生死记硬背了所有课后习题的答案，但没理解原理，考试题目一变就不会了。
• 如何识别？最根本的方法是数据集划分。将数据随机分为训练集（如70%）和测试集（如30%）。用训练集数据拟合模型，然后分别在训练集和测试集上计算评价指标（如R^2,RMSE）。如果训练集指标远好于测试集指标（例如训练集R^2=0.99，测试集R^2=0.70），那就是明显的过拟合。
• 如何避免？
  1. 简化模型：优先选择更简洁、参数更少的模型。能用线性就别用二次，能用二次就别用五次。奥卡姆剃刀原理在这里非常适用。
  2. 交叉验证：将数据分成k份（如5份），轮流将其中一份作为测试集，其余作为训练集，重复k次，最后取k次测试结果的平均值作为模型泛化能力的评估。这比单次划分更稳定。
  3. 正则化：在目标函数RSS中加入一个对参数大小的惩罚项。例如岭回归在RSS后加上λΣ(参数^2)，拉索回归加上λΣ|参数|。这样可以在拟合数据的同时，迫使参数值不要太大，从而抑制模型的复杂度，缓解过拟合。参数λ控制了惩罚的力度。

踩坑实录：我曾处理过一组仅含7个数据点的实验数据，试图用6次多项式去拟合（6个参数）。结果R^2高达0.999，曲线完美穿过所有点。但当我用这个公式去预测一个新的x值时，得到的y值荒谬至极。这就是典型的过拟合——模型自由度（参数个数）几乎和数据点一样多，它已经不是在找规律，而是在“连接点”了。最终，我通过观察数据散点图，结合物理背景，选择了一个简单的指数衰减模型（2个参数），虽然训练集R^2只有0.95，但对新数据的预测能力却可靠得多。

5. 从理论到代码：Python/MATLAB实战指南

理解了原理，最终要落到实操上。这里以最常用的 Python 科学计算栈和 MATLAB 为例，展示完整的拟合流程。

5.1 Python 实现（使用 NumPy 和 SciPy）假设我们有一组数据，疑似符合指数衰减关系y = a * exp(-b*x) + c。

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit # 1. 准备数据 x_data = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) y_data = np.array([2.1, 1.2, 0.8, 0.5, 0.3, 0.2, 0.15, 0.1, 0.08, 0.05]) # 2. 定义待拟合的模型函数 def exponential_decay(x, a, b, c): return a * np.exp(-b * x) + c # 3. 执行拟合。p0是参数的初始猜测值，这对非线性拟合至关重要！ initial_guess = (2.0, 0.5, 0.0) # 根据数据趋势粗略估计：a约2，b约0.5，c约0 params, params_covariance = curve_fit(exponential_decay, x_data, y_data, p0=initial_guess) # 4. 获取拟合参数 a_fit, b_fit, c_fit = params print(f"拟合参数: a = {a_fit:.4f}, b = {b_fit:.4f}, c = {c_fit:.4f}") # 5. 计算拟合值及评价指标 y_fit = exponential_decay(x_data, a_fit, b_fit, c_fit) # 计算R^2 ss_res = np.sum((y_data - y_fit) ** 2) ss_tot = np.sum((y_data - np.mean(y_data)) ** 2) r_squared = 1 - (ss_res / ss_tot) print(f"R^2 = {r_squared:.6f}") # 计算RMSE rmse = np.sqrt(np.mean((y_data - y_fit) ** 2)) print(f"RMSE = {rmse:.6f}") # 6. 可视化 plt.figure(figsize=(10, 6)) plt.scatter(x_data, y_data, label='原始数据', color='blue', s=50) plt.plot(x_data, y_fit, label=f'拟合曲线: y={a_fit:.2f}*exp(-{b_fit:.2f}x)+{c_fit:.2f}', color='red', linewidth=2) plt.xlabel('X') plt.ylabel('Y') plt.title('非线性最小二乘拟合示例') plt.legend() plt.grid(True, linestyle='--', alpha=0.7) plt.show()

5.2 MATLAB 实现MATLAB 的语法更加直观，特别适合矩阵运算和快速建模。

% 1. 准备数据 x_data = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]; y_data = [2.1, 1.2, 0.8, 0.5, 0.3, 0.2, 0.15, 0.1, 0.08, 0.05]; % 2. 定义模型函数句柄 model = @(params, x) params(1) * exp(-params(2) * x) + params(3); % 3. 执行拟合 (使用 lsqcurvefit) initial_guess = [2.0, 0.5, 0.0]; % 初始猜测值 options = optimoptions('lsqcurvefit', 'Display', 'iter'); % 显示迭代过程 [params, resnorm, residual, exitflag, output] = lsqcurvefit(model, initial_guess, x_data, y_data); % 4. 获取拟合参数 a_fit = params(1); b_fit = params(2); c_fit = params(3); fprintf('拟合参数: a = %.4f, b = %.4f, c = %.4f\n', a_fit, b_fit, c_fit); % 5. 计算拟合值及评价指标 y_fit = model(params, x_data); % 计算R^2 ss_res = sum((y_data - y_fit).^2); ss_tot = sum((y_data - mean(y_data)).^2); r_squared = 1 - (ss_res / ss_tot); fprintf('R^2 = %.6f\n', r_squared); % 计算RMSE rmse = sqrt(mean((y_data - y_fit).^2)); fprintf('RMSE = %.6f\n', rmse); % 6. 可视化 figure; scatter(x_data, y_data, 50, 'b', 'filled'); hold on; plot(x_data, y_fit, 'r-', 'LineWidth', 2); xlabel('X'); ylabel('Y'); title('非线性最小二乘拟合示例 (MATLAB)'); legend('原始数据', sprintf('拟合曲线: y=%.2f*exp(-%.2fx)+%.2f', a_fit, b_fit, c_fit)); grid on; hold off;

关键技巧：无论是 Python 还是 MATLAB，非线性拟合的成功率极大依赖于initial_guess（初始猜测值）。一个实用的方法是先绘制数据散点图，根据图形趋势手动估算大致参数。例如，对于指数衰减y = a*exp(-b*x)，当x=0时y≈a，衰减速度由b控制。也可以先对数据取对数，将指数拟合转化为线性拟合来获得粗略的初始值。

6. 进阶话题：稳健回归与异常值处理

最小二乘法有一个潜在的软肋：它对异常值非常敏感。因为它的目标是最小化平方误差，一个偏离很远的异常点会产生巨大的平方误差，为了“讨好”这个异常点，拟合曲线可能会被强行“拉偏”，导致对大多数正常数据点的拟合变差。

6.1 异常值的影响想象一下，你在用尺子画一条穿过一堆点的直线，大部分点都密集在一条狭长地带，但有一个点离得非常远。最小二乘法就像一根有弹性的橡皮筋，它会被那个远处的点狠狠拉扯，使得直线为了靠近那个点而偏离了主要点群的中心。这显然不是我们想要的。

6.2 稳健回归方法为了解决这个问题，统计学家提出了稳健回归。其核心思想是降低异常点在目标函数中的权重。常见的方法有：

• 最小一乘法：将目标函数从最小化平方和改为最小化绝对值和RSS = Σ |y_i - f(x_i)|。绝对值对极端值的惩罚小于平方，因此更稳健。但绝对值函数在零点不可导，求解需要用线性规划或其他迭代方法，计算比最小二乘复杂。
• M-估计法：这是更一般的框架。它用一个增长慢于平方函数的ρ函数来代替平方。例如 Huber 损失函数，它在误差较小时是二次的（保持效率），在误差较大时是线性的（降低异常值影响）。目标函数为Σ ρ(y_i - f(x_i))。
• RANSAC：这是一种完全不同的思路。它通过迭代随机采样子集来拟合模型。具体步骤：
  1. 随机选择能确定模型的最小数据点集（如线性拟合选2个点）。
  2. 用这些点拟合一个模型。
  3. 找出所有数据中，符合该模型（误差小于某个阈值）的点，这些点称为“内点”。
  4. 用所有内点重新拟合模型。
  5. 重复以上过程多次，最终选择内点最多（或内点误差最小）的那个模型。 RANSAC 能很好地从包含大量异常值的数据中估计出模型参数，在计算机视觉中应用极广。

6.3 实践建议在数学建模中，处理异常值的标准流程应该是：

1. 可视化：首先绘制散点图，直观检查是否有明显远离群体的点。
2. 诊断：用普通最小二乘拟合后，绘制残差图（残差 vs. 拟合值或 vs. 自变量）。如果残差图呈现随机分布，说明模型基本合适；如果发现有规律的模式或个别点残差极大，则可能存在模型误设或异常值。
3. 处理：
   • 谨慎剔除：如果确认某个点是记录错误或实验失误导致，可以剔除。但必须有充分理由，不能为了追求高R^2而随意删点。
   • 使用稳健方法：当怀疑存在异常值但无法或不宜剔除时，应采用稳健回归方法重新拟合，并对比其结果与普通最小二乘结果的差异。
   • 数据变换：有时对因变量y进行变换（如取对数、开方）可以稳定方差，减少异常值的影响。

我个人在处理一组关于城市交通流量的数据时，就曾遇到一个点：因为当天有大型活动，流量激增，成为明显的异常值。使用普通最小二乘拟合的直线被它明显拉高。我首先分析了原因，确认是特殊事件后，如果建模目的是研究日常流量规律，则可以选择剔除该点；如果建模需要包含此类事件，则保留该点，但同时在论文中明确指出其影响，并考虑使用稳健回归或引入事件哑变量来改进模型。