预处理共轭梯度法(1)
1.共轭梯度法的收敛性分析
定理:设A为n x n对称正定矩阵,其最大与最小特征值分别为λ1\lambda _{1}λ1和λn\lambda _{n}λn,Ax=bAx = bAx=b的精确解为x∗x^{*}x∗,则对任意初始x(0)x^{(0)}x(0),求解Ax=bAx = bAx=b的共轭梯度法有
∥x(k)−x∗∥A⩽2(cond(A)−1cond(A)+1)k∥x(0)−x∗∥A\left \|x^{(k)}-x^{*}\right \|_{A}\leqslant 2(\frac{\sqrt{cond(A)}-1}{\sqrt{cond(A)}+1})^{k}\left \|x^{(0)}-x^*\right \|_{A}∥∥∥x(k)−x∗∥∥∥A⩽2(cond(A)+1cond(A)−1)k∥∥∥x(0)−x∗∥∥∥A
其中cond(A)2=λ1(A)λn(A)cond(A)_2=\frac{\lambda_1(A)}{\lambda_n(A)}cond(A)2=λn(A)λ1(A)。
我们观察(cond(A)−1cond(A)+1)k(\frac{\sqrt{cond(A)}-1}{\sqrt{cond(A)}+1})^{k}(cond(A)+1cond(A)−1)k可以知道
当λ1>>λn\lambda_1>>\lambda_nλ1>>λn时,或者cond(A)2cond(A)_2cond(A)2较大时,共轭梯度法的收敛效率会变得比较低。
2.预处理的基本思想
预处理被称为PCG方法(preconditioned conjugated gradient method)
既然共轭梯度法的收敛速度取决于系数矩阵的特征值,那么我们可以将Ax=bAx = bAx=b转化为等价的A~x=b~\widetilde{A}x = \widetilde{b}Ax=b。
使得A~x=b~\widetilde{A}x = \widetilde{b}Ax=b在与Ax=bAx = bAx=b同解的前提下,而A~\widetilde{A}A的最大、最小的特征值之比远小于A的最大、最小的特征值之比。
从而再次运用共轭梯度法求解方程组能够达到提高收敛速度的效果。
3.预处理共轭梯度法
求解Ax=bAx = bAx=b
其中A为n阶系数的对称正定矩阵。现寻找非奇异的n阶矩阵C,使得A‾=C−1A(C−1)T\overline{A} = C^{-1}A (C^{-1})^TA=C−1A(C−1)T的条件数比A的条件数小。而A‾\overline{A}A也是对称正定矩阵。
我们接着令x‾=CTx,b=C−1b\overline{x} = C^{T}x,b = C^{-1}bx=CTx,b=C−1b,最终将问题转化为求解A‾x‾=b‾\overline{A}\overline{x} = \overline{b}Ax=b
然后通过x=(CT)−1x‾x = (C^{T})^{-1}\overline{x}x=(CT)−1x求出原问题的解。
4.预处理矩阵
已知在通过迭代法求解线性方程组的过程中,我们需要清楚每次下降的方向和下降的步长。
(1)下降的方向
d‾(k+1)=d‾(r+1)+β‾kd‾(k)\overline{d}^{(k+1)}=\overline{d}^{(r+1)}+\overline{\beta}_{k}\overline{d}^{(k)}d(k+1)=d(r+1)+βkd(k)
通过已知的变换,我们可以得出d‾(k)=CTd(k)\overline{d}^{(k)} = C^Td^{(k)}d(k)=CTd(k)(与xxx的变换是一致的)和r‾(k)=C−1r(k)\overline{r}^{(k)}=C^{-1}r^{(k)}r(k)=C−1r(k)。
我们有:d(k+1)=C−Td‾(k+1)=C−T(r‾(k+1)+β‾kd‾(k))d^{(k+1)}=C^{-T}\overline{d}^{(k+1)}=C^{-T}(\overline{r}^{(k+1)}+\overline{\beta}_{k}\overline{d}^{(k)})d(k+1)=C−Td(k+1)=C−T(r(k+1)+βkd(k))
其中,β‾k=(r‾(k+1),r‾(k+1))(r‾(k),r‾(k))\overline{\beta}_{k}=\frac{(\overline{r}^{(k+1)},\overline{r}^{(k+1)})}{(\overline{r}^{(k)},\overline{r}^{(k)})}βk=(r(k),r(k))(r(k+1),r(k+1))
带入求解后:d(k+1)=(CCT)−1r(k+1)+β‾kd(k)d^{(k+1)}=(CC^T)^{-1}r^{(k+1)}+\overline{\beta}_kd^{(k)}d(k+1)=(CCT)−1r(k+1)+βkd(k)
我们令M=CCTM = CC^{T}M=CCT,并称其为预处理矩阵
同时我们可以得到:β‾k=(r‾(k+1),r‾(k+1))(r‾(k),r‾(k))=(z(k+1),r(k+1))(z(k),r(k))\overline{\beta}_{k}=\frac{(\overline{r}^{(k+1)},\overline{r}^{(k+1)})}{(\overline{r}^{(k)},\overline{r}^{(k)})}=\frac{({z}^{(k+1)},{r}^{(k+1)})}{({z}^{(k)},{r}^{(k)})}βk=(r(k),r(k))(r(k+1),r(k+1))=(z(k),r(k))(z(k+1),r(k+1))
(2)下降的步长
α‾k=(r‾(k),r‾(k))(d(k),Ad(k))=(z(k),r(k+1))(Ad(k),d(k))\overline{\alpha}_{k}=\frac{(\overline{r}^{(k)},\overline{r}^{(k)})}{({d}^{(k)},{Ad}^{(k)})}=\frac{({z}^{(k)},{r}^{(k+1)})}{(Ad^{(k)},{d}^{(k)})}αk=(d(k),Ad(k))(r(k),r(k))=(Ad(k),d(k))(z(k),r(k+1))
万事俱备后,我么就可以开始将迭代法进行下去了。
关于预处理矩阵的选取会放在后面讲
5.PCG算法
我们先构造预处理矩阵M(对称正定)
算法:
于是我们接下来的问题就是如何选取预处理矩阵M了