随着计算技术的发展, 计算机的存储量日益增大, 计算速度也迅速提高. 可是, 需要求解的问题的规模也越来越大, 对计算存储和速度的要求也越来越大. 譬如, 从地球物理和大气预测等问题中来的线性方程组的未知量有十几亿 (10 ) 个, 甚至更多. 用直接法求解这样的问题需要很大的内存空间, 而且计算量也非常大, 计算时间通常需要几天, 甚至几个月.

因此, 在现有的计算条件下, 针对应用问题设计一个高效稳定的算法便显得至关重要. 在实际应用中, 迭代法是目前求解大规模稀疏线性方程组的主要方法之一. 它能够充分利用系数矩阵的稀疏性质, 占用内存少, 运算方便, 计算程序也较为简单.

目前, 迭代法可以分为定常迭代法和不定常迭代法两大类. 定常迭代法的迭代矩阵通常保持不变, 包括雅可比 (Jacobi) 迭代法、高斯-赛德尔 (Gauss-Seidel) 迭代法和超松弛 (SOR) 迭代法等. 不定常迭代法一般指基于变分方法来最小化线性方程组的残量 $(\boldsymbol{r} = \boldsymbol{b} - \boldsymbol{A}\boldsymbol{x})$ 的一大类迭代方法, 通常没有明显的迭代矩阵, 包括求解对称正定线性方程组的共轭梯度法和求解不对称线性方程组的广义极小残量 (GMRES) 法等.

本章主要讨论这些基本的迭代方法, 包括如何建立迭代格式, 迭代收敛的条件, 误差估计和收敛速度等课题.

6.1 范数和条件数#

线性方程组的解是一个向量, 称为解向量. 近似解向量与精确解向量之差成为近似解的误差向量. 为了估计误差的大小, 需要引入衡量向量和矩阵大小的度量概念——范数.

6.1.1 向量范数和矩阵范数#

首先给出向量范数的定义.

定义 6.1.1　对任意 n 维向量 x, 若对应非负实数 $\|x\|$ , 且满足

(1) $\|x\|\geqslant0$ , 当且仅当 x=0 时等号成立;

(2) 对任意实数 $\alpha, \|\alpha x\| = |\alpha| \cdot \|\boldsymbol{x}\|;$

(3) 对任意的 $n$ 维向量 $\pmb{x}$ 和 $\pmb{y}, \| \pmb{x} + \pmb{y} \| \leqslant \| \pmb{x} \| + \| \pmb{y} \|$ ,

则称 $\| x\|$ 为向量 $x$ 的范数.

设 $\boldsymbol{x}=(x_{1},x_{2},\cdots,x_{n})^{\mathrm{T}}$ ，定义

为向量的1 范数、2 范数和无穷范数.

易证, 这三个范数都满足上述范数的三条性质.

譬如, 对于向量的 2 范数 $\| \boldsymbol{x} \|_2$ , 前两条性质显然成立.

对于性质 3, 使用柯西不等式, 有

即性质 3 成立, 所以 $\|x\|_{2}$ 是一个向量范数.

类似于向量范数的定义, 下面给出矩阵范数的定义.

定义 6.1.2　对任意 n 阶方阵 A, 若对应一个非负实数 $\|A\|$ , 满足:

(1) $\| \pmb{A}\| \geqslant 0$ ，当且仅当 $A = 0$ 时等号成立；

(2) 对任意实数 $\alpha, \|\alpha A\| = |\alpha| \cdot \|A\|;$

(3) 对任意两个 n 阶方阵 A 和 B, $\|A + B\| \leqslant \|A\| + \|B\|$ ;

(4) 对任意两个 $n$ 阶方阵 $\mathbf{A}$ 和 $\mathbf{B}, \| \mathbf{A} \mathbf{B} \| \leqslant \| \mathbf{A} \| \cdot \| \mathbf{B} \|$ ,

则称 $\| A\|$ 为方阵 $A$ 的矩阵范数.

与向量范数定义相比较, 前三个性质只是向量范数定义的推广, 而第四个性质则是矩阵乘法的要求.

对于 n 阶方阵 A, 定义

分别为矩阵的1 范数、2 范数、无穷范数和F 范数.

注意, 矩阵的 F 范数是向量 2 范数的直接推广, 矩阵的 2 范数的计算是 $A^{T}A$ 的谱半径的开方, 所以又称为谱范数. 可以证明, 这些范数都满足矩阵范数的定义要求.

在误差估计中, 我们希望不等式

成立. 但是, 对于任意的向量范数和矩阵范数, 上述不等式不一定成立.

因此, 我们特别把满足上述关系的向量范数和矩阵范数称为相容的. 在使用范数时, 注意保持矩阵范数与向量范数相容是很重要的.

就常用范数而言, 有下列一些相容关系:

6.1.2 扰动分析和条件数#

线性方程组

的解是由它的系数矩阵 A 和它的右端项 b 决定的. 而在实际问题中, 数据 A 和 b 由于各种各样的原因都带有一些误差. 譬如观测时引入的测量误差, 数学模型时引入的离散误差, 存储到计算机上运算时的舍入误差等.

通常, 通过观察或者计算引入 A 和 b 的误差相对于精确数据都是微小的, 我们称之为小的扰动, 分别计为 $\delta A$ 和 $\delta b$ .

那么, 自然要问: A 和 b 的微小扰动将对计算线性方程组的解有何影响? 很自然的想法认为, 既然 A, b 受到的扰动是微小的, 那么对应的方程组的解 x 的变化也应该是微小的. 然而, 实际上这是不对的. 请看下例.

考察线性方程组

它的解为 $\boldsymbol{x}=(1,1)^{\mathrm{T}}$ ，若方程组右端有扰动 $\delta\boldsymbol{b}=(2\times10^{-4},-2\times10^{-4})^{T}$ ，则原方程组变为

其解为 $\boldsymbol{x}=(1.5,0.5)^{\mathrm{T}}$ . 这样, 就有

即解的相对误差是右端项相对误差的 10 000 倍.

这个例子表明, 确实有一些线性方程组, 其系数的微小变化会引起解的巨大变化. 下面就一般的非奇异线性方程组 Ax = b 来讨论其扰动后的性态问题.

首先, 如果方程组右端 b 有一个小扰动 $\delta b$ , 则解 x 产生一个扰动 $\delta x$ , 即

于是有

两边取范数, 可得

另一方面, 在方程组 (6.1) 两边取范数可得

将式 (6.2) 和式 (6.3) 的两边分别相乘, 不等式依然成立, 经整理可得

如果系数矩阵 A 有一个小扰动 $\delta A$ ，则解 x 相应产生一个扰动 $\delta x$ ，即

于是有

两边取范数, 可得

对 (6.5) 进行整理, 可得

观察 (6.4) 和 (6.6) 可以得出, 无论方程组 (6.1) 中的系数矩阵 $\mathbf{A}$ 有扰动, 还是右端 $\mathbf{b}$ 有扰动, 解 $\mathbf{x}$ 的相对误差除了受相应扰动的相对误差以外, 还与 $\| \mathbf{A} \| \| \mathbf{A}^{-1} \|$ 的大小有关. 实际上, 这个数 $\| \mathbf{A} \| \| \mathbf{A}^{-1} \|$ 在问题中起到了放大倍数的作用, 对方程组的解的相对误差起了关键的控制作用. 于是, 有如下定义.

条件数的大小与所取的范数有关, 最常用的是

如果 A 对称正定, 则有

其中, $\lambda_{1}(\boldsymbol{A})$ 和 $\lambda_{n}(\boldsymbol{A})$ 分别是最大特征值和最小特征值.

条件数有如下四条性质.

(1) 对于任意的 $n$ 阶非奇异矩阵 $\mathbf{A}$ , $\operatorname{cond}(\mathbf{A}) \geqslant 1$ 成立;

(2) 对于任意的 n 阶非奇异矩阵 A 及任意非零常数 c, $\text{cond}(cA) = \text{cond}(A)$ 成立;

(3) 对于任意的正交矩阵 A, 有 $\mathrm{cond}_{2}(A)=1;$

(4) 对于任意的 $n$ 阶非奇异矩阵 $\mathbf{A}$ 及任意 $n$ 阶正交矩阵 $\mathbf{P}$ , 有 $\mathrm{cond}_2(\mathbf{PA}) = \mathrm{cond}_2(\mathbf{AP}) = \mathrm{cond}_2(\mathbf{A}).$

在实际计算中, 由于观测误差和运算过程的舍入误差, 线性方程组的系数矩阵和右端项总会带来一些扰动. 因此, 实际求解的对象总是有扰动的方程组, 得到的解是一个近似方程组的解. 这个解的好坏与条件数有极其密切的关系. 因而, 系数矩阵的条件数刻画了线性方程组的性态.

对于一个确定的线性方程组, 若系数矩阵的条件数相对地小, 就称这个方程组是良态的, 矩阵为良态矩阵; 反之, 条件数相对地大, 就称方程组是病态的, 矩阵为病态矩阵. 用一个稳定的方法去解一个良态的方程组, 必然得到较为准确的结果; 用一个稳定的方法去解一个病态的方程组, 结果就有可能很差.

一类十分典型的病态矩阵是希尔伯特 (Hilbert) 矩阵, 其定义为

其条件数 $\mathrm{cond}_{2}(\boldsymbol{H}_{n})\approx\mathrm{e}^{3.5n}$ ，随着 n 的增大，条件数非常迅速地增加。因而，希尔伯特矩阵阶数越高，病态程度就越为严重。

6.2 基本迭代法#

设给定一个线性方程组

其中 $A \in R^{n \times n}$ 和 $b \in R^{n}$ 已知, $x \in R^{n}$ 未知. 现在, 我们来考虑如何构造迭代法求 (6.7) 的解.

一般来说, 假定 A 有分解

其中 M 是非奇异方阵. 那么, 方程组 (6.7) 可以改写成

或者

其中 $B = M^{-1}N$ 且 $g = M^{-1}b$ .

从而, 我们建立迭代公式

给定一个初始向量 $\boldsymbol{x}^{(0)}$ ，按公式 (6.9) 进行迭代，就可以得到一个向量序列 $\{\boldsymbol{x}^{(k)}\}$ 。若 $\boldsymbol{x}^{(k)}$ 收敛于确定的向量 $x^{*}$ ，则 $x^{*} = Bx^{*} + g$ ，亦即 $Ax^{*} = b$ ，则 $x^{*}$ 就是方程组 (6.7) 的解。这种方法就是解线性方程组的基本迭代解法。

在公式 (6.9) 中, 为了避免 $M^{-1}N$ 和 $M^{-1}b$ 的计算, 我们可以按如下方式进行迭代

但这样一来, 每次迭代需要求解一个系数矩阵为 $M$ 的线性方程组. 因此, 我们自然希望 $M$ 应具有某种特殊性质, 使得这样的方程组易于求解. 譬如, $M$ 是对角矩阵或者上三角矩阵等. 令

其中

分别为对角矩阵、严格下三角矩阵和严格上三角矩阵.

下面, 我们介绍三种基本迭代解法: 雅可比迭代法、高斯-赛德尔迭代法和SOR迭代法, 并对它们的适用性、收敛性质和收敛速度作简要的分析.

6.2.1 雅可比迭代法#

在迭代公式 (6.10) 中, 取 $M = D$ 和 $N = L + U$ , 则得到雅可比迭代法. 其迭代格式为

算法 6.2.1 (雅可比迭代算法)

(1) 选定初值 $x_0$ ;

(2) 对 $k = 0, 1, 2, \cdots$ ，计算:

(3) $x_{i}^{(k + 1)} = \frac{1}{a_{ii}}\Big(b^{i} - \sum_{j = 1,j\neq i}^{n}a_{ij}x_{j}^{(k)}\Big);$

(4) 如果近似解达到收敛条件, 退出; 否则, 继续第 (2)~(3) 步的计算.

此时, 迭代矩阵为

Matlab 程序 jacobi.m 如下:

function [x,iter] = jacobi(A,b,tol)
D = diag(diag(A));
L = D-tril(A);
U = D-triu(A);
x = zeros(size(b));
for iter=1:500
    x = D$(b+L*x+U*x)$;

error = norm( b-A*x ) / norm(b);
if ( error < tol )
    break;
end
end

6.2.2 高斯-赛德尔迭代法#

如果雅可比迭代法是收敛的, 那么在雅可比迭代算法的第三步中, 将计算出来的 $\boldsymbol{x}^{(k+1)}$ 的分量马上投入到下一个迭代方程中使用, 这样可能会收到更好的效果.

按这样的方式建立起来的迭代格式成为高斯-赛德尔迭代法, 简称 GS 迭代, 其迭代格式为

或者

算法 6.2.2 (高斯-赛德尔迭代算法)

(1) 选定初值 $\boldsymbol{x}^{(0)}$ ;

(2) 对 $k = 0, 1, 2, \cdots$ ，计算:

(4) 如果近似解达到收敛条件, 退出; 否则, 继续第 (2)~(3) 步的计算.

Matlab 程序 gs.m 如下:

function [x,iter] = gs(A,b,tol)
D = diag(diag(A));
L = D-tril(A);
U = D-triu(A);
x = zeros(size(b));
for iter=1:500
    x = (D-L)$b+U*x$;
    error = norm(b-A*x) / norm(b);
    if (error < tol)
    break;
    end
end
此时，迭代矩阵为

6.2.3 超松弛 (SOR) 迭代法#

GS 迭代格式可以改写成

为了加快迭代的收敛速度, 我们将上式等号右端的第二项看成是修正量, 为了获得更快的收敛效果, 在修正量前乘以一个参数 $\omega$ , 即得所谓的逐次超松弛迭代法, 简称SOR迭代法, 其中, $\omega$ 称为松弛因子. SOR迭代法的迭代格式为

即

算法 6.2.3 (SOR 迭代算法)

(1) 选定初值 $\boldsymbol{x}^{(0)}$ ;

(2) 对 $k = 0, 1, 2, \cdots$ ，计算:

(3) $x_{i}^{(k + 1)} = x_{i}^{(k)} + \frac{\omega}{a_{ii}} (b_{i} - \sum_{j = 1}^{i - 1}a_{ij}x_{j}^{(k + 1)} - \sum_{j = i}^{n}a_{ij}x_{j}^{(k)});$

(4) 如果近似解达到收敛条件, 退出; 否则, 继续第 (2)~(3) 步的计算.

Matlab 程序 sor.m 如下:

function [x,iter] = sor(A,b,omega,tol)
D = diag(diag(A));
L = D-tril(A);
U = D-triu(A);
x = zeros(size(b));
for iter=1:500
    x = (D-omega*L)$b+(1-\omega)*D*x+omega*U*x$;
    error = norm(b-A*x) / norm(b);
    if (error < tol)
    break;
    end
end
此时，迭代矩阵为

注意到, 当 $\omega = 1$ 时, SOR 迭代法其实就是 GS 迭代法.

下面来看一个数值试验的例子.

例 6.2.4 取初值为 $\boldsymbol{x}^{(0)}=(0,0,0)^{\mathrm{T}}$ ，试用雅可比迭代、GS 迭代以及 SOR 迭代分别求解线性方程组

迭代过程保留 5 位有效数字 (精确解 $\boldsymbol{x}^{*} = (2, 3, -1)^{\mathrm{T}}$ ).

解：建立雅可比迭代格式

取 $x_{1}^{(0)} = x_{2}^{(0)} = x_{3}^{(0)} = 0,$ 对 $k = 0, 1, 2, \cdots$ 计算可得

这里, 雅可比方法迭代 21 步就得到了保留 5 位有效数字的近似解.

建立 GS 迭代格式

取 $x_{1}^{(0)} = x_{2}^{(0)} = x_{3}^{(0)} = 0,$ 对 $k = 0, 1, 2, \cdots$ 计算可得

这里, GS 方法迭代 9 步就得到了保留 5 位有效数字的近似解. 对于这个例子而言, GS 迭代的收敛速度大约是雅可比迭代的 2 倍.

选取 $\omega = 1.1$ ，建立SOR迭代格式

取 $x_{1}^{(0)} = x_{2}^{(0)} = x_{3}^{(0)} = 0,$ 对 $k = 0, 1, 2, \cdots$ 计算可得

这里, SOR 方法迭代 7 步就得到了保留 5 位有效数字的近似解. 对于这个例子而言, $\omega = 1.1$ 时, SOR 迭代法的收敛速度比 GS 迭代法要快.

虽然上述数值实验中所使用的三种基本迭代法都是收敛的, 但对于一般的系数矩阵而言, 这三种基本迭代的收敛性是无法得到保证的.

6.2.4 迭代的收敛性分析和误差估计#

讨论线性方程组的迭代解法的收敛性问题, 要用到与迭代解法相关的一些矩阵的基本概念和性质, 首先引进两个重要的概念.

一个排列矩阵左乘或右乘一个矩阵, 相当于对这个矩阵进行行或者列的相互调换. 设线性方程组 Ax = b 的系数矩阵 A 是一个可约矩阵, 利用排列阵 P 将方程组变换为等价的线性方程组

这里的系数矩阵 $P^{\mathrm{T}}AP$ 有分块形式 (6.21), 若令 $\pmb{y} = \pmb{P}^{\mathrm{T}}\pmb{x}, \pmb{f} = \pmb{P}^{\mathrm{T}}\pmb{b}$ , 将式 (6.21) 按式 (6.20) 写为分块形式, 则有

先从第二个方程中解出 $y_{2}$ ，再从第一个方程中解出 $y_{1}$ ，最后按排列阵 P 对 y 的分量重新排列即可得到原方程组的解。这样，就将求解方程组的问题化为求解两个阶数更低的方程组的问题。在实际计算中，后者往往比前者容易处理，计算工作量也大大减少。

引理 6.2.7　设 n 阶矩阵 A 是严格对角占优或不可约弱对角占优矩阵, 则 A 是非奇异矩阵. 证: 用反证法, 仅对第一个结论给出证明. 假定 A 是奇异矩阵, 则存在非零向量 v, 使 Av = 0. 不失一般性, 可设 $\|v\|_{\infty} = 1$ . 若令 $|v_{r}| = 1$ , 则 $|v_{j}| \leqslant |v_{r}| = 1, j = 1, \cdots, n$ , 由 Av = 0 的第 r 个方程立即得到

若 A 是严格对角占优的, 则式 (6.23) 不可能成立, 因此, A 一定是非奇异矩阵.

利用矩阵的 Jordan 标准型, 可以得出下述引理.

任意一个矩阵的谱半径与其范数是有关系的.

证：按照谱半径的定义 $\rho(\boldsymbol{A}) = \max_{1 \leqslant i \leqslant n} |\lambda_i|$ ，其中， $\lambda_i (i = 1, 2, \cdots, n)$ 为 n 阶方阵 A 的特征值。设 x 为对应于 $\lambda$ 的 A 的特征向量，则有 $\lambda x = Ax$ ，两边取范数可得

因为 x 为非零向量, $\|x\| \neq 0$ , 故有

上述这个不等式对 A 的任何特征值均成立, 即可得

定理 6.2.10　对于迭代格式 (6.9), 给定任意的初值 $x^{(0)}$ , 有下列收敛结果和误差估计:

(1) 迭代格式 (6.9) 收敛的充要条件为谱半径 $\rho(\boldsymbol{B}) < 1$ ;

(2) 如果 $\| B\| < 1$ ，则有估计

和

其中, $x^{*}$ 为(6.1)的真解.

证：(1) 迭代格式 (6.9) 为 $x^{(k+1)} = B x^{(k)} + g$ ，而真解 $x^*$ 满足 $x^* = B x^* + g$ ，两式相减可得

进而递推得

也就是

必要性设迭代收敛, 则有 $\boldsymbol{x}^{(k)} - \boldsymbol{x}^{*} \rightarrow 0 (k \rightarrow \infty)$ . 在式 (6.28) 两边取极限可得

由 $\pmb{x}^{(0)}$ 的任意性可得

即 $B^{k}$ 的任意元素都随着 $k \to \infty$ 趋向于 0. 由引理 6.2.5 可得 $\rho(B) < 1$ .

充分性设 $\rho(B) < 1$ ，则由引理6.2.5可得 $B^k \to 0 (k \to +\infty)$ ，再由式(6.28)，可得 $\pmb{x}^{(k)} - \pmb{x}^* \to 0 (k \to +\infty)$ ，即迭代收敛.

(2) 由 $\| B\| < 1$ 及定理6.2.6, 有

于是

由式 (6.28), 且在式 (6.29) 两边取范数可得

即式 (6.25) 成立.

由式 (6.27) 可得

从而，

即式 (6.26) 亦成立.

证：(1) 对于雅可比迭代, 迭代矩阵 $B_{J} = D^{-1}(L + U)$ . 使用反证法, 假设 $B_{J}$ 有一特征值 $|\lambda| \geqslant 1$ , 则有

则有

由于 A = D - L - U, 如果 A 是严格对角占优或不可约弱对角占优, 则 $\lambda D - L - U$ 仍然是严格对角占优或不可约弱对角占优矩阵. 根据引理 6.2.7, $\lambda D - L - U$ 一定是非奇异矩阵, 这显然同式 (6.30) 相矛盾. 所以 $B_{J}$ 的任一特征值按模必须小于 1, 也就是, $\rho(B_{J}) < 1$ , 故雅可比迭代收敛.

(2) 对于 GS 迭代, 迭代矩阵 $\boldsymbol{B}_{G} = (\boldsymbol{D} - \boldsymbol{L})^{-1}\boldsymbol{U}$ . 使用反证法, 假定 $B_{G}$ 有一特征值 $|\lambda| \geqslant 1$ , 则有

即有

由于 A = D - L - U 是严格对角占优或不可约弱对角占优，则 $\lambda(D - L) - U$ 仍是严格对角占优或不可约弱对角占优矩阵。由引理 6.2.7, $\lambda(D - L) - U$ 一定是非奇异矩阵。这显然同式 (6.31) 相矛盾。所以 $B_{J}$ 的任一特征值按模必须小于 1，也就是， $\rho(B_{G}) < 1$ ，故 GS 迭代收敛。定理 6.2.12 若 A 是对称正定矩阵，则雅可比迭代收敛的充要条件是 2D - A 也是对称正定矩阵。

证：由于 A 是对称正定矩阵，则 $\alpha_{ii} > 0, i = 1, 2, \cdots, n.$

其中， $D^{\frac{1}{2}}=\mathrm{diag}(\sqrt{a_{11}},\sqrt{a_{22}},\cdots,\sqrt{a_{nn}})$ .

显然 $I - D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ 对称, 从而由 $B_{J}$ 相似于对称矩阵 $I - D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ 可知, $B_{J}$ 的 n 个特征值均为实数.

必要性：设雅可比迭代收敛，则有

于是, $D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ 的任一特征值 $\mu$ 均满足

即 $0 < \mu < 2$ .

注意到

且 $2\boldsymbol{I} - \boldsymbol{D}^{-\frac{1}{2}}\boldsymbol{A}\boldsymbol{D}^{-\frac{1}{2}}$ 的特征值 $2 - \mu \in (0,2)$ ，故 $2\boldsymbol{I} - \boldsymbol{D}^{-\frac{1}{2}}\boldsymbol{A}\boldsymbol{D}^{-\frac{1}{2}}$ 对称正定.由式(6.33)，进而可得 $2D - A$ 对称正定.

充分性设 2D - A 对称正定. 一方面, 由 A 对称正定知, $D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ 对称正定, 其任一特征值 $\mu > 0$ , 于是 $I - D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ 的特征值 $1 - \mu < 1$ , 由式 (6.32) 可知 $B_{J}$ 的任一特征值 $\lambda(B_{J}) < 1$ .

另一方面, 由 $2D - A$ 对称正定可得 $2I - D^{-\frac{1}{2}}AD^{-\frac{1}{2}} = D^{-\frac{1}{2}}(2D - A)D^{-\frac{1}{2}}$ 也对称正定. 由于

因此, 可得 $I-D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ 的特征值全大于 -1, 由式 (6.32), 可知 $B_{J}$ 的任一特征值 $\lambda(B_{J}) > -1$ . 于是 $\rho(B_{J}) < 1$ , 雅可比迭代收敛.

证： $B_{S} = (D - \omega L)^{-1}[(1 - \omega)D + \omega U]$ ，设 $\lambda_{1}, \lambda_{2}, \cdots, \lambda_{n}$ 是 $B_{S}$ 的 n 个特征值，利用根与系数的关系，有

进而

如果SOR迭代收敛,则 $\rho (B_S) < 1$ ，由式(6.34)可得松弛因子 $\omega$ 必须满足 $0 < \omega < 2.$ 定理6.2.14 设系数矩阵 $\mathbf{A}$ 对称正定，则 $0 < \omega < 2$ 时SOR迭代收敛.

证：只需证明迭代矩阵 $B_{S}$ 的任一特征值 $\lambda$ 均按模小于 1.

设 $\lambda$ 为 $B_{S}$ 的任一特征值, x 为相应的特征向量, 则

由于 $x \in C^{n}$ , $x^{H}$ 是 x 的复共轭转置向量. 在酉空间 $C^{n}$ 中, 对式 (6.35) 两边作内积可得

如果记 $x^{H}Lx = \alpha + i\beta$ ，则 $x^{H}Ux = x^{H}L^{T}x = (x^{H}Lx)^{H} = \alpha - i\beta$ ，记 $x^{H}Dx = \sigma$ ，于是

在式 (6.36) 两边取模的平方

式 (6.37) 中分子与分母之差为

由于 A 对称正定, 所以 $a_{ii} > 0 (i = 1, 2, \cdots, n)$ , 可得 $\sigma = x^{H} D x = \sum_{i=1}^{n} \alpha_{ii} |x_i|^2 > 0$ . 同样利用 A 的正定性可得

利用式 (6.37)、式 (6.38) 和式 (6.39)，以及 $\sigma > 0$ 和 $0 < \omega < 2$ ，可得 $|\lambda| < 1$ ，从而谱半径 $\rho(B_S) < 1$ ，即SOR迭代收敛.

6.3 不定常迭代法#

本节将介绍两类最基本的不定常迭代方法：一类是求解对称正定线性方程组的最速下降法和共轭梯度法；另一类是求解不对称线性方程组的广义极小残量法.

最速下降法和共轭梯度法本质上来说是一种变分方法, 对应于求一个二次函数的极值. 所以也说它是一种极小化方法. 共轭梯度法一开始被当作一个直接方法出现在 20 世纪 50 年代. 八九十年代以来, 预处理共轭梯度法得到很大的发展并被广泛运用到数值计算的各个领域, 使得这类方法成为求解大型稀疏对称正定方程组的最有效方法.

6.3.1 最速下降法#

设 A 对称正定, 线性方程组为

其中 $\boldsymbol{A}=(a_{ij})\in\mathbb{R}^{n\times n},\boldsymbol{x}=(x_{1},\cdots,x_{n})^{\mathrm{T}},\boldsymbol{b}=(b_{1},\cdots,b_{n})^{\mathrm{T}}.$

首先, 我们介绍一下与求解线性方程组 (6.40) 等价的变分问题.

任取 x, 对于公式 (6.40) 中给定的 A 和 b, 定义 n 元二次函数 $\varphi: R^{n} \rightarrow R$ 为

其中, $(\boldsymbol{x},\boldsymbol{y})=\sum_{i=1}^{n}x_{i}y_{i}.$

该二次函数具有以下三条性质.

(1) 对于一切 $x \in R^{n}$ ，通过直接计算

知函数 $\varphi(\boldsymbol{x})$ 的梯度为

(2) 对于一切 $x, y \in R^{n}, a \in R,$

(3) $x^{*} = A^{-1}b$ 为方程组(6.40)的解, 则有 $\varphi(x^{*}) = -\frac{1}{2}(x^{*}, b) = -\frac{1}{2}(x^{*}, Ax^{*})$ , 且对于一切 $x \in R^{n}$ 有

以上性质可以通过直接运算验证, 其中要用到 A 的对称性.

定理6.3.1 设 $A$ 对称正定, $\pmb{x}^{*}$ 是方程组 $Ax = b$ 的解的充要条件是 $\pmb{x}^{*}$ 为二次函数 $\varphi (\pmb{x})$ 的极小值点, 即

证：必要性设 $Ax^{*}=b$ ，则由 A 的对称正定性有

上式对于任意的 $x \in R^{n}$ 都成立, 即 $\varphi(x^{*}) = \min_{x \in \mathbb{R}^{n}} \varphi(x)$ .

充分性设 $\varphi (\pmb{x}^{*}) = \min_{\pmb{x}\in \mathbb{R}^{n}}\varphi (\pmb{x})$ ，则对任意的 $\pmb{y}\in \mathbb{R}^n$ ，有

即

由

可知, 对任意的向量 $y \in R^{n}$ 都有 $(y, Ax^{*} - b) = 0$ . 由此可得 $Ax^{*} = b$ .

为了找到 $\varphi(\boldsymbol{x})$ 的极小点 $\boldsymbol{x}^*$ , 可以从任一点 $\boldsymbol{x}^{(k)}$ 出发, 沿某一指定的方向 $\boldsymbol{y}^{(k)} \in \mathbb{R}^n$ , 搜索下一个近似点 $\boldsymbol{x}^{(k+1)} = \boldsymbol{x}^{(k)} + \alpha_k \boldsymbol{y}^{(k)}$ , 使得 $\varphi(\boldsymbol{x}^{(k+1)})$ 在该方向上达到极小值. 选择 $\boldsymbol{y}^{(k)}$ 的方式不同时, 将会得到不同的算法.

令 $\boldsymbol{y}^{(k)}$ 为某一搜索方向, $\boldsymbol{r}^{(k)} = \boldsymbol{b} - \boldsymbol{A}\boldsymbol{x}^{(k)}$ 为 $\boldsymbol{x}^{(k)}$ 对应的残量, 则有函数

由

可知

同时, 由 A 的正定性有

所以 $a_{k}$ 是 $\varphi(\boldsymbol{x}^{(k)} + a\boldsymbol{y}^{(k)})$ 的极小点.

将 (6.47) 代入 (6.46), 可得

当 $(\pmb{r}^{(k)},\pmb{y}^{(k)})\neq 0,$ 即 $\pmb{y}^{(k)}$ 不与 $\pmb{r}^{(k)}$ 正交时， $\varphi (\pmb{x}^{(k + 1)}) < \varphi (\pmb{x}^{(k)})$ 成立.

相对于 $\varphi(\boldsymbol{x}^{(k)})$ ，由式 (6.48) 知 $\varphi(\boldsymbol{x}^{(k+1)})$ 的下降量取决于 $\boldsymbol{y}^{(k)}$ 的方向，而与 $\boldsymbol{y}^{(k)}$ 的大小无关。函数 $\varphi(\boldsymbol{x})$ 在点 $\boldsymbol{x}^{(k)}$ 处下降最快的方向应该是在该点的负梯度方向。

从而, 取 $\boldsymbol{y}^{(k)} = \boldsymbol{r}^{(k)}$ , 并构造 $\boldsymbol{x}^{(k+1)} = \boldsymbol{x}^{(k)} + \alpha_{k} \boldsymbol{r}^{(k)}$ , 可得求 $\varphi(\boldsymbol{x})$ 的极小点的最速下降法.

算法 6.3.2 (最速下降法)#

(1) 选定初值 $\boldsymbol{x}^{(0)}$ ;

(2) 对 $k = 0, 1, 2, \cdots$ 直到收敛, 计算:

(3) $\pmb{r}^{(k)} = \pmb{b} - \pmb{A}\pmb{x}^{(k)}$ ;

(4) $\alpha^{(k)} = \frac{(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})}{(\boldsymbol{A}\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})};$

(5) $\pmb{x}^{(k + 1)} = \pmb{x}^{(k)} + \alpha_k\pmb{r}^{(k)}$

下面讨论最速下降法的收敛性.

引理6.3.3 设 $\pmb{A}$ 是 $n$ 阶实对称正定矩阵, $f(t)$ 是 $m$ 次实系数多项式, 则对任意 $\pmb{x} \in \mathbb{R}^n$ 有

其中, $\lambda_{i}$ 是 $\pmb{A}$ 的特征值, $\| \pmb{x}\|_{\pmb{A}} = \sqrt{\pmb{x}^{\mathrm{T}}\pmb{A}\pmb{x}}$ 是 $\mathbb{R}^n$ 中的向量范数.

证：由于 A 是对称正定矩阵, 可设 A 的对应于特征值 $\lambda_{1}, \lambda_{2}, \cdots, \lambda_{n}$ 的特征向量 $y_{1}, y_{2}, \cdots, y_{n}$ 标准正交, 则对任意的 $x \in R^{n}$ , 有 $x = \sum_{i=1}^{n} c_{i} y_{i}$ , 由此可得

上式两边开方即得式 (6.49).

证：对任意实数 $\alpha$ 有 $\varphi (\pmb{x}^{(k)})\leqslant \varphi (\pmb{x}^{(k - 1)} + \alpha \pmb{r}^{(k - 1)})$ 成立,所以

利用式 (6.45) 及引理 6.3.3 可得

特别地, 取 $\alpha = 2 / (\lambda_1 + \lambda_n)$ 时,

由此可得

由误差估计式 (6.50) 可得最速下降法的收敛性. 不过, 当 $\lambda_1 \gg \lambda_n$ 时, $\frac{\lambda_1 - \lambda_n}{\lambda_1 + \lambda_n} \approx 1$ , 这时, 最速下降法的收敛速度将会很慢.

6.3.2 共轭梯度法#

对上述最速下降法作一简单分析, 就会发现, 负梯度方向虽从局部来看是最佳的搜索方向, 但从整体来看并非最优. 这就促使人们去寻找更好的搜索方向, 当然, 希望每一步确定新的搜索方向时付出的代价也不要太大.

共轭梯度法就是根据这一思想设计的．仍设 A 对称正定，我们还是采用一维极小搜索的概念．但是我们不再沿负梯度方向 $\boldsymbol{r}^{(0)}$ , $\boldsymbol{r}^{(1)}$ , $\cdots$ , $\boldsymbol{r}^{(k)}$ 搜索，而要另找一组方向 $\boldsymbol{p}^{(0)}$ , $\boldsymbol{p}^{(1)}$ , $\cdots$ , $\boldsymbol{p}^{(k)}$ ，使得进行 k 次一维搜索后，求得近似解 $\boldsymbol{x}^{(k)}$ .

对一维极小问题 $\min_{\alpha}\varphi(\boldsymbol{x}^{(k)}+\alpha\boldsymbol{p}^{(k)})$ ，我们令

可得

从而, 下一个近似解和对应的残量分别为

为了讨论方便并不失一般性, 设 $\boldsymbol{x}^{(0)} = \boldsymbol{0}$ , 可知

从而， $\boldsymbol{x}\in\operatorname{span}\{\boldsymbol{p}^{(0)},\boldsymbol{p}^{(1)},\cdots,\boldsymbol{p}^{(k)}\}$ ，其中 $\operatorname{span}\{\boldsymbol{p}^{(0)},\boldsymbol{p}^{(1)},\cdots,\boldsymbol{p}^{(k)}\}$ 是指由 $\boldsymbol{p}^{(0)},\boldsymbol{p}^{(1)},\cdots,\boldsymbol{p}^{(k)}$ 生成的子空间。记 $\boldsymbol{x}=\boldsymbol{y}+\alpha\boldsymbol{p}^{(k)}$ ， $\boldsymbol{y}\in\operatorname{span}\{\boldsymbol{p}^{(0)},\boldsymbol{p}^{(1)},\cdots,\boldsymbol{p}^{(k-1)}\}$ 。

现在考虑 $\boldsymbol{p}^{(0)},\boldsymbol{p}^{(1)},\cdots,\boldsymbol{p}^{(k)}$ 取什么方向使得二次函数 $\varphi(\boldsymbol{x}^{(k)}+\alpha\boldsymbol{p}^{(k)})$ 下降最快.

一开始, 可以令 $\boldsymbol{p}^{(0)} = \boldsymbol{r}^{(0)}$ . 当 $k \geqslant 1$ 时, 我们不但希望满足

而且希望 $\{\pmb{p}^{(k)}\}$ 的选择满足

将 $\varphi(\boldsymbol{x})$ 对 y 展开, 有

上式出现交叉项 $(\pmb{A}\pmb{y},\pmb{p}^{(k)})$ ，使求 $\varphi (\pmb{x})$ 极小问题复杂化了.

为了把上述问题分离为分别对 y 和对 $\alpha$ 求极小, 我们令

也就是

对于每一个 $k = 1, 2, \cdots, n$ ，我们都要选择 $\boldsymbol{p}^{(k)}$ 使得满足上述条件.

显然, 当 $l < n$ 时, 不含零向量的 $A$ -共轭向量组线性无关. 当 $\mathbf{A} = \mathbf{I}$ 时, $A$ -共轭性质就是一般的正交性.

若向量组 $\{\boldsymbol{p}^{(0)},\boldsymbol{p}^{(1)},\cdots,\boldsymbol{p}^{(l)}\}$ 是 A- 共轭的, 可知

不妨设 $\boldsymbol{p}^{(k)} = \boldsymbol{r}^{(k)} + \beta_{k-1}\boldsymbol{p}^{(k-1)}$ ，利用 $(\boldsymbol{p}^{(k)}, \boldsymbol{A}\boldsymbol{p}^{(k-1)}) = 0$ ，可以求出

这样的得到的 $\boldsymbol{p}^{(k)}, \boldsymbol{p}^{(k-1)}$ 是 A- 共轭的.

算法 6.3.6 (共轭梯度法)

(1) 选定初值 $x^{(0)}$ , 设 $p^{(0)} = r^{(0)} = b - Ax^{(0)}$ ;

(2) 对 $k = 0, 1, 2, \cdots$ 直到收敛, 计算:

$x^{(k+1)} = x^{(k)} + \alpha_k p^{(k)};$ $r^{(k+1)} = r^{(k)} - \alpha A p^{(k)};$
如果 $\|r^{(k+1)}\|$ 足够小，
退出算法;
否则,
$\beta_k = \|r^{(k+1)}\|_2^2 / \|r^{(k)}\|_2^2;$ $p^{(k+1)} = r^{(k+1)} + \beta_k p^{(k)}.$
Matlab 程序 cg.m 如下:
function [x,iter] = cg(A,b,tol)
x = zeros(size(b));
r = b-A*x;
for iter = 1:500
    rho = r'*r;
    if (iter == 1)
    p = r;
    else
    beta = rho / rho_1;
    p = r + beta * p;
    end
    q = A*p;
    alpha = rho / (p'*q);
    x = x + alpha * p;
    r = r - alpha * q;
    rho_1 = rho;
    error = norm(r) / norm(b);
    if (error < tol)
    break;
    end
end
下面不加证明地给出一些理论结果.
**定理 6.3.7**　由共轭梯度法得到的向量组 $\{r^{(i)}\}$ 和 $\{p^{(i)}\}$ 具有如下性质:
(1) $(p^{(i)}, r^{(j)}) = 0, \quad 0 \leqslant i < j \leqslant k;$
(2) $(r^{(i)}, r^{(j)}) = 0, \quad 0 \leqslant i, j \leqslant k; \quad i \neq j;$
(3) $(p^{(i)}, A p^{(j)}) = 0, \quad 0 \leqslant i, j \leqslant k; \quad i \neq j.$

(3) $(\pmb{p}^{(i)},\pmb{A}\pmb{p}^{(j)}) = 0,\quad 0\leqslant i,j\leqslant k;i\neq j.$

由误差估计式 (6.51) 可得共轭梯度法的收敛性. 比较式 (6.50) 和式 (6.51), 可得

由于

所以共轭梯度法比最速下降法的收敛性好得多.

由于 A 是对称正定的, 故可记 $\mu = \text{cond}_{2}(A) = \lambda_{1}/\lambda_{n}$ , 这时误差估计式 (6.51) 可写成

从式 (6.52) 可以看出, 当系数矩阵 A 的条件数很大时, 共轭梯度法的收敛速度可能很慢. 条件数较小时, 收敛很快.

理论上, 如果所有计算都没有舍入误差的话, 利用定理 6.3.7 可以证明共轭梯度法最多经过 $n$ 次迭代即可得到精确解 $\pmb{x}^*$ . 但是, 由于计算过程中不可避免要有舍入误差, 所以实际上是将共轭梯度法作为一种迭代法使用.

例 6.3.9 试用共轭梯度法求解下述线性方程组,

其中初值为 $\boldsymbol{x}^{(0)}=(0,0,0)^{\mathrm{T}}$ ，使得最终迭代误差 $\boldsymbol{r}^{(k)}=\boldsymbol{b}-\boldsymbol{A}\boldsymbol{x}^{(k)}$ 达到 $\|\boldsymbol{r}^{(k)}\|_{2}/\|\boldsymbol{r}^{(0)}\|_{2}<1\times10^{-4}$ ，最大迭代步数设为 10.

解：设 $\pmb{x}^{(0)} = (0,0,0)^{\mathrm{T}}$ ，则

第一步迭代如下:

我们发现 $\| \pmb{r}^{(1)}\| _2 / \| \pmb{r}^{(0)}\| _2 = 0.3094 > 1\times 10^{-4}$ 因此进行第二步迭代如下：

此时 $\| \pmb{r}^{(2)}\| _2 / \| \pmb{r}^{(0)}\| _2 = 7.04\times 10^{-5} < 1\times 10^{-4}$ ，满足迭代收敛条件，因此迭代终止.综上，总共需要两步迭代达到 $\| \pmb{r}^{(k)}\| _2 / \| \pmb{r}^{(0)}\| _2 < 1\times 10^{-4}$ ，其中迭代解为(1.0001,1.0001,1.0001)T.

6.3.3 广义极小残量法#

本节中我们介绍求解非对称线性方程组的一类算法: 广义极小残量法 (Gerneral Minimal RESidual). 这个算法自 20 世纪 80 年代以来, 在许多研究者的努力下得到很大的完善, 已经成为当前求解大型稀疏非对称线性方程组的主要手段. 本节中的范数 ||·|| 均为 2- 范数.

设所求线性方程组为

取 $\boldsymbol{x}^{(0)}\in\mathbb{R}^{n}$ 为任一向量, 令 $\boldsymbol{x}=\boldsymbol{x}^{(0)}+\boldsymbol{z}$ , 则上式等价于

其中 $\boldsymbol{r}^{(0)} = \boldsymbol{b} - \boldsymbol{A}\boldsymbol{x}^{(0)}$ . 下面我们讨论方程 (6.54) 的求解问题.

我们从 $\boldsymbol{r}^{(0)}$ 开始, 构造一组相互正交且范数为 1 的向量 $\boldsymbol{v}^{(1)}, \boldsymbol{v}^{(2)}, \cdots, \boldsymbol{v}^{(m)}$ 如下.

首先, 为了满足 $\| \pmb{v}^{(1)} \| = 1$ , 令 $\beta = \| \pmb{r}^{(0)} \|$ , 则 $\pmb{v}^{(1)} = \pmb{r}^{(0)} / \beta$ .

其次, 计算 $\boldsymbol{w} = \boldsymbol{A}\boldsymbol{v}^{(1)}$ . 为了构造 $\boldsymbol{v}^{(2)}$ 与 $\boldsymbol{v}^{(1)}$ 正交, 令 $h_{11} = (\boldsymbol{w}, \boldsymbol{v}^{(1)})$ 且 $\widetilde{\boldsymbol{w}} = \boldsymbol{w} - h_{11}\boldsymbol{v}^{(1)}$ , 则

令 $h_{21} = \|\widetilde{w}\|$ ，则 $\boldsymbol{v}^{(2)} = \widetilde{\boldsymbol{w}} / h_{21}$ 满足 $(\boldsymbol{v}^{(2)}, \boldsymbol{v}^{(1)}) = 0$ 且 $\|\boldsymbol{v}^{(2)}\| = 1$ .

类似地, 计算 $\boldsymbol{v}^{(3)}$ 的过程如下:

(1) 计算 $w = A v^{(2)}$ ,

(2) 令 $h_{12} = (\pmb{w}, \pmb{v}^{(1)})$ , $h_{22} = (\pmb{w}, \pmb{v}^{(2)})$ ,

(3) 计算 $\widetilde{w} = w - h_{12}v^{(1)} - h_{22}v^{(2)}$ ,

(4) 令 $h_{32} = \|\widetilde{\pmb{w}}\|$ ,

(5) 计算 $\boldsymbol{v}^{(3)} = \widetilde{\boldsymbol{w}} / h_{32}$ ,

像这样继续下去, 经过 m 步, 我们就可以得到 $\boldsymbol{V}_{m} = (\boldsymbol{v}^{(1)}, \boldsymbol{v}^{(2)}, \cdots, \boldsymbol{v}^{(m)})$ . 回顾上述过程, 我们发现

记 $V_{m}=(\boldsymbol{v}^{(1)},\boldsymbol{v}^{(2)},\cdots,\boldsymbol{v}^{(m)})$ ，上述性质可以写成

其中，

为一个上海森伯格矩阵.

或者

其中，

假设存在一个 $\pmb{y} \in \mathbb{R}^n$ 满足 $z = V_{m}\pmb{y}$ , 我们可以看到

其中, $e_{1}$ 是第一个分量为 1, 其余分量为 0 的 $m+1$ 维列向量.

由于 $V_{m+1}^{T}V_{m+1}=I$ ，所以

于是, 极小化 $\|r_{0}-Az\|_{2}$ 就相当于极小化 $\|\beta e_{1}-\overline{H}_{m}y\|_{2}$ , 而后者通过 QR 分解很容易求解.

我们可以把广义极小残量法详细叙述如下:

算法 6.3.10 (广义极小残量法)

(1) 选定初值 $x^{(0)}$ , 设 $r^{(0)} = b - Ax, \beta = \| r^{(0)} \|_2$ 和 $v^{(1)} = r^{(0)} / \beta$ ;

(2) 对 $j = 1, 2, \cdots$ 直到收敛, 计算步骤 (3):

(3) $\pmb{w}_j = \pmb{A}\pmb{v}^{(j)}$ ;

对 $i=1,2,\cdots,j$ ，计算:

如果 $h_{j + 1,j}$ 足够小，

退出算法;否则
 $\boldsymbol{v}^{(j+1)}=\boldsymbol{w}^{(j)}/h_{j+1,j};$

(4) 求解最小二乘问题 $\min_{y\in\mathbb{R}^{k}}\|\beta e_{1}-\overline{H}_{k}y\|_{2}$ ，得到的解为 $y_{k}$ ;

(5) 计算解为 $\boldsymbol{x}_{k} = \boldsymbol{x}^{(0)} + \boldsymbol{V}_{k} \boldsymbol{y}_{k}$ .

理论上讲, 如果 $\{A^i r^{(0)}\}_{i=0}^{n-1}$ 线性无关, 当 $m=n$ 时, 在精确计算的条件下, GMRES 算法应当给出方程组 (6.53) 的准确解. Matlab 程序 gmres.m 如下:

function [x,iter] = gmres(A,b,tol)
    bnrm = norm( b );
    x = zeros(size(b));
    [n,n] = size(A);
    V(1:n,1:n+1) = zeros(n,n+1);
    H(1:n+1,1:n) = zeros(n+1,n);
    cs(1:n) = zeros(n,1);
    sn(1:n) = zeros(n,1);
    e1    = zeros(n,1);
    e1(1) = 1.0;
    for iter = 1:500
    r = b-A*x;
    V(:,1) = r / norm( r );
    s = norm( r )*e1;
    for i = 1:n,
    w = A*V(:,i);
    for k = 1:i,
    H(k,i) = w'*V(:,k);
    w = w - H(k,i)*V(:,k);
    end
    H(i+1,i) = norm( w );
    V(:,i+1) = w / H(i+1,i);
    for k = 1:i-1,
    temp    = cs(k)*H(k,i) + sn(k)*H(k+1,i);
    H(k+1,i) = -sn(k)*H(k,i) + cs(k)*H(k+1,i);
    H(k,i)    = temp;
    end
    cs(i)   = H(i,i)/\sqrt{H(i,i)^2+ H(i+1,i)^2};
    sn(i)   = H(i+1,i)/\sqrt{H(i,i)^2+ H(i+1,i)^2};
    temp    = cs(i)*s(i);
    s(i+1) = -sn(i)*s(i);
    s(i)    = temp;

H(i,i) = cs(i)*H(i,i) + sn(i)*H(i+1,i);
H(i+1,i) = 0.0;
error = abs(s(i+1)) / bnrm;
if (error < tol),
    y = H(1:i,1:i) \ s(1:i);
    x = x + V(:,1:i)*y;
    break;
end
end
y = H(1:m,1:m) \ s(1:m);
x = x + V(:,1:m)*y;

但是, 不难看出, 当 $m$ 很大时, 计算中需要保存所有的 $\{\pmb{v}^{(i)}\}_{i=1}^{m}$ . 对于大规模问题, 这对存储空间的要求比较高, 而且计算量也随之大大增加. 因此, 在实际使用中往往采用重新启动的 GMRES 算法. 采用重新启动的好处是可以节省计算量和存储空间, 但是理论上却不能保证算法的收敛性.

下面我们不加证明地给出广义极小残量法的收敛性定理.

例 6.3.12 试用广义极小残量法法求解下述线性方程组

其中, 初值为 $\boldsymbol{x}^{(0)}=(0,0,0)^{\mathrm{T}}$ , 使得迭代第 k 步的 $h_{k+1,k}$ 满足 $|h_{k+1,k}|<1\times10^{-4}$ , 最大迭代步数设为 10.

解：设初始解为 $\boldsymbol{x}^{(0)}=(0,0,0)^{\mathrm{T}}$ ，计算 $\boldsymbol{r}^{(0)}=(2.8,2,3.2)^{\mathrm{T}},\beta=4.6989,\boldsymbol{v}^{(1)}=(0.5959,0.4256,0.6810)^{\mathrm{T}}$ .

第一步迭代如下:

我们发现 $h_{21}=0.932\ 5>1\times10^{-4}$ . 因此进行第二步迭代, 如下:

我们发现 $h_{32} = 0.5194 > 1 \times 10^{-4}$ . 因此进行第三步迭代, 如下:

此时 $h_{43}=0.000\times10^{-16}<1\times10^{-4}$ ，满足迭代收敛条件，迭代终止。此时， $4\times3$ 阶海森伯格矩阵为

求解 $\min_{y\in R^3}\| \beta e^{(1)} - \overline{H}_{43}y\| _2$ 相当于求解一个上三角矩阵，我们得到 $\pmb{y} = (1.7025, - 0.3159,$ $0.0404)^{\mathrm{T}}$ 因此，最终迭代近似解为 $\pmb{x} = 1.7025\pmb{v}_1 - 0.3159\pmb{v}_2 + 0.0404\pmb{v}_3 = (1.0000,1.0000,$ $1.0000)^{\mathrm{T}}$

6.3.4 预处理技术#

由于存在浮点运算的误差, 共轭梯度法和广义极小残量法计算得到的向量会逐渐失去正交性, 因而都不能在 $n$ 步之内得到原方程的精确解. 况且, 遇到求解大规模的线性代数方程组, 即使能够在 $n$ 步收敛的话, 这个收敛速度也不能令人满意. 预处理技术能有效地改善收敛性质并加快收敛速度, 因而在实际使用中应用得非常广泛.

预处理技术从广义上来说可以指对原方程组进行的任何显示的或者隐式的修正, 使得该方程组通过迭代法更容易求解.

简单地说, 花比较小的代价找到一个矩阵 M, 然后用迭代法求解如下的同解线性方程组

或

新得到的算法分别称为左预处理或者右预处理的迭代方法.

特别地, 如果存在矩阵 L 使得 $M = LL^{T}$ , 那么我们可以求解如下的同解线性方程组:

这称为对称预处理方法, 特别在预处理共轭梯度方法时经常使用.

这里, $M$ 称为预处理矩阵. 一个好的预处理矩阵至少能够满足如下的条件:

(1) 构造 M 的代价很小;

(2) M 跟 A 足够接近;

(3) 关于 $M^{-1}$ 的线性方程组很容易求解.

特别地, 当系数矩阵 $\mathbf{A}$ 的对角元非零时, 取 $M = \mathrm{diag}(A)$ , 我们可以得到一个最常用的预处理矩阵.

这里, 我们给出左预处理共轭梯度法的具体计算过程.

算法 6.3.13 (左预处理共轭梯度法)

(1) 选定初值 $x^{(0)}$ , 设 $\pmb{p}^{(0)} = \pmb{r}^{(0)} = M^{-1}(\pmb{b} - \pmb{A}\pmb{x})$ ;

(2) 对 $k = 0, 1, 2, \cdots$ 直到收敛, 计算:

如果 $\| r^{(k + 1)}\|$ 足够小，

退出算法;

否则，

一般认为, 如果预处理后的系数矩阵 $M^{-1}A$ 的特征值更加聚集的话, 不管是用共轭梯度法还是用广义极小残量法都会得到更好的收敛效果.

评 注#

本章主要讨论了解线性方程组的各种迭代解法。迭代法是计算技术中常用的一种有效方法。使用各种迭代格式时，最主要的就是判断它的收敛性以及了解收敛速度。当然，在实际计算中，对一种迭代格式，不必事先判断了收敛性才敢使用，它完全可以在计算过程中判断是否收敛，雅可比迭代与GS迭代的收敛域并不互相包含，所以不能相互代替，但当两者皆收敛时，一般来说GS迭代比雅可比迭代的收敛速度快。实用中更多的是使用SOR迭代，选择松弛因子有赖于实际经验，对于一类具有相容次序的矩阵，可以得到最佳的松弛因子。共轭梯度法(简称CG法)是求解系数矩阵为对称正定的线性方程组的非常有效的方法。近20年来有关的研究得到了前所未有的发展，目前有关的方法和理论已经相当成熟，并且已经成为求解大型系数线性方程组最受欢迎的一类方法. 当 Ax = b 为病态方程组时, $\text{cond}(A)$ 很大, 共轭梯度法收敛缓慢. 这时, 可以使用预条件共轭梯度法来计算, 往往其收敛速度大大提高.

习题六#

1. 用迭代法求解下述线性方程组

(1) 分别写出雅可比迭代、GS 迭代、SOR 迭代 (取 $\omega = 1.35$ ) 的迭代格式；

(2) 判别上述三个迭代格式的收敛性, 并说明理由;

(3) 用收敛的迭代格式分别计算方程组的解, 要求满足

2. 分别用雅可比迭代和 GS 迭代求解下述线性方程组:

取初值 $\boldsymbol{x}^{(0)}=(0,0,0)^{\mathrm{T}}$ ，精确到小数点后四位，并在理论上判别这两个迭代格式的收敛性.

3. 对下述线性方程组:

分别讨论用雅可比迭代和 GS 迭代的收敛性.

4. 试证对于 $n$ 维向量 $x$ 有如下关系式成立:

5. 证明：对于矩阵范数, 如果 $\|A\| < 1$ , 则

6. 设线性方程组 Ax = b, 其中, A 为 n 阶对称正定矩阵 (设 A 的特征值满足 $0 < \alpha \leqslant \lambda(A) \leqslant \beta$ ), 建立如下迭代格式:

证明：当 $0 < \omega < \frac{2}{\beta}$ 时，上述迭代法收敛.

7. 试证当 $-0.5 < \alpha < 1$ 时矩阵

是正定的. 当 $-0.5 < \alpha < 0.5$ 时, 用雅可比迭代求解 Ax = b 是收敛的.

8. 证明对 GS 迭代有

其中， $\mu=\max_{i}(\sum_{j=1}^{n}|b_{ij}|/(1-\sum_{j=1}^{i-1}|b_{ij}|))$ .

9. 试用共轭梯度法求解下述线性方程组,

其中，初值为 $\boldsymbol{x}^{(0)} = (0, 0, 0, 0, 0, 0)^{\mathrm{T}}$ ，使得最终迭代误差 $\boldsymbol{r}^{(k)} = \boldsymbol{b} - \boldsymbol{A}\boldsymbol{x}^{(k)}$ 达到 $\|r^{(k)}\|_{2}/\|r^{(0)}\|_{2} < 1 \times 10^{-4}$ ，最大迭代步数设为 10.

10. 试用广义极小残量法法求解线性方程组

其中, 初值为 $\boldsymbol{x}^{(0)} = [0, 0, 0, 0, 0, 0]^{\mathrm{T}}$ , 使得最终迭代误差 $\boldsymbol{r}^{(k)} = \boldsymbol{b} - \boldsymbol{A}\boldsymbol{x}^{(k)}$ 达到 $\|r^{(k)}\|_{2}/\|r^{(0)}\|_{2} < 1 \times 10^{-4}$ , 最大迭代步数设为 10.

数值实验六#

试用SOR迭代计算线性方程组

取 $\boldsymbol{x}^{(0)}=(0,0,0)^{\mathrm{T}}$ ，松弛因子分别选取为 $\omega=0.1t,1\leqslant t\leqslant19$ ，要求达到精度 $\|\boldsymbol{x}^{(k+1)}-\boldsymbol{x}^{(k)}\|\leqslant10^{-4}$ 。试通过数值计算得出不同的松弛因子所需要的迭代次数和收敛最快的松弛因子，并指出哪些松弛因子使得迭代发散。

第 7 章 非线性方程求根#