线性方程组的迭代法

不定常迭代法

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

不定常迭代法

不定常迭代法

不定常迭代包括2类：

  • 求解对称正定阵的最速下降法(steepest descent method)和共轭梯度法(cg)
  • 求解不对称阵的广义极小残量法(Gerneral Minimal RESidual，简写成GMRES)

最速下降法和共轭梯度法本质上来说是一种变分方法，对应于求一个二次函数的极值，是一种极小化方法。

共轭梯度法在最初是当作直接法出现在20世纪50年代。到了八九十年代，预处理共轭梯度法得到很大的发展，成为求解大型稀疏对称正定矩阵的最有效方法。

最速下降法

设$n$阶方阵$A=(a_{ij})$对称正定，$x$, $b$为$n$阶向量，对于线性方程组$Ax=b$， 先给出等价的变分问题。定义$n$元二次函数$\phi:\mathbb{R}^n\to\mathbb{R}$,

\[\phi(x)=\frac12(x,Ax)-(x,b) \]

其中$(x,y)=\displaystyle\sum_{i=1}^nx_iy_i$。

设$x=(x_1,\cdots,x_n)^T$, $b=(b_1,\cdots,b_n)^T$，则

\[\phi(x)=\frac12\sum_{i=1}^n\sum_{j=1}^na_{ij}x_ix_j-\sum_{j=1}^nb_jx_j \]

二次函数$\phi(x)$具有如下特性：

（1） 可以得到

\[\frac{\partial \phi}{\partial x_i}=a_{i1}x_1+\cdots+a_{in}x_n-b_i , i=1,2,\cdots,n \]

也就是说，函数$\phi$的梯度为

\[\nabla \phi(x)=Ax-b \]

（2） 对于$x,y\in\mathbb{R}^n$, $\alpha\in\mathbb{R}$，有

\[\phi(x+\alpha y)=\phi(x)+\alpha(y,Ax-b)+\frac{\alpha^2}2(y,Ay) \]

（3） 若$x^*=A^{-1}b$为方程的解，则有

\[\begin{split} \phi(x^*)&=\frac12(x^*,b)-(x^*,b) \\ &=-\frac12(x^*,b)=-\frac12(x^*,Ax^*) \end{split} \]

且，对于$x\in\mathbb{R}^n$，有

\[\begin{split} \phi(x)-\phi(x^*)& =\frac12(x,Ax)-(x,b)+\frac12(x^*,b)\\ &=\frac12(x,Ax)-(x,Ax^*)+\frac12(x^*,Ax^*) \\ &=\frac12(x-x^*, A(x-x^*)) \end{split} \]

定理 1.
$A$对称正定，$x^*$是方程$Ax=b$的解的充要条件是$x^*$为二次函数$\phi(x)$的极小值点，即

\[\phi(x^*)=\min_{x\in\mathbb{R}^n}\phi(x) \]

证明. 必要性： 设$Ax^*=b$，则由性质（3）

\[\phi(x)-\phi(x^*)=\frac12((x-x^*),A(x-x^*)) \]

上式为$(x-x^*)$的二次型，由$A$的正定性知，$\phi(x)-\phi(x^*)\geq 0$，即

\[\phi(x^*)=\min_{x\in\mathbb{R}^n}\phi(x) \]

充要性： 设$\phi(x^*)=\displaystyle\min_{x\in\mathbb{R}^n}\phi(x)$，则有

\[\frac{d\phi(x^*+\alpha y)}{d\alpha}|_{\alpha=0}=0 \]

即

\[\lim_{\alpha\to0}\frac{\phi(x^*+\alpha y)-\phi(x^*)}{\alpha}=0 \]

由性质（2）知，

\[\frac{\phi(x^*+\alpha y)-\phi(x^*)}{\alpha}=(y,Ax^*-b)+\frac12\alpha(y,Ay) \]

所以，$(y,Ax^*-b)=0, \forall y\in\mathbb{R}$。即有$Ax^*-b=0$

如何找到$\phi(x)$的极小值点$x^*$呢？

  • 由点$x^{(k)}$出发，沿指定的方向$y^{(k)}$搜索点$x^{(k+1)}=x^{(k)}+\alpha y^{(k)}$，使$\phi(x^{(k+1)}$最小。
  • 选择$y^{(k)}$的方式不同，算法也不同

记$r^{(k)}=b-Ax^{(k)}$(称为$x^{(k)}$对应的残量)，由性质（2）

\[\phi(x^{(k)}+a y^{(k)})=\phi(x^{(k)})-a(y^{(k)}, r^{(k)})+\frac12 a^2(y^{(k)},Ay^{(k)}) \]

这是一个关于$a$的二次多项式，

由$A$的正定性知$(y^{(k)},Ay^{(k)})>0$，由二次多项式的性质知，当

\[a=a_k=\frac{(y^{(k)},r^{(k)})}{(y^{(k)}, Ay^{(k)})} \]

时，$\phi(x^{(k)}+a y^{(k)})$取到最小值。

且，有

\[\begin{split} \phi((x^{(k)}+a_k y^{(k)})-\phi(x^{(k)})=-\frac12 a_k (y^{(k)}, r^{(k)}) \\ =-\frac12\frac{(y^{(k)},r^{(k)})^2}{(y^{(k)}, Ay^{(k)})} \end{split} \]

当$(y^{(k)},r^{(k)})\neq 0$时，$\phi(x^{(k+1)})<\phi(x^{(k)})$

注意到对于任意非$0$实数$\lambda$，有

\[\frac12\frac{(\lambda y^{(k)},r^{(k)})^2}{(\lambda y^{(k)}, A\lambda y^{(k)})} =\frac12\frac{\lambda^2(y^{(k)},r^{(k)})^2}{\lambda^2(y^{(k)}, Ay^{(k)})} \]

  • $\phi(x^{(k)}$的下降量只与$y^{(k)}$的方向有关，而与长度无关。
  • 函数$\phi(x)$在$x^{(k)}$处下降最快的方向是该点的负梯度方向，由性质（1）知，取$y^{(k)}=r^{(k)}$

这样，得到整个最速下降法算法为

  1. 选定初值$x^0$
  2. 对$k=0,1,2,\cdots$，计算：
    1. $r^{k}=b-Ax^k$
    2. $\alpha_k=\frac{(r^k,r^k)}{(r^k, Ar^k)}$
    3. $x^{k+1}=x^k+\alpha_k r^k$

定理 2.
$A$对称正定，$\lambda_1$，$\lambda_n$分别是$A$的最大和最小的特征值。则由最速下降法得到的序列$x^k$满足误差估计

\[\|x^k-x^*\|_A\leq(\frac{\lambda_1-\lambda_n}{\lambda_1+\lambda_n})^k\|x^0-x^*\|_A \]

其中$x^*$是方程$Ax=b$的根，$\|x\|_A=\sqrt{x^TAx}$是$\mathbb{R}^n$中的向量范数。

注. 可以算出，当$\lambda_1>>\lambda_n$时，收敛会很慢

共轭梯度法

  • 最速下降法中，取负梯度方向从局部来说，是最佳的搜索方向。但从整体来看，并非最优。
  • 在每一步使用小代价得到新的搜索方向，使得整体比最速下降法要快

现在，搜索方向不再是$r^0, r^1, \cdots, r^k$，假定新的搜索方向是$p^0,p^1,\cdots,p^k$。

则由前知

\[\alpha=\alpha_k=\frac{(r^k, p^k)}{(Ap^k,p^k)} \]

时，$\phi(x^k+\alpha p^k)$取到最小。这样，有

\[\begin{split} x^{k+1}=x^k+\alpha_k p^k \\ r^{k+1}= b-Ax^{k+1}=r^k-\alpha_k Ap^k \end{split} \]

可以得到

\[\begin{array}{rl} x^{k+1}&=x^k+\alpha_k p^k \\ &=x^{k-1}+\alpha_{k-1} p^{k-1}+\alpha_k p^k \\ &=x^0+\alpha_0 p^0+\cdots+\alpha_k p^k \end{array} \]

不失一般性，令$x^0=0$，则

\[x^{k+1}=\alpha_0 p^0+\cdots+\alpha_k p^k \]

现在，希望$\phi(x^{k+1}$不仅对$\alpha_k$最小，同时对$p^{i}, i=0,1,\cdots,k$也最小。

记$y=\alpha_0 p^0+\cdots+\alpha_{k-1} p^{k-1}$，求$\phi(y+\alpha p^{k})$对$\alpha$和$y$的最小值。

可以得到

\[\phi(y+\alpha p^{k})=\phi(y)+\alpha \phi(Ay, p^{k})-\alpha(b, p^{k})+\frac{\alpha^2}2(Ap^k, p^k) \]

上式的交叉项$\alpha(Ay,p^k)$不好处理，为此，令

\[(Ay,p^k)=0, \forall y=\alpha_0 p^0+\cdots+\alpha_{k-1} p^{k-1} \]

即，对$k=1,2,\cdots,n$，成立

\[(A p^i, p^k)=0, \forall i=0,\cdots,k-1 \]

定义 1.
$A$对称正定，若$\mathbb{R}^n$中向量组$\{p^0,p^1,\cdots,p^l\}$满足

\[(Ap^i,p^j)=0, \forall i\neq j \]

则称向量组为$\mathbb{R}^n$中的一个$A$-共轭向量组，或$A$-正交向量组，或$A$共轭的

当$\{p^0,\cdots, p^k\}$是$A$共轭时，有

\[\min_{y,\alpha}\phi(y+\alpha p^k)=\min_y\phi(y)+\min_{\alpha}(-\alpha(b, p^{k})+\frac{\alpha^2}2(Ap^k, p^k)) \]

设$p^k=r^k+\beta_{k-1}p^{k-1}$，利用A正交性，可以得到

\[\beta_{k-1}=-\frac{(r^k,Ap^{k-1})}{(p^{k-1}, Ap^{k-1})} \]

共轭梯度法算法

  1. 选定初值$x^0$，计算$p^0=r^0=b-Ax^0$
  2. 对$k=0,1,2,\cdots$，计算：
    1. $\alpha_k=\frac{\|r^k\|_2^2}{(Ap^k, p^k)}$
    2. $x^{k+1}=x^k+\alpha_kp^k$
    3. $r^{k+1}=r^k-\alpha_kAp^k$
    4. 若$\|r^{k+1}\|$足够小，得到解$x^{k+1}$
    5. 否则
      1. $\beta_k=\frac{\|r^{k+1}\|_2^2}{\|r^k\|_2^2}$
      2. $p^{k+1}=r^{k+1}+\beta_k p^k$

定理 3.
由共轭梯度法得到的向量组$\{r^i\}$和$\{p^i\}$满足

  1. $(p^i, r^j)=0$, $\forall 0\leq i<j\leq k$
  2. $(r^i, r^j)=0$, $\forall 0\leq i,j\leq k$, $i\neq j$
  3. $(p^i, Ap^j)=0$, $\forall 0\leq i,j\leq k$, $i\neq j$

注. 由此，可以证明，共轭梯度法最多经过$n$步迭代后，可以得到精确解。

定理 4.
$A$对称正定，$\lambda_1$，$\lambda_n$分别是$A$的最大和最小的特征值。则由共轭梯度法得到的序列$x^k$满足误差估计

\[\|x^k-x^*\|_A\leq 2(\frac{\sqrt{\lambda_1}-\sqrt{\lambda_n}}{\sqrt{\lambda_1}+\sqrt{\lambda_n}})^k\|x^0-x^*\|_A \]

其中$x^*$是方程$Ax=b$的根，$\|x\|_A=\sqrt{x^TAx}$是$\mathbb{R}^n$中的向量范数。

注. 由

\[\begin{split} \frac{\lambda_1-\lambda_n}{\lambda_1+\lambda_n} =\frac{\sqrt{\lambda_1}-\sqrt{\lambda_n}}{\sqrt{\lambda_1}+\sqrt{\lambda_n}} \frac{(\sqrt{\lambda_1}+\sqrt{\lambda_n})^2}{\lambda_1+\lambda_n} \\ =2\frac{\sqrt{\lambda_1}-\sqrt{\lambda_n}}{\sqrt{\lambda_1}+\sqrt{\lambda_n}} \frac12(1+\frac{2\sqrt{\lambda_1\lambda_n}}{\lambda_1+\lambda_n}) \\ \end{split} \]

知，共轭梯度法比最速下降法的收敛性要好很多

广义极小残量法

求解非对称线性方程组的算法：广义极小残量法（General Minimal RESidual）。现在已经成为当前求解大型稀疏非对称线性方程组的主要手段。

预处理技术

由于浮点运算的误差，共轭梯度法和广义极小残量法在计算中得到的向量会逐渐失去正交性。况且在大型稀疏矩阵中，$n$步收敛仍然不能令人满意。 预处理技术用来加速收敛。

简单地说，以较小的代价找到矩阵$M$，然后解方程

\[M^{-1}Ax=M^{-1}b \]

或

\[AM^{-1}y=b, x=M^{-1}y \]

得到解$Ax=b$的算法。分别称为左预处理或右预处理的迭代方法。若存在$L$使得$M=LL^T$，则可以解

\[L^{-1}AL^{-T}y=L^{-1}b, x=L^{-T}y \]

称为对称预处理方法。对称预处理方法在共轭梯度法中经常使用。称$M$为预处理矩阵

预处理矩阵$M$需要满足：

  1. 构造$M$的代价很小
  2. $M$跟$A$足够接近
  3. 关于$M^{-1}$的线性方程组很容易求解

如，取$M=diag(A)$，就是一个常用的预处理矩阵。

  • 预处理后的矩阵$M^{-1}A$(或$AM^{-1}$,$L^{-1}AL^{-T}$)的特征值越集中，收敛效果越好。

目录

本节读完

例 1.

[#ex9-1-0]. 参考书： 现代数值计算（第2版），同济大学计算数学教研室，人民邮电出版社

P171 , 最速下降法误差估计