非线性最小二乘（Levenberg-Marquardt）

一、直观理解：参数“弯着”进模型时怎么拟合

当参数以非线性方式进入模型时，例如

$$y=a\,e^{-bx}+c,\quad y=\frac{A}{1+e^{-k(x-x_0)}}$$

残差平方和不再是 $\beta$ 的二次函数，法方程不能一步解出。非线性最小二乘通过迭代：从初值出发，反复线性化模型、更新参数，直到残差足够小。

Levenberg-Marquardt（L-M）算法在高斯—牛顿法与梯度下降之间自适应插值：信任靠近解时采用近似二阶（快），远离时则更保守（稳）。它是曲线拟合中最常用的工作马之一。

二、适用场景与优缺点

类型	说明
适用场景	机理模型参数非线性；指数、幂律、Logistic、酶动力学等；需要比黑箱机器学习更可解释的参数。
核心优点	能拟合真正的非线性结构；L-M 对中等规模问题稳健高效；可结合边界约束；参数有物理意义。
局限与缺点	依赖初值；可能收敛到局部极小；参数不可识别时会漂移；需要雅可比或差商；噪声大时要正则/约束。

三、核心原理：高斯—牛顿与阻尼

目标

$$\min_{\boldsymbol\theta} S(\boldsymbol\theta)=\frac12\sum_{i=1}^n r_i(\boldsymbol\theta)^2,\quad r_i=y_i-f(x_i;\boldsymbol\theta)$$

令雅可比 $J_{ij}=\partial r_i/\partial\theta_j$。高斯—牛顿步满足

$$(J^\mathsf{T}J)\delta=-J^\mathsf{T}r$$

L-M 引入阻尼 $\lambda>0$：

$$(J^\mathsf{T}J+\lambda D)\delta=-J^\mathsf{T}r$$

其中 $D$ 常取 $\mathrm{diag}(J^\mathsf{T}J)$ 或单位阵。若试探步使 $S$ 下降，则接受并减小 $\lambda$；否则增大 $\lambda$ 重试。

收敛与识别

• 多组 $\boldsymbol\theta$ 可能给出相近曲线（不可识别）；
• 应检查参数置信区间、轮廓 SSE、相关矩阵；
• 多初值重复优化，报告最佳与稳定性。

四、关键实践决策

1. 模型形式必须有机理或强经验支持；
2. 初值用线性化近似、图形估读或网格搜索；
3. 缩放参数使量级接近，改善条件数；
4. 边界约束（lsqcurvefit 的 lb/ub）避免非物理值；
5. 权重：异方差时用加权最小二乘；
6. 验证：残差图、样本外点、极限行为（$x\to\infty$）。

五、MATLAB 示例

%% 指数衰减 + 常数：y = a*exp(-b*x) + c
rng(6);
xdata = linspace(0, 5, 40)';
aTrue = 3; bTrue = 1.2; cTrue = 0.5;
ydata = aTrue*exp(-bTrue*xdata) + cTrue + 0.08*randn(size(xdata));

model = @(theta, x) theta(1).*exp(-theta(2).*x) + theta(3);
theta0 = [1, 1, 0];   % 初值

opts = optimoptions("lsqcurvefit", ...
    "Algorithm", "levenberg-marquardt", ...
    "Display", "off");

% L-M 在 lsqcurvefit 中不支持边界；无边界时可用
thetaHat = lsqcurvefit(model, theta0, xdata, ydata, [], [], opts);

yhat = model(thetaHat, xdata);
RMSE = sqrt(mean((ydata - yhat).^2));

fprintf("估计 a,b,c = %.4f, %.4f, %.4f\n", thetaHat);
fprintf("真值 a,b,c = %.4f, %.4f, %.4f\n", [aTrue bTrue cTrue]);
fprintf("RMSE = %.4f\n", RMSE);

xq = linspace(0, 5, 200)';
figure;
scatter(xdata, ydata, 40, "filled"); hold on;
plot(xq, model(thetaHat, xq), "r-", "LineWidth", 1.5);
plot(xq, model([aTrue bTrue cTrue], xq), "k--");
xlabel("x"); ylabel("y");
legend("数据", "L-M 拟合", "真值曲线", "Location", "best");
title("Levenberg-Marquardt 非线性拟合");
grid on;

Logistic 曲线与多初值：

%% Logistic: y = A / (1+exp(-k*(x-x0)))
rng(9);
x = linspace(-8, 8, 35)';
thetaTrue = [2.5, 1.0, 0.5];  % A,k,x0
f = @(th,x) th(1)./(1+exp(-th(2)*(x-th(3))));
y = f(thetaTrue, x) + 0.05*randn(size(x));

starts = [
    1, 0.5, 0;
    3, 2.0, 1;
    2, 0.2, -1
];
opts = optimoptions("lsqcurvefit", "Algorithm","levenberg-marquardt", "Display","off");
bestS = inf; bestTh = starts(1,:);
for i = 1:size(starts,1)
    th = lsqcurvefit(f, starts(i,:), x, y, [], [], opts);
    S = sum((y - f(th,x)).^2);
    if S < bestS
        bestS = S; bestTh = th;
    end
end
fprintf("多初值最佳参数 = [%.3f %.3f %.3f], SSE = %.4f\n", bestTh, bestS);

若必须加边界，可改用 trust-region-reflective 算法：

optsTR = optimoptions("lsqcurvefit", "Display","off");
lb = [0, 0, -10]; ub = [10, 10, 10];
thB = lsqcurvefit(f, [1 1 0], x, y, lb, ub, optsTR);

六、如何解释结果

应包含：

1. 非线性模型公式与参数物理意义；
2. 优化算法（L-M）、初值策略与是否多起点；
3. 参数估计值及不确定度（可用残差方差与近似协方差 $\sigma^2(J^\mathsf{T}J)^{-1}$）；
4. RMSE/残差图；
5. 与线性化近似或多项式基线的比较；
6. 参数相关性：若两参数高度相关，解释要谨慎。

七、常见误区

1. 初值随意导致“收敛成功但参数荒谬”

始终检查曲线是否贴合、参数是否在物理范围。

2. 模型不可识别仍强行拟合

如 $ae^{-(b+c)x}$ 中 $b,c$ 不能分开，需要重新参数化。

3. 把迭代收敛当成全局最优

非线性最小二乘是局部方法，需多初值。

4. 忽视权重与异方差

方差随 $y$ 变化时，普通 SSE 可能偏向高方差区。

5. 混淆插值与拟合

非线性最小二乘不要求过点，过点应使用插值或增参到插值极限（不推荐）。

附：建模流程中的落地清单（非线性最小二乘）

要把方法真正写进论文，建议补齐以下模块，而不是只贴一个最终数字。

1. 问题对齐与假设

说明题目中的输入/输出如何映射到本方法，并显式列出关键假设。假设越清楚，后面的敏感性分析越有对象。

2. 关键关系

最小化 $\sum r_i(\theta)^2$，L-M 在 $J^\top J+\lambda D$ 上求步长，兼顾 Gauss-Newton 与梯度下降。

3. 补充实现片段

x=linspace(0,3,40)'; y=2*exp(-1.3*x)+0.4+0.05*randn(size(x));
model=@(t,x)t(1)*exp(-t(2)*x)+t(3); th=lsqcurvefit(model,[1 1 0],x,y); disp(th);

4. 报告要素

模型公式、初值策略、参数估计、RMSE、多起点结果。

5. 不适用情形

不可识别参数化；单起点收敛就当全局最优。

6. 与验证章节的衔接

至少做一项：参数扰动、与基线方法对照、极限情形检查。把“方法正文”和“全局方法论（第15门类）”连起来，结论才站得住。

八、小结

非线性最小二乘处理参数非线性的机理型曲线拟合，Levenberg-Marquardt 通过阻尼在速度与稳定性之间取得平衡。成功的关键是：模型可识别、初值有依据、诊断要充分、不确定度要报告。在数学建模中，它往往比黑箱高阶多项式更能讲清“参数意味着什么”。