矩阵计算与分析-求解方程组-预处理共轭梯度法

(k ) −T ( k ) x = S u ,有 把这组公式换回原来的变量，令 ~ ( k ) = g − Fu( k ) = S −1 (b − AS −T S T x ( k ) ) = S −1r ( k ) r
再令 p( k ) = S −T ~ p ( k ) , p( 0 ) = S −T S −1 p( 0 ) , 再引入
( k +1 )
=z
( k +1 )
+ βk p
(k )
一般 z
( k +1 )
=M r
−1 ( k + 1 )
用解 Mz
( k +1 )
=r
( k +1 )
来实现 .
下面讨论 M 和 S 的选择．首先可以考虑 M 的 Cholesky分解： M = LLT，即 S = L，而且，当 LLT ≈ A时,有F ≈ L−1(LLT ) L−T = I, 有cond(F )≈1, 这就改善了条件数．一般可以考虑A的一种分裂： A = M − N，其中M对称正定，M = LLT，而N尽可 能小，这称为A的不完全Cholesky分解． M 的一种最简单的选择是取为对角阵，这里仍 设 A =D −L −U，严格下三角阵L =UT，设M = D, 即Jacobi迭代法的预处理矩阵，就有 S = S T = diag ( a11 , " , ann ) 这种方法虽然简单，但往往不能有效地降低条件数.
z ( k ) = S −T S −1 r ( k ) = M −1 r ( k ) ,
得到预处理共轭梯度法( PCG方法)为 (1) ∀x(0)∈R n，
( 2) r ( 0 ) = b − Ax ( 0 ) , z ( 0 ) = M −1r ( 0 ) , p( 0 ) = z ( 0 )
(3) 对k =0 ,1 , …
M = SS T (7.19) x = S −T u g= S b,
−1
M是对称正定阵，把 Ax = b改写为等价的方程组
S −1 AS −T u = S −1b , F = S AS
−1 −T
新方程组的系数矩阵保持了对称正定．令
,
新方程源自文库写成
Fu = g (7.20)
对(7.20)用CG法计算，即
可以证明，经过这样的预处理，F =S−1AS−T 的条 件数大约是A的条件数的平方根．特别是ω = 1情 形，即对称的GS迭代预处理会有好的效果．
( z(k ) , r (k ) ) αk = (k ) ( p , Ap( k ) ) x ( k +1 ) = x ( k ) + α k p ( k ) r ( k +1) = r ( k ) − α k Ap( k )
z ( k + 1 ) = M −1 r ( k + 1 )
( z ( k +1 ) , r ( k +1 ) ) βk = ( z(k ) , r (k ) ) p
另外一种方法是取M 为SSOR迭代法的预处理 矩阵，即 1 M SSOR = ( D − ω L ) D −1 ( D − ω U ) ω( 2 − ω) 它是一个对称正定阵，令M =SST，其中
S = [ω( 2 − ω)]−1 / 2 ( D − ωL) D −1 / 2 S T = [ω( 2 − ω)]−1 / 2 D −1 / 2 ( D − ωLT )
9.4.5 预处理共轭梯度法
由 x
(k )
⎛ K − 1 ⎞ (0) − x * A ≤ 2⎜ ⎟ x − x* A ⎝ K + 1⎠
k
可以看到，当A病态时有K>>1， CG法收敛将很 慢．为了改善收敛性，可以设法先降低矩阵的条 件数，这就是预处理方法．
当 A对称正定时，希望预处理后的方程组仍保 持对称正定，为此，设 S 可逆，
(1) ∀u(0)∈R n， ~ ( 0 ) = g − Fu( 0 ) , ~ ~ (0) , ( 2) 令 r p(0) = r (3)对k =0 ,1 , … , ~(k ) , r ~(k ) ) r ( ~ = α k (~ p ( k ) , S −1 AS −T ~ p(k ) )
(k ) ~ ~ u( k +1 ) = u( k ) + α p k (k ) ~ ( k +1 ) = r ~(k ) − α ~ S −1 AS −T ~ r p k ( k +1 ) ~ ( k +1 ) ~ r ) ( ,r ~ βk = ~ ( k ) ~ ( k ) (r , r ) ~ ~(k ) ( k +1 ) ( k +1 ) ~ ~ p =r + βk p