数值分析读书笔记（3）求解线性代数方程组的迭代法

数值分析读书笔记（3）求解线性代数方程组的迭代法

1.基本迭代法及其构造

考虑方程组Ax=b，其中A属于n*n维的矩阵空间，b和x属于n维向量空间，一般来说我们需要从这个隐式的方程组转变成显示的等价方程，一般具有形式

x=ϕ(x) $$x=\varphi (x)$$

,这样的方程为不动点方程，我们可以通过不断迭代，计算出等式右端然后赋值给变量x

对于Ax=b而言，如果我们简单取A=I-B，可以得到等价的x=Bx+b,从而构造迭代格式

x(k+1)=Bx(k)+b $$x^{(k+1)}=B{x}^{(k)}+b$$

需要注意的是，迭代法不一定会是收敛的，也就是说x不一定会收敛到某个值，这样并不是我们所希望的，故后面会讨论一下迭代法的收敛性，我们先来谈谈迭代法的构造，从迭代格式中可以看到，我们对矩阵A进行了一次分裂，分裂的格式显然不是唯一的，我们将采取A=M-N这样一种分裂方法，可以得到线性的不动点方程

x=Bx+f,B=M−1N=I−M−1A,f=M−1b $$x=Bx+f,B=M^{-1}N=I-M^{-1}A,f=M^{-1}b$$

不同的分裂法将使我们的B和f不一致，我们通常将A分裂成三个部分A=D-L-U，其中D为对角矩阵，L和U分别是下三角和上三角矩阵，这里需要注意一下符号的不同，L和U的元素都是原矩阵上三角和下三角元素符号的相反数

下面通过这样的一种分裂方式，我们介绍几种迭代形式，首先是Jacobi迭代法（同步迭代）

注意到线性迭代中的M和N，在Jacobi迭代中，我们令M=D, N=L+U，构造出迭代格式，即

AxDxxx(k+1)x(k+1)=((D)−(L+U))x=b=(L+U)x+b=D−1(L+U)x+D−1b=D−1(L+U)x(k)+D−1b=BJx(k)+fJ(1)(2)(3)(4)(5) $$\begin{array}{}\text{(1)}& Ax& =((D)-(L+U))x=b\text{(2)}& Dx& =(L+U)x+b\text{(3)}& x& =D^{-1}(L+U)x+D^{-1}b\text{(4)}& x^{(k+1)}& =D^{-1}(L+U)x^{(k)}+D^{-1}b\text{(5)}& x^{(k+1)}& =B_J{x}^{(k)}+f_J\end{array}$$

直观上来看Jacobi迭代，就是把方程n行对应的x保留，其余维度的x移到方程的左端，用这n维的左端的式子来迭代更新n个维度的x

那么这样看就可以理解Jacobi迭代为什么是同步迭代了，因为所有的维度的x必须全部更新完成之后才可以进行下一步的迭代，由此我们介绍下一种迭代，Gauss-Seidel迭代（异步迭代）

在Gauss-Seidel迭代中，我们和上面对A进行同样的分裂方式，只不过在M，N的选取上做出了改变，我们令M=D-L，N=U，构造出迭代格式，即

Ax(D−L)xxx(k+1)x(k+1)=((D−L)−(U))x=b=Ux+b=(D−L)−1Ux+(D−L)−1b=(D−L)−1Ux(k)+(D−L)−1b=BGx(k)+fG(6)(7)(8)(9)(10) $$\begin{array}{}\text{(6)}& Ax& =((D-L)-(U))x=b\text{(7)}& (D-L)x& =Ux+b\text{(8)}& x& =(D-L{)}^{-1}Ux+(D-L{)}^{-1}b\text{(9)}& x^{(k+1)}& =(D-L{)}^{-1}U{x}^{(k)}+(D-L{)}^{-1}b\text{(10)}& x^{(k+1)}& =B_G{x}^{(k)}+f_G\end{array}$$

直观来看Gauss-Seidel迭代，和Jacobi一样就是把方程n行对应的第n个x保留，其余的x移到方程的左端，只不过在我们更新第k个的时候会利用前面迭代更新完成了的前k-1个x进行带入计算后面的n-k个x仍然使用初始值，也就是一种异步的思想

在实际中，我们使用Jacobi迭代或者是Gauss-Seidel迭代都可能会出现不收敛或者收敛速度比较慢这样的情况，我们是不是可以试着去构造一种带参数的迭代方法，这样我们就能够利用参数的调整来实现对整个迭代方法调整已增加其收敛速度等，下面介绍SOR(Successive over relaxation method)迭代

SOR迭代的基本思想是，在以求出第k个x值的基础上用Gauss-Seidel解出第k+1个x值，然后利用这k+1个x的值和第k个x的值的线性组合来的出更好的近似解

同样的我们利用和Jacobi，Gauss-Seidel迭代方法一样的分裂方法，对A进行分裂，注意到一种矩阵的记法

Dx(k+1)=(1−ω)Dx(k)+ω(b+Lx(k+1)+Ux(k))=(1−ω)Dx(k)+ωb+ωLx(k+1)+ωUx(k)(11)(12) $$\begin{array}{}\text{(11)}& D{x}^{(k+1)}& =(1-\omega )D{x}^{(k)}+\omega (b+L{x}^{(k+1)}+U{x}^{(k)})\text{(12)}& & =(1-\omega )D{x}^{(k)}+\omega b+\omega L{x}^{(k+1)}+\omega U{x}^{(k)}\end{array}$$

注意到等式右端的第一部分，即 (1−ω)Dx(k) $(1-\omega )D{x}^{(k)}$,这一部分其实就是之前的上一步的值

等式的第二部分就是类似Gauss-Seidel迭代得到的解然后乘上一个松弛因子（ ω $\omega$）

通过对上式进行化简，我们可以简单的得到

(D−ωL)x(k+1)=(1−ω)Dx(k)+ωb+ωUx(k)=((1−ωD)−ωU)x(k)+ωb(13)(14) $$\begin{array}{}\text{(13)}& (D-\omega L)x^{(k+1)}& =(1-\omega )D{x}^{(k)}+\omega b+\omega U{x}^{(k)}\text{(14)}& & =((1-\omega D)-\omega U)x^{(k)}+\omega b\end{array}$$

对于 D−ωL $D-\omega L$这个矩阵我们知道，当 ω $\omega$取任意值的时候，总为非奇异矩阵，也就是所总存在逆矩阵，所以上式可以继续化简（证明过程只需将矩阵的形式画出来，由于对角线元素不等于0，即可得证）

得到 x(k+1)=(D−ωL)−1((1−ωD)−ωU)x(k)+ω(D−ωL)−1b $x^{(k+1)}=(D-\omega L{)}^{-1}((1-\omega D)-\omega U)x^{(k)}+\omega (D-\omega L{)}^{-1}b$

所以我们可以根据上式建立起来SOR的迭代格式

x(k+1)=Bωx(k)+fω $$x^{(k+1)}=B_{\omega }{x}^{(k)}+f_{\omega }$$

其中， Bω=(D−ωL)−1((1−ωD)−ωU),fω=ω(D−ωL)−1b $B_{\omega }=(D-\omega L{)}^{-1}((1-\omega D)-\omega U),f_{\omega }=\omega (D-\omega L{)}^{-1}b$

注意到SOR中的Guass-Seidel迭代也有区分向前或者向后Gauss-Seidel迭代，由此可以引申出SSOR（Sysmetrical Successive Over Relaxation method），此外还有Richardson迭代以及块迭代的方法，这里暂时不做出说明，以后补充

2.迭代法的收敛性分析

我们利用上面的知识可以构造迭代格式来进行迭代，但是这里就出现一个问题，所有的迭代方法都是可行的吗？这里指的可行可能是收不收敛，收敛的快慢等等，所以我们有必要对其进行规范化，系统化的分析

先来看看一般的迭代格式

x(k+1)=Bx(k)+f $$x^{(k+1)}=B{x}^{(k)}+f$$

我们可以从n=0一直写到n=k+1，这样来看就有

x(k+1)=B(k+1)x(0)+(Bk+⋯+B+I)f $$x^{(k+1)}=B^{(k+1)}{x}^{(0)}+(B^k+\cdots +B+I)f$$

我们必须研究矩阵 B $B$是否会变的越来越大，或者趋近于一个常数矩阵 B∗ $B^{\ast }$,以及 x $x$是否去趋近于 x∗ $x^{\ast }$

这里引入误差向量

e(k+1)=x(k+1)−x∗ $$e^{(k+1)}=x^{(k+1)}-x^{\ast }$$

只要误差向量趋近于0即可收敛，这里做个简单的计算

e(k+1)=x(k+1)−x∗=(Bx(k)+f)−(Bx∗+f)=B(x(k)−x∗)=Be(k)=B(k+1)(x(0)−x∗) $$e^{(k+1)}=x^{(k+1)}-x^{\ast }=(B{x}^{(k)}+f)-(B{x}^{\ast }+f)=B(x^{(k)}-x^{\ast })=B{e}^{(k)}=B^{(k+1)}(x^{(0)}-x^{\ast })$$

下面给出收敛定理

设 B=(bij)n×n $B=(b_{ij}{)}_{n\times n}$, 则 Bk→O(k→∞) $B^k\to O(k\to \mathrm{\infty })$的充要条件是 B $B$的谱半径满足

ρ(B)<1 $$\rho (B)<1$$

这里直接给出这个收敛定理，并没有给出证明，具体的证明思路为利用Jordan标准型以及极限关系可以证明，后续等待补充

继续不加证明地给出一个定理

设方程组 Ax=b $Ax=b$的不动点方程组为 x=Bx+f $x=Bx+f$,则对于任意初始近似向量 x(0) $x^{(0)}$与任意常数向量 f $f$,求解 Ax=b $Ax=b$的基本迭代法 x(k+1)=Bx(k)+f $x^{(k+1)}=B{x}^{(k)}+f$收敛的充要条件为

ρ(B)<1 $$\rho (B)<1$$

一般来说我们通过上述几个定理就可以通过计算矩阵的谱半径，也就是先计算最大的特征值，然后判断其是否大于1来判别一个基本迭代格式是否是收敛的，但是对于维度等级十分高的矩阵，我们可能得寻求一种较为简单的判别方法，下面我们不加证明地给出利用范数来判别的一个定理

求解 Ax=b $Ax=b$的基本迭代法 x(k+1)=Bx(k)+f $x^{(k+1)}=B{x}^{(k)}+f$收敛的充分条件为

∥B∥<1 $$\Vert \begin{array}{c}B\end{array}\Vert <1$$

其中 ∥⋅∥ $\Vert \begin{array}{c}\cdot \end{array}\Vert$为任意一种矩阵范数

3.误差估计

对于迭代格式的收敛性我们已经讨论过了，下面给出误差的估计，主要是用来计算相应到达误差范围相应迭代次数的值，下面给出一个定理

设求解 Ax=b $Ax=b$ 的基本迭代法为 x(k+1)=Bx(k)+f $x^{(k+1)}=B{x}^{(k)}+f$, x(0) $x^{(0)}$为初始迭代向量,且迭代矩阵的一矩阵范数 ∥B∥=q<1 $\Vert \begin{array}{c}B\end{array}\Vert =q<1$, 则

∥x(k)−x∗∥≤∥B∥1−∥B∥∥x(k)−x(k−1)∥ $$\Vert \begin{array}{c}{x}^{(k)}-x^{\ast }\end{array}\Vert \le \frac{\Vert \begin{array}{c}B\end{array}\Vert }{1-\Vert \begin{array}{c}B\end{array}\Vert }\Vert \begin{array}{c}{x}^{(k)}-x^{(k-1)}\end{array}\Vert$$

∥x(k)−x∗∥≤∥B∥k1−∥B∥∥x(1)−x(0)∥ $$\Vert \begin{array}{c}{x}^{(k)}-x^{\ast }\end{array}\Vert \le \frac{{\Vert \begin{array}{c}B\end{array}\Vert }^k}{1-\Vert \begin{array}{c}B\end{array}\Vert }\Vert \begin{array}{c}{x}^{(1)}-x^{(0)}\end{array}\Vert$$

该定理证明可以利用之前所介绍的Banach引理来证明

用上面的式子，可以求解出来精度 ∥e(k)∥≤ϵ $\Vert \begin{array}{c}{e}^{(k)}\end{array}\Vert \le ϵ$的迭代步数，令步数为k，B的范数为q，则有

k≥ϵ(1−q)∥x(0)−x(1)∥|ln(q)| $$k\ge \frac{ϵ(1-q)}{\Vert \begin{array}{c}{x}^{(0)}-x^{(1)}\end{array}\Vert \left|\begin{array}{c}ln(q)\end{array}\right|}$$

4.关于几类特殊矩阵进行迭代法的收敛性

首先我们来看看对角占优矩阵进行基本迭代法的收敛性

对角占优矩阵可以简单的划分成严格对角占优和弱对角占优

严格对角占优是指对角线上的元素的绝对值比相同行其他元素的绝对值的和都大，这里不存在等号的条件
弱对角占优是严格对角占优的基础上添加等号的条件,也就是说对角线上的元素的绝对值大于等于相同行其他元素的绝对值的和

我们直接不加证明的给出一个定理：

对于严格对角占优矩阵和弱对角占优矩阵，我们使用任意的初始向量，构造Jacobi迭代格式或者Gauss-Seidel迭代格式，结果均收敛

接下来，我们看一下另外一种特殊矩阵，即对称正定矩阵

给出一个定理

设A对称，且对角元素为正，则方程组Ax=b的Jacobi迭代格式收敛的充分必要条件就是A和2D-A均正定，其中D为A生成的对角矩阵

同样的，我们对于Gauss-Seidel迭代也给出一个定理

设A对称正定，则方程组Ax=b的Gauss-Seidel迭代收敛

对SOR的迭代格式，比较复杂，对于某些特定的矩阵有着一定规律，给出一个定理

对任意的 A∈Rn×n $A\in R^{n\times n}$,如果它的对角元素皆非0，则SOR迭代的迭代矩阵 Bω $B_{\omega }$的谱半径 ρ(Bω) $\rho (B_{\omega })$与松弛系数 ω $\omega$有着下列关系

ρ(Bω)≥|(1−ω)| $$\rho (B_{\omega })\ge \left|\begin{array}{c}(1-\omega )\end{array}\right|$$

由上述定理可以推出，方程组使用SOR方法收敛的一个必要条件

0<ω<2 $$0<\omega <2$$

反过来，也有一个定理

设 A∈Rn×n $A\in R^{n\times n}$，且 A $A$对称正定，如果 0<ω<2 $0<\omega <2$, 则求解 Ax=b $Ax=b$的SOR迭代格式收敛