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)这个方程被称为高斯牛顿方程。解出Δβ后我们就更新参数β_{k1} β_k Δβ然后用新的β_{k1}计算新的残差和雅可比矩阵重复上述过程直到满足停止条件例如Δβ的范数很小或者目标函数S(β)下降不明显。注意这里有一个关键的细节。雅可比矩阵J的元素是J_ij ∂r_i/∂β_j -∂f(x_i, β)/∂β_j。在实际编程中我们需要提供这个导数信息。有时我们可以推导出解析的导数公式这样最精确高效如果不行也可以用数值差分如有限差分来近似但这会引入误差并增加计算量。2.3 优势与天生的缺陷高斯牛顿法的优势很明显它通常比最速下降法收敛快得多。因为它利用了目标函数平方和的二阶信息近似海森矩阵J^T J迭代方向更智能。在许多问题中它表现出接近二阶收敛的速度。但它也有两个著名的“阿喀琉斯之踵”初始值依赖性强因为它基于局部线性近似如果初始猜测β_0离真实解太远线性近似可能非常糟糕导致算法收敛到错误的局部极小点甚至发散。矩阵病态问题方程中的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 第一步准备数据与模型函数首先我们生成一组模拟数据。这里我们设定真实参数为a2.5,b0.8,c0.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到550个点列向量 % 计算无噪声的理论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)是否小于容忍度一个更全面的收敛判据应该同时考虑参数变化norm(delta_p) tol_p函数值变化abs(current_cost - previous_cost) tol_cost梯度范数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 内置函数如何选择选择手写实现的情况教学与理解为了彻底理解高斯牛顿/L-M算法的每一个步骤自己实现一遍是无价的。高度定制化需求内置函数可能无法满足某些特殊需求比如你需要使用特定的线性方程组求解器、嵌入特殊的参数约束逻辑、或者与自定义的数值模拟代码深度耦合。轻量级依赖在不方便安装优化工具箱的环境下一个自己写的、功能明确的简单实现可能更合适。性能极限调优对于超大规模或特定结构的问题你可能有比内置函数更高效的雅可比矩阵计算或线性求解方法。选择内置函数的情况生产环境与科研对于绝大多数实际问题内置函数lsqcurvefit或lsqnonlin是首选。它们经过严格测试鲁棒性强功能丰富支持多种算法、边界约束、提供详细的输出信息。快速原型当你需要快速验证一个想法时内置函数能让你在几分钟内完成拟合。避免重复造轮子内置函数已经处理了数值稳定性、算法切换、收敛判断等复杂细节比自己从头实现更可靠。我个人在数学建模竞赛或快速分析数据时几乎总是先用lsqcurvefit尝试。只有当它出现问题比如不收敛、结果不合理或者我需要向学生/队友解释算法原理时才会去深入底层自己编写优化循环。自己实现的最大收获是当内置函数给出警告或错误时你能明白背后可能的原因并知道如何去调整选项或修改问题表述。最后无论用哪种方法可视化都是不可或缺的一环。永远要绘制拟合曲线与原始数据的对比图、残差图、以及参数收敛过程图。图形能最直观地告诉你拟合得好不好算法行为是否正常这是任何数值指标都无法替代的。