博客：Python 实现 LM 算法（Levenberg-Marquardt）
目录
- 弁言
- 什么是 Levenberg-Marquardt (LM) 算法？
- LM 算法的应用场景
- LM 算法的优点与局限性
- LM 算法的原理
- LM 算法的基本思想
- LM 算法的数学推导
- 与高斯牛顿法和梯度下降法的比较
- Python 实现 LM 算法
- LM 算法应用实例：非线性曲线拟合
- LM 算法的改进与扩展
- LM 算法中的参数调节
- LM 算法的改进与其他变种
- 总结
- LM 算法的实用场景
- 何时选择 LM 算法
- 与其他算法的对比
1. 弁言
什么是 Levenberg-Marquardt (LM) 算法？
Levenberg-Marquardt（LM）算法是求解非线性最小二乘问题的著名方法。它结合了高斯牛顿法和梯度下降法的优点，广泛用于非线性曲线拟合、机器学习、最优参数估计等领域。
LM 算法的一个主要特点是可以大概在问题接近线性时，像高斯牛顿法一样快速收敛；而在问题严峻非线性时，它通过梯度下降法的机制制止了求解发散。它是一个调和局部曲率信息与梯度信息的肴杂算法。
LM 算法的应用场景
- 非线性回归：通过 LM 算法来找到模子参数，使得拟合曲线与数据误差最小。
- 机器学习中的参数优化：在训练神经网络等模子时，LM 算法可以用于权重更新。
- 图像处置惩罚中的形变分析：用 LM 算法来计算形变参数，以最小化图像配准的误差。
LM 算法的优点与局限性
优点：
- 快速收敛：当模子接近线性时，LM 算法会像高斯牛顿法一样快速收敛。
- 制止局部极小：通过引入正则化，LM 算法可以大概更好地处置惩罚复杂的优化问题，制止陷入局部极小值。
局限性：
- 计算复杂度较高：每次迭代都须要计算雅可比矩阵和更新参数，复杂度较高。
- 对初始推测敏感：尽管比牛顿法鲁棒，但如果初始参数推测得过于禁绝确，收敛速度会减慢。
2. LM 算法的原理
LM 算法的基本思想
LM 算法是一种在高斯牛顿法和梯度下降法之间动态调整的算法。当误差函数的二次近似足够好时，算法运动类似高斯牛顿法；当近似不好时，算法则趋向梯度下降法。这通过一个阻尼因子 (\lambda) 来调节，该因子控制了每次迭代时的步长和方向。
LM 算法的数学推导
LM 算法的目标是最小化非线性最小二乘问题：
min θ ∑ i = 1 n ( y i − f ( x i , θ ) ) 2 \min_{\theta} \sum_{i=1}^{n} \left( y_i - f(x_i, \theta) \right)^2 θmini=1∑n(yi−f(xi,θ))2
迭代的更新公式为：
Δ θ = − ( J T J + λ I ) − 1 J T r ( θ ) \Delta \theta = -(J^T J + \lambda I)^{-1} J^T r(\theta) Δθ=−(JTJ+λI)−1JTr(θ)
此中：
- J J J 是误差函数 r ( θ ) r(\theta) r(θ) 的雅可比矩阵；
- λ \lambda λ 是调节参数；
- I I I 是单位矩阵；
- r ( θ ) r(\theta) r(θ) 是残差向量。
λ \lambda λ 决定了算法的特性。当 λ \lambda λ 较大时， λ I \lambda I λI 的影响增大，算法运动趋向于梯度下降法；当 λ \lambda λ 较小时，算法运动接近高斯牛顿法。
与高斯牛顿法和梯度下降法的比较
- 与高斯牛顿法的关系：高斯牛顿法不使用正则化项 λ = 0 \lambda = 0 λ=0，因此当误差模子线性化的假设不建立时，大概会导致收敛迟钝或不稳定。
- 与梯度下降法的关系：梯度下降法通过减少步长来制止发散，而 LM 算法通过动态调节 λ \lambda λ 更灵活地控制步长。
3. Python 实现 LM 算法
面向对象的设计思绪
为了实现 LM 算法，首先须要构建一个代表非线性模子的类，用于计算误差和雅可比矩阵；其次，设计一个 LM 算法的类，用于执行参数优化。通过这种面向对象的方式，我们可以清楚地封装模子和优化算法的功能。
代码实现
- import numpy as np
- class NonlinearModel:
- """表示非线性模型的类，包含残差和Jacobian矩阵的计算。"""
- def __init__(self, func, jacobian):
- """
- :param func: 非线性模型函数
- :param jacobian: Jacobian矩阵的计算函数
- """
- self.func = func
- self.jacobian = jacobian
-
- def residuals(self, x_data, y_data, theta):
- """计算残差向量"""
- return y_data - self.func(x_data, theta)
-
- def jacobian_matrix(self, x_data, theta):
- """计算给定参数下的Jacobian矩阵"""
- return self.jacobian(x_data, theta)
- class LevenbergMarquardt:
- """Levenberg-Marquardt算法的实现类。"""
- def __init__(self, model, tolerance=1e-6, max_iters=100, lambda_init=0.01):
- """
- :param model: 待拟合的非线性模型对象
- :param tolerance: 收敛阈值
- :param max_iters: 最大迭代次数
- :param lambda_init: 初始阻尼因子lambda
- """
- self.model = model
- self.tolerance = tolerance
- self.max_iters = max_iters
- self.lambda_init = lambda_init
-
- def fit(self, x_data, y_data, initial_theta):
- """使用Levenberg-Marquardt算法拟合模型参数"""
- theta = initial_theta
- lambda_factor = self.lambda_init
-
- for i in range(self.max_iters):
- residuals = self.model.residuals(x_data, y_data, theta)
- jacobian = self.model.jacobian_matrix(x_data, theta)
-
- # 计算Hessian近似
- H = jacobian.T @ jacobian
-
- # 更新公式中增加lambda项
- delta_theta = np.linalg.inv(H + lambda_factor * np.eye(H.shape[0])) @ jacobian.T @ residuals
-
- # 更新参数
- theta_new = theta + delta_theta
-
- # 计算新的残差
- residuals_new = self.model.residuals(x_data, y_data, theta_new)
-
- # 判断是否收敛
- if np.linalg.norm(delta_theta) < self.tolerance:
- print(f"迭代收敛，共迭代 {i+1} 次")
- return theta_new
-
- # 动态调整lambda
- if np.linalg.norm(residuals_new) < np.linalg.norm(residuals):
- lambda_factor /= 10 # 减少lambda
- theta = theta_new
- else:
- lambda_factor *= 10 # 增加lambda
-
- print("达到最大迭代次数，未能完全收敛。")
- return theta
- # 使用示例
- if __name__ == "__main__":
- # 定义非线性模型 y = a * exp(b * x)
- def func(x, theta):
- return theta[0] * np.exp(theta[1] * x)
-
- # 定义Jacobian矩阵
- def jacobian(x, theta):
- J = np.zeros((len(x), len(theta)))
- J[:, 0] = np.exp(theta[1] * x)
- J[:, 1] = theta[0] * x * np.exp(theta[1] * x)
- return J
-
- # 创建模型和LM算法实例
- model = NonlinearModel(func, jacobian)
- lm_solver = LevenbergMarquardt(model)
-
- # 生成数据
- x_data = np.linspace(0, 1,
- 10)
- y_data = 2 * np.exp(3 * x_data) + np.random.normal(0, 0.1, size=x_data.shape)
-
- # 初始参数猜测
- initial_theta = np.array([1, 1])
-
- # 执行拟合
- optimal_theta = lm_solver.fit(x_data, y_data, initial_theta)
- print("最优参数:", optimal_theta)
4. LM 算法应用实例：非线性曲线拟合
场景形貌
我们将使用 Levenberg-Marquardt 算法来拟合非线性模子 y = a ⋅ e b ⋅ x y = a \cdot e^{b \cdot x} y=a⋅eb⋅x。生成一组具有噪声的样本数据，并通过算法找到最优的 a a a 和 b b b 参数，使得模子曲线与数据点的误差最小。
结果分析与可视化
通过 matplotlib 可视化拟合结果：
- import matplotlib.pyplot as plt
- # 绘制拟合曲线与真实数据点
- y_fit = func(x_data, optimal_theta)
- plt.scatter(x_data, y_data, label="数据点")
- plt.plot(x_data, y_fit, label="拟合曲线", color='r')
- plt.xlabel("x")
- plt.ylabel("y")
- plt.legend()
- plt.show()
5. LM 算法的改进与扩展
LM 算法中的参数调节
LM 算法中的阻尼因子 (\lambda) 是一个至关重要的调节参数。初始值的选择以及动态调整策略会影响算法的收敛速度和结果。通常的做法是根据误差的变化动态增大或减小 (\lambda) 的值，从而在梯度下降和高斯牛顿法之间做出均衡。
LM 算法的改进与其他变种
- Levenberg-Marquardt-Bayesian：结合贝叶斯理论的 LM 算法，使得模子的参数更新不仅仅依靠于当前的误差，还考虑参数的先验分布。
- Trust-region-reflective：一种在约束条件下的 LM 算法变种，实用于有约束的最小二乘问题。
6. 总结
在本篇博客中，我们详细介绍了 Levenberg-Marquardt 算法的基本思想、数学原理以及其 Python 的面向对象实现。通过非线性曲线拟合的示例，我们展示了该算法在实际应用中的有效性。LM 算法以其快速收敛和较好的鲁棒性，成为处置惩罚非线性最小二乘问题的重要工具。
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。 |
