Jacobi 预处理共轭梯度法¶

约 610 个字 57 行代码预计阅读时间 3 分钟

PCG 的基本概念¶

PCG 即 预处理共轭梯度法（Preconditioned Conjugate Gradient）。
它是在共轭梯度法（CG）基础上加入了 预处理算子，用以加快对大型稀疏对称正定线性方程组 $\mathbf{A}\mathbf{x} = \mathbf{b}$ 的收敛速度。

预处理的作用¶

预处理 的目标是把原方程组通过“变换”，使得新的线性系统的 谱性质更好（条件数更小），从而使迭代方法收敛更快、更稳定。
实质是引入一个“预处理矩阵” $\mathbf{M}$（一般对称正定且易于求逆），将原系统左乘 $\mathbf{M}^{-1}$，实际迭代时只需解 $\mathbf{M}z = r$ 这样的小方程。

PCG 的算法流程（分量表达）¶

假设有 $\mathbf{A}\mathbf{x} = \mathbf{b}$，$\mathbf{M}$ 为预处理矩阵：

初始化 $x^0$，$r^0 = b - Ax^0$
解预处理方程 $\mathbf{M}z^0 = r^0$，令 $p^0 = z^0$
对 $k=0,1,2,\ldots$，迭代如下：

\[ \begin{aligned} & \alpha_k = \frac{(r^k, z^k)}{(p^k, Ap^k)} \\ & x^{k+1} = x^k + \alpha_k p^k \\ & r^{k+1} = r^k - \alpha_k A p^k \\ & \text{解}~ \mathbf{M}z^{k+1} = r^{k+1} \\ & \beta_k = \frac{(r^{k+1}, z^{k+1})}{(r^k, z^k)} \\ & p^{k+1} = z^{k+1} + \beta_k p^k \end{aligned} \]

直到 $\|r^{k+1}\| < \epsilon$

本系列文章主要讲述较为常见的Jacobi预处理器和Cholesky预处理器，本次推文主要分享Jacobi预处理器。

Jacobi 预处理 PCG 的分量形式¶

Jacobi预处理： 预处理矩阵 $\mathbf{M}$ 取 $\mathbf{A}$ 的对角元组成的对角阵
$\mathbf{M} = \mathrm{diag}(a_{11}, a_{22}, \ldots, a_{nn})$

每步的PCG主要分量公式如下：

预处理向量 $z$ 的计算（同步更新）： $$ z_i = \frac{r_i}{a_{ii}} $$
步长 $\alpha_k$ 计算：

\[ \alpha_k = \frac{\sum_{i=1}^n r_i^{(k)} z_i^{(k)}}{\sum_{i=1}^n p_i^{(k)} (A p^{(k)})_i} \]

新解向量更新： $$ x_i^{(k+1)} = x_i^{(k)} + \alpha_k p_i^{(k)} $$
残差向量更新： $$ r_i^{(k+1)} = r_i^{(k)} - \alpha_k (A p^{(k)})_i $$
新的预处理向量： $$ z_i^{(k+1)} = \frac{r_i^{(k+1)}}{a_{ii}} $$
$\beta_k$ 的计算：

\[ \beta_k = \frac{\sum_{i=1}^n r_i^{(k+1)} z_i^{(k+1)}}{\sum_{i=1}^n r_i^{(k)} z_i^{(k)}} \]

新的搜索方向： $$ p_i^{(k+1)} = z_i^{(k+1)} + \beta_k p_i^{(k)} $$

subroutine pcg_j(a, b, x, x0, n)
    ! 预处理共轭梯度法（PCG，Jacobi前置/斜前置）
    integer, intent(in) :: n
    real*8, intent(in) :: a(n,n), b(n), x0(n)
    real*8, intent(out) :: x(n)

    real*8 :: x1(n), x2(n), r0(n), r1(n)
    real*8 :: z0(n), z1(n), p0(n), p1(n)
    real*8 :: m(n)    ! 预处理向量，对角元的倒数
    integer :: i, k, iter_max = 1000
    real*8 :: eps = 1.0d-7
    real*8 :: tmp1, tmp2, tmp3, afa, beta, dr_s

    write(*,501)

    ! 构造预处理向量 M = diag(A) 的倒数
    do i = 1, n
        m(i) = 1.0d0 / a(i, i)
    end do

    x1 = x0
    r0 = b - matmul(a, x1)
    z0 = r0 * m          ! 逐元素相乘，z0 = M r0
    p0 = z0

    do k = 1, iter_max
        tmp1 = dot_product(r0, z0)
        tmp2 = dot_product(p0, matmul(a, p0))
        afa = tmp1 / tmp2

        x2 = x1 + afa * p0

        ! 收敛判据：残差欧氏范数
        dr_s = sqrt(dot_product(r0, r0))
        if (dr_s < eps) exit

        r1 = r0 - afa * matmul(a, p0)
        z1 = r1 * m        ! 预处理，z1 = M r1

        tmp3 = dot_product(r1, z1)
        beta = tmp3 / tmp1

        p1 = z1 + beta * p0

        write(*,503) k, (x2(i), i=1,n)

        ! 全部更新
        r0 = r1
        z0 = z1
        p0 = p1
        x1 = x2
    end do
    x = x2

501 format(//,18x,"Preconditioned Conjugate Gradient Iterative",//)
503 format(i5, 4(3x, f12.8))
end subroutine pcg_j

在进行有限元计算时，常常会遇到较大型的稀疏矩阵，这时需要我们在上面程序的基础上稍作修改，适应常见的稀疏矩阵存储方式，如COO、CSR、CSC等，后续推文都会讲到，感兴趣的读者，可持续关注，谢谢。