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