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

文章目录

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

PCG共轭梯度法

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

共轭梯度法（CG）：
共轭梯度法是求解对称正定线性系统的一种迭代方法。它的基本思想是，在每一步中找到一个搜索方向，使得搜索方向与之前所有搜索方向共轭，从而在有限的步数内达到精确解。
预条件共轭梯度法（PCG）：
预条件共轭梯度法是在共轭梯度法的基础上，引入了预条件矩阵 $M$ ，目的是改善原始矩阵 $A$ 的条件数，从而加速收敛。预条件矩阵 $M$ 应该是容易求逆的，并且 $M^{-1}A$ 的条件数较小。

PCG算法步骤

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

初始化：设初始解 $x_0$ ，计算初始残差 $r_0 = b - Ax_0$ ，设置初始搜索方向 $p_0 = M^{-1}r_0$ 。
迭代步骤：
- 计算 $\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$ 。
返回最终解 $x_k$ 。

示例程序实现

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

import numpy as npdef 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 * nx = x0r = b - A @ xz = M_inv @ rp = zrz_old = np.dot(r, z)for i in range(max_iter):Ap = A @ palpha = rz_old / np.dot(p, Ap)x = x + alpha * pr = r - alpha * Apif np.linalg.norm(r) < tol:breakz = M_inv @ rrz_new = np.dot(r, z)beta = rz_new / rz_oldp = z + beta * prz_old = rz_newreturn 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优势

适用于大规模稀疏线性系统：
- PCG特别适用于求解大型稀疏对称正定矩阵的线性系统，能够高效处理大规模数据。
内存使用效率高：
- 由于PCG在每次迭代中只需要存储少量的向量，因此对于内存的使用非常高效，适合内存受限的环境。
渐近最优性：
- PCG在理论上可以在有限步内收敛到精确解（对于理想的预条件情况下），适用于一些特殊的工程和科学计算问题。
易于并行化：
- 由于PCG的矩阵-向量乘法和预条件应用可以并行化，PCG算法非常适合在多核处理器和分布式计算环境中实现。

PCG劣势

依赖预条件：
- PCG的性能高度依赖于预条件矩阵的选择。如果预条件选择不当，可能会导致收敛速度很慢甚至不收敛。
适用范围有限：
- PCG主要适用于对称正定矩阵，对于非对称或非正定矩阵，PCG不能直接应用。

LM优势

解决非线性最小二乘问题的标准方法：
- LM算法是求解非线性最小二乘问题的标准方法，特别适合于曲线拟合、参数估计和图像处理中的许多问题。
鲁棒性强：
- LM算法结合了高斯-牛顿法和梯度下降法的优点，具有较强的鲁棒性和稳定性，能够处理一些具有噪声和非线性特性的实际问题。
良好的全局收敛性：
- LM算法通过调节阻尼参数，能够在全局范围内有效收敛，适用于各种初始条件。

LM劣势

计算复杂度高：
- LM算法在每次迭代中需要计算雅可比矩阵并求解线性系统，对于大规模问题，计算和存储成本非常高。
内存占用大：
- 由于需要存储和操作较大的雅可比矩阵和海森矩阵，LM算法对内存的需求较大，不适合非常大规模的问题。
依赖初始猜测：
- LM算法的收敛速度和效果依赖于初始参数的选择，对于一些初始猜测不佳的情况，可能会收敛到局部最优解。

综合比较

适用场景：
- PCG适合于大规模稀疏对称正定线性系统。
- LM适合于中小规模的非线性最小二乘优化问题。
计算性能：
- PCG在处理大规模稀疏问题时具有高效性和内存使用效率。
- LM在求解复杂非线性问题时具有更好的鲁棒性和全局收敛性，但计算和内存成本较高。
预条件依赖性：
- PCG的性能依赖于预条件矩阵的选择。
- LM不需要预条件，但依赖于合理的初始参数选择。

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

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

$\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 = paramsreturn a * np.exp(b * x) - y# 定义雅可比矩阵
def jacobian(params, x, y):a, b = paramsJ = 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在处理非线性问题时需要线性化近似。