迭代法01--定常迭代法

1. 简介

用于求解线性方程的迭代法可分为两类：定常迭代法（stationary iterative method）和_Krylov_法。

1. 定常迭代法包括：

• Jacobi 迭代
• Gauss-Seidel 迭代
• Successive Over-Relaxation
• ……

2. Krylov法包括：

• Conjugate Gradient
• GMRES
• BiCGStab
• ……

3. 区别：

• 定常迭代法相对古老，容易了解与实现，但效果通常不好。
• Krylov法相对年轻，虽然较不易理解分析，但效果普遍优异。

本文主要介绍定常迭代法：通式、收敛性分析以及常用定常迭代方法，Krylov法以后有机会再介绍。

2. 定常迭代法通式及其收敛性

2.1 定常迭代法的通式

考虑线性方程$Ax=b$，其中$A$是一$n\times n$阶系数矩阵，$b$是一$n$维常数向量。分解$A=M-N$，其中$M$可逆。另$B=M^{-1}N$，$c=M^{-1}b$，定常迭代法的通式如下：
$$x^{(k)}=B{x}^{(k-1)}+c,k=1,2,...$$
称$B$为迭代矩阵(iteration matrix),之所以称这种形式的迭代方法为定常迭代法(stationary iterative methods)是因为从$x_k$向$x_{k+1}$过渡不依赖于迭代历史。
##2.2 收敛性分析
定常迭代法的收敛性完全由迭代矩阵$B$的特征值决定。令$\sigma (B)$为$B$的所有相异特征值所形成的集合，称为矩阵谱（spectrum），并令$\rho (B)$为$B$的绝对值最大的特征值，即，$\rho (B)=ma{x}_{\lambda \in \sigma (B)}|\lambda |$，称为普半径。
###2.2.1 定理一

定理一：若$A$是一$n\times n$阶矩阵，则以下称述是等价的：

• (1) $\rho (A)<1$
• (2) ${\lim }_{k\to \infty }{A}^k=0$
• (3) ${\sum }_{k=0}^{\infty }{A}^k$，且$(I-A{)}^{-1}={\sum }_{k=0}^{\infty }{A}^k$

本文只需要用到$(1)\Rightarrow (2)$和$(2)\Rightarrow (3)$，下面给出证明：

1. $(1)\Rightarrow (2)$

若$A$是可对角化矩阵，即，$A=S\Lambda S^{-1}$，$\Lambda$为特征值构成的对角矩阵，则$A^k=S{\Lambda }^k{S}^{-1}$，其中，
\begin{equation}
{\left[ {\begin{array}{*{20}c}
{\lambda _1^k } & {} & {} & {} & {} \
{} & {\lambda _2 ^k} & {} & {} & {}\
{} & {} & \ddots & {} & {}\
{} & {} & {} & {\lambda _n^k }
\end{array} }
\right ]}
\end{equation}
若每一特征值都满足$|{\lambda }_i|<1$，当$k\to \infty$时，${\lambda }_i^k\to 0$，即知${\Lambda }^k\to 0$,也就有$A^k=S{\Lambda }^k{S}^{-1}\to 0$

若$A$是不可对角化矩阵，此性质仍然成立。本文不再证明。

2. $(2)\Rightarrow (3)$
　考虑，
$$(I-A)(I+A+A^2+\cdots +A^{k-1})=I-A^k$$
当$k\to \infty$，若$A^k\to 0$，则$(I-A{)}^{-1}={\sum }_{k=0}^{\infty }{A}^k$

2.2.2 收敛性说明

定理二：若$\rho (B)<1$，则，

• (1) $B$ 为一收敛矩阵;
• (2) $A^{-1}$ 存在;
• (3) 对于任何一个初始向量$x^{(0)}$，
  $$\underset{k\to \infty }{\lim }{x}^{(k)}=x=A^{-1}b$$

证明：

1. (1)已由定理一证明。

2. 由定理一知，若$\rho (B)<1$，则，$(I-B{)}^{-1}={\sum }_{k=0}^{\infty }{B}^k$，写出$$A=M-N=M(I-B)$$，可知$A$是一可逆矩阵。即证(2).

3. 逐次迭代可得：
   $$x^{(k)}=B^k{x}^{(0)}+(I+B+B^2+\cdots +B^{k-1})c$$
   由定理一知，若$\rho (B)<1$，则${\lim }_{k\to \infty }{B}^k=0$，则，
   $$\underset{k\to \infty }{\lim }{x}^{(k)}=(I-B{)}^{-1}c=(I-B)-1M^{-1}b=A^{-1}b=x$$

收敛速度：
令${ϵ}^{(k)}=x^{(k)}-x$表示第$k$次迭代后的估计误差，则，
$$|\frac{{ϵ}^k}{{ϵ}^{k-1}}|=\rho (B)$$
故每次迭代的精确度该善良是$-lo{g}_{10}\rho (B)$。因此，我们定义渐进收敛速率(asymptotic rate of convergence)为$-ln\rho (B)$。

2.2.3 总结

综上：定常迭代法的设计要领是在容易计算$B=M^{-1}N$和$c=M^{-1}b$的前提下，找出迭代矩阵$B$使其普半径$\rho (B)$越小越好。

3. 常用定常迭代法

以下各种方法的区别仅在于对系数矩阵$A$分解(splitting)的不同。

3.1 Jacobi 法

分解$A=D-N$，其中，$D$是$A$的主对角部分，$-N$是$A$除去$D$后剩余的部分。假设每一$-a_{ii}̸=0$，即，$D$可逆，则，Jacobi法的迭代公式为：
$$x^{(k)}=D^{-1}N{x}^{(k-1)}+D^{-1}b$$
显然，$B=D^{-1}N$和$c=D^{-1}b$仅耗费很少量的计算。
分量迭代公式为：
$$x_i^{(k)}=(b_i-\sum _{j̸=i}{a}_{ij}{x}_j^{(k-1)})/a_{ii}$$

收敛条件：
若$A$为对角占优矩阵，则对于任意的$x^{(0)}$和$b$，Jacobi法必定收敛。

3.2 Gauss-Seidel 法

分解$A=(D-L)-U$，其中$D$是$A$的主对角部分，$-L$和$-U$分别是$A$的严格下三角矩阵和严格上三角矩阵，则，Gauss-Seidel 的迭代公式为：
$$x^{(k)}=(D-L{)}^{-1}U{x}^{(k-1)}+(D-L{)}^{-1}b$$
即，$B=(D-L{)}^{-1}U$，$c=(D_L{)}^{-1}b$。

注意，在实际求解中，求解矩阵的逆$(D-L{)}^{-1}$仍然十分耗时，因此会使用如下公式简化求解：
$$(D-L)x^{(k)}=U{x}^{(k-1)}+b$$

分量迭代公式为：
$$x_i^{(k)}=(b_i-\sum _{j<i}{a}_{ij}{x}_j^{(k)}-\sum _{j>i}{a}_{ij}{x}_j^{(k-1)})/a_{ii},i=1,2,\cdots ,n$$

收敛条件：
若A为对角占优矩阵，则对于任意的$x^{(0)}$和$b$，Gauss-Seidel法必定收敛。
证明略。

3.3 Backward Gauss-Seidel & Symmetric Gauss-Seidel

1. Backward Gauss-Seidel
   　与Jacobi迭代不同，Gauss-Seidel法依赖于未知量的顺序。Backward Gauss-Seidel首先更新第$N$个未知量，而不是第一个，即，分解$A=(D-U)-L$
   得到迭代矩阵为：
   $$B_{BGS}=(D-U{)}^{-1}L$$
2. Symmetric Gauss-Seidel
   先进行一次forward Gauss-Seidel iteration，再进行一次backward Gauss-Seidel iteration就称为Symmetirc Gauss-Seidel iteration. 其迭代矩阵为：
   $$B_{SGS}=B_{BGS}{B}_{GS}=(D-U{)}^{-1}L(D-L{)}^{-1}U$$
   这两种方法的迭代公式很容易推出，这里不推导了。

3.４ SOR(Successive Over-Relaxtion)法

改写Gauss-Seidel的分量形式为：
$$x_i^{(k)}=(1-\omega )x_i^{(k-1)}+\omega (b_i-\sum _{j<i}{a}_{ij}{x}_j^{(k)}-\sum _{j>i}{a}_{ij}{x}_j^{(k-1)})/a_{ii},i=1,2,\cdots ,n$$
称为SOR(successive overrelaxation)。其中，$\omega$称为松弛因子。显然当$\omega =1$时，即为Gauss-Seidel法。

以矩阵表达，分解$A=M-N$，其中，$M={\omega }^{-1}D-L,N=({\omega }^{-1}-1)D+U$。如前，$D$是$A$的主对角部分，$-L$和$-U$分别是$A$的严格下三角矩阵，和严格上三角矩阵。因此迭代矩阵，
$$B_{\omega }=M^{-1}N=(I-\omega D^{-1}L{)}^{-1}((1-\omega )I+\omega D^{-1}U)$$
则SOR法的矩阵表达式如下：
$$x^{(k)}=B_{\omega }{x}^{(k-1)}+\omega (I-\omega D^{-1}L{)}^{-1}{D}^{-1}b$$

收敛性：
１. 若$\rho (B_{\omega })<1$，则，$0<\omega <2$。即说明了$\omega$的取值范围。证明略。
２. 若，$\omega >1$，称为过松弛（overrelaxation）;若$\omega <1$，则称为欠松弛。显然松弛因子$\omega$影响SOR法的收敛性，但寻找使$\rho (B_{\omega })$最小化的松弛参数通常是一个相当困难的工作。

4. 总结

1. 本文介绍了定常迭代法的通式，并分析了其收敛性（即，迭代矩阵的普半径$\rho (B)<1$时，算法收敛），并给出了常用定常迭代方法的迭代公式。
2. Jacobi适合并行，收敛速度慢；
3. Gauss-Seidel与迭代的顺序相关，并行性差（针对不同问题，可以实现部分并行，见参考文献3），收敛速度快。
4. 一般而言，除了除少数特例，寻找使 $\rho (B_{\omega })$ 最小化的松弛因子 $\omega$ 是一个相当困难的工作。为了近似最佳的 $\omega$，我们必须另外引进适应性程序 (adaptive) 计算程序，即便如此，SOR 法的表现也常令人失望。1950年诞生的 SOR 法虽然比其他定常迭代法往前迈出了一大步，但终究还是逐渐被二十世纪后半期发展起来的 Krylov 法所取代。

5. 参考文献

1. 线代启示录：
   https://ccjou.wordpress.com/2013/08/22/
2. Iterative Methods for Linear and Nonlinear Equations
   C. T. Kelley North Carolina State University Philadelphia 1995
   https://www.siam.org/books/textbooks/fr16_book.pdf
3. 华东师范大学 线性方程组的迭代解法
   http://math.ecnu.edu.cn/~jypan/Teaching/MatrixComp/ch06_iter_s.pdf