线性方程组的迭代法

经典迭代

张瑞
中国科学技术大学数学科学学院

迭代

在实际应用中，经常会遇上大型的稀疏矩阵。直接法会破坏矩阵的稀疏结构。这样，直接法在应用中，就需要$n^2$的存贮和$n^3$的计算量。 这对于大型的稀疏矩阵非常不适用。

迭代法

与非线性方程的迭代法类似，需要构造一个收敛到解的向量数组。首先，将$Ax=b$变为等价的$x=Gx+g$的形式。

取$A=M-N$，其中$M$可逆。则

\[\begin{aligned} Ax=b \Rightarrow & (M-N)x=b \Rightarrow Mx=Nx+b \\ \Rightarrow & x=M^{-1}Nx+M^{-1}b \end{aligned} \]

然后，取$x^{(k+1)}=Gx^{(k)}+g$得到向量数组。

设线性方程组的根为$x^*$，即

\[Ax^*=b, \quad x^*=Gx^*+g \]

则

\[x^{(k+1)}-x^*=Gx^{(k)}+g-Gx^*-g=G(x^{(k)}-x^*) \]

这样

\[x^{(k)}-x^*=G^k(x^{(0)}-x^*) \]

可以看到，当$G^k\to0$时，$x^{(k)}\to x^*$。

定理 1.
对任意初值序列$x^{(k)}$收敛，当且仅当$G^k\to0,k\to\infty$。

问题. 若$\rho(G)>1$，点列$x^{(k)}$有可能收敛吗？

定义 1.
若矩阵$G$满足$G^k\to0, k\to \infty$，则称$G$为收敛矩阵

定理 2.
矩阵$G$为收敛矩阵，当且仅当$\rho(G)<1$

推论 1.
对任意初值$x^{(0)}$迭代格式$x^{(k+1)}=Gx^{(k)}+g$都收敛的充要条件是$\rho(G)<1$。

推论 2.
若有$\|G\|<1$，则迭代格式$x^{(k+1)}=Gx^{(k)}+g$对任意初值$x^{(0)}$都收敛。

Jacobi 迭代

对于线性方程组

\[\begin{cases} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1 \\ a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n=b_2 \\ \vdots \\ a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n=b_n \\ \end{cases} \]

变换为等价形式

\[\begin{cases} x_1=-\frac1{a_{11}}(+a_{12}x_2+\cdots+a_{1n}x_n-b_1) \\ x_2=-\frac1{a_{22}}( a_{21}x_1+a_{23}x_3+\cdots+a_{2n}x_n-b_2) \\ \vdots \\ x_n=-\frac1{a_{nn}}(a_{n1}x_1+a_{n2}x_2+\cdots+a_{n,n-1}x_{n-1}-b_n) \\ \end{cases} \]

产生迭代格式

\[\begin{cases} x_1^{(k+1)}=-\frac1{a_{11}}(\quad +a_{12}x_2^{(k)}+\cdots+a_{1n}x_n^{(k)}-b_1) \\ x_2^{(k+1)}=-\frac1{a_{22}}( a_{21}x_1^{(k)}+a_{23}x_3^{(k)}+\cdots+a_{2n}x_n^{(k)}-b_2) \\ \vdots \\ x_n^{(k+1)}=-\frac1{a_{nn}}(a_{n1}x_1^{(k)}+\cdots+a_{n,n-1}x_{n-1}^{(k)}-b_n) \\ \end{cases} \]

写成分量形式为

\[x_i^{(k+1)}=-\frac1{a_{ii}}(\sum_{j=1}^{i-1}a_{ij}x_j^{(k)}+\sum_{j=i+1}^na_{ij}x_j^{(k)}-b_i) \]
\[\forall i=1,2,\cdots,n \]

将格式写成矩阵形式

\[\begin{pmatrix} x_1^{(k+1)} \\ x_2^{(k+1)} \\ \vdots \\ x_n^{(k+1)} \end{pmatrix} =\begin{pmatrix} 0 & -\frac{a_{12}}{a_{11}} & -\frac{a_{13}}{a_{11}} & \cdots & -\frac{a_{1n}}{a_{11}} \\ -\frac{a_{21}}{a_{22}} & 0 & -\frac{a_{23}}{a_{22}} &\cdots & -\frac{a_{2n}}{a_{22}} \\ -\frac{a_{31}}{a_{33}} & -\frac{a_{32}}{a_{33}} & 0 & \cdots & -\frac{a_{3n}}{a_{33}} \\ \vdots & \vdots & \vdots &\ddots &\vdots \\ -\frac{a_{n1}}{a_{nn}} & -\frac{a_{n2}}{a_{nn}} & -\frac{a_{n3}}{a_{nn}} & \cdots & 0 \\ \end{pmatrix} \begin{pmatrix} x_1^{(k)} \\ x_2^{(k)} \\ \vdots \\ x_n^{(k)} \end{pmatrix} +g \]

迭代矩阵可以写成

\[-\begin{pmatrix} \frac{1}{a_{11}} \\ & \frac{1}{a_{22}} \\ & & \frac{1}{a_{33}} \\ & & & \ddots \\ & & & & \frac{1}{a_{nn}} \end{pmatrix} \begin{pmatrix} 0 & {a_{12}} & {a_{13}} & \cdots & {a_{1n}} \\ {a_{21}} & 0 & {a_{23}} &\cdots & {a_{2n}} \\ {a_{31}} & {a_{32}} & 0 & \cdots & {a_{3n}} \\ \vdots & \vdots & \vdots &\ddots &\vdots \\ {a_{n1}} & {a_{n2}} & {a_{n3}} & \cdots & 0 \\ \end{pmatrix} \]

若记

\[D=\begin{pmatrix} {a_{11}} \\ & {a_{22}} \\ & & {a_{33}} \\ & & & \ddots \\ & & & & {a_{nn}} \end{pmatrix} \]

则Jacoib迭代的迭代矩阵可以写为

\[G_J=-D^{-1}(A-D)=I-D^{-1}A \]

Jacobi迭代可以写成

\[x^{(k+1)}=-D^{-1}(A-D)x^{(k)}+D^{-1}b \]

迭代收敛的充要条件是$\rho(G_J)<1$。这个条件对于对一个任意形态的$n$阶矩阵使用。 可以给出一些方便使用的充分条件

定理 3.
矩阵$A$满足如下条件之一，则Jacobi迭代收敛。

  1. $A$是严格行或列对角占优阵
  2. $A$满足
    \[\sum_{i\neq j}\frac{|a_{ij}|}{|a_{ii}|}<1 , \forall j=1,2,\cdots,n \]

证明. 2. 利用迭代矩阵的$\| G_J\|_1<1$，可以证明。

1. 用迭代矩阵的$\| G_J\|_\infty<1$，可以证明$A$严格行对角占优时，Jacoib迭代收敛。

当$A$严格列对角占优时，

\[\|I-AD^{-1}\|_1=\max_j\sum_{i\neq j}\frac{|a_{ij}|}{|a_{jj}|}<1 \]

由

\[\rho(I-D^{-1}A)=\rho(I-AD^{-1}) \]

即证

注. 若有$x\neq0$，满足$(I-D^{-1}A)x=\lambda x$，则有

\[\begin{aligned} A(I-D^{-1}A)x &=\lambda Ax \\ Ax-AD^{-1}Ax &=\lambda Ax \\ (I-AD^{-1})y &=\lambda y, y=Ax \end{aligned} \]

即$I-D^{-1}A$与$I-AD^{-1}$有相同的特征值

注. $A^T$为严格行对角占优阵，则$A^T$的Jacobi迭代收敛。 由

\[\rho(I-D^{-1}A)=\rho(I-D^{-1}A^T) \]

可证明Jacobi迭代收敛。

算法 

1、输入系数矩阵A和向量b，和误差控制eps
2、x1={0,0,…..,0} , x2={1,1,…..,1} //赋初值
3、while( ||x1-x2||>eps)  {
            x1=x2; // x1 是第k步，x2是第k+1步
            for(i=0;i<n;i++)  {
                  x2[i]=0;
                  for(j=0;j<i;j++)  {
                        x2[i] += A[i][j]*x1[j]
                  }
                  for(j=i+1;j<n;j++)  {
                        x2[i] += A[i][j]*x1[j]
                  }
                  x2[i]=-(x2[i]-b[i])/A[i][i]
            }
     }
4、输出解x2

例 1. 计算如下矩阵的Jacobi迭代的迭代矩阵，并判定对任意初值的收敛性

\[A=\begin{pmatrix} 1. & 2. & -2. \\ 1. & 1. & 1. \\ 2. & 2. & 1. \\ \end{pmatrix} \]
\[B=\begin{pmatrix} 2. & -1. & 1. \\ 1. & 1. & 1. \\ 1. & 1. & -2. \\ \end{pmatrix} \]

解. 矩阵$A$的Jacobi 迭代矩阵为

\[-\begin{pmatrix} 0. & 2.& -2. \\ 1. & 0. & 1. \\ 2. & 2. & 0. \\ \end{pmatrix} \]

特征多项式$\lambda^3$，特征值$0$

矩阵$B$的Jacibo迭代矩阵

\[G=-\begin{pmatrix} 0 & 0.5 & -0.5 \\ -1 & 0 & -1 \\ 0.5 & 0.5 & 0 \\ \end{pmatrix} \]

特征多项式为$|\lambda I-G|=\lambda^3+\frac{5}4\lambda$，可得

\[\lambda_1=0, \lambda_{2,3}=\pm\frac{\sqrt 5}{2}i \]

Gauss-Seidel迭代

在Jacibo迭代中，使用新计算的分量来代替前一步的分量，得到格式

\[\begin{cases} x_1^{(k+1)}=-\frac1{a_{11}}(\quad +a_{12}x_2^{(k)}+\cdots+a_{1n}x_n^{(k)}-b_1) \\ x_2^{(k+1)}=-\frac1{a_{22}}( a_{21}x_1^{(k+1)}+a_{23}x_3^{(k)}+\cdots+a_{2n}x_n^{(k)}-b_2) \\ \vdots \\ x_n^{(k+1)}=-\frac1{a_{nn}}(a_{n1}x_1^{(k+1)}+\cdots+a_{n,n-1}x_{n-1}^{(k+1)}-b_n) \\ \end{cases} \]

写成分量形式为

\[x_i^{(k+1)}=-\frac1{a_{ii}}(\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}+\sum_{j=i+1}^na_{ij}x_j^{(k)}-b_i) \]
\[\forall i=1,2,\cdots,n \]

算法 

1、输入系数矩阵A和向量b，和误差控制eps
2、x1={0,...,0}, x2={1,1,…..,1} //赋初值
3、while( ||x1-x2||>eps)  {
            x1=x2; //x1 是第k步，x2是第k+1步
            for(i=0;i<n;i++)  {
                  t=0.0
                  for(j=0;j<i;j++)  {
                        t += A[i][j]*x2[j]
                  }
                  for(j=i+1;j<n;j++)  {
                        t += A[i][j]*x2[j]
                  }
                  x2[i]=-(t-b[i])/A[i][i]
            }
     }
4、输出解x2

注. 由表达式可以看到，若迭代收敛，则收敛到的向量满足方程组。

记

\[L=\begin{pmatrix} 0 \\ a_{21} & 0 \\ \vdots & & \ddots \\ a_{n1} & a_{n2} & \cdots & 0 \end{pmatrix} U=\begin{pmatrix} 0 & u_{12} & \cdots & u_{1n} \\ & 0 &\cdots & u_{2n} \\ & & \ddots \\ & & & 0 \\ \end{pmatrix} \]

则Gauss-Seidel迭代的可以写成

\[x^{(k+1)}=-D^{-1}(Lx^{(k+1)}+Ux^{(k)}-b) \]

进一步有

\[(D+L)x^{(k+1)}=-Ux^{(k)}+b \]

Gauss-Seidel迭代的迭代矩阵为

\[G_{GS}=-(D+L)^{-1}U \]

注. Gauss-Seidel迭代收敛的充要条件是$\rho(G_{GS})<1$

定理 4.
若$A$满足如下条件之一，则Gauss-Seidel迭代收敛：

  1. $A$为严格行或列对角占优阵
  2. $A$为对称正定阵

证明. 若$A$为行或列对角占优阵，则它的Gauss-Seidel迭代的迭代矩阵的特征多项式为

\[\begin{aligned} P_{GS}(\lambda) & =|\lambda I-G_{GS}|=|\lambda I-(D+L)^{-1}U| \\ &=|\lambda(D+L)^{-1}(D+L)+(D+L)^{-1}U| \\ &=|(D+L)^{-1}|\cdot|\lambda(D+L)+U| \end{aligned} \]

可以看到，当$|\lambda|\geq 1$时，$\lambda(D+L)+U$仍然是一个对角占优阵。 这样，就有$P_{GS}(\lambda)\neq 0$。

也就是说，若$P_{GS}(\lambda)=0$，则有$|\lambda|<1$。

即$\rho(G_{GS})<1$

例 2. 计算如下矩阵的Gauss-Seidel迭代的迭代矩阵，并判定对任意初值的收敛性

\[A=\begin{pmatrix} 1. & 2. & -2. \\ 1. & 1. & 1. \\ 2. & 2. & 1. \\ \end{pmatrix} \]
\[B=\begin{pmatrix} 2. & -1. & 1. \\ 1. & 1. & 1. \\ 1. & 1. & -2. \\ \end{pmatrix} \]

解. 矩阵$A$的 Gauss-Seidel 迭代矩阵为

\[-\begin{pmatrix} 0.& -2.& 2. \\ 0.& 2.& -3. \\ 0.& 0.& 2. \\ \end{pmatrix} \]

特征多项式$\lambda(\lambda-2)^2$，特征值$0,2$

矩阵$B$的Jacibo迭代矩阵

\[G=-\begin{pmatrix} 0 & 0.5 & -0.5 \\ 0 & -0.5 & -0.5 \\ 0 & 0 & -0.5 \\ \end{pmatrix} \]

特征多项式为$|\lambda I-G|=\lambda(\lambda+0.5)^2$，特征值$0,-0.5$

例 3. 给出如下矩阵的Jacobi和Gauss-Seidel迭代的迭代矩阵

\[A=\begin{pmatrix} 4. & 0. & 3. \\ 1. & 3. & 0. \\ 0. & 1. & 3. \\ \end{pmatrix} \]

解. Jacobi迭代矩阵

\[-\begin{pmatrix} 0 & 0& \frac34 \\ \frac13 & 0 & 0 \\ 0 & \frac13 & 0 \\ \end{pmatrix} \]

特征多项式$\lambda^3+\frac1{12}$，谱半径为$\frac1{\sqrt[3]{12}}$

Gauss-Seidel迭代矩阵

\[\begin{pmatrix} 0 & 0 & -\frac34 \\ 0 & 0 & \frac14 \\ 0 & 0 & -\frac1{12} \\ \end{pmatrix} \]

谱半径为$\frac1{12}$

可以看到Gauss-Seidel迭代的谱半径比Jacobi迭代的谱半径小，因此收敛更快。

Gauss-Seidel迭代的迭代矩阵，可以直接由公式$-(D+L)^{-1}U$来计算。或者用如下方法：

Gauss-Seidel迭代的分量形式为

\[\begin{cases} x^{(k+1)}_1= -\frac14(3x^{(k)}_3-b_1) \\ x^{(k+1)}_2= -\frac13(x^{(k+1)}_1-b_2)\\ x^{(k+1)}_3= -\frac13(x^{(k+1)}_2-b_3)\\ \end{cases} \]

将上式中$x_1^{(k+1)}$的表达式代入其它的式子，

\[\begin{cases} x^{(k+1)}_1= -\frac14(3x^{(k)}_3-b_1) \\ x^{(k+1)}_2= -\frac13(-\frac14(3x^{(k)}_3-b_1)-b_2)\\ x^{(k+1)}_3= -\frac13(x^{(k+1)}_2-b_3)\\ \end{cases} \]

将上式中$x_1^{(k+1)}$, $x_2^{(k+1)}$的表达式代入其它的式子，

\[\begin{cases} x^{(k+1)}_1= -\frac14(3x^{(k)}_3-b_1) \\ x^{(k+1)}_2= \frac14(x^{(k)}_3-\tilde b_2)\\ x^{(k+1)}_3= -\frac13(\frac14(x^{(k)}_3-\tilde b_2)-b_3)\\ \end{cases} \]

最终，可以得到

\[\begin{cases} x^{(k+1)}_1= -\frac34 x^{(k)}_3 +\bar b_1 \\ x^{(k+1)}_2= \frac14 x^{(k)}_3 +\bar b_2\\ x^{(k+1)}_3= -\frac1{12} x^{(k)}_3+\bar b_3\\ \end{cases} \]

松弛迭代

记

\[\Delta x^{(k)}=x^{(k+1)}-x^{(k)} \]

为迭代格式的修正量，将迭代格式看成是前一步值加上一个修正量的过程。

在修正量$\Delta x^{(k)}$上乘上一个控制因子$\omega$，可以更好的控制收敛的速度。 这样，就可以得到新的迭代格式，这个格式更灵活

\[x^{(k+1)}=x^{(k)}+\omega\Delta x^{(k)}=x^{(k)}+\omega(\tilde x^{(k+1)}-x^{(k)}) \]

其中的$\tilde x^{(k+1)}$是某个格式得到的。

取$\tilde x^{(k+1)}$为Jacobi迭代格式，则有

\[\begin{aligned} x^{(k+1)}=& (1-\omega)x^{(k)}+\omega((I-D^{-1}A)x^{(k)}+D^{-1}b) \\ =& (I-\omega D^{-1}A)x^{(k)}+\omega D^{-1}b \end{aligned} \]

称为Damped Jacobi Method。

取$\tilde x^{(k+1)}$为Gauss-Seidel迭代格式，则有

\[\begin{aligned} x^{(k+1)}=&(1-\omega)x^{(k)}+\omega(-D^{-1}(L x^{(k+1)}+D x^{(k)}-b))\\ =&(D+\omega L)^{-1}((1-\omega)D-\omega U)x^{(k)}+\omega(D+\omega L)^{-1}b \end{aligned} \]

称为松弛迭代 (Successive over-relaxation method)。

  • $\omega=1$，即为Gauss-Seidel迭代。
  • 调节因子$\omega$，可以调节收敛的速度。

算法 

1、输入系数矩阵A,向量b,松弛因子omega,和误差控制eps
2、x1={0,...,0}, x2={1,1,…..,1} //赋初值
3、while( ||x1-x2||>eps)  {
            x1=x2; //x1 是第k步，x2是第k+1步
            for(i=0;i<n;i++)  {
                  t=0.0
                  for(j=0;j<i;j++)  {
                        t += A[i][j]*x2[j]
                  }
                  for(j=i+1;j<n;j++)  {
                        t += A[i][j]*x2[j]
                  }
                  t=-(t-b[i])/A[i][i]
                  x2[i]=(1-omega)*x2[i]+omega*t
            }
     }
4、输出解x2

定理 5.
若SOR方法收敛，则$0<\omega<2$

证明. 当SOR方法收敛，则$\rho(G_\omega)<1$，所以有

\[|det(G_\omega)|=|\lambda_1\lambda_2\cdots\lambda_n|<1 \]

另外，

\[\begin{aligned} det(G_\omega)&=|(D+\omega L)^{-1}((1-\omega)D-\omega U)| \\ &=|(I+\omega D^{-1}L)^{-1}D^{-1}((1-\omega)D-\omega U)| \\ &=|(I+\omega D^{-1}L)^{-1}|\cdot|(1-\omega)I-\omega D^{-1}U| \\ &=(1-\omega)^n \end{aligned} \]

即有

\[|1-\omega|<1 , \quad 0<\omega<2 \]

定理 6.
若$A$是对称正定阵，则SOR方法当$0<\omega<2$时收敛。

证明. 设$\lambda$是$G_\omega$的任一特征值，$v$是对应的特征向量，则

\[((1-\omega)D-\omega U)v=\lambda (D+\omega L)v \]

两边对$v$求内积

\[(1-\omega)(Dv,v)-\omega(Uv,v)=\lambda((Dv,v)+\omega(Lv,v)) \]

可得

\[\lambda=\frac{(1-\omega)(Dv,v)-\omega(Uv,v)}{(Dv,v)+\omega(Lv,v)} \]

由$A=D+L+U$是正定的，则$a_{ii}>0$，即$D$是正定的，且$L=U^T$。 记$(Lv,v)=\alpha+i\beta$，则有

\[(Dv,v)=\sigma>0 \]
\[(Uv,v)=(v,Lv)=\overline{(Lv,v)}=\alpha-i\beta \]
\[0<(Av,v)=(Dv,v)+(Lv,v)+(Uv,v)=\sigma+2\alpha \]

则

\[\lambda=\frac{(1-\omega)\sigma-\omega(\alpha-i\beta)}{\sigma+\omega(\alpha+i\beta)} \]
\[|\lambda|^2=\frac{((1-\omega)\sigma-\omega\alpha)^2+\omega^2\beta^2}{(\sigma+\omega\alpha)^2+\omega^2\beta^2} \]

当$0<\omega<2$时，

\[\begin{aligned} (\sigma(1-\omega)-\omega\alpha)^2-(\sigma+\omega\alpha)^2 &=(-2\omega\alpha-\sigma \omega)(2\sigma-\sigma\omega) \\ &=\omega(-2\alpha-\sigma)(2-\omega)\sigma<0 \end{aligned} \]

即$|\lambda|^2<1$。因此SOR方法收敛。




注. 类似，可以证明，若$A$对称正定，Jacobi迭代收敛$\Leftrightarrow$ $2D-A$正定

程序作业

编写Gauss-Seidel迭代与松弛迭代的程序，并且计算

\[\begin{pmatrix} 31 & -13 & 0 & 0 & 0 & -10 \\ -13 & 35 & -9 & 0 & -11 \\ & -9 & 31 & -10 \\ & & -1 & 79 & -30 & 0 & 0 & 0 & -9 \\ & & &-30 & 57 & -7 & 0 & -5 \\ & & & & -7 & 47 & -30 \\ & & & & & -30 & 41 \\ & & & & -5 & & & 27 & -2 \\ & & & -9 & & & &-2 & 29 \\ \end{pmatrix} x= \begin{pmatrix} -15 \\ 27 \\ -23 \\ 0 \\ -20 \\ -12 \\ -7 \\ 7 \\ 10 \end{pmatrix} \]

给出迭代步数。

  1. 给出Gauss-Seidel迭代需要的迭代步数
  2. 对$\omega=\frac{i}{50}, i=1,2,\cdots,99$，给出SOR迭代需要的迭代步数

输出： 

初值：
0.0
0.0
...

Gauss-Seidel迭代的迭代步数： 100
SOR的迭代步数：
0.02, 300
0.04, 320
...

谢谢

例 4. 对线性方程组

\[\begin{pmatrix} 1 & 2 & -2 \\ 1 & 1 & 1 \\ 2 & 2 & 1 \\ \end{pmatrix} x= \begin{pmatrix} 1 \\ 2 \\ 3 \end{pmatrix} \]
  1. 写出Jacobi迭代和Gauss-Seidel迭代的分量形式
  2. 判断（1）中两种方法的收敛性，如果均收敛，说明哪一种收敛更快
  3. 取初值$(1,1,1)$，用Gauss-Seidel迭代法计算2步

目录

本节读完

例 5.

[#ex9-1-0].