简介：基于Levenberg-Marquardt算法的非线性最小二乘拟合Matlab实现，面向需要处理曲线拟合、参数估计与数据建模的科研与工程人员。该压缩包共5个文件，包括3个.m源文件与2张拟合效果图：源码覆盖LM核心迭代逻辑、角度相关数据示例及主调用脚本，图片直观展示拟合前后的残差收敛情况，便于对照调试。已有1320人学习，具有较强的实用参考价值。文件中通过angleLoveForYou示例演示了带角度约束数据的非线性拟合流程，帮助理解LM算法在梯度下降法和高斯-牛顿法之间的自适应调节机制，以及阻尼参数对收敛速度与稳定性的影响。压缩包仅46KB，轻量易读，适合正在学习最优化理论、数值算法或Matlab仿真的中高级用户，既能作为教学辅助材料，也可扩展用于物理、化学、生物统计等领域的自定义模型拟合任务。

1. LM算法与非线性最小二乘：先回答“拟合效果是什么意思”

数据拿回来，第一件事永远是拟合。海浪谱、潮位时序、传感器校准曲线，说到底都是给一组散点找一个函数，然后让残差尽可能小。很多人开口就问“拟合效果是什么意思”，其实就是问残差压下去了没有、参数靠不靠谱。这个问题的标准解法是非线性最小二乘，而工程上最常用的求解器就是LM算法，也叫Levenberg-Marquardt方法。它不需要你保证函数可导到三阶，不需要手工调学习率，只要给出残差表达式和初值，就能在几步迭代内收敛到局部最优点。这篇围绕LM算法展开：先讲清它的内部机制，再用手写代码和海洋水文波浪拟合函数两个案例把它跑通，最后把参数设置和踩坑记录列出来。适合刚接触数据拟合的工程师，也适合被初值问题、局部极值折磨过一阵子的老手。

2. 拆解LM算法的三层内部逻辑：阻尼因子、雅可比矩阵与步长折中

2.1 从梯度下降到高斯牛顿：为什么工程上最终常用LM

非线性最小二乘的目标函数是残差平方和，残差向量记作 (r(x))，目标就是让 (F(x)=\frac{1}{2}|r(x)|^2) 尽量小。梯度下降的思路最简单：沿着负梯度方向走一步，(x_{k+1}=x_k-\alpha J^T r)，但这只用了梯度信息，没有考虑曲率，在病态问题上经常走一步晃三下，收敛慢得让人着急。

高斯牛顿法向前走了一步，把目标函数的二阶近似拉进来，用 (J^T J \delta=-J^T r) 求步长 (\delta)。这里的思路是把残差函数线性化，再用线性最小二乘解出最优步长，所以在残差接近零时收敛速度非常快，二阶收敛不是说着玩的。但代价是它太自信了：一旦初始点离最优解远，或者 (J^T J) 近似得不好，一步就可能迈到十万八千里之外，目标函数不降反升，直接翻车。

LM算法的高明之处在于把两者的优势接在一起，在正规方程里加了一个阻尼项：

[ (J^T J + \lambda D)\delta = -J^T r ]

其中 (D) 通常取 (\mathrm{diag}(J^T J))，是为了让不同参数的量纲差异不至于把方程搞病态。当 (\lambda) 很小时，方程趋近于高斯牛顿，迭代后期快速收敛；当 (\lambda) 很大时，步长方向会被强行掰向负梯度方向，等于退回最速下降。也就是说，LM在迭代过程中自动在“大胆逼近”和“小心试探”之间切换，这就是它作为非线性最小二乘默认求解器的根本原因。

2.2 λ阻尼因子如何控制步长：增益比ρ的调节逻辑

阻尼因子 (\lambda) 不是拍脑袋设一个固定值就完事，它需要在每轮迭代里被动态调整。调整的依据是“实际下降量”和“模型预测下降量”的比值，工程上叫增益比 (\rho)：

[ \rho = \frac{F(x_k)-F(x_{k+1})}{F(x_k)-L(\delta)} ]

分子是目标函数真实的减少量，分母是根据线性化模型推算出来的预测减少量。如果 (\rho) 接近1，说明二阶模型把地形描述得很准，那就可以放心地把 (\lambda) 调小，让算法往高斯牛顿方向靠，加快收敛；如果 (\rho) 明显小于某个阈值，说明模型预测和实际结果对不上，这步走得不靠谱，算法会把 (\lambda) 调大，让下一步更保守。

我常用的阈值是一组经验值：(\rho > 0.75) 就把 (\lambda) 乘以0.3，(\rho < 0.25) 就把 (\lambda) 乘以2到10，收到0.0001到1e10之间。初始 (\lambda) 一般取1e-3到1之间，视残差尺度而定。这个机制的优点是不需要人盯着迭代过程，算法自己会判断该激进还是该保守，这也是LM少有的几个“不用玄学调参”的环节之一。

对应到落地工具里，(xtol)、(ftol)、(max_nfev) 这些参数控制的是收敛判断和迭代上限，而不是阻尼因子本身。阻尼因子的变化细节被封装在求解器内部了。你要理解的是：如果看到迭代很快停在半路，多半是 (xtol) 设得太松，或者初值差到连LM的自动调节都救不回来。

2.3 数值雅可比与解析雅可比：精度、速度与实现成本

LM每轮迭代都要用到雅可比矩阵 (J)，它描述的是每个残差分量对每个参数的偏导数。如果残差函数足够简单，解析雅可比自然是最好的选择，精度高、计算快，还不用担心差分步长选不好。但在真实工程里，模型经常带分支条件，波浪谱的峰形因子就是典型例子：频率小于等于峰值频率时 (\sigma=0.07)，大于峰值频率时 (\sigma=0.09)，这就在峰值频率处形成了一个分段点。手推解析导数容易漏掉分段点处的贡献，写出来一堆代码还不敢保证正确。

常规做法是用中心差分近似数值雅可比：

def numerical_jacobian(residual, x, r, eps=1e-6): n = len(x) m = len(r) J = np.zeros((m, n)) for i in range(n): x_plus = x.copy() x_minus = x.copy() x_plus[i] += eps x_minus[i] -= eps J[:, i] = (residual(x_plus) - residual(x_minus)) / (2 * eps) return J

中心差分比前向差分多算一次残差函数，但截断误差小一个数量级，在LM这种对梯度质量敏感的算法里很值得。当参数个数不超过10个时，我不建议手推解析雅可比，写错一次排错的时间够跑几百次数值差分。如果参数超过几十个，或者性能有硬性要求，再考虑用自动微分或者手推稀疏雅可比。

3. 手写一个LM拟合器：求解骨架与指数衰减参数估计

3.1 残差函数与雅可比矩阵：先定义问题再写迭代

用LM做拟合，第一步不是调库，而是把残差函数写清楚。所谓残差函数，就是“观测值减去模型预测值”组成的向量。以指数衰减模型 (y=a e^{-b x}+c) 为例，残差函数长这样：

import numpy as np def model(theta, x): a, b, c = theta return a * np.exp(-b * x) + c def residual(theta, x, y): return model(theta, x) - y

注意这里返回的是向量而不是标量，LM的核心运算 (J^T r) 和 (J^T J) 都依赖这个向量结构。如果只返回总误差，那就丢掉了每个数据点的信息，没法构造雅可比矩阵。所以任何LM拟合脚本的第一步，永远是先把“残差向量”定义出来，这一步写不对，后面再漂亮的代码也是空转。

3.2 正规方程的稳定求解：QR分解、SVD与直接求逆的差别

每一轮迭代都要解一个线性方程组 ((J^T J+\lambda D)\delta=-J^T r)。最容易踩坑的写法是直接求逆：

# 不推荐：数值上容易丢精度 delta = -np.linalg.inv(A_reg) @ g

直接求逆的问题在于，它把矩阵求逆过程中的舍入误差放大了。尤其当参数量级差异很大时，(J^T J) 的条件数可能达到1e12以上，求逆结果可能已经完全失真。更稳妥的办法是用np.linalg.solve：

delta = np.linalg.solve(A_reg, -g)

np.linalg.solve走的是LAPACK的LU分解，比显式求逆稳定得多，而且在多数场景下速度差距可以忽略。如果要进一步保险，可以把正规方程改写成增广最小二乘形式：把 (\sqrt{\lambda D^{1/2}}) 作为额外的行拼到 (J) 下面，用np.linalg.lstsq求解。这一步相当于用QR分解或SVD直接处理带正则化的最小二乘问题，避免先构造 (J^T J) 再求逆带来的平方精度损失。

我用下来的感受是：对小规模拟合问题，np.linalg.solve已经够用，但如果你发现同一组数据换一种初值就报“奇异矩阵”，不妨改成np.linalg.lstsq版本，多数情况下能多撑一阵子。

3.3 完整实现：用NumPy跑通一个带噪声的指数衰减拟合

把上面的思路拼起来，一轮完整的LM拟合器可以写得很紧凑。下面是我常用的骨架实现，包含数值雅可比、阻尼因子更新和收敛判断：

def lm_fit(residual_func, x0, max_iter=200, xtol=1e-8, ftol=1e-8, lam0=1e-3): x = np.array(x0, dtype=float) lam = lam0 r = residual_func(x) cost = 0.5 * np.dot(r, r) for _ in range(max_iter): J = numerical_jacobian(residual_func, x, r, eps=1e-7) A = J.T @ J g = J.T @ r D = np.diag(np.diag(A)) # 只保留对角，避免量纲差异 A_reg = A + lam * D delta = np.linalg.solve(A_reg, -g) # 预测下降量：-0.5 * g^T delta pred = -0.5 * np.dot(g, delta) if pred <= 0: lam *= 10 continue x_new = x + delta r_new = residual_func(x_new) cost_new = 0.5 * np.dot(r_new, r_new) rho = (cost - cost_new) / pred if rho > 1e-4: # 这一步有效，接受新点并减小阻尼 x, r, cost = x_new, r_new, cost_new lam = max(lam * 0.3, 1e-12) # 收敛判断：步长足够小 if np.linalg.norm(delta) < xtol * (np.linalg.norm(x) + xtol): break else: # 这一步无效，增大阻尼，重新尝试 lam = min(lam * 10, 1e12) return x, cost

这段代码每一行都有具体的工程含义。pred <= 0说明当前步长方向根本不会让线性模型下降，这时候直接放大阻尼换方向；增益比rho > 1e-4表示实际下降虽然小，但不是系统误差导致的，可以接受；阻尼增大上限设到1e12，是为了防止数值溢出，同时也是给算法一个退出机制——如果阻尼大到上限仍然无法下降，基本可以判断问题出在初值或者模型本身。

用一组伪造的带噪声数据来验证：

np.random.seed(0) x_data = np.linspace(0, 5, 80) y_true = 3.0 * np.exp(-1.2 * x_data) + 0.8 y_data = y_true + np.random.normal(0, 0.05, x_data.size) theta0 = np.array([2.0, 1.0, 0.5]) theta_fit, final_cost = lm_fit( lambda th: residual(th, x_data, y_data), theta0 ) print(theta_fit) # 大约输出 [3.01, 1.21, 0.79]，和真实参数很接近

这里初值故意给得偏离真实值，但LM仍然在几十轮迭代内收敛到了正确区域。这就是它在工程里的价值：不需要太精确的初值，只要别偏到另一个吸引域里，它就能自己摸到谷底。注意每迭代一次，数值雅可比要算 (2n) 次残差函数，对80个数据点、3个参数来说毫无压力，如果参数上升到20个以上，性能才需要认真考虑。

4. 海洋水文里的波浪拟合函数：用LM给JONSWAP谱估参数

4.1 把JONSWAP谱写成残差函数：从谱密度公式到可优化代码

海洋水文里的波浪拟合函数，最典型的就是JONSWAP谱。它描述的是充分成长的风浪谱密度随频率的分布，表达式是：

[ S(f)=\frac{\alpha g^2}{(2\pi)^4 f^5}\exp\left[-1.25\left(\frac{f}{f_p}\right)^{-4}\right]\gamma^{\exp\left[-\frac{(f-f_p)^2}{2\sigma^2 f_p^2}\right]} ]

其中 (\alpha) 是能量尺度参数，(f_p) 是谱峰频率，(\gamma) 是峰升因子，(\sigma) 在 (f\le f_p) 时取0.07，否则取0.09。拟合目标是根据实测谱密度序列估计 (\alpha)、(f_p)、(\gamma) 这三个非线性参数。写成残差函数时有一个关键技巧：直接对谱密度做残差，低频段的大能量值会主导整个拟合过程，高频段的细节全被丢掉。我一般会对谱密度做对数变换后再计算残差：

def jonswap_spectrum(f, alpha, fp, gamma): sigma = np.where(f <= fp, 0.07, 0.09) exponent = np.exp(-(f - fp) ** 2 / (2 * sigma**2 * fp**2)) base = alpha * 9.81**2 / (2 * np.pi**4 * f**5) fetch = np.exp(-1.25 * (f / fp) ** (-4)) return base * fetch * gamma ** exponent def jonswap_residual(theta, f, S_obs): model = jonswap_spectrum(f, *theta) return np.log1p(S_obs) - np.log1p(model)

这里np.log1p作用是在对数变换前加1，避免谱密度接近0时对数发散。对数残差本质上给低频和高频分配了更均衡的权重。波浪谱跨好几个数量级，不这么做，LM的雅可比矩阵里高频段的列几乎全是0，拟合结果必然被低频主导。

4.2 初值从哪来：峰值频率、谱峰参量的工程估算

LM算法不是万能的，初值给得太离谱，阻尼因子再会调节也救不回来。JONSWAP谱三个参数里，(f_p) 是最敏感的非线性参数，初值偏差超过20%，迭代就可能跑到别的峰上去。工程上最简单可靠的初值策略，是从实测谱的峰值频率和峰值大小反推：

def initial_guess(f, S_obs): idx = np.argmax(S_obs) fp0 = f[idx] peak = S_obs[idx] # 忽略gamma项，从峰值处反推alpha的近似值 alpha0 = peak * (2 * np.pi**4 * fp0**5) / 9.81**2 * np.exp(1.25) gamma0 = 1.0 return np.array([alpha0, fp0, gamma0])

这个初值估算的思路是：在峰值频率处，指数项 (e^{-1.25}) 相对固定，gamma项贡献未知就先忽略，反算出来的 (\alpha) 初始值可能偏大，但LM一般能在前几次迭代里快速修正。(\gamma) 初值取1.0是稳妥做法，因为 (\gamma=1) 时JONSWAP谱退化为Pierson-Moskowitz谱，相当于从一个已知的基准形态开始调整。

如果实测谱峰值不尖锐，低频还有涌浪能量干扰，np.argmax找出来的可能不是风浪主峰。这种时候我更建议先对谱做滑动平均平滑，再找峰值。峰找错了，后续拟合全部白费。

4.3 拟合效果量化：RMSE、R²与残差自相关三个指标怎么用

有了初值和残差函数，调库求解就是几行代码的事：

from scipy.optimize import least_squares theta0 = initial_guess(f, S_obs) res = least_squares( jonswap_residual, theta0, method='lm', xtol=1e-10, ftol=1e-10, max_nfev=2000, args=(f, S_obs) ) S_fit = jonswap_spectrum(f, *res.x) rmse = np.sqrt(np.mean((S_fit - S_obs) ** 2)) ss_res = np.sum((S_obs - S_fit) ** 2) ss_tot = np.sum((S_obs - S_obs.mean()) ** 2) r2 = 1 - ss_res / ss_tot print(f"RMSE={rmse:.4f}, R2={r2:.4f}")

这里method='lm'直接调用Levenberg-Marquardt算法，xtol和ftol设到1e-10是为了让收敛更充分，max_nfev=2000是安全阀，防止某些难缠的初值导致无限迭代。

拿到RMSE和R²之后，再回答“拟合效果是什么意思”这个问题：R²接近1只能说明曲线形状大体对上了，不能说明模型已经完全捕捉了数据中的结构。真正要看的还有残差自相关。如果拟合残差在相邻频点之间呈现系统性相关，说明模型缺了一个峰或者某个参数被硬凑出来了。我习惯顺手算一个一阶自相关系数：

resid = jonswap_residual(res.x, f, S_obs) r_acf = np.corrcoef(resid[1:], resid[:-1])[0, 1]

如果这个系数绝对值超过0.5，我就会怀疑单峰JONSWAP谱本身不足以描述当前数据，而不是继续调参数。这个判断习惯帮我少走了很多弯路。

5. LM拟合避坑指南：初值敏感、矩阵奇异与“假收敛”的排查

5.1 现象：换一组初值，参数估计结果差一个量级

这是LM拟合里最常见的翻车现场。同一组波浪谱数据，(f_p) 初值取0.10拟合结果是0.15，取0.12拟合结果却是0.12，两个结果都能收敛，但只有后者落在真实值附近。

原因不是LM算法坏了，而是非线性最小二乘的目标函数非凸，存在多个局部极小值。JONSWAP谱的参数 (\alpha)、(\gamma) 之间存在强耦合：把 (\gamma) 调大可以把峰拉高，把 (\alpha) 调大也可以整体抬升谱密度，两个方向都能降低残差，于是目标函数等高线图里出现了脊线。

解决思路就一条：多起点搜索。把 (f_p) 在可能范围内均匀取5到10个初值，每个初值跑一次LM，然后对比最终残差。LM收敛快，跑10次也就几百毫秒，这点成本远低于事后发现结果不合理返工的成本。另一个实用技巧是把拟合分两阶段：第一阶段固定 (\gamma=1.0)，只拟合 (\alpha) 和 (f_p)，拿到稳定结果后再放开 (\gamma) 做精细调整。分步拟合在波浪谱这类多参数强耦合问题里几乎每次都能提高成功率。

5.2 现象：迭代全部收敛，但拟合曲线明显偏离数据

还有一种情况更气人：LM报告收敛成功，残差也在下降，但画出来的拟合曲线明显偏离实测谱的峰值。

原因通常是目标函数里不同频段的权重失衡。如果直接在原始谱密度域做残差，峰值附近的大数值会主导梯度，两侧高频段的细节完全被忽略。LM只能保证让加权后的整体误差最小，但你的加权方式不对，拟合结果就会“狠心”牺牲一部分区域去保另一部分区域。

解决方法是改用对数残差，或者在残差里乘以一个与频率相关的权重向量。JONSWAP谱拟合我基本都是用np.log1p变换，就是为了均衡高低频段的贡献。如果你已经在用对数残差还是出现系统性偏离，检查一下模型本身：单峰JONSWAP谱没法描述双峰谱，这是模型表达能力的问题，不是算法问题。

5.3 现象：提示矩阵奇异或步长不下降，LM卡死

LM在迭代中途报“奇异矩阵”或者步长半天不下降，我遇到过好几次。最典型的场景是三个参数里有两个高度线性相关，比如在某个频段上 (\alpha) 和 (\gamma) 几乎可以互相补偿，(J^T J) 的一个特征值接近0，正规方程解出来步长巨大。

原因分析起来不复杂，(J^T J) 奇异意味着当前参数空间里存在一个方向，沿这个方向移动残差几乎不变。LM的阻尼项 (D=\mathrm{diag}(J^T J)) 只能缓解量纲差异，不能消除本质上的参数冗余。如果 (J^T J) 本身秩亏，阻尼再大也救不了。

解决的第一步是减少参数自由度：固定 (\gamma) 或固定 (\sigma)，先拟合两个参数，再逐层放开。第二步是把正规方程求解改成np.linalg.lstsq，用SVD处理秩亏问题，至少能让迭代继续走，不至于一上来就崩溃。第三步才是调大初始阻尼 (lam0)，这能缓解症状，但没有解决参数冗余这个病根。

5.4 现象：曲线拟合得很好，参数方差却大得离谱

拟合曲线贴着数据走，残差很小，R²高达0.98，但打印出来的参数标准差让人不敢用——某个参数的相对不确定性超过100%。

这是典型的“过参数化”信号。一个参数对残差的影响可以被另一个参数替换，数据本身提供的信息不足以同时确定两者的精确值。JONSWAP谱里的 (\alpha) 和 (\gamma) 就经常掉进这个坑。R²高只能说明模型在给定参数下能很好地复现数据，不能说明参数本身被数据牢牢约束住了。

我在交付拟合结果前，一定会用协方差矩阵估算参数标准差，如果发现参数标准差过大，就考虑固定次要参数或者增加约束。另外一个实用验证方法是bootstrap重采样：把数据按残差分布加扰动，重新拟合20次，看参数分布的离散程度。如果20次拟合结果发散得厉害，说明这组数据对这个参数的约束力不足。

5.5 现象：波浪谱低频段残差出现系统性条带

用单峰JONSWAP谱拟合实测波浪谱，低频段残差经常出现连续的、一正一负的条带结构。第一次看到这种残差图，我以为是LM收敛不到位，把xtol从1e-8调到1e-12，结果条带纹丝不动。

原因非常明确：实测谱里往往同时含有风浪和涌浪两种成分，涌浪能量集中在比风浪更低的频段，频谱呈现双峰甚至多峰形态。单峰JONSWAP模型只能描述主峰，低频的第二个峰自然就变成残差里的系统结构。把这类数据强塞给单峰模型，LM再努力也只是在“折中”。

解决方法是改用双峰谱叠加模型：一个JONSWAP风浪谱加上一个高斯形的涌浪峰，比如 (S(f)=S_{JONSWAP}(f)+A_g\exp(-(f-f_g)^2/(2w_g^2)))，参数从4个增加到7个。LM完全能处理这种扩展模型的拟合，前提是初值要把两个峰的位置都识别出来。如果实测数据本身只有单峰，那就别提前给患者开双峰的药方，低频条带也可能是 (\sigma) 固定值不精确导致的局部失配，先试试把 (\sigma) 也放开拟合。

6. 让LM结果可信：协方差、稳健权重与残差白噪声验证

6.1 参数不确定性：从协方差矩阵读标准差

拟合做完，不能只输出三个参数就交差。LM求解器给出的雅可比矩阵 (J) 可以近似评估参数的不确定性，公式是 (\mathrm{Cov}=s^2(J^T J)^{-1})，其中 (s^2) 是残差方差的估计：

J = res.jac sigma2 = 2 * res.cost / (len(f) - len(theta_fit)) cov = sigma2 * np.linalg.inv(J.T @ J) std = np.sqrt(np.diag(cov)) print(f"fp = {theta_fit[1]:.4f} +/- {std[1]:.4f}")

这个协方差矩阵是在“残差独立同分布”的假设下推导的，数据如果有明显的异方差或自相关，标准差会被低估。所以它的正确用法是相对比较：同一组数据用不同模型拟合，协方差大的模型参数可信度低。遇到标准差超过参数值一半的情况，基本可以判定参数冗余或数据信息不足。

6.2 Huber损失在LM里的落地：压制野值的两种写法

实测谱里偶尔会出现单个频点的强干扰，比如仪器噪声尖峰。普通最小二乘会被这个点拖着走，把整个拟合结果带偏。如果新版SciPy的least_squares支持在LM模式下用稳健损失，可以直接写：

res_robust = least_squares( jonswap_residual, theta0, method='lm', loss='huber', f_scale=0.2, args=(f, S_obs) )

f_scale的含义是“残差超过这个阈值时开始降权”，数值需要按对数残差的实际尺度去试。如果求解器不支持loss参数在LM模式下生效，退路是自己构造加权残差：先跑一轮标准LM得到残差分布，然后用Huber权函数重新定义残差向量再跑一轮。这种做法等于把稳健统计的思想手动接进了LM。注意稳健拟合不是万能药，如果干扰点太多，Huber损失也救不回来，应该先从数据质量入手。

6.3 残差白噪声检验：判断模型是否已经够用

判断拟合是否到头，我最后一步永远是画残差图并检验相关性。视觉上扫一遍是最快的，没有系统性的峰和条带，再算一阶自相关：

resid = jonswap_residual(theta_fit, f, S_obs) r_acf = np.corrcoef(resid[1:], resid[:-1])[0, 1] print(f"lag-1 auto correlation: {r_acf:.3f}")

绝对值小于0.2可以认为是白噪声，模型基本把数据里的可解释结构都吃干净了；大于0.5说明还有结构没被模型覆盖，继续调参数没用，要去升级模型结构。

我以前吃过不看残差的亏：有一组潮位数据拟合R²到了0.99，自以为万事大吉，交付后被现场人员指出特定时段有明显周期性偏差。后来养成了一个习惯，每次LM拟合收敛后必须把残差序列画出来，看有没有颜色、有没有周期性、有没有单点尖峰。这个习惯帮我发现的模型缺陷，比任何自动指标都多。LM是把非线性最小二乘问题做扎实的利器，但它不会替你判断模型本身对不对，最后那一眼残差图，永远是自己的事。希望帮到你。

本文还有配套的精品资源，点击获取