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

PCG共轭梯度法

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

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

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

PCG算法步骤

A A A是一个对称正定矩阵, M M M是预条件矩阵,求解 A x = b Ax = b Ax=b

  1. 初始化:设初始解 x 0 x_0 x0,计算初始残差 r 0 = b − A x 0 r_0 = b - Ax_0 r0=bAx0,设置初始搜索方向 p 0 = M − 1 r 0 p_0 = M^{-1}r_0 p0=M1r0
  2. 迭代步骤:
    • 计算 α k = r k T M − 1 r k p k T A p k \alpha_k = \frac{r_k^T M^{-1} r_k}{p_k^T A p_k} αk=pkTApkrkTM1rk
    • 更新解: x k + 1 = x k + α k p k x_{k+1} = x_k + \alpha_k p_k xk+1=xk+αkpk
    • 更新残差: r k + 1 = r k − α k A p k r_{k+1} = r_k - \alpha_k A p_k rk+1=rkαkApk
    • 检查收敛:如果 r k + 1 r_{k+1} rk+1的范数小于给定的阈值,则停止迭代。
    • 计算新的搜索方向系数: β k = r k + 1 T M − 1 r k + 1 r k T M − 1 r k \beta_k = \frac{r_{k+1}^T M^{-1} r_{k+1}}{r_k^T M^{-1} r_k} βk=rkTM1rkrk+1TM1rk+1
    • 更新搜索方向: p k + 1 = M − 1 r k + 1 + β k p k p_{k+1} = M^{-1}r_{k+1} + \beta_k p_k pk+1=M1rk+1+βkpk
  3. 返回最终解 x k x_k xk

示例程序实现

以下是用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 ) = ∑ i = 1 n ( y i − ( a ⋅ e b ⋅ x i ) ) 2 f(x) = \sum_{i=1}^n (y_i - (a \cdot e^{b \cdot x_i}))^2 f(x)=i=1n(yi(aebxi))2

其中,参数 a a a b b 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参考链接

  • 13
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值