高斯牛顿法原理与Matlab实现：非线性最小二乘拟合实战指南

1. 从“猜”到“算”：非线性参数拟合的工程困境

在科研和工程领域，我们常常会遇到一个经典问题：手里有一个描述现象的数学模型，比如描述化学反应速率的阿伦尼乌斯方程、描述生物种群增长的逻辑斯蒂方程，或者描述传感器输出与温度关系的校准曲线。这个模型的数学形式（解析式）是已知的，但里面有几个关键的“魔法数字”——也就是未知参数——我们不知道。我们的任务，就是利用一组实际观测到的数据，把这些参数给“揪”出来。

这听起来像是一个“猜谜”游戏，但数学家们把它变成了一个“计算”问题，这就是参数拟合或参数估计。对于线性模型，我们有最小二乘法这种优雅且直接的解法。但现实世界往往是弯曲的、非线性的。当模型关于参数是非线性的时候，比如指数函数、对数函数或者更复杂的组合，问题就变得棘手了。你没法直接套用公式得到一个闭式解。

这时候，迭代优化算法就登场了。它们的基本思路是：先“蒙”一组参数初始值，然后计算当前参数下模型预测值与实际观测值的差距（残差），接着根据某种规则，朝着让差距缩小的方向“微调”参数，如此反复，直到差距小到我们认为可以接受为止。高斯牛顿法，就是这类算法中一位久经沙场、效率突出的“老将”。它特别擅长处理最小二乘形式的问题，也就是目标是最小化残差平方和的情况。在数学建模竞赛和许多工程实践中，当你需要在Matlab环境下快速、可靠地解决一个中小规模的非线性最小二乘问题时，高斯牛顿法往往是工具箱里的首选利器。

2. 高斯牛顿法：在“局部”寻找最优路径

在深入代码之前，我们必须先理解高斯牛顿法到底在做什么。它不是一个黑箱，理解了其内核，你才能用好它，并在它“卡壳”时知道如何调整。

2.1 核心思想：用线性近似解决非线性问题

高斯牛顿法的聪明之处在于“以直代曲”。对于一个非线性函数f(x, β)，其中x是自变量（可能是一个向量），β是我们要求解的未知参数向量。假设我们有m个观测数据点(x_i, y_i)，我们的目标是找到参数β，使得模型预测值f(x_i, β)尽可能接近y_i。

我们定义残差r_i(β) = y_i - f(x_i, β)。目标函数（损失函数）就是所有残差的平方和：S(β) = Σ [r_i(β)]^2。最小化S(β)就是我们的任务。

现在，假设我们有一个参数猜测值β_k（第k次迭代的值）。在β_k这个点附近，我们可以把非线性的残差函数r_i(β)用一阶泰勒展开来近似，也就是把它线性化：r_i(β) ≈ r_i(β_k) + J_i(β_k) * (β - β_k)其中，J_i(β_k)是残差r_i在β_k处关于参数向量β的梯度（雅可比矩阵的行）。注意，这里是对参数β求导，而不是对自变量x。J_i(β_k) = -∂f(x_i, β)/∂β |_{β=β_k}。因为f是非线性的，所以这个导数通常依赖于β。

将所有的残差堆叠成一个向量r(β) = [r_1, r_2, ..., r_m]^T，所有的雅可比行堆叠成矩阵J(β_k)，那么线性化后的残差向量可以写成：r(β) ≈ r(β_k) + J(β_k) * (β - β_k)

我们的目标函数S(β)就近似为：S(β) ≈ || r(β_k) + J(β_k) * (β - β_k) ||^2这里||·||表示向量的2-范数（平方和）。你看，现在这个问题变成了一个关于Δβ = β - β_k的线性最小二乘问题！因为r(β_k)和J(β_k)在本次迭代中都是已知的常数。

2.2 迭代步骤与“正规方程”

对于线性最小二乘问题min ||b - A*x||^2，其解析解可以通过求解正规方程(A^T A) x = A^T b得到。对应到我们的近似问题：

• A对应J(β_k)
• b对应-r(β_k)
• x对应参数增量Δβ

因此，我们每一步迭代需要求解的方程是：[J(β_k)^T J(β_k)] Δβ = -J(β_k)^T r(β_k)这个方程被称为高斯牛顿方程。解出Δβ后，我们就更新参数：β_{k+1} = β_k + Δβ

然后，用新的β_{k+1}计算新的残差和雅可比矩阵，重复上述过程，直到满足停止条件（例如Δβ的范数很小，或者目标函数S(β)下降不明显）。

注意：这里有一个关键的细节。雅可比矩阵J的元素是J_ij = ∂r_i/∂β_j = -∂f(x_i, β)/∂β_j。在实际编程中，我们需要提供这个导数信息。有时我们可以推导出解析的导数公式，这样最精确高效；如果不行，也可以用数值差分（如有限差分）来近似，但这会引入误差并增加计算量。

2.3 优势与天生的缺陷

高斯牛顿法的优势很明显：它通常比最速下降法收敛快得多。因为它利用了目标函数（平方和）的二阶信息（近似海森矩阵J^T J），迭代方向更智能。在许多问题中，它表现出接近二阶收敛的速度。

但它也有两个著名的“阿喀琉斯之踵”：

1. 初始值依赖性强：因为它基于局部线性近似，如果初始猜测β_0离真实解太远，线性近似可能非常糟糕，导致算法收敛到错误的局部极小点，甚至发散。
2. 矩阵病态问题：方程中的H = J^T J矩阵要求是正定的才能求解。如果J是病态的（即列近似线性相关，意味着参数之间存在强耦合或某个参数对输出影响甚微），那么H可能奇异或病态，导致Δβ的计算极不稳定，数值误差巨大。

为了解决第二个问题，实践中更常用的是Levenberg-Marquardt (L-M) 算法。你可以把L-M算法理解为高斯牛顿法和最速下降法的“平滑切换”。它在高斯牛顿方程中加了一个阻尼因子λ：[J^T J + λ I] Δβ = -J^T r当λ很大时，方程近似为λ I Δβ = -J^T r，即Δβ方向接近最速下降方向（步长很小），适合在远离解时使用；当λ很小时，方程退化为标准高斯牛顿方程，适合在接近解时快速收敛。L-M算法会根据每次迭代的效果动态调整λ，鲁棒性更强。我们实现的高斯牛顿法可以看作是L-M算法在λ=0时的一个特例，理解它对于掌握更高级的L-M算法至关重要。

3. 手把手实现：一个完整的Matlab案例

理论说得再多，不如动手实现一遍。我们通过一个具体的例子，将上述过程转化为Matlab代码。假设我们要拟合一个常见的非线性模型——指数衰减模型：y = a * exp(-b * x) + c其中，a,b,c是待求参数。我们有了一组模拟的带噪声数据。

3.1 第一步：准备数据与模型函数

首先，我们生成一组模拟数据。这里我们设定真实参数为a=2.5,b=0.8,c=0.5，并加上一点随机噪声。

% 1. 生成模拟数据 clear; clc; rng(2024); % 固定随机种子，确保结果可复现 % 真实参数 a_true = 2.5; b_true = 0.8; c_true = 0.5; % 生成自变量x x_data = linspace(0, 5, 50)'; % 从0到5，50个点，列向量 % 计算无噪声的理论y值 y_true = a_true * exp(-b_true * x_data) + c_true; % 添加高斯白噪声 noise_level = 0.1; y_data = y_true + noise_level * randn(size(x_data)); % 绘制数据点 figure(1); scatter(x_data, y_data, 40, 'b', 'filled', 'DisplayName', '观测数据'); hold on; plot(x_data, y_true, 'r-', 'LineWidth', 2, 'DisplayName', '真实模型'); xlabel('x'); ylabel('y'); legend('Location', 'best'); title('待拟合的指数衰减数据'); grid on;

接下来，定义我们的模型函数和残差函数。在高斯牛顿法中，我们需要计算残差向量和雅可比矩阵。

% 2. 定义模型函数（正向预测） % 输入：参数向量p = [a; b; c]， 自变量x % 输出：模型预测值y model_func = @(p, x) p(1) * exp(-p(2) * x) + p(3); % 3. 定义用于高斯牛顿法的残差函数和雅可比计算函数 % 残差函数：计算所有数据点的残差向量 r = y_data - y_pred % 输入：当前参数p % 输出：残差向量 (m x 1) residual_func = @(p) y_data - model_func(p, x_data); % 雅可比函数：计算残差关于参数的雅可比矩阵J % J的第i行：第i个数据点处，残差对每个参数的偏导: [dr/da, dr/db, dr/dc] % 对于我们的模型： y = a*exp(-b*x) + c % r_i = y_i - (a*exp(-b*x_i) + c) % 因此： % dr/da = -exp(-b*x_i) % dr/db = a * x_i * exp(-b*x_i) （注意：对b求导，指数函数求导会产生一个负号，再与r定义中的负号抵消？这里要小心！） % dr/dc = -1 % 我们来仔细推导一下，避免符号错误，这是算法成败的关键。 % r_i = y_i - f_i, 其中 f_i = a*exp(-b*x_i) + c % 所以 ∂r_i/∂a = -∂f_i/∂a = -exp(-b*x_i) % ∂r_i/∂b = -∂f_i/∂b = -[a * (-x_i) * exp(-b*x_i)] = a * x_i * exp(-b*x_i) % ∂r_i/∂c = -∂f_i/∂c = -1 % 因此，雅可比矩阵J的第i行为：[-exp(-b*x_i), a*x_i*exp(-b*x_i), -1] jacobian_func = @(p) [ -exp(-p(2) * x_data), % 对a的偏导列 p(1) * x_data .* exp(-p(2) * x_data), % 对b的偏导列 -ones(size(x_data)) % 对c的偏导列 ];

3.2 第二步：实现高斯牛顿法迭代核心

现在，我们编写高斯牛顿法的主循环。我们需要设置初始猜测、最大迭代次数、容忍度等。

% 4. 高斯牛顿法实现 % 初始参数猜测（故意给得离真值远一点，增加挑战性） p0 = [1.0; 0.3; 1.0]; % [a; b; c] current_p = p0; % 算法参数 max_iter = 100; % 最大迭代次数 tolerance = 1e-6; % 参数更新量的容忍度 history.p = []; % 记录参数迭代历史 history.cost = []; % 记录损失函数历史 fprintf('开始高斯牛顿法迭代...\n'); fprintf('迭代 | 损失函数值 | 参数更新范数\n'); fprintf('-----|------------|----------------\n'); for iter = 1:max_iter % 计算当前参数下的残差和雅可比 r = residual_func(current_p); J = jacobian_func(current_p); % 计算当前损失（残差平方和） current_cost = sum(r.^2); % 记录历史 history.p = [history.p, current_p]; history.cost = [history.cost, current_cost]; % 构建高斯牛顿方程： (J^T * J) * delta_p = -J^T * r % 在Matlab中，我们可以用反斜杠运算符直接求解线性最小二乘问题，它更稳定。 % 即求解： J * delta_p ≈ -r % 这等价于求解正规方程，但数值上更优。 delta_p = J \ (-r); % 这是求解 min ||J*delta_p + r||^2 % 更新参数 new_p = current_p + delta_p; % 计算参数更新量范数 delta_norm = norm(delta_p); fprintf('%4d | %10.6e | %12.6e\n', iter, current_cost, delta_norm); % 检查收敛条件 if delta_norm < tolerance fprintf('在 %d 次迭代后收敛。\n', iter); break; end % 为下一次迭代准备 current_p = new_p; if iter == max_iter fprintf('达到最大迭代次数 %d，可能未完全收敛。\n', max_iter); end end fitted_p = current_p; fprintf('\n拟合结果：\n'); fprintf('参数 a: 真实值 = %.4f, 拟合值 = %.4f, 误差 = %.4f\n', a_true, fitted_p(1), fitted_p(1)-a_true); fprintf('参数 b: 真实值 = %.4f, 拟合值 = %.4f, 误差 = %.4f\n', b_true, fitted_p(2), fitted_p(2)-b_true); fprintf('参数 c: 真实值 = %.4f, 拟合值 = %.4f, 误差 = %.4f\n', c_true, fitted_p(3), fitted_p(3)-c_true);

3.3 第三步：结果可视化与算法分析

运行完迭代，我们需要看看拟合效果如何，并分析算法的行为。

% 5. 结果可视化 % 绘制拟合曲线 y_fitted = model_func(fitted_p, x_data); figure(2); subplot(2,1,1); scatter(x_data, y_data, 40, 'b', 'filled', 'DisplayName', '观测数据'); hold on; plot(x_data, y_true, 'r-', 'LineWidth', 2, 'DisplayName', '真实模型'); plot(x_data, y_fitted, 'g--', 'LineWidth', 2, 'DisplayName', '高斯牛顿法拟合'); xlabel('x'); ylabel('y'); legend('Location', 'best'); title('模型拟合结果对比'); grid on; % 绘制残差图 residuals = y_data - y_fitted; subplot(2,1,2); scatter(x_data, residuals, 40, 'k', 'filled'); hold on; plot([min(x_data), max(x_data)], [0,0], 'r-', 'LineWidth', 1); % 零线 xlabel('x'); ylabel('残差'); title('拟合残差分布'); grid on; % 绘制损失函数下降曲线 figure(3); plot(1:length(history.cost), history.cost, 'bo-', 'LineWidth', 1.5, 'MarkerFaceColor', 'b'); xlabel('迭代次数'); ylabel('损失函数 (残差平方和)'); title('高斯牛顿法损失函数收敛过程'); set(gca, 'YScale', 'log'); % 使用对数坐标更清晰地观察下降 grid on; % 绘制参数迭代路径（以a和b为例） figure(4); plot(history.p(1,:), history.p(2,:), 'bd-', 'LineWidth', 1.5, 'MarkerFaceColor', 'b', 'MarkerSize', 8); hold on; plot(a_true, b_true, 'rp', 'MarkerSize', 20, 'LineWidth', 3, 'DisplayName', '真实值'); plot(p0(1), p0(2), 'gs', 'MarkerSize', 15, 'LineWidth', 2, 'DisplayName', '初始猜测'); xlabel('参数 a'); ylabel('参数 b'); title('参数空间迭代路径 (a-b平面)'); legend('迭代路径', '真实值', '初始猜测', 'Location', 'best'); grid on;

运行这段完整的代码，你应该能看到算法在几十次迭代内收敛，拟合曲线与真实曲线基本重合，残差随机分布，损失函数单调下降。参数迭代路径图能直观地展示参数是如何从初始猜测点“走”向真实值点的。

4. 关键实现细节与“避坑”指南

自己动手实现一遍后，你会发现有几个细节决定了算法的成败和效率。这些是教科书上不一定强调，但实践中必须注意的。

4.1 雅可比矩阵的计算：解析法 vs 数值法

在上面的例子中，我们幸运地推导出了雅可比矩阵的解析表达式。这提供了最高的精度和计算效率。然而，很多复杂的模型，其导数可能很难甚至无法手动推导。这时就需要使用数值微分来近似雅可比矩阵。

最常见的方法是前向有限差分：J_ij ≈ [r_i(p + ε*e_j) - r_i(p)] / ε其中e_j是第j个分量为1的单位向量，ε是一个很小的正数（如1e-7）。

在Matlab中，你可以这样实现一个通用的数值雅可比计算函数：

function J = numerical_jacobian(residual_func, p, epsilon) % residual_func: 函数句柄，输入参数向量p，输出残差向量r % p: 当前参数向量 (n x 1) % epsilon: 差分步长，可选，默认1e-7 if nargin < 3 epsilon = 1e-7; end n = length(p); % 参数个数 r0 = residual_func(p); % 当前残差 m = length(r0); % 数据点个数 J = zeros(m, n); % 初始化雅可比矩阵 for j = 1:n p_perturbed = p; p_perturbed(j) = p_perturbed(j) + epsilon; r_perturbed = residual_func(p_perturbed); J(:, j) = (r_perturbed - r0) / epsilon; end end

然后在主循环中，将J = jacobian_func(current_p);替换为J = numerical_jacobian(residual_func, current_p);即可。

注意：数值微分有几个坑。第一，步长ε的选择是个权衡：太小会放大舍入误差，太大则截断误差大。通常取ε = sqrt(eps)，其中eps是Matlab的浮点精度。第二，计算量是解析法的n倍（n为参数个数），对于参数多或残差计算昂贵的问题，这可能成为瓶颈。第三，对于具有不连续或剧烈变化的函数，数值微分可能不准确。因此，只要可能，尽量使用解析导数。

4.2 线性方程组的求解：慎用inv(J'*J)

在理论推导中，我们得到了正规方程(J^T J) Δp = -J^T r。一个危险的诱惑是直接计算H = J'*J和g = J'*r，然后求逆得到Δp = -inv(H) * g。千万不要这样做！

原因有二：1) 计算H再求逆，在数值上比直接求解原方程更不稳定，会放大J的条件数（条件数平方）。2) 效率更低。Matlab中的反斜杠运算符\在求解J \ (-r)时，内部会采用QR分解或SVD等稳定的数值方法，自动处理秩亏或病态问题，远比显式求逆来得稳健。所以，记住这个黄金法则：在Matlab中，永远用A \ b来求解线性方程组A*x = b或最小二乘问题，而不是inv(A)*b。

4.3 初始值的选择：决定收敛的“起跑线”

高斯牛顿法对初始值非常敏感。如果初始值离全局最优解太远，很容易收敛到局部极小点，甚至发散。在实践中，选择初始值没有万能公式，但有一些策略：

• 物理意义猜测：如果参数有物理意义（如衰减率、振幅、基线），可以根据对问题的理解给出一个合理的数量级估计。
• 网格搜索：对于2-3个参数的问题，可以在一个合理的范围内进行粗略的网格搜索，选择使损失函数最小的点作为初始值。
• 线性化近似：对于一些可线性化的模型（如指数模型取对数），可以先通过线性回归得到一个粗略估计，再作为非线性拟合的初值。对于我们例子中的y = a*exp(-b*x)+c，当c已知或可估计时，取对数可化为线性问题。
• 随机多起点：从多个随机初始点开始运行算法，选择最终损失函数最小的结果。这在一定程度上可以缓解局部极小问题。

在我们的代码示例中，初始值[1.0; 0.3; 1.0]虽然离真值[2.5; 0.8; 0.5]有距离，但仍在“吸引盆”内，所以能成功收敛。你可以尝试将p0改为[10; 0.1; 5]看看，算法可能就会发散或收敛到一个很差的解。

4.4 收敛判据与迭代控制

除了检查参数增量norm(delta_p)是否小于容忍度，一个更全面的收敛判据应该同时考虑：

1. 参数变化：norm(delta_p) < tol_p
2. 函数值变化：abs(current_cost - previous_cost) < tol_cost
3. 梯度范数：norm(J'*r) < tol_grad（在最优解处梯度应为零）

通常将1和2结合使用。此外，必须设置最大迭代次数max_iter以防止无限循环。在迭代过程中，打印或记录每次迭代的损失函数值和参数更新量（如我们代码中所做），对于调试和监控算法行为至关重要。如果看到损失函数在若干次迭代后不再下降甚至上升，或者参数更新量出现振荡，那可能就是算法遇到问题了。

5. 当高斯牛顿法“失灵”时：诊断与进阶策略

即使你小心翼翼地实现了代码，高斯牛顿法有时还是会“罢工”。常见的症状包括：迭代不收敛（损失函数震荡或爆炸）、收敛速度极慢、或者得到明显错误的解。这时你需要成为一名“算法医生”进行诊断。

5.1 问题一：矩阵J^T J奇异或病态

这是高斯牛顿法最常见的问题。在迭代中，如果J的列向量近似线性相关（即参数之间存在强共线性），或者某个参数在当前点对残差的敏感性极低（雅可比矩阵对应列接近零向量），那么J^T J就会接近奇异矩阵。在Matlab中，用反斜杠求解J \ (-r)时，可能会收到“矩阵接近奇异或缩放错误”的警告，计算出的Δp可能异常巨大，导致参数更新步长爆炸，算法发散。

诊断方法：在每次迭代中，计算雅可比矩阵J的条件数cond(J)。如果条件数非常大（比如大于1e10），就说明矩阵病态。

解决方案：

• Levenberg-Marquardt 方法：如前所述，这是最直接有效的改进。通过添加阻尼项λI，确保系数矩阵总是正定的。Matlab优化工具箱中的lsqnonlin函数默认使用L-M算法或其变种。你可以自己实现一个简单的L-M：在求解delta_p时，不是用J \ (-r)，而是用(J'*J + lambda*eye(n)) \ (-J'*r)，并根据本次更新是否降低了损失函数来动态调整lambda（降低则减小lambda，接受更新；升高则增大lambda，拒绝更新并重试）。
• 参数缩放：如果各个参数的数量级相差巨大（例如a约等于1000，而b约等于0.001），这本身就会导致雅可比矩阵的列尺度差异大，从而引起病态。可以对参数进行缩放，令其处于同一数量级（例如令b' = 1000 * b），在优化完成后，再将结果缩放回来。
• 奇异值分解（SVD）：对于病态方程，可以使用SVD来求解最小二乘问题，并可以截断小的奇异值来获得一个稳定的解（正则化）。Matlab中可以用[U,S,V] = svd(J, 'econ');然后手动处理奇异值。

5.2 问题二：残差函数或模型函数存在数值问题

有时问题不出在算法，而出在模型本身。例如，在计算exp(-b*x)时，如果b*x很大，可能导致下溢（结果为零）；如果b为负且x很大，可能导致上溢（结果为无穷大Inf）。这些都会导致雅可比矩阵计算出错。

诊断方法：在残差函数和雅可比函数中加入调试输出，检查是否有NaN或Inf值出现。使用isfinite()函数进行检查。

解决方案：

• 参数约束：如果参数有物理意义（如衰减率b应为正数），可以在迭代过程中加入约束。简单的方法是在参数更新后，将其截断到合理范围内（例如b = max(b, 1e-6)）。更严谨的方法是使用带约束的优化算法，如Matlab的lsqnonlin可以设置lb和ub参数。
• 模型重参数化：有时改变参数的表达方式可以改善数值稳定性。例如，对于指数衰减，如果担心b为负，可以令b = exp(θ)，然后优化θ，这样无论θ取何值，b总是正的。
• 稳健的代码实现：在计算指数、对数等函数时，考虑使用expm1,log1p等函数来提高小参数值附近的精度。

5.3 问题三：陷入局部极小点

非线性最小二乘问题通常是非凸的，可能存在多个局部极小点。高斯牛顿法只能保证收敛到初始点附近的局部极小点，而不一定是全局最小点。

诊断方法：从多个不同的、分散的初始点运行算法。如果总是收敛到同一个点，那很可能是全局最优。如果从不同起点收敛到不同的点，且损失函数值差异较大，那就存在局部极小问题。

解决方案：

• 多起点优化：如前所述，这是最实用的方法。结合随机初始化和确定性网格搜索。
• 全局优化算法：对于困难的全局优化问题，可以考虑使用模拟退火、遗传算法、粒子群优化等全局搜索方法先进行粗搜索，将其结果作为高斯牛顿法的初始值进行精炼。Matlab的全局优化工具箱提供了这些功能。
• 修改损失函数：有时使用更稳健的损失函数（如Huber损失、Cauchy损失）代替平方损失，可以减少对异常值的敏感性，并可能改变优化问题的景观，使得全局最优点更容易被找到。但这本质上已经不再是标准的高斯牛顿法了。

6. 与Matlab内置函数的对比与实践建议

我们费了很大功夫自己实现了高斯牛顿法，但在实际工作中，我们更常使用Matlab优化工具箱中成熟的函数，比如lsqnonlin或lsqcurvefit。了解它们与我们的手写实现有何异同，以及何时该用哪个，非常重要。

6.1 使用lsqcurvefit一键拟合

对于我们的例子，用lsqcurvefit可以极其简洁地完成：

% 定义模型函数 (注意：lsqcurvefit要求函数形式为 f(p, x)，输出预测值y) model_for_lsq = @(p, x) p(1) * exp(-p(2) * x) + p(3); % 设置选项：采用Levenberg-Marquardt算法，显示迭代过程 options = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', 'Display', 'iter'); % 调用lsqcurvefit [p_fit_lsq, resnorm, residual, exitflag, output] = lsqcurvefit(model_for_lsq, p0, x_data, y_data, [], [], options); fprintf('\nlsqcurvefit 拟合结果：\n'); disp(p_fit_lsq); fprintf('残差平方和： %.6e\n', resnorm); fprintf('迭代次数： %d\n', output.iterations);

lsqcurvefit内部使用的就是L-M算法，它自动处理了阻尼因子的调整、雅可比矩阵的数值估计（如果你不提供的话）、以及各种收敛判断。它比我们写的基础版高斯牛顿法要稳健得多。

6.2 使用lsqnonlin处理更一般的残差形式

如果你的问题不是标准的曲线拟合（yvsx），而是更一般的最小化残差平方和问题，lsqnonlin更合适。它要求你提供一个返回残差向量的函数。

% 定义残差函数 (输入参数p，输出残差向量) residual_for_lsq = @(p) y_data - (p(1) * exp(-p(2) * x_data) + p(3)); options = optimoptions('lsqnonlin', 'Algorithm', 'levenberg-marquardt', 'Display', 'iter'); [p_fit_nonlin, resnorm, residual, exitflag, output] = lsqnonlin(residual_for_lsq, p0, [], [], options);

6.3 手写实现 vs 内置函数：如何选择？

• 选择手写实现的情况：

  1. 教学与理解：为了彻底理解高斯牛顿/L-M算法的每一个步骤，自己实现一遍是无价的。
  2. 高度定制化需求：内置函数可能无法满足某些特殊需求，比如你需要使用特定的线性方程组求解器、嵌入特殊的参数约束逻辑、或者与自定义的数值模拟代码深度耦合。
  3. 轻量级依赖：在不方便安装优化工具箱的环境下，一个自己写的、功能明确的简单实现可能更合适。
  4. 性能极限调优：对于超大规模或特定结构的问题，你可能有比内置函数更高效的雅可比矩阵计算或线性求解方法。

• 选择内置函数的情况：

  1. 生产环境与科研：对于绝大多数实际问题，内置函数lsqcurvefit或lsqnonlin是首选。它们经过严格测试，鲁棒性强，功能丰富（支持多种算法、边界约束、提供详细的输出信息）。
  2. 快速原型：当你需要快速验证一个想法时，内置函数能让你在几分钟内完成拟合。
  3. 避免重复造轮子：内置函数已经处理了数值稳定性、算法切换、收敛判断等复杂细节，比自己从头实现更可靠。

我个人在数学建模竞赛或快速分析数据时，几乎总是先用lsqcurvefit尝试。只有当它出现问题（比如不收敛、结果不合理），或者我需要向学生/队友解释算法原理时，才会去深入底层，自己编写优化循环。自己实现的最大收获，是当内置函数给出警告或错误时，你能明白背后可能的原因，并知道如何去调整选项或修改问题表述。

最后，无论用哪种方法，可视化都是不可或缺的一环。永远要绘制拟合曲线与原始数据的对比图、残差图、以及参数收敛过程图。图形能最直观地告诉你拟合得好不好，算法行为是否正常，这是任何数值指标都无法替代的。