第六讲线性方程组基本迭代法

考虑线性方程组

A x = b, \quad A \in \mathbb {R} ^ {n \times n}, b \in \mathbb {R} ^ {n}.

目前, 求解线性方程组的方法有:

直接法: PLU 分解, LDL ${}^{\mathrm{T}}$ 分解,Cholesky 分解等
迭代法：

$\triangleright$ 基本迭代法: Jacobi, Gauss-Seidel, SOR, SSOR, Richardson, ADI 等
$\triangleright$ Krylov子空间迭代法：CG, MINRES, GMRES, BiCGStab等

快速方法:

$\triangleright$ 基于各种快速变换的求解方法, 如 FFT, DCT, DST 等; 代数多重网格法 (Algebraic multigrid); 快速多极子法 (Fast multipole), 等等.

快速方法通常只适用于某些具有特殊结构的方程组.
实际应用中, 这些方法常常结合使用, 如混合 (hybrid) 方法, 预处理 (preconditioning) 方法等.

直接法的优点是稳定可靠, 能在有限步内得到近似解 (如果不考虑误差, 则得到精确解), 而且所需存储量和运算量都是可知的. 缺点是所需运算量约为 $\mathcal{O}(n^3)$ , 这对于大规模线性方程组来说是非常巨大的. 而且在实际应用中, 很多问题中需要求解的大规模线性方程组都是稀疏的, 如偏微分方程的有限差分/有限元离散, 但直接法很难有效地利用问题的稀疏性来降低总运算量, 而迭代法则可以很好地利用问题的稀疏性, 大大降低运算量.

从历史上看, 最早的迭代法可以追溯到十九世纪 Gauss, Jacobi, Seidel 和 Nekrasov 等的工作 [13, 76]. 但是针对迭代法的系统研究主要还是在计算机出现以后, 大约是从二十世纪五十年代开始. 在开始阶段主要研究的是基本迭代法 [6, 30, 139] (也称经典迭代法 [57]), 典型代表有 Jacobi, GS, SOR, SSOR, ADI, Chebyshev 迭代法等. 在这期间, 有两本非常有名的经典著作, 一本是 Varga 的 “Matrix Iterative Analysis” (1962) [139], 另一本是 Yong 的 “Iterative Solution of Large Linear Systems” (1971) [141]. 基本迭代法的收敛性有着非常完美的理论分析, 但在实际使用中却存在着许多不足, 比如收敛速度较慢, 最优参数估计困难, 等等.

从二十世纪七十年代中期开始, 研究重点慢慢转向 Krylov 子空间迭代法. 事实上, 早在 1952 年, Lanczos [85] 和 Hestenes & Stiefel [68] 就同时独立地提出了求解对称正定线性方程组的共轭梯度法 (CG). 对于一个 $n$ 阶的线性方程组, 如果不考虑舍入误差的影响, 则共轭梯度法在 $n$ 步后就一定会得到精确解. 因此共轭梯度法一开始被认作是直接法. 但在实际使用中发现, 由于舍入误差的影响, 迭代步数可能会超过 $n$ , 特别是对于坏条件问题. 而对于条件数较小的线性方程组, 在给定精度下, 所需迭代步数则可能远远小于 $n$ , 这使得共轭梯度法具有一定的吸引力. 但由于种种原因 [55], 共轭梯度法提出后并没有受到重视, 在其出现后的近二十年里, 主流方法仍然是 Gauss

消去法, SOR 迭代法和 Chebyshev 迭代法.

1971年, Reid[106]指出, 对于好条件的大规模稀疏线性方程组, 共轭梯度法能在很少的迭代步数内得到一个很好的近似解(事实上, Engeli等[41]在1959年就发现了该现象, 但并没有引起关注). 特别是预处理方法的引入[94], 大幅提升了共轭梯度法的收敛速度, 这极大地促发了大家对共轭梯度法的研究兴趣, 包括各种改进和推广, 如求解对称不定线性方程组的MINRES迭代法和SYMMLQ迭代法[99], 求解一般线性方程组的GMRES迭代法[111], QMR迭代法[48], BiCGSTAB迭代法[129], 等等. 目前, 带预处理的Krylov子空间迭代法已成为求解大规模稀疏线性方程组的主流方法.

本讲介绍常用的基本迭代法, 关于 Krylov 子空间迭代法, 我们将在下一讲介绍.

关于线性方程组基本迭代法的相关参考资料

G. H. Golub and C. F. Van Loan, Matrix Computations, 2013. [57]
R. S. Varga, Matrix Iterative Analysis, 2nd edition, 2000. [139]
D. M. Young, Iterative Solution of Large Linear Systems, 1971. [141]
$\triangleright$ O. Axelsson, Iterative Solution Methods, 1994. [6]
R. Barrett, et.al, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 1994. [11]
$\triangleright$ 徐树方, 矩阵计算的理论与方法, 1995. [150]
$\triangleright$ Y. Saad and H. A. van der Vorst, Iterative solution of linear systems in the 20th century, 2000. [114]

随着矩阵规模的增大, 直接法的运算量也随之快速增长. 对于大规模的线性方程组, 由于运算量太大, 往往会采用迭代法.

当直接求解 $Ax = b$ 比较困难时, 我们可以求解一个比较容易求解的近似等价方程组 $Mx = b$ , 其中 $M$ 可以看作是 $A$ 在某种意义下的近似. 设 $Mx = b$ 的解为 $x^{(1)}$ . 易知它与原方程组的解 $x_{*} = A^{-1}b$ 之间的差距满足

A \left(x _ {*} - x ^ {(1)}\right) = b - A x ^ {(1)}.

如果 $x^{(1)}$ 已经满足精度要求, 即非常接近真解 $x_{*}$ , 则可以停止计算, 否则需要修正. 记 $\Delta x \triangleq x_{*} - x^{(1)}$ , 则 $\Delta x$ 满足方程 $A\Delta x = b - Ax^{(1)}$ . 但由于直接求解该方程组比较困难 (与求解原方程组一样困难), 因此我们还是通过求解近似方程组

M \Delta x ^ {(1)} = b - A x ^ {(1)},

得到一个修正量 $\Delta x^{(1)}$ . 于是修正后的近似解为

x ^ {(2)} = x ^ {(1)} + \Delta x ^ {(1)} = x ^ {(1)} + M ^ {- 1} (b - A x ^ {(1)}).

如果 $x^{(2)}$ 已经满足精度要求, 则停止计算, 否则继续按以上的方式进行修正, 即求解 $M\Delta x^{(2)} = b - Ax^{(2)}$ 得到修正量 $\Delta x^{(2)}$ , 然后加到 $x^{(2)}$ 上得到 $x^{(3)}$ :

x ^ {(3)} = x ^ {(2)} + \Delta x ^ {(2)} = x ^ {(2)} + M ^ {- 1} (b - A x ^ {(2)}).

不断重复以上步骤, 于是, 我们就得到一个序列

x ^ {(1)}, x ^ {(2)}, \ldots ,, x ^ {(k)}, \ldots .

满足以下递推关系

x ^ {(k + 1)} = x ^ {(k)} + M ^ {- 1} (b - A x ^ {(k)}), \quad k = 1, 2, \dots .

这就构成了一个迭代法

注记

通常, 构造一个好的迭代法, 需要考虑以下两点:

(1) 以 $M$ 为系数矩阵的线性方程组要比原线性方程组更容易求解；
(2) $M$ 应该是 $A$ 的一个很好的近似, 而且迭代序列 $\{x^{(k)}\}$ 要收敛.

常用的基本迭代法主要包括：

Richardson 迭代法
Jacobi 迭代法
Gauss-Seidel (G-S) 迭代法
超松弛 (SOR, Successive Over-Relaxation) 迭代法
对称超松弛 (SSOR, Symmetric SOR) 迭代法
加速超松弛 (AOR, Accelerated Over-Relaxation) 迭代法
Chebyshev (加速) 迭代法
交替方向 (ADI, Alternating Direction Implicit) 迭代法

README

第六讲 线性方程组基本迭代法

关于线性方程组基本迭代法的相关参考资料

注记

第六讲线性方程组基本迭代法