共轭梯度法¶

约 275 个字 72 行代码预计阅读时间 2 分钟

计算原理¶

共轭梯度法是一种高效求解 对称正定矩阵 线性方程组 $\mathbf{A}\mathbf{x} = \mathbf{b}$ 的迭代算法。与最速下降法相比，共轭梯度法每次沿着共轭方向搜索，极大提升了收敛速度。

算法步骤¶

设：

$\mathbf{A}$ 为 $n \times n$ 对称正定系数矩阵
$\mathbf{b}$ 为已知向量
$x^0$ 为初始猜测
$TOL$ 为误差容限
$IMAX$ 为最大允许迭代次数

具体步骤如下：

初始化 $$ \mathbf{x}^0 = x_0 $$ $$ \mathbf{r}^0 = \mathbf{b} - \mathbf{A} x^0 $$ $$ \mathbf{p}^0 = \mathbf{r}^0 $$
迭代循环 对 $k = 0, 1, 2, \ldots, IMAX-1$，做以下步骤：
1. 计算步长因子 $\alpha_k = \frac{(\mathbf{r}^k, \mathbf{r}^k)}{(\mathbf{p}^k, \mathbf{A} \mathbf{p}^k)}$
2. 更新解向量 $\mathbf{x}^{k+1} = \mathbf{x}^k + \alpha_k \mathbf{p}^k$
3. 更新残差 $\mathbf{r}^{k+1} = \mathbf{r}^k - \alpha_k \mathbf{A} \mathbf{p}^k$
4. 收敛判据（停止条件） $\|\mathbf{r}^{k+1}\| < TOL$ 其中 $\|\cdot\|$ 为向量欧氏范数。如果满足则终止迭代。
5. 计算共轭系数 $\beta_k = \frac{(\mathbf{r}^{k+1}, \mathbf{r}^{k+1})}{(\mathbf{r}^k, \mathbf{r}^k)}$
6. 更新搜索方向 $\mathbf{p}^{k+1} = \mathbf{r}^{k+1} + \beta_k \mathbf{p}^k$
输出结果 $\mathbf{x} = \mathbf{x}^{k+1}$

Fortran数值实现¶

subroutine conjugate_gradient(a, b, x, x0, n)
! 共轭梯度法
    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),p0(N),p1(N)
    integer :: i, k, iter_max = 1000
    real*8 :: eps = 1.0d-7

    write(*,501)
    x1 = x0
    r0 = b - matmul(a, x1)
    p0 = r0
    do k = 1, iter_max
        tmp1 = dot_product(r0, r0)
        tmp2 = dot_product(p0, matmul(A, p0))

        afa = tmp1 / tmp2

        x2 = x1 + afa * p0

        ! 如果残差的欧氏范数小于 tol，退出
        dr_s = sqrt(dot_product(r0, r0))
        if (dr_s < eps) exit

        r1 = r0 - afa * matmul(A, p0)

        tmp3 = dot_product(r1, r1)

        beta = tmp3 / tmp1

        p1 = r1 + beta * p0

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

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

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

数值案例¶

\[ { \left( \begin{array} { l l l l } { 5 } & { 7 } & { 6 } & { 5 } \\ { 7 } & { 1 0 } & { 8 } & { 7 } \\ { 6 } & { 8 } & { 1 0 } & { 9 } \\ { 5 } & { 7 } & { 9 } & { 1 0 } \end{array} \right) } \mathbf { x } = { \left( \begin{array} { l } { 6 2 } \\ { 8 7 } \\ { 9 1 } \\ { 9 0 } \end{array} \right) } \]

主程序：

program main
    implicit none

    integer, parameter :: n = 4
    real*8 :: a(n,n), b(n), x(n), x0(n)

    a = reshape([ &
            5, 7, 6, 5, &
            7, 10, 8, 7, &
            6, 8, 10, 9, &
            5, 7, 9, 10], [4,4], order=[2,1])

    b = [62.0, 87.0, 91.0, 90.0]

    x0 = 0.0d0

    call conjugate_gradient(a, b, x, x0, n)

end program main

输出：

                  Conjugate Gradient Iterative


    1     2.04792193     2.87369690     3.00582089     2.97278989
    2     1.56058255     2.51989545     2.58844050     4.12469697
    3     1.53689971     3.27970581     1.11628376     4.93103148
    4     2.00000000     3.00000000     1.00000000     5.00000000

注意：子程序中默认：implicit real*8(a-z)

参考文献¶

宋叶志,茅永兴,赵秀杰.Fortran 95/2003科学计算与工程[M].清华大学出版社,2011.