7.5 预处理方法 Nothing will be more central to computational science in the next century than the art of transforming a problem that appears intractable into another whose solution can be approximated rapidly. For Krylov subspace matrix iterations, this is preconditioning.
Trefethen & Bau III [125, page 319], 1997.
为了改善 krylov 迭代方法的收敛性与可靠性, 我们通常需要运用预处理技术 (preconditioning, 也称为预条件). 通俗地说, 预处理就是将原来难以求解的问题转化成一个等价的但比较容易求解的新问题. 预处理技术的研究是目前科学计算领域中的重要研究课题之一.
早在1948年，著名数学家Turing(也被称为计算机科学之父和人工智能之父)就提出了预处理这个概念，用于改善线性方程组的性态.
The experiences of Fox, Huskey, and Wilkinson prompted Turing to write a remarkable paper "Rounding-off errors in matrix processes" [127, 1948]. In this paper, Turing made several important contributions. He formulated the LU (actually, the LDU) factorization of a matrix, .... He used the word "preconditioning" to mean improving the condition of a system of linear equations (a term that did not come into popular use until the 1970s).
— Nicholas J. Higham [70, page 184], 2002.
随后, Evans 在 1968 年提出了类似的想法, 用于降低线性方程组的条件数.
To obtain improved convergence rates for the methods of successive displacement we require the coefficient matrix to have a P P P
D.J.Evans [42],1968.
关于预处理的主要参考文献有[13, 131].
7.5.1 预处理方法介绍 对于线性方程组而言, 预处理就是对系数矩阵进行适当的线性转换, 转换为一个新矩阵. 考虑线性方程组
A x = b , A ∈ R n × n 非 奇 异 , b ∈ R n . ( 7.34 ) A x = b, \qquad A \in \mathbb {R} ^ {n \times n} \text {非 奇 异}, \quad b \in \mathbb {R} ^ {n}. \qquad \qquad (7. 3 4) A x = b , A ∈ R n × n 非 奇 异 , b ∈ R n . ( 7.34 ) 如果在方程组的两边同时左乘一个非奇异矩阵 P ∈ R n × n P \in \mathbb{R}^{n \times n} P ∈ R n × n
P − 1 A x = P − 1 b . (7.35) P ^ {- 1} A x = P ^ {- 1} b. \tag {7.35} P − 1 A x = P − 1 b . ( 7.35 ) 这个新方程组就是预处理后的方程组, P P P
理论上讲, 任何一个非奇异矩阵都可以作为预处理子. 但一个好的预处理子 P P P
(1) 预处理后的线性方程组更容易求解;P P P
(1) P − 1 A P^{-1}A P − 1 A P z = r Pz = r P z = r
其中第一个要求是为了确保预处理后的线性方程组更容易求解, 也就是说, 选取的预处理方法是有效的, 通常要求 P − 1 A P^{-1}A P − 1 A P P P
预处理方式 (7.34) 称为左预处理. 相应地, 我们可以在方程组两边同时右乘 P − 1 P^{-1} P − 1
A P − 1 u = b , x = P − 1 u . (7.36) A P ^ {- 1} u = b, \quad x = P ^ {- 1} u. \tag {7.36} A P − 1 u = b , x = P − 1 u . ( 7.36 ) 另外, 我们也可以将 P P P P = L R P = LR P = L R
L − 1 A R − 1 u = b , x = R − 1 u . (7.37) L ^ {- 1} A R ^ {- 1} u = b, \quad x = R ^ {- 1} u. \tag {7.37} L − 1 A R − 1 u = b , x = R − 1 u . ( 7.37 ) 这就是两边预处理 以上是三种常用的预处理方式。这三种方式预处理后的系数矩阵分别为 P − 1 A P^{-1}A P − 1 A A P − 1 AP^{-1} A P − 1 L − 1 A R − 1 L^{-1}AR^{-1} L − 1 A R − 1 A A A
在实际使用中, 该选取哪种预处理方式, 需要根据问题本身和所用的方法来确定. 如对于对称正定线性方程组的 CG 方法, 三种方式都可以, 而对于 GMRES 方法, 则选取右预处理比较合适. 一方面是实际使用时, 得到的残量范数与原方程组的残量范数是一样的, 另一方面是, 右预处理极小化的是原始残量范数, 而左预处理极小化的是预处理后的残量.
这里需要指出的是, 在实际求解预处理后的方程组时, 我们并不会显式地计算 P − 1 P^{-1} P − 1 P − 1 P^{-1} P − 1 P − 1 A P^{-1}A P − 1 A
7.5.2 预处理CG方法 设 A ∈ R n × n A \in \mathbb{R}^{n \times n} A ∈ R n × n P P P P P P
P = L L ⊺ . P = L L ^ {\intercal}. P = L L ⊺ . 于是我们得到下面的预处理方程组
L − 1 A ( L T ) − 1 u = L − 1 b , x = ( L T ) − 1 u . (7.38) L ^ {- 1} A \left(L ^ {\mathsf {T}}\right) ^ {- 1} u = L ^ {- 1} b, \quad x = \left(L ^ {\mathsf {T}}\right) ^ {- 1} u. \tag {7.38} L − 1 A ( L T ) − 1 u = L − 1 b , x = ( L T ) − 1 u . ( 7.38 ) 用CG方法求解上述方程组,迭代 k k k u ( k ) u^{(k)} u ( k ) r ~ k ≜ L − 1 b − L − 1 A ( L T ) − 1 u ( k ) \tilde{r}_k\triangleq L^{-1}b - L^{-1}A(L^{\mathsf{T}})^{-1}u^{(k)} r ~ k ≜ L − 1 b − L − 1 A ( L T ) − 1 u ( k )
算法7.8.两边预处理CG方法 1: 给定初值 x ( 0 ) x^{(0)} x ( 0 ) r 0 = b − A x ( 0 ) r_0 = b - Ax^{(0)} r 0 = b − A x ( 0 ) r ~ 0 = L − 1 r 0 , p ~ 1 = r ~ 0 \tilde{r}_0 = L^{-1}r_0, \tilde{p}_1 = \tilde{r}_0 r ~ 0 = L − 1 r 0 , p ~ 1 = r ~ 0 k = 1 , 2 , … k = 1,2,\ldots k = 1 , 2 , …
5: ξ k = ( r ~ k − 1 , r ~ k − 1 ) ( L − 1 A ( L ⊤ ) − 1 p ~ k , p ~ k ) \xi_{k} = \frac{(\tilde{r}_{k - 1},\tilde{r}_{k - 1})}{(L^{-1}A(L^{\top})^{-1}\tilde{p}_{k},\tilde{p}_{k})} ξ k = ( L − 1 A ( L ⊤ ) − 1 p ~ k , p ~ k ) ( r ~ k − 1 , r ~ k − 1 )
6: u ( k ) = u ( k − 1 ) + ξ k p ~ k u^{(k)} = u^{(k - 1)} + \xi_k\tilde{p}_k u ( k ) = u ( k − 1 ) + ξ k p ~ k r ~ k = r ~ k − 1 − ξ k L − 1 A ( L ⊤ ) − 1 p ~ k \tilde{r}_k = \tilde{r}_{k - 1} - \xi_kL^{-1}A(L^\top)^{-1}\tilde{p}_k r ~ k = r ~ k − 1 − ξ k L − 1 A ( L ⊤ ) − 1 p ~ k μ k = ( r ~ k , r ~ k ) ( r ~ k − 1 , r ~ k − 1 ) \mu_{k} = \frac{(\tilde{r}_{k},\tilde{r}_{k})}{(\tilde{r}_{k - 1},\tilde{r}_{k - 1})} μ k = ( r ~ k − 1 , r ~ k − 1 ) ( r ~ k , r ~ k ) p ~ k + 1 = r ~ k + μ k p ~ k \tilde{p}_{k + 1} = \tilde{r}_k + \mu_k\tilde{p}_k p ~ k + 1 = r ~ k + μ k p ~ k
10: end for
由于算法7.8中得到的是预处理后的近似解 u ( k ) u^{(k)} u ( k ) r ~ k \tilde{r}_k r ~ k
我们对上述算法中的迭代公式进行适当的改写. 首先引入一个辅助变量 p k p_k p k
p k ≜ ( L T ) − 1 p ~ k , k = 1 , 2 , … . p _ {k} \triangleq (L ^ {\mathsf {T}}) ^ {- 1} \tilde {p} _ {k}, \quad k = 1, 2, \ldots . p k ≜ ( L T ) − 1 p ~ k , k = 1 , 2 , … . 于是有
p k + 1 = ( L T ) − 1 p k + 1 = ( L T ) − 1 r ~ k + μ k ( L T ) − 1 p ~ k = ( L T ) − 1 r ~ k + μ k p k . p _ {k + 1} = (L ^ {\mathsf {T}}) ^ {- 1} p _ {k + 1} = (L ^ {\mathsf {T}}) ^ {- 1} \tilde {r} _ {k} + \mu_ {k} (L ^ {\mathsf {T}}) ^ {- 1} \tilde {p} _ {k} = (L ^ {\mathsf {T}}) ^ {- 1} \tilde {r} _ {k} + \mu_ {k} p _ {k}. p k + 1 = ( L T ) − 1 p k + 1 = ( L T ) − 1 r ~ k + μ k ( L T ) − 1 p ~ k = ( L T ) − 1 r ~ k + μ k p k . 又由 r ~ k \tilde{r}_k r ~ k r k = L r ~ k r_k = L\tilde{r}_k r k = L r ~ k
r k = L r ~ k = L ( r ~ k − 1 − ξ k L − 1 A ( L ⊤ ) − 1 p ~ k ) = r k − 1 − ξ k A p k , r _ {k} = L \tilde {r} _ {k} = L \left(\tilde {r} _ {k - 1} - \xi_ {k} L ^ {- 1} A \left(L ^ {\top}\right) ^ {- 1} \tilde {p} _ {k}\right) = r _ {k - 1} - \xi_ {k} A p _ {k}, r k = L r ~ k = L ( r ~ k − 1 − ξ k L − 1 A ( L ⊤ ) − 1 p ~ k ) = r k − 1 − ξ k A p k , p k + 1 = ( L T ) − 1 r ~ k + μ k p k = ( L T ) − 1 L − 1 r k + μ k p k = P − 1 r k + μ k p k , p _ {k + 1} = \left(L ^ {\mathsf {T}}\right) ^ {- 1} \tilde {r} _ {k} + \mu_ {k} p _ {k} = \left(L ^ {\mathsf {T}}\right) ^ {- 1} L ^ {- 1} r _ {k} + \mu_ {k} p _ {k} = P ^ {- 1} r _ {k} + \mu_ {k} p _ {k}, p k + 1 = ( L T ) − 1 r ~ k + μ k p k = ( L T ) − 1 L − 1 r k + μ k p k = P − 1 r k + μ k p k , x ( k ) = ( L T ) − 1 u ( k ) = ( L T ) − 1 ( u ( k − 1 ) + ξ k p ~ k ) = x ( k − 1 ) + ξ k p k , x ^ {(k)} = \left(L ^ {\mathsf {T}}\right) ^ {- 1} u ^ {(k)} = \left(L ^ {\mathsf {T}}\right) ^ {- 1} \left(u ^ {(k - 1)} + \xi_ {k} \tilde {p} _ {k}\right) = x ^ {(k - 1)} + \xi_ {k} p _ {k}, x ( k ) = ( L T ) − 1 u ( k ) = ( L T ) − 1 ( u ( k − 1 ) + ξ k p ~ k ) = x ( k − 1 ) + ξ k p k , 其中
ξ k = ( L − 1 r k − 1 , L − 1 r k − 1 ) ( L − 1 A ( L T ) − 1 p ~ k , p ~ k ) = ( r k − 1 , ( L T ) − 1 L − 1 r k − 1 ) ( A ( L T ) − 1 p ~ k , ( L T ) − 1 p ~ k ) = ( r k − 1 , P − 1 r k − 1 ) ( A p k , p k ) , \xi_ {k} = \frac {(L ^ {- 1} r _ {k - 1} , L ^ {- 1} r _ {k - 1})}{(L ^ {- 1} A (L ^ {\mathsf {T}}) ^ {- 1} \tilde {p} _ {k} , \tilde {p} _ {k})} = \frac {(r _ {k - 1} , (L ^ {\mathsf {T}}) ^ {- 1} L ^ {- 1} r _ {k - 1})}{(A (L ^ {\mathsf {T}}) ^ {- 1} \tilde {p} _ {k} , (L ^ {\mathsf {T}}) ^ {- 1} \tilde {p} _ {k})} = \frac {(r _ {k - 1} , P ^ {- 1} r _ {k - 1})}{(A p _ {k} , p _ {k})}, ξ k = ( L − 1 A ( L T ) − 1 p ~ k , p ~ k ) ( L − 1 r k − 1 , L − 1 r k − 1 ) = ( A ( L T ) − 1 p ~ k , ( L T ) − 1 p ~ k ) ( r k − 1 , ( L T ) − 1 L − 1 r k − 1 ) = ( A p k , p k ) ( r k − 1 , P − 1 r k − 1 ) , μ k = ( L − 1 r k , L − 1 r k ) ( L − 1 r k − 1 , L − 1 r k − 1 ) = ( r k , P − 1 r k ) ( r k − 1 , P − 1 r k − 1 ) . \mu_ {k} = \frac {(L ^ {- 1} r _ {k} , L ^ {- 1} r _ {k})}{(L ^ {- 1} r _ {k - 1} , L ^ {- 1} r _ {k - 1})} = \frac {(r _ {k} , P ^ {- 1} r _ {k})}{(r _ {k - 1} , P ^ {- 1} r _ {k - 1})}. μ k = ( L − 1 r k − 1 , L − 1 r k − 1 ) ( L − 1 r k , L − 1 r k ) = ( r k − 1 , P − 1 r k − 1 ) ( r k , P − 1 r k ) . 记 z k ≜ P − 1 r k z_{k}\triangleq P^{-1}r_{k} z k ≜ P − 1 r k
算法7.9.PCG：预处理CG方法 1: 给定初值 x ( 0 ) x^{(0)} x ( 0 ) ε > 0 \varepsilon > 0 ε > 0 r 0 = b − A x ( 0 ) r_0 = b - Ax^{(0)} r 0 = b − A x ( 0 ) β = ∥ r 0 ∥ 2 \beta = \| r_0\| _2 β = ∥ r 0 ∥ 2
3: 令 z 0 = P − 1 r 0 , p 1 = z 0 z_0 = P^{-1}r_0, p_1 = z_0 z 0 = P − 1 r 0 , p 1 = z 0 ρ = r 0 ⊤ z 0 \rho = r_0^\top z_0 ρ = r 0 ⊤ z 0 k = 1 k = 1 k = 1 q k = A p k q_{k} = Ap_{k} q k = A p k ξ k = ρ / ( p k T q k ) \xi_{k} = \rho /(p_{k}^{\mathsf{T}}q_{k}) ξ k = ρ / ( p k T q k ) x ( k ) = x ( k − 1 ) + ξ k p k x^{(k)} = x^{(k - 1)} + \xi_{k}p_{k} x ( k ) = x ( k − 1 ) + ξ k p k r k = r k − 1 − ξ k q k r_k = r_{k-1} - \xi_k q_k r k = r k − 1 − ξ k q k < ε < \varepsilon < ε ρ 0 = ρ \rho_0 = \rho ρ 0 = ρ z k = P − 1 r k z_{k} = P^{-1}r_{k} z k = P − 1 r k ρ = r k T z k \rho = r_k^{\mathsf{T}}z_k ρ = r k T z k μ k = ρ / ρ 0 \mu_{k} = \rho /\rho_{0} μ k = ρ / ρ 0 p k + 1 = z k + μ k p k p_{k + 1} = z_k + \mu_k p_k p k + 1 = z k + μ k p k < ε < \varepsilon < ε x ( k ) x^{(k)} x ( k )
我们注意到, 矩阵 L L L P P P z k = P − 1 r k z_{k} = P^{-1}r_{k} z k = P − 1 r k P z k = r k Pz_{k} = r_{k} P z k = r k P − 1 P^{-1} P − 1
事实上, PCG 方法 7.9 也可以从左预处理方式导出. 考虑左预处理方程组
P − 1 A x = P − 1 b . (7.39) P ^ {- 1} A x = P ^ {- 1} b. \tag {7.39} P − 1 A x = P − 1 b . ( 7.39 ) 易知 P − 1 A P^{-1}A P − 1 A P P P
( x , y ) P ≜ ( P x , y ) = ( x , P y ) . (7.40) (x, y) _ {P} \triangleq (P x, y) = (x, P y). \tag {7.40} ( x , y ) P ≜ ( P x , y ) = ( x , P y ) . ( 7.40 ) 由于 P P P
( P − 1 A x , y ) P = ( A x , y ) = ( x , A y ) = ( x , P ( P − 1 A y ) ) = ( x , P − 1 A y ) P , (P ^ {- 1} A x, y) _ {P} = (A x, y) = (x, A y) = (x, P (P ^ {- 1} A y)) = (x, P ^ {- 1} A y) _ {P}, ( P − 1 A x , y ) P = ( A x , y ) = ( x , A y ) = ( x , P ( P − 1 A y )) = ( x , P − 1 A y ) P , 即 P − 1 A P^{-1}A P − 1 A P P P P P P P − 1 A P^{-1}A P − 1 A P P P
引入辅助变量 z k ≜ P − 1 r k z_{k}\triangleq P^{-1}r_{k} z k ≜ P − 1 r k
算法7.10.基于 P P P 1: 给定初值 x ( 0 ) x^{(0)} x ( 0 ) ε > 0 \varepsilon > 0 ε > 0 r 0 = b − A x ( 0 ) r_0 = b - Ax^{(0)} r 0 = b − A x ( 0 ) β = ∥ r 0 ∥ 2 \beta = \| r_0\| _2 β = ∥ r 0 ∥ 2 z 0 = P − 1 r 0 , p 1 = z 0 z_0 = P^{-1}r_0, p_1 = z_0 z 0 = P − 1 r 0 , p 1 = z 0 k = 1 k = 1 k = 1
5 : ξ k = ( z k − 1 , z k − 1 ) P ( P − 1 A p k , p k ) P = ( r k − 1 , z k − 1 ) ( A p k , p k ) 6 : x ( k ) = x ( k − 1 ) + ξ k p k 7 : r k = r k − 1 − ξ k A p k 8 : relres = ∥ r k ∥ 2 / β 9 : i f relres < ε t h e n \begin{array}{l} 5: \quad \xi_ {k} = \frac {\left(z _ {k - 1} , z _ {k - 1}\right) _ {P}}{\left(P ^ {- 1} A p _ {k} , p _ {k}\right) _ {P}} = \frac {\left(r _ {k - 1} , z _ {k - 1}\right)}{\left(A p _ {k} , p _ {k}\right)} \\ 6: \quad x ^ {(k)} = x ^ {(k - 1)} + \xi_ {k} p _ {k} \\ 7: \quad r _ {k} = r _ {k - 1} - \xi_ {k} A p _ {k} \\ 8: \quad \operatorname {r e l r e s} = \| r _ {k} \| _ {2} / \beta \\ 9: \quad \text {i f} \operatorname {r e l r e s} < \varepsilon \text {t h e n} \\ \end{array} 5 : ξ k = ( P − 1 A p k , p k ) P ( z k − 1 , z k − 1 ) P = ( A p k , p k ) ( r k − 1 , z k − 1 ) 6 : x ( k ) = x ( k − 1 ) + ξ k p k 7 : r k = r k − 1 − ξ k A p k 8 : relres = ∥ r k ∥ 2 / β 9 : i f relres < ε t h e n 10: break
11: end if
z k = z k − 1 − α k P − 1 A p k = P − 1 r k μ k = ( z k , z k ) P ( z k − 1 , z k − 1 ) P = ( r k , z k ) ( r k − 1 , z k − 1 ) 14 : p k + 1 = z k + μ k p k 15 : e n d f o r \begin{array}{l} z _ {k} = z _ {k - 1} - \alpha_ {k} P ^ {- 1} A p _ {k} = P ^ {- 1} r _ {k} \\ \mu_ {k} = \frac {\left(z _ {k} , z _ {k}\right) _ {P}}{\left(z _ {k - 1} , z _ {k - 1}\right) _ {P}} = \frac {\left(r _ {k} , z _ {k}\right)}{\left(r _ {k - 1} , z _ {k - 1}\right)} \\ 1 4: \quad p _ {k + 1} = z _ {k} + \mu_ {k} p _ {k} \\ 1 5: \quad e n d \quad f o r \\ \end{array} z k = z k − 1 − α k P − 1 A p k = P − 1 r k μ k = ( z k − 1 , z k − 1 ) P ( z k , z k ) P = ( r k − 1 , z k − 1 ) ( r k , z k ) 14 : p k + 1 = z k + μ k p k 15 : e n d f or 不难看出, 算法 7.10 和算法 7.9 是完全一样的.
7.5.3 预处理GMRES方法 对于非对称 Krylov 子空间方法, 也有三种预处理方式. 但与对称正定情形不同的是, 这三种方式并不等价, 而且有时效果会相差很大.
如果采用左预处理方式，则原方程组转化为
P − 1 A x = P − 1 b . (7.41) P ^ {- 1} A x = P ^ {- 1} b. \tag {7.41} P − 1 A x = P − 1 b . ( 7.41 ) 记 r k r_k r k r k = b − A x ( k ) r_k = b - Ax^{(k)} r k = b − A x ( k ) r ~ k = P − 1 r k \tilde{r}_k = P^{-1}r_k r ~ k = P − 1 r k r ~ k \tilde{r}_k r ~ k r k r_k r k
设 P P P
A P − 1 u = b , x = P − 1 u . (7.42) A P ^ {- 1} u = b, \quad x = P ^ {- 1} u. \tag {7.42} A P − 1 u = b , x = P − 1 u . ( 7.42 ) 给定迭代初值 x ( 0 ) x^{(0)} x ( 0 ) u ( 0 ) = P x ( 0 ) u^{(0)} = Px^{(0)} u ( 0 ) = P x ( 0 )
r ~ 0 = b − A P − 1 u ( 0 ) = b − A x ( 0 ) = r ( 0 ) , \tilde {r} _ {0} = b - A P ^ {- 1} u ^ {(0)} = b - A x ^ {(0)} = r ^ {(0)}, r ~ 0 = b − A P − 1 u ( 0 ) = b − A x ( 0 ) = r ( 0 ) , 即预处理残量与原方程组残量是一样的. 因此无需计算 u ( 0 ) u^{(0)} u ( 0 ) u ( m ) = u ( 0 ) + V m y m u^{(m)} = u^{(0)} + V_{m}y_{m} u ( m ) = u ( 0 ) + V m y m m m m
后的近似解, 则对应的原方程组近似解 x ( m ) x^{(m)} x ( m )
x ( m ) = P − 1 u ( m ) = P − 1 ( u ( 0 ) + V m y m ) = x ( 0 ) + P − 1 V m y m , x ^ {(m)} = P ^ {- 1} u ^ {(m)} = P ^ {- 1} (u ^ {(0)} + V _ {m} y _ {m}) = x ^ {(0)} + P ^ {- 1} V _ {m} y _ {m}, x ( m ) = P − 1 u ( m ) = P − 1 ( u ( 0 ) + V m y m ) = x ( 0 ) + P − 1 V m y m , 原方程组的残量为
r m = b − A x ( m ) = r 0 − A P − 1 V m y m = β v 1 − V m + 1 H m + 1 , m y m = V m + 1 ( β e 1 − H m + 1 , m y m ) . r _ {m} = b - A x ^ {(m)} = r _ {0} - A P ^ {- 1} V _ {m} y _ {m} = \beta v _ {1} - V _ {m + 1} H _ {m + 1, m} y _ {m} = V _ {m + 1} (\beta e _ {1} - H _ {m + 1, m} y _ {m}). r m = b − A x ( m ) = r 0 − A P − 1 V m y m = β v 1 − V m + 1 H m + 1 , m y m = V m + 1 ( β e 1 − H m + 1 , m y m ) . 由于 y m y_{m} y m
min y ∈ R m ∥ β e 1 − H m + 1 , m y ∥ 2 \min _ {y \in \mathbb {R} ^ {m}} \| \beta e _ {1} - H _ {m + 1, m} y \| _ {2} y ∈ R m min ∥ β e 1 − H m + 1 , m y ∥ 2 的解, 因此
∥ r m ∥ 2 = ∥ β e 1 − H m + 1 , m y m ∥ 2 = β ⋅ ∣ q 1 ( m + 1 ) ∣ . \left\| r _ {m} \right\| _ {2} = \left\| \beta e _ {1} - H _ {m + 1, m} y _ {m} \right\| _ {2} = \beta \cdot | q _ {1} (m + 1) |. ∥ r m ∥ 2 = ∥ β e 1 − H m + 1 , m y m ∥ 2 = β ⋅ ∣ q 1 ( m + 1 ) ∣. 这与不带预处理的 GMRES 方法是一样的. 因此在实际求解过程中我们可以直接得到原方程组残量的模, 而无需计算 u ( m ) u^{(m)} u ( m ) x ( m ) x^{(m)} x ( m )
右预处理GMRES方法具体描述如下：
算法7.11.右预处理GMRES方法 1: 给定初值 x ( 0 ) x^{(0)} x ( 0 ) ε > 0 \varepsilon > 0 ε > 0 r 0 = b − A x ( 0 ) r_0 = b - Ax^{(0)} r 0 = b − A x ( 0 ) β = ∥ r 0 ∥ 2 \beta = \| r_0\| _2 β = ∥ r 0 ∥ 2 v 1 = r ~ 0 / β , ξ = β e 1 v_{1} = \tilde{r}_{0} / \beta, \xi = \beta e_{1} v 1 = r ~ 0 / β , ξ = β e 1 j = 1 , 2 , … j = 1,2,\ldots j = 1 , 2 , … w ~ j = P − 1 v j \tilde{w}_j = P^{-1}v_j w ~ j = P − 1 v j w j = A w ~ j w_{j} = A\tilde{w}_{j} w j = A w ~ j i = 1 , 2 , … , j i = 1,2,\ldots ,j i = 1 , 2 , … , j h i j = ( w j , v i ) h_{ij} = (w_j,v_i) h ij = ( w j , v i ) w j = w j − h i j v i w_{j} = w_{j} - h_{ij}v_{i} w j = w j − h ij v i
10: end for 11: h j + 1 , j = ∥ w j ∥ 2 h_{j + 1,j} = \| w_j\| _2 h j + 1 , j = ∥ w j ∥ 2 i = 1 , 2 , … , j − 1 i = 1,2,\ldots ,j - 1 i = 1 , 2 , … , j − 1 % \% % G j − 1 , … , G 1 G_{j - 1},\dots,G_1 G j − 1 , … , G 1 H j + 1 , j H_{j + 1,j} H j + 1 , j
[ h i , j h i + 1 , j ] = [ c i s i − s i c i ] [ h i , j h i + 1 , j ] (13:) \left[ \begin{array}{c} h _ {i, j} \\ h _ {i + 1, j} \end{array} \right] = \left[ \begin{array}{c c} c _ {i} & s _ {i} \\ - s _ {i} & c _ {i} \end{array} \right] \left[ \begin{array}{c} h _ {i, j} \\ h _ {i + 1, j} \end{array} \right] \tag {13:} [ h i , j h i + 1 , j ] = [ c i − s i s i c i ] [ h i , j h i + 1 , j ] ( 13: ) 14: end for 15: if h j + 1 , j = 0 h_{j + 1,j} = 0 h j + 1 , j = 0
16: set m = j m = j m = j
17: end if 18: v j + 1 = w j / h j + 1 , j v_{j + 1} = w_j / h_{j + 1,j} v j + 1 = w j / h j + 1 , j ∣ h j , j ∣ > ∣ h j + 1 , j ∣ |h_{j,j}| > |h_{j+1,j}| ∣ h j , j ∣ > ∣ h j + 1 , j ∣ G j G_j G j c j = 1 1 + τ 2 , s j = c j τ c_{j} = \frac{1}{\sqrt{1 + \tau^{2}}}, s_{j} = c_{j}\tau c j = 1 + τ 2 1 , s j = c j τ τ = h j + 1 , j h j , j \tau = \frac{h_{j+1,j}}{h_{j,j}} τ = h j , j h j + 1 , j
22: s j = 1 1 + τ 2 , c j = s j τ s_j = \frac{1}{\sqrt{1 + \tau^2}}, c_j = s_j \tau s j = 1 + τ 2 1 , c j = s j τ τ = h j , j h j + 1 , j \tau = \frac{h_{j,j}}{h_{j+1,j}} τ = h j + 1 , j h j , j h j , j = c j h j , j + s j h j + 1 , j , h j + 1 , j = 0 h_{j,j} = c_j h_{j,j} + s_j h_{j+1,j}, h_{j+1,j} = 0 h j , j = c j h j , j + s j h j + 1 , j , h j + 1 , j = 0 % Apply G j to last column of H j + 1 , j \% \text{Apply } G_j \text{ to last column of } H_{j+1,j} % Apply G j to last column of H j + 1 , j [ ξ j ξ j + 1 ] = [ c j s j − s j c j ] [ ξ j 0 ] \begin{bmatrix} \xi_j \\ \xi_{j+1} \end{bmatrix} = \begin{bmatrix} c_j & s_j \\ -s_j & c_j \end{bmatrix} \begin{bmatrix} \xi_j \\ 0 \end{bmatrix} [ ξ j ξ j + 1 ] = [ c j − s j s j c j ] [ ξ j 0 ] G j G_j G j ∣ ξ j + 1 ∣ / β < ε |\xi_{j + 1}| / \beta < \varepsilon ∣ ξ j + 1 ∣/ β < ε % \% % ∥ r j ∥ 2 = ∣ ξ j + 1 ∣ \| r_j\| _2 = |\xi_{j + 1}| ∥ r j ∥ 2 = ∣ ξ j + 1 ∣ m = j m = j m = j y ( m ) = R m − 1 ξ ( m ) y^{(m)} = R_m^{-1}\xi^{(m)} y ( m ) = R m − 1 ξ ( m ) R m = H ( 1 : m , 1 : m ) R_m = H(1:m,1:m) R m = H ( 1 : m , 1 : m ) ξ ( m ) = ξ ( 1 : m ) \xi^{(m)} = \xi(1:m) ξ ( m ) = ξ ( 1 : m ) x ( m ) = x ( 0 ) + P − 1 V m y ( m ) x^{(m)} = x^{(0)} + P^{-1}V_{m}y^{(m)} x ( m ) = x ( 0 ) + P − 1 V m y ( m )
关于左、右预处理GMRES,我们有下面的最优性质
定理7.16 设 x L ( k ) x_{L}^{(k)} x L ( k ) x R ( k ) x_{R}^{(k)} x R ( k ) k k k x ( 0 ) x^{(0)} x ( 0 ) x L ( k ) x_{L}^{(k)} x L ( k ) ∥ P − 1 ( b − A x ) ∥ 2 \| P^{-1}(b - Ax)\|_{2} ∥ P − 1 ( b − A x ) ∥ 2 x ( 0 ) + K k ( P − 1 A , P − 1 r 0 ) x^{(0)} + \mathcal{K}_k(P^{-1}A,P^{-1}r_0) x ( 0 ) + K k ( P − 1 A , P − 1 r 0 ) x R ( k ) x_{R}^{(k)} x R ( k ) ∥ b − A x ∥ 2 \| b - Ax\|_{2} ∥ b − A x ∥ 2
7.5.4 预处理子构造 This process of "preconditioning" is essential to most successful applications of iterative methods.
Trefethen & Bau III [125, page 313], 1997.
Finding a good preconditioner to solve a given sparse linear system is often viewed as a combination of art and science. Theoretical results are rare and some methods work surprisingly well, often despite expectations.
— Saad [113], 2003.
预处理能否取得成功的关键就是能否找到一个好的预处理子。一个公认的事实是，好的预处理子通常是与问题本身的特征密切相关的。
A既要有很好的通用性, 又要有很好的加速效果的预处理子往往是可遇而不可求的.
关于预处理技术的理论分析很少, 大多数情况下只能根据经验来构造. 尽管如此, 在实际应用中, 这些根据经验构造出来的预处理子往往能取得很好的数值效果, 有时甚至会大大出乎人们的意料.
预处理方法主要是用来提升Krylov子空间迭代法的收敛速度，因此，对Krylov子空间迭代法的收敛性质的深刻理解对构造好的预处理方法有很大的帮助.
A preconditioner M M M M − 1 A M^{-1}A M − 1 A
Trefethen & Bau III [125, page 314], 1997.
一般来说，预处理子可以分为两大类
(a) 代数预处理子 (Algebraic Preconditioner), 即仅仅根据所给的矩阵构造出来的预处理子.
我们这里只考虑代数预处理子, 即仅仅根据所给的系数矩阵来构造预处理方法. 这种预处理方法具有较好的通用性. 这类预处理子的一般来源有:
定常迭代法: 设有矩阵分裂 A = M − N A = M - N A = M − N M M M
直接法: 如不完全 LU 分解, 不完全 Cholesky 分解等;
近似逆: 即选取矩阵 P P P P − 1 ≈ A − 1 P^{-1} \approx A^{-1} P − 1 ≈ A − 1
块状结构：如基于Schur补的块对角矩阵，块三角矩阵等。
基于矩阵分裂的预处理子 考虑线性方程组
A x = b , A ∈ R n × n , A x = b, \quad A \in \mathbb {R} ^ {n \times n}, A x = b , A ∈ R n × n , 对 A A A
A = M − N (7.43) A = M - N \tag {7.43} A = M − N ( 7.43 ) 其中 M M M
x ( k + 1 ) = M − 1 N x ( k ) + M − 1 b = x ( k ) + ( M − 1 b − M − 1 A x ( k ) ) . k = 0 , 1 , 2 , … . x ^ {(k + 1)} = M ^ {- 1} N x ^ {(k)} + M ^ {- 1} b = x ^ {(k)} + \left(M ^ {- 1} b - M ^ {- 1} A x ^ {(k)}\right). \quad k = 0, 1, 2, \dots . x ( k + 1 ) = M − 1 N x ( k ) + M − 1 b = x ( k ) + ( M − 1 b − M − 1 A x ( k ) ) . k = 0 , 1 , 2 , … . 这等价于求解下面的方程组
M − 1 A x = M − 1 b . (7.44) M ^ {- 1} A x = M ^ {- 1} b. \tag {7.44} M − 1 A x = M − 1 b . ( 7.44 ) 这就是与矩阵分裂 (7.42) 相对应的预处理线性方程组。将 Krylov 子空间方法用于求解方程组 (7.43)，就得到预处理 Krylov 子空间方法。矩阵 M M M
理论上讲, 任何一个矩阵分裂都可以定义一个预处理子. 但为了使得预处理子能有很好的预处理效果, 往往需要其在一定意义下与 A A A
设 A = D − L − U A = D - L - U A = D − L − U D , − L , − U D, -L, -U D , − L , − U A A A D D D
Jacobi预处理子, 即取 A A A
G-S预处理子, 即取 A A A
P G S = D − L . P _ {\mathrm {G S}} = D - L. P GS = D − L . P S O R = 1 ω ( ω D − L ) . P _ {\mathrm {S O R}} = \frac {1}{\omega} (\omega D - L). P SOR = ω 1 ( ω D − L ) . P S S O R = 1 ω ( 2 − ω ) [ ( 1 − ω ) D + ω L ] D − 1 [ ( 1 − ω ) D + ω U ] . P _ {\text {S S O R}} = \frac {1}{\omega (2 - \omega)} \left[ (1 - \omega) D + \omega L \right] D ^ {- 1} \left[ (1 - \omega) D + \omega U \right]. P S S O R = ω ( 2 − ω ) 1 [ ( 1 − ω ) D + ω L ] D − 1 [ ( 1 − ω ) D + ω U ] . 由于SSOR对参数 ω \omega ω ω = 1 \omega = 1 ω = 1
P S G S = ( D − L ) D − 1 ( D − U ) . P _ {S G S} = (D - L) D ^ {- 1} (D - U). P SGS = ( D − L ) D − 1 ( D − U ) . HSS预处理子，即
P H S S = 1 2 α ( α I + H ) ( α I + S ) , P _ {\mathrm {H S S}} = \frac {1}{2 \alpha} (\alpha I + H) (\alpha I + S), P HSS = 2 α 1 ( α I + H ) ( α I + S ) , 其中 H = 1 2 ( A + A T ) , S = 1 2 ( A − A T ) . H = \frac{1}{2} (A + A^{\mathsf{T}}),S = \frac{1}{2} (A - A^{\mathsf{T}}). H = 2 1 ( A + A T ) , S = 2 1 ( A − A T ) .
不完全LU分解 对于大规模稀疏线性方程组, 有一类比较常用的预处理方法是不完全分解. 不完全分解首先由 Buleev [18, 19] 于上个世纪五十年代提出, 当时主要用来求解五点差分方程 (2D) 和七点差分方程 (3D) [77], 取得了非常好的效果, 但没有任何理论结果. 之后有更多学者对不完全分解进行了研究, 如 Varga (1960) [138], Oliphant (1962) [97], 等等 (参见 [77]). 但真正的突破是 1977 年 Meijerink 和 vander Vorst [94] 将不完全 Cholesky 分解用作共轭梯度法的预处理子, 提出了不完全 Cholesky 分解共轭梯度法 (ICCG). 它的出现, 才使得共轭梯度法成为求解大规模对称正定线性方程组的首选方法. 这时离共轭梯度法的提出已经过去了整整 25 年. 同时, 这也标志着 Krylov 子空间方法逐渐替代定常迭代方法, 成为求解大规模稀疏线性方程组的主流方法. 从那以后, 关于不完全矩阵分解的研究工作越来越多, 直至今天, 仍然是科学与工程计算领域的一个核心研究课题.
不完全 LU 分解的基本思想就是在对系数矩阵 A A A L L L U U U A A A A A A ( i , j ) ( i ≠ j ) (i,j) (i \neq j) ( i , j ) ( i = j ) L ( i , j ) L(i,j) L ( i , j ) U ( i , j ) U(i,j) U ( i , j )
A = L U + R , A = L U + R, A = LU + R , 其中 R R R I L U ( 0 ) ILU(0) I LU ( 0 )
为了增加 L U LU LU A A A L L L U U U L L L U U U L U LU LU