直线拟合三大核心方法:最小二乘法、梯度下降与高斯-牛顿法详解

📅 2026/8/1 9:29:36 👁️ 阅读次数 📝 编程学习
直线拟合三大核心方法:最小二乘法、梯度下降与高斯-牛顿法详解

1. 从“画一条线”到“找一条线”:直线拟合的本质是什么?

在数据分析、图像处理、机器人定位、金融建模等无数领域,我们常常会遇到一个看似简单却至关重要的任务:给定一组离散的数据点,如何找到一条最能代表它们整体趋势的直线?这就是直线拟合。新手可能会觉得,这不就是凭感觉画一条线吗?但当你面对成千上万个点,或者这条线的斜率直接决定了某个物理参数、预测了明天的股价、校准了传感器的精度时,“凭感觉”就完全不可靠了。直线拟合,就是从数学和计算的角度,把“画线”这个主观行为,变成一个客观、可量化、可复现的求解过程。

今天,我们就来深入聊聊三种最核心、最实用的直线拟合方法:最小二乘法、梯度下降法以及高斯-牛顿法(及其变种列文伯格-马夸尔特算法)。这三种方法并非简单的并列关系,它们背后代表了不同的数学思想和适用场景。最小二乘法是“一步到位”的解析解,优雅但有其局限;梯度下降法是“步步为营”的迭代通用解,稳健但可能缓慢;高斯-牛顿法则是针对特定问题(非线性最小二乘)的“加速专家”。理解它们,你就能在面对“找直线”这个问题时,不再迷茫,而是能根据数据特点、精度要求和计算资源,做出最合适的选择。

2. 最小二乘法:经典解析解的优雅与局限

最小二乘法是直线拟合领域当之无愧的“老大哥”,它的核心思想直观而强大:寻找一条直线,使得所有数据点到这条直线的垂直距离(即残差)的平方和最小。为什么是平方和?而不是直接求距离和?这主要是为了数学上的便利:平方操作让距离函数处处可导(避免了绝对值的不可导点),并且对大误差给予更大的惩罚,使得拟合结果对异常值不那么敏感(虽然仍有一定敏感性)。

2.1 数学推导与“闭式解”

假设我们有n个数据点(x_i, y_i),要拟合的直线方程为y = kx + b。我们的目标是找到参数k(斜率) 和b(截距),使得损失函数L最小:L(k, b) = Σ(y_i - (k*x_i + b))^2

这是一个关于kb的二元二次函数。为了找到最小值,我们分别对kb求偏导数,并令其等于零:∂L/∂k = -2 * Σ[x_i * (y_i - k*x_i - b)] = 0∂L/∂b = -2 * Σ(y_i - k*x_i - b) = 0

整理后,我们得到一个关于kb的线性方程组,这就是著名的正规方程

Σ(x_i^2) * k + Σ(x_i) * b = Σ(x_i * y_i) Σ(x_i) * k + n * b = Σ(y_i)

这个方程组可以直接求解,得到kb的解析表达式(闭式解):

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

这个解是唯一的,并且可以通过一次矩阵运算(求解正规方程)或直接套用上述公式得到,计算效率非常高。

注意:在计算时,分母(n * Σ(x_i^2) - (Σ(x_i))^2)可能接近于零。这通常发生在所有x_i值都相同或非常接近时,意味着数据点在x方向上几乎没有变化,此时直线斜率趋于无穷大(垂直线),最小二乘法的这种形式失效。在实际编程中,需要加入判断以避免除以零的错误。

2.2 Python实现与实战

用Python的NumPy库实现最小二乘法拟合非常简洁,我们可以用两种方式:一种是基于上述公式手动计算,另一种是利用NumPy的线性代数功能。

import numpy as np import matplotlib.pyplot as plt # 生成示例数据 np.random.seed(42) x = np.linspace(0, 10, 50) true_k, true_b = 2.5, 1.0 y = true_k * x + true_b + np.random.randn(50) * 2 # 添加噪声 # 方法1:手动套用公式 def linear_regression_manual(x, y): n = len(x) sum_x = np.sum(x) sum_y = np.sum(y) sum_xy = np.sum(x * y) sum_x2 = np.sum(x ** 2) denominator = n * sum_x2 - sum_x ** 2 if abs(denominator) < 1e-10: raise ValueError("数据点在x方向上无变化,无法计算斜率。") k = (n * sum_xy - sum_x * sum_y) / denominator b = (sum_y - k * sum_x) / n return k, b k_manual, b_manual = linear_regression_manual(x, y) print(f"手动计算: 斜率 k = {k_manual:.4f}, 截距 b = {b_manual:.4f}") # 方法2:使用NumPy的polyfit(1次多项式拟合) coefficients = np.polyfit(x, y, 1) # deg=1 表示一次多项式,即直线 k_np, b_np = coefficients print(f"NumPy polyfit: 斜率 k = {k_np:.4f}, 截距 b = {b_np:.4f}") # 方法3:使用正规方程矩阵求解 (更通用的形式,便于扩展到多元) # 构造设计矩阵 X,增加一列1用于截距 X = np.vstack([x, np.ones_like(x)]).T # 正规方程解: theta = (X^T * X)^(-1) * X^T * y theta = np.linalg.inv(X.T @ X) @ X.T @ y k_matrix, b_matrix = theta print(f"矩阵求解: 斜率 k = {k_matrix:.4f}, 截距 b = {b_matrix:.4f}") # 可视化 plt.scatter(x, y, alpha=0.6, label='原始数据') plt.plot(x, k_manual * x + b_manual, 'r-', linewidth=2, label=f'拟合直线: y={k_manual:.2f}x+{b_manual:.2f}') plt.xlabel('X') plt.ylabel('Y') plt.legend() plt.grid(True, linestyle='--', alpha=0.5) plt.title('最小二乘法直线拟合') plt.show()

三种方法的结果在数值上应该完全一致(忽略微小的浮点数误差)。np.polyfit是最方便快捷的;手动公式有助于理解原理;矩阵形式则是理解更复杂回归模型(如多元线性回归)的基础。

2.3 优势、局限与常见陷阱

优势

  1. 解析解,计算高效:一次计算即可得到全局最优解,速度极快,尤其适合数据量不大或需要频繁拟合的场景。
  2. 理论完备:有严格的统计解释。在误差满足独立同分布且服从正态分布的假设下,最小二乘估计是参数的最佳线性无偏估计。
  3. 实现简单:公式清晰,几乎所有编程语言和数据分析库都有现成实现。

局限与陷阱

  1. 对异常值敏感:由于损失函数使用误差平方,个别远离群体的异常点(Outliers)会对拟合结果产生巨大的“拉扯”效应,导致直线严重偏离大多数数据点。在实际应用中,拟合前进行异常值检测或使用鲁棒性更强的损失函数(如Huber损失)是常见做法。
  2. 仅限于线性参数:这里指的是参数kb以线性形式出现在方程y = kx + b中。最小二乘法本身解决的是线性参数估计问题。如果模型本身是非线性的(例如y = a * exp(b*x)),则需要通过变换或使用非线性最小二乘法(这正是高斯-牛顿法的用武之地)。
  3. 矩阵求逆的数值稳定性:在矩阵求解形式(X^T X)^(-1)中,当X^T X矩阵接近奇异(即列向量之间存在近似线性关系,称为多重共线性)时,求逆运算会变得非常不稳定,导致结果误差极大。对于这种情况,通常采用岭回归(Ridge Regression)或使用更稳定的数值算法(如奇异值分解SVD)来求解。

3. 梯度下降法:通用迭代求解的“慢工细活”

当问题变得复杂,无法像最小二乘法那样直接求出解析解时,梯度下降法就登场了。它的思想源于最朴素的直觉:如果你想最快地下到山谷底部,就沿着当前最陡峭的方向往下走。在优化问题中,“山谷”就是我们的损失函数曲面,“最陡峭的方向”就是该点损失函数的负梯度方向。

3.1 核心思想与迭代过程

对于直线拟合,损失函数依然是L(k, b) = Σ(y_i - (k*x_i + b))^2。梯度下降法通过以下步骤迭代更新参数kb

  1. 初始化:随机猜测一组参数值,例如k=0,b=0
  2. 计算梯度:计算损失函数在当前参数(k, b)处的梯度。梯度是一个向量,其两个分量分别是Lkb的偏导数:∂L/∂k = -2 * Σ[x_i * (y_i - k*x_i - b)]∂L/∂b = -2 * Σ(y_i - k*x_i - b)这个梯度的方向指向了损失函数在当前点上升最快的方向。
  3. 沿负梯度方向更新:为了让损失函数减小,我们沿着梯度的反方向(即负梯度方向)迈出一步。更新公式为:k_new = k_old - α * (∂L/∂k)b_new = b_old - α * (∂L/∂b)其中α是一个关键的超参数,称为学习率。它决定了每一步迈多大。
  4. 重复迭代:用新的(k_new, b_new)替换旧的参数,回到第2步,直到满足停止条件(例如梯度变得非常小、损失函数变化很小,或达到预设的迭代次数)。

3.2 学习率的选择与挑战

学习率α是梯度下降法的“命门”。选择不当会导致严重问题:

  • 学习率过大:更新步伐太大,可能会在最小值点附近来回震荡,甚至直接“飞越”最小值,导致算法无法收敛,甚至发散(损失函数越来越大)。
  • 学习率过小:更新步伐太小,收敛速度会非常缓慢,需要大量的迭代步数才能接近最优解,计算成本高。

在实际操作中,我通常会从一个较小的值开始尝试(如0.001或0.01),观察损失函数在迭代过程中的下降曲线。一个健康的下降曲线应该是初期快速下降,后期平缓收敛。如果曲线震荡,就调小学习率;如果下降太慢,可以适当调大。更高级的策略是使用自适应学习率的优化器,如Adam、Adagrad等,它们能根据历史梯度信息动态调整每个参数的学习率,在实践中(尤其是深度学习领域)几乎成为标配。

3.3 Python实现:从零开始与优化器对比

让我们手动实现一个基础的批量梯度下降(Batch Gradient Descent),并与使用PyTorch内置优化器的版本进行对比。

import numpy as np import matplotlib.pyplot as plt import torch import torch.optim as optim # 使用相同的数据 np.random.seed(42) x_np = np.linspace(0, 10, 50) true_k, true_b = 2.5, 1.0 y_np = true_k * x_np + true_b + np.random.randn(50) * 2 # 转换为PyTorch张量,便于后续使用优化器 x_tensor = torch.from_numpy(x_np).float() y_tensor = torch.from_numpy(y_np).float() # --- 方法1:手动实现梯度下降 --- def gradient_descent_manual(x, y, lr=0.01, epochs=1000): """ 手动实现批量梯度下降 """ n = len(x) k, b = 0.0, 0.0 # 初始化参数 history = {'loss': [], 'k': [k], 'b': [b]} # 记录历史 for epoch in range(epochs): # 计算预测值 y_pred = k * x + b # 计算损失 (均方误差) loss = np.mean((y - y_pred) ** 2) history['loss'].append(loss) # 计算梯度 dk = (-2/n) * np.sum(x * (y - y_pred)) db = (-2/n) * np.sum(y - y_pred) # 更新参数 k = k - lr * dk b = b - lr * db history['k'].append(k) history['b'].append(b) # 简单停止条件:梯度很小 if epoch % 200 == 0: print(f'Epoch {epoch}: loss={loss:.4f}, k={k:.4f}, b={b:.4f}') if np.sqrt(dk**2 + db**2) < 1e-5: print(f'在 epoch {epoch} 提前收敛') break return k, b, history k_gd, b_gd, history_gd = gradient_descent_manual(x_np, y_np, lr=0.02, epochs=2000) print(f"\n手动梯度下降结果: k={k_gd:.4f}, b={b_gd:.4f}") # --- 方法2:使用PyTorch和Adam优化器 --- def gradient_descent_torch(x, y, lr=0.1, epochs=1000): """ 使用PyTorch和Adam优化器 """ # 定义需要优化的参数 k = torch.tensor(0.0, requires_grad=True) b = torch.tensor(0.0, requires_grad=True) # 选择优化器,这里使用Adam optimizer = optim.Adam([k, b], lr=lr) history_torch = {'loss': [], 'k': [k.item()], 'b': [b.item()]} for epoch in range(epochs): # 前向传播:计算预测和损失 y_pred = k * x + b loss = torch.mean((y - y_pred) ** 2) # 反向传播:计算梯度 optimizer.zero_grad() # 清除旧梯度 loss.backward() # 计算新梯度 # 更新参数 optimizer.step() history_torch['loss'].append(loss.item()) history_torch['k'].append(k.item()) history_torch['b'].append(b.item()) if epoch % 200 == 0: print(f'Epoch {epoch}: loss={loss.item():.4f}, k={k.item():.4f}, b={b.item():.4f}') if epoch > 10 and abs(history_torch['loss'][-1] - history_torch['loss'][-2]) < 1e-7: print(f'在 epoch {epoch} 提前收敛') break return k.item(), b.item(), history_torch k_torch, b_torch, history_torch = gradient_descent_torch(x_tensor, y_tensor, lr=0.1, epochs=1000) print(f"\nPyTorch Adam优化器结果: k={k_torch:.4f}, b={b_torch:.4f}") # 与最小二乘法结果对比 k_ls, b_ls = np.polyfit(x_np, y_np, 1) print(f"\n最小二乘法结果: k={k_ls:.4f}, b={b_ls:.4f}") print(f"差异: Δk = {abs(k_gd - k_ls):.6f}, Δb = {abs(b_gd - b_ls):.6f}") # 可视化拟合结果和损失下降曲线 fig, axes = plt.subplots(1, 2, figsize=(12, 4)) # 左图:拟合直线对比 axes[0].scatter(x_np, y_np, alpha=0.5, label='数据') axes[0].plot(x_np, k_ls*x_np + b_ls, 'r-', label=f'最小二乘 (基准)') axes[0].plot(x_np, k_gd*x_np + b_gd, 'g--', label=f'手动GD') axes[0].plot(x_np, k_torch*x_np + b_torch, 'b:', linewidth=2, label=f'PyTorch Adam') axes[0].set_xlabel('X') axes[0].set_ylabel('Y') axes[0].legend() axes[0].grid(True, linestyle='--', alpha=0.5) axes[0].set_title('不同方法拟合直线对比') # 右图:损失下降曲线 axes[1].plot(history_gd['loss'], label='手动GD (lr=0.02)') axes[1].plot(history_torch['loss'], label='PyTorch Adam (lr=0.1)') axes[1].set_xlabel('迭代次数') axes[1].set_ylabel('损失 (MSE)') axes[1].set_yscale('log') # 使用对数坐标更清晰地观察下降 axes[1].legend() axes[1].grid(True, linestyle='--', alpha=0.5) axes[1].set_title('损失函数下降曲线 (对数坐标)') plt.tight_layout() plt.show()

运行这段代码,你会发现几个关键点:

  1. 收敛性:手动实现的梯度下降和PyTorch的Adam优化器最终都能收敛到与最小二乘法非常接近的结果(差异在可接受的浮点误差范围内)。
  2. 收敛速度:Adam优化器通常比固定学习率的手动梯度下降收敛得更快、更稳定,这得益于其自适应学习率机制。
  3. 学习率敏感度:手动梯度下降对学习率lr非常敏感。如果我把lr从0.02改为0.05,可能会看到震荡;改为0.001,则需要更多迭代次数。而Adam在lr=0.1时依然能稳定收敛,显示了其鲁棒性。

3.4 适用场景与心得

梯度下降法的最大优势在于其通用性。它不仅适用于线性模型的参数求解,更是训练神经网络、逻辑回归、支持向量机等几乎所有复杂机器学习模型的基石。当模型没有解析解,或者数据量太大以至于无法一次性加载计算(此时可以使用随机梯度下降SGD小批量梯度下降)时,梯度下降法是唯一可行的选择。

我的几点实操心得

  • 监控是关键:始终绘制损失函数随迭代次数的变化曲线。这是诊断学习率是否合适、算法是否收敛的最直观工具。
  • 初始化很重要:虽然对于凸问题(如线性回归)梯度下降最终能收敛到全局最优,但好的初始化可以大大减少迭代次数。通常可以用最小二乘法的解作为梯度下降的初始值,这是一个非常有效的“热启动”策略。
  • 试试Adam:对于大多数不太极端的问题,使用Adam优化器(默认参数lr=0.001)通常是一个安全且高效的选择,它能省去大量手动调参的麻烦。

4. 高斯-牛顿法与列文伯格-马夸尔特算法:非线性最小二乘的利器

前面讨论的最小二乘法和梯度下降法,主要针对的是线性模型y = kx + b。但在现实中,大量关系是非线性的,例如指数衰减y = a * exp(-b*x)、幂律关系y = a * x^b等。对于这类问题,我们通常将其转化为非线性最小二乘问题:寻找一组参数θ,使得残差平方和Σ [y_i - f(x_i; θ)]^2最小,其中f是非线性函数。

高斯-牛顿法和列文伯格-马夸尔特算法就是专门为解决非线性最小二乘问题而设计的优化算法。它们可以看作是梯度下降法的“升级版”和“稳健版”。

4.1 高斯-牛顿法:利用局部线性化的快速收敛

高斯-牛顿法的核心思想是迭代重加权线性最小二乘。在每一次迭代中,它都对非线性函数f(x; θ)在当前参数估计值θ_k处进行一阶泰勒展开(即线性化):f(x; θ) ≈ f(x; θ_k) + J(θ_k) * (θ - θ_k)其中J(θ_k)是函数f关于参数θθ_k处的雅可比矩阵(一阶偏导数矩阵)。

将线性化后的近似代入损失函数,原来的非线性最小二乘问题就变成了一个关于参数增量Δθ = θ - θ_k线性最小二乘问题。求解这个线性问题,得到参数增量Δθ,然后更新参数:θ_{k+1} = θ_k + Δθ。重复这个过程直到收敛。

优势:当初始猜测接近真实解,且残差较小时,高斯-牛顿法具有二次收敛速度,比梯度下降法快得多。致命缺点:它要求近似的海森矩阵J^T J是良态的(可逆且条件数好)。如果J^T J接近奇异,或者初始猜测离解太远导致线性化近似很差,算法可能根本不收敛,甚至发散。

4.2 列文伯格-马夸尔特算法:自适应信赖域的稳健策略

列文伯格-马夸尔特算法(简称L-M算法)可以看作是高斯-牛顿法和梯度下降法之间的一个自适应桥梁。它通过引入一个阻尼因子λ来巧妙地解决高斯-牛顿法的不稳定问题。

L-M算法的参数更新公式为:(J^T J + λ * I) * Δθ = J^T * r其中r是残差向量,I是单位矩阵。

这个公式非常精妙:

  • λ很大时,λ * I占主导,方程近似为λ * I * Δθ ≈ J^T * r,即Δθ ≈ (1/λ) * J^T * r。这其实就是梯度下降法的方向(J^T * r是梯度的负方向),只是步长受λ控制。此时算法行为像梯度下降,稳健但收敛慢,适用于离解较远的情况。
  • λ很小时,J^T J占主导,方程退化为高斯-牛顿法的正规方程。此时算法行为像高斯-牛顿法,收敛速度快,适用于接近解的情况。

L-M算法的智能之处在于,它在每次迭代中动态调整λ

  1. 计算试探步长Δθ
  2. 用新参数θ_new = θ + Δθ计算实际损失减少量。
  3. 与基于线性模型预测的损失减少量进行比较。
  4. 如果实际减少量符合预期(甚至更好),则接受这一步,并减小λ(信任模型,下次更接近高斯-牛顿法)。
  5. 如果实际减少量不符合预期,则拒绝这一步,增大λ(不信任模型,下次更接近梯度下降法,步长更小更谨慎)。

这种机制使得L-M算法既能拥有高斯-牛顿法在接近解时的快速收敛性,又具备梯度下降法的全局稳健性。

4.3 Python实战:拟合指数衰减曲线

让我们用一个具体的例子来演示L-M算法的威力。假设我们有一组数据,它遵循指数衰减模型y = a * exp(-b * x) + c,我们要拟合参数[a, b, c]。这是一个典型的非线性最小二乘问题。

我们将使用SciPy库中的curve_fit函数,它内部默认使用的就是L-M算法。

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit, least_squares # 1. 生成模拟的非线性数据 (指数衰减 + 基线) np.random.seed(123) x_data = np.linspace(0, 5, 50) a_true, b_true, c_true = 5.0, 1.2, 0.5 y_true = a_true * np.exp(-b_true * x_data) + c_true # 添加噪声 noise = np.random.randn(len(x_data)) * 0.2 y_data = y_true + noise # 2. 定义要拟合的非线性模型函数 def exp_decay(x, a, b, c): """指数衰减模型:y = a * exp(-b*x) + c""" return a * np.exp(-b * x) + c # 3. 使用SciPy的curve_fit进行拟合(默认使用L-M算法) # 提供参数的初始猜测值,这对非线性拟合很重要 initial_guess = [1.0, 0.5, 0.0] # 猜测 [a, b, c] popt, pcov = curve_fit(exp_decay, x_data, y_data, p0=initial_guess) a_fit, b_fit, c_fit = popt print(f"真实参数: a={a_true:.3f}, b={b_true:.3f}, c={c_true:.3f}") print(f"拟合参数: a={a_fit:.3f}, b={b_fit:.3f}, c={c_fit:.3f}") print(f"参数协方差矩阵的对角线(方差):\n{np.diag(pcov)}") # 4. 计算拟合优度 R-squared residuals = y_data - exp_decay(x_data, *popt) ss_res = np.sum(residuals**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:.4f}") # 5. 可视化 x_fine = np.linspace(0, 5, 200) y_fine_fit = exp_decay(x_fine, *popt) plt.figure(figsize=(10, 6)) plt.scatter(x_data, y_data, alpha=0.7, label='带噪声数据', color='blue') plt.plot(x_data, y_true, 'k--', linewidth=2, label='真实模型', alpha=0.8) plt.plot(x_fine, y_fine_fit, 'r-', linewidth=2, label=f'L-M算法拟合: y={a_fit:.2f}*exp(-{b_fit:.2f}x)+{c_fit:.2f}') plt.fill_between(x_fine, exp_decay(x_fine, *(popt - 1.96*np.sqrt(np.diag(pcov)))), exp_decay(x_fine, *(popt + 1.96*np.sqrt(np.diag(pcov)))), color='red', alpha=0.2, label='95%置信区间') plt.xlabel('X') plt.ylabel('Y') plt.legend() plt.grid(True, linestyle='--', alpha=0.5) plt.title('列文伯格-马夸尔特算法拟合指数衰减曲线') plt.show() # 6. 对比:如果使用不合适的初始猜测 print("\n--- 测试不同初始猜测的影响 ---") bad_guess = [10.0, 0.1, 2.0] # 一个很差的初始值 try: popt_bad, _ = curve_fit(exp_decay, x_data, y_data, p0=bad_guess, maxfev=5000) # 增加最大迭代次数 print(f"差初始值拟合结果: {popt_bad}") except RuntimeError as e: print(f"差初始值导致拟合失败: {e}") # 可以尝试使用更鲁棒的方法,比如差分进化算法提供初始值 from scipy.optimize import differential_evolution # 定义参数边界 bounds = [(0, 10), (0, 5), (-2, 2)] def sum_of_squares(params): a, b, c = params y_pred = exp_decay(x_data, a, b, c) return np.sum((y_data - y_pred) ** 2) result = differential_evolution(sum_of_squares, bounds, maxiter=100, seed=42) print(f"差分进化算法提供的初始值: {result.x}") # 再用L-M算法精细优化 popt_refined, _ = curve_fit(exp_decay, x_data, y_data, p0=result.x) print(f"经L-M算法精细优化后: {popt_refined}")

4.4 关键要点与选择建议

通过这个例子,我们可以总结出关于高斯-牛顿法和L-M算法的几个关键点:

  1. 初始值至关重要:对于非线性问题,损失函数可能存在多个局部极小值。算法的收敛结果严重依赖于初始猜测。一个糟糕的初始值可能导致算法收敛到错误的局部最优,甚至发散。在实践中,通常需要:

    • 基于物理意义或经验给出初始值。
    • 使用全局优化算法(如差分进化、模拟退火)先进行粗略搜索,再用L-M算法进行精细优化。
    • 多次尝试不同的随机初始值,选择损失最小的结果。
  2. L-M算法是实际首选:由于L-M算法集成了梯度下降的稳健性和高斯-牛顿的快速收敛性,它已成为解决非线性最小二乘问题的事实标准。SciPy的curve_fitleast_squares,MATLAB的lsqnonlin,以及许多其他科学计算库的默认算法都是L-M算法。

  3. 理解输出信息curve_fit返回的pcov是参数估计的协方差矩阵,其对角线元素的平方根给出了参数的标准误差,可用于计算置信区间。R^2值可以量化拟合效果,但要注意对于非线性模型,R^2的解释与线性模型略有不同。

如何在这三种方法中选择?

  • 如果你的模型关于参数是线性的(如y = kx + b),且数据量不大、没有异常值困扰,最小二乘法是你的首选,因为它简单、快速、精确。
  • 如果你的模型关于参数是线性的,但数据量极大(无法一次性计算),或者你正在训练一个更复杂的模型(如神经网络),那么梯度下降法(及其变种SGD, Adam)是必由之路。
  • 如果你的模型关于参数是非线性的(如指数、对数、幂函数等),那么列文伯格-马夸尔特算法是你应该首先尝试的工具。它高效、稳健,并且有成熟的库支持。

直线拟合的旅程,从一条可以直接写出的公式,到需要迭代探索的优化路径,再到处理更复杂非线性关系的稳健策略,背后是数学工具不断适应现实问题复杂性的演进。掌握这三种方法,你就拥有了从处理简单线性关系到攻克复杂非线性拟合问题的全套工具箱。下次当你面对一堆散点图时,你将清楚地知道,该用哪把“钥匙”去解开数据背后的趋势之谜。