LM(Levenberg-Marquardt) & PCG(预条件共轭梯度法)

文章目录

• • PCG共轭梯度法
  • • PCG算法步骤
    • 示例程序实现
    • 运行结果
    • PCG vs LM
    • • PCG优势
      • PCG劣势
      • LM优势
      • LM劣势
      • 综合比较
      • 对比示例代码
      • 结果与对比
      • 总结

PCG共轭梯度法

预条件共轭梯度法（PCG）是一种用于求解大规模稀疏线性系统$Ax=b$的迭代方法。它在求解对称正定矩阵时特别高效。PCG 通过结合预条件技术，可以加速共轭梯度法的收敛速度。

1. 共轭梯度法（CG）：
   共轭梯度法是求解对称正定线性系统的一种迭代方法。它的基本思想是，在每一步中找到一个搜索方向，使得搜索方向与之前所有搜索方向共轭，从而在有限的步数内达到精确解。

2. 预条件共轭梯度法（PCG）：
   预条件共轭梯度法是在共轭梯度法的基础上，引入了预条件矩阵$M$，目的是改善原始矩阵$A$的条件数，从而加速收敛。预条件矩阵$M$应该是容易求逆的，并且$M^{-1}A$的条件数较小。

PCG算法步骤

设$A$是一个对称正定矩阵，$M$是预条件矩阵，求解$Ax=b$：

1. 初始化：设初始解$x_0$，计算初始残差$r_0=b-A{x}_0$，设置初始搜索方向$p_0=M^{-1}{r}_0$。
2. 迭代步骤：
   • 计算${\alpha }_k=\frac{r_k^T{M}^{-1}{r}_k}{p_k^{T}A{p}_k}$。
   • 更新解：$x_{k+1}=x_k+{\alpha }_k{p}_k$。
   • 更新残差：$r_{k+1}=r_k-{\alpha }_{k}A{p}_k$。
   • 检查收敛：如果$r_{k+1}$的范数小于给定的阈值，则停止迭代。
   • 计算新的搜索方向系数：${\beta }_k=\frac{r_{k+1}^T{M}^{-1}{r}_{k+1}}{r_k^T{M}^{-1}{r}_k}$。
   • 更新搜索方向：$p_{k+1}=M^{-1}{r}_{k+1}+{\beta }_k{p}_k$。
3. 返回最终解$x_k$。

示例程序实现

以下是用Python实现的PCG共轭梯度法：

import numpy as np  def pcg(A, b, M_inv, x0=None, tol=1e-8, max_iter=None):     n = len(b)     if x0 is None:         x0 = np.zeros(n)     if max_iter is None:         max_iter = 2 * n          x = x0     r = b - A @ x     z = M_inv @ r     p = z     rz_old = np.dot(r, z)          for i in range(max_iter):         Ap = A @ p         alpha = rz_old / np.dot(p, Ap)         x = x + alpha * p         r = r - alpha * Ap                  if np.linalg.norm(r) < tol:             break                  z = M_inv @ r         rz_new = np.dot(r, z)         beta = rz_new / rz_old         p = z + beta * p         rz_old = rz_new          return x  # 示例 if __name__ == "__main__":     A = np.array([[4, 1], [1, 3]], dtype=float)     b = np.array([1, 2], dtype=float)     M_inv = np.linalg.inv(np.array([[4, 0], [0, 3]], dtype=float))  # 预条件矩阵的逆     x0 = np.zeros_like(b)          x = pcg(A, b, M_inv, x0)     print("解 x =", x)     print("验证 Ax =", A @ x)     print("原向量 b =", b) 

运行结果

解 x = [0.09090909 0.63636364] 验证 Ax = [1. 2.] 原向量 b = [1. 2.] 

• A 是待求解的对称正定矩阵。
• b 是等式右侧的向量。
• M_inv 是预条件矩阵的逆。
• x0 是初始解。

通过预条件共轭梯度法，我们可以快速求解大规模稀疏线性系统。在实际应用中，选择适当的预条件矩阵是关键，它直接影响算法的收敛速度和计算效率。

预条件共轭梯度法（PCG）和Levenberg-Marquardt（LM）算法都是解决非线性优化问题的有效方法，但它们在原理、适用场景和性能方面有明显的不同。以下是PCG和LM算法的优缺点比较：

PCG vs LM

PCG优势

1. 适用于大规模稀疏线性系统：

   • PCG特别适用于求解大型稀疏对称正定矩阵的线性系统，能够高效处理大规模数据。

2. 内存使用效率高：

   • 由于PCG在每次迭代中只需要存储少量的向量，因此对于内存的使用非常高效，适合内存受限的环境。

3. 渐近最优性：

   • PCG在理论上可以在有限步内收敛到精确解（对于理想的预条件情况下），适用于一些特殊的工程和科学计算问题。

4. 易于并行化：

   • 由于PCG的矩阵-向量乘法和预条件应用可以并行化，PCG算法非常适合在多核处理器和分布式计算环境中实现。

PCG劣势

1. 依赖预条件：

   • PCG的性能高度依赖于预条件矩阵的选择。如果预条件选择不当，可能会导致收敛速度很慢甚至不收敛。

2. 适用范围有限：

   • PCG主要适用于对称正定矩阵，对于非对称或非正定矩阵，PCG不能直接应用。

LM优势

1. 解决非线性最小二乘问题的标准方法：

   • LM算法是求解非线性最小二乘问题的标准方法，特别适合于曲线拟合、参数估计和图像处理中的许多问题。

2. 鲁棒性强：

   • LM算法结合了高斯-牛顿法和梯度下降法的优点，具有较强的鲁棒性和稳定性，能够处理一些具有噪声和非线性特性的实际问题。

3. 良好的全局收敛性：

   • LM算法通过调节阻尼参数，能够在全局范围内有效收敛，适用于各种初始条件。

LM劣势

1. 计算复杂度高：

   • LM算法在每次迭代中需要计算雅可比矩阵并求解线性系统，对于大规模问题，计算和存储成本非常高。

2. 内存占用大：

   • 由于需要存储和操作较大的雅可比矩阵和海森矩阵，LM算法对内存的需求较大，不适合非常大规模的问题。

3. 依赖初始猜测：

   • LM算法的收敛速度和效果依赖于初始参数的选择，对于一些初始猜测不佳的情况，可能会收敛到局部最优解。

综合比较

• 适用场景：

  • PCG适合于大规模稀疏对称正定线性系统。
  • LM适合于中小规模的非线性最小二乘优化问题。

• 计算性能：

  • PCG在处理大规模稀疏问题时具有高效性和内存使用效率。
  • LM在求解复杂非线性问题时具有更好的鲁棒性和全局收敛性，但计算和内存成本较高。

• 预条件依赖性：

  • PCG的性能依赖于预条件矩阵的选择。
  • LM不需要预条件，但依赖于合理的初始参数选择。

PCG和LM各有其优势和劣势，选择哪种算法取决于具体问题的性质和需求。如果需要解决大规模稀疏线性系统，尤其是对称正定矩阵，PCG是更好的选择。而对于复杂的非线性最小二乘问题，特别是中小规模的数据集，LM则更为适用。

将使用一个简单的非线性最小二乘问题来进行对比，定义目标函数为以下形式：

$$f(x)=\sum _{i=1}^n(y_i-(a\cdot e^{b\cdot x_i}){)}^2$$

其中，参数$a$和$b$是待优化的变量。

对比示例代码

我们将使用Python和常用的科学计算库（NumPy和SciPy）来实现PCG和LM算法，并进行对比。

import numpy as np from scipy.optimize import least_squares from scipy.sparse import diags from scipy.sparse.linalg import cg  # 生成示例数据 np.random.seed(0) n = 100 x = np.linspace(0, 10, n) a_true, b_true = 2.5, 1.3 y = a_true * np.exp(b_true * x) + 0.1 * np.random.randn(n)  # 定义目标函数 def residuals(params, x, y):     a, b = params     return a * np.exp(b * x) - y  # 定义雅可比矩阵 def jacobian(params, x, y):     a, b = params     J = np.zeros((x.size, 2))     J[:, 0] = np.exp(b * x)     J[:, 1] = a * x * np.exp(b * x)     return J  # 预条件共轭梯度法求解 def pcg_solve(A, b, tol=1e-8, max_iter=None):     M_inv = diags([1.0 / A.diagonal()], [0])  # 使用对角预条件     x, info = cg(A, b, tol=tol, maxiter=max_iter, M=M_inv)     return x  # Levenberg-Marquardt算法求解 def lm_solve(x, y, initial_params):     result = least_squares(residuals, initial_params, jac=jacobian, args=(x, y), method='lm')     return result.x  # 初始猜测 initial_params = [1.0, 0.5]  # LM算法求解 params_lm = lm_solve(x, y, initial_params) print("LM求解参数:", params_lm)  # 生成线性系统Ax = b用于PCG J = jacobian(initial_params, x, y) A = J.T @ J b = -J.T @ residuals(initial_params, x, y)  # PCG算法求解 params_pcg = pcg_solve(A, b) print("PCG求解参数:", params_pcg)  # 对比结果 print("真实参数:", [a_true, b_true]) 

结果与对比

• 真实参数：[2.5, 1.3]

• LM求解参数：接近真实参数；

  • 适用于非线性最小二乘问题，直接使用SciPy中的least_squares实现。
  • 结果较为准确，接近真实参数。

• PCG求解参数：由于PCG本质上是用于线性系统求解，这里的实现是通过将非线性问题线性化，可能与真实参数有一定偏差；

  • 通常用于解决线性系统，这里通过构造Jacobian矩阵近似求解非线性问题。
  • 结果可能与真实参数有一定偏差，因为PCG在处理非线性问题时需要线性化近似。

总结

• LM算法在处理非线性最小二乘问题时表现更好，能够更准确地求解参数。
• PCG算法虽然主要用于线性系统求解，但在引入预条件和线性化近似后，也能在一定程度上解决非线性问题。

对于大规模稀疏线性系统，PCG有显著优势；而对于非线性优化问题，LM算法则表现出色。

LM参考链接