7.1 Krylov子空间 7.1.1 Arnoldi过程与Lanczos过程 设 A ∈ R n × n A \in \mathbb{R}^{n \times n} A ∈ R n × n r ∈ R n r \in \mathbb{R}^n r ∈ R n A A A r r r
K m ( A , r ) = span { r , A r , A 2 r , … , A m − 1 r } , m ≤ n , \mathcal {K} _ {m} (A, r) = \operatorname {s p a n} \{r, A r, A ^ {2} r, \dots , A ^ {m - 1} r \}, \quad m \leq n, K m ( A , r ) = span { r , A r , A 2 r , … , A m − 1 r } , m ≤ n , 通常简记为 κ m \kappa_{m} κ m
设解析解在 K m \mathcal{K}_m K m x ( m ) x^{(m)} x ( m ) r , A r , A 2 r , … , A m − 1 r r, Ar, A^2 r, \ldots, A^{m-1}r r , A r , A 2 r , … , A m − 1 r dim K m = m \dim \mathcal{K}_m = m dim K m = m v 1 , v 2 , … , v m v_1, v_2, \ldots, v_m v 1 , v 2 , … , v m K m \mathcal{K}_m K m K m \mathcal{K}_m K m x x x
x = y 1 v 1 + y 2 v 2 + ⋯ + y m v m = V m y , x = y _ {1} v _ {1} + y _ {2} v _ {2} + \dots + y _ {m} v _ {m} = V _ {m} y, x = y 1 v 1 + y 2 v 2 + ⋯ + y m v m = V m y , 其中 y = [ y 1 , y 2 , … , y m ] ⊤ y = [y_{1},y_{2},\ldots ,y_{m}]^{\top} y = [ y 1 , y 2 , … , y m ] ⊤ V m = [ v 1 , v 2 , … , v m ] V_{m} = [v_{1},v_{2},\dots ,v_{m}] V m = [ v 1 , v 2 , … , v m ] x ( m ) x^{(m)} x ( m )
就转化为下面两个子问题:
(1) 寻找一组合适的基 v 1 , v 2 , … , v m v_{1}, v_{2}, \ldots, v_{m} v 1 , v 2 , … , v m x ( m ) x^{(m)} x ( m ) y ( m ) = [ y 1 ( m ) , y 2 ( m ) , … , y m ( m ) ] ⊤ y^{(m)} = \left[y_1^{(m)}, y_2^{(m)}, \ldots, y_m^{(m)}\right]^\top y ( m ) = [ y 1 ( m ) , y 2 ( m ) , … , y m ( m ) ] ⊤
Arnoldi过程 首先考虑基的选取. 假定 r , A r , A 2 r , … , A m − 1 r r, Ar, A^2 r, \ldots, A^{m-1}r r , A r , A 2 r , … , A m − 1 r K m \mathcal{K}_m K m { r , A r , A 2 r , … , A m − 1 r } \{r, Ar, A^2 r, \ldots, A^{m-1}r\} { r , A r , A 2 r , … , A m − 1 r }
算法7.1.Arnoldi过程(MGS)v 1 = r / ∥ r ∥ 2 v_{1} = r / \| r \|_{2} v 1 = r /∥ r ∥ 2 j = 1 j = 1 j = 1 m − 1 m - 1 m − 1 z = A v j z = A v_{j} z = A v j i = 1 i = 1 i = 1 j j j h i , j = ( v i , z ) h_{i,j} = (v_{i}, z) h i , j = ( v i , z ) z = z − h i , j v i z = z - h_{i,j} v_{i} z = z − h i , j v i h j + 1 , j = ∥ z ∥ 2 h_{j+1,j} = \| z \|_{2} h j + 1 , j = ∥ z ∥ 2 h j + 1 , j = 0 h_{j+1,j} = 0 h j + 1 , j = 0 v j + 1 = z / h j + 1 , j v_{j+1} = z / h_{j+1,j} v j + 1 = z / h j + 1 , j
可以证明, Arnoldi 过程生成的向量 v 1 , v 2 , … , v m v_{1}, v_{2}, \ldots, v_{m} v 1 , v 2 , … , v m K m \mathcal{K}_{m} K m
引理7.1 如果Arnoldi过程不中断, 则
K m = span { v 1 , v 2 , … , v m } . \mathcal {K} _ {m} = \operatorname {s p a n} \left\{v _ {1}, v _ {2}, \dots , v _ {m} \right\}. K m = span { v 1 , v 2 , … , v m } . (留作练习, 归纳法)
记 V m = [ v 1 , v 2 , … , v m ] , V_{m} = [v_{1},v_{2},\ldots ,v_{m}], V m = [ v 1 , v 2 , … , v m ] ,
H m + 1 , m = [ h 1 , 1 h 1 , 2 h 1 , 3 … h 1 , m h 2 , 1 h 2 , 2 h 2 , 3 … h 2 , m h 3 , 2 h 3 , 3 … h 3 , m ⋱ ⋱ ⋮ ⋱ h m , m h m + 1 , m ] ∈ R ( m + 1 ) × m , H _ {m + 1, m} = \left[ \begin{array}{c c c c c} h _ {1, 1} & h _ {1, 2} & h _ {1, 3} & \dots & h _ {1, m} \\ h _ {2, 1} & h _ {2, 2} & h _ {2, 3} & \dots & h _ {2, m} \\ & h _ {3, 2} & h _ {3, 3} & \dots & h _ {3, m} \\ & & \ddots & \ddots & \vdots \\ & & & \ddots & h _ {m, m} \\ & & & & h _ {m + 1, m} \end{array} \right] \in \mathbb {R} ^ {(m + 1) \times m}, H m + 1 , m = h 1 , 1 h 2 , 1 h 1 , 2 h 2 , 2 h 3 , 2 h 1 , 3 h 2 , 3 h 3 , 3 ⋱ … … … ⋱ ⋱ h 1 , m h 2 , m h 3 , m ⋮ h m , m h m + 1 , m ∈ R ( m + 1 ) × m , 则由Arnoldi过程可知
h j + 1 , j v j + 1 = A v j − h 1 , j v 1 − h 2 , j v 2 − ⋯ − h j , j v j , h _ {j + 1, j} v _ {j + 1} = A v _ {j} - h _ {1, j} v _ {1} - h _ {2, j} v _ {2} - \dots - h _ {j, j} v _ {j}, h j + 1 , j v j + 1 = A v j − h 1 , j v 1 − h 2 , j v 2 − ⋯ − h j , j v j , 即
A v j = ∑ i = 1 j + 1 h i , j v i = V m + 1 [ h 1 , j ⋮ h j + 1 , j 0 ⋮ 0 ] = V m + 1 H m + 1 , m ( : , j ) . A v _ {j} = \sum_ {i = 1} ^ {j + 1} h _ {i, j} v _ {i} = V _ {m + 1} \left[ \begin{array}{c} h _ {1, j} \\ \vdots \\ h _ {j + 1, j} \\ 0 \\ \vdots \\ 0 \end{array} \right] = V _ {m + 1} H _ {m + 1, m} (:, j). A v j = i = 1 ∑ j + 1 h i , j v i = V m + 1 h 1 , j ⋮ h j + 1 , j 0 ⋮ 0 = V m + 1 H m + 1 , m ( : , j ) . 所以有
A V m = V m + 1 H m + 1 , m = V m H m + h m + 1 , m v m + 1 e m ⊤ , (7.1) A V _ {m} = V _ {m + 1} H _ {m + 1, m} = V _ {m} H _ {m} + h _ {m + 1, m} v _ {m + 1} e _ {m} ^ {\top}, \tag {7.1} A V m = V m + 1 H m + 1 , m = V m H m + h m + 1 , m v m + 1 e m ⊤ , ( 7.1 ) 其中 H m H_{m} H m H m + 1 , m H_{m + 1,m} H m + 1 , m m m m H m = H m + 1 , m ( 1 : m , 1 : m ) , e m = [ 0 , … , 0 , 1 ] T ∈ H_{m} = H_{m + 1,m}(1:m,1:m),e_{m} = [0,\dots ,0,1]^{\mathsf{T}}\in H m = H m + 1 , m ( 1 : m , 1 : m ) , e m = [ 0 , … , 0 , 1 ] T ∈ R m \mathbb{R}^m R m V m V_{m} V m V m T V_{m}^{\mathsf{T}} V m T
V m ⊤ A V m = H m . (7.2) V _ {m} ^ {\top} A V _ {m} = H _ {m}. \tag {7.2} V m ⊤ A V m = H m . ( 7.2 ) 等式 (7.1) 和 (7.2) 是 Arnoldi 过程的两个重要性质. 这两个性质也可以通过下图来表示.
A V m = V m + 1 H m + 1 , m = V m H m + A \boxed {V _ {m}} = V _ {m + 1} \boxed {H _ {m + 1, m}} = V _ {m} \boxed {H _ {m}} + A V m = V m + 1 H m + 1 , m = V m H m + 需要指出的是, 如果 r , A r , A 2 r , … , A m − 1 r r, Ar, A^2 r, \ldots, A^{m-1} r r , A r , A 2 r , … , A m − 1 r
定理7.2 如果Arnoldi过程在第 k k k h k + 1 , k = 0 h_{k+1,k} = 0 h k + 1 , k = 0 k < m k < m k < m A V k = V k H k AV_k = V_kH_k A V k = V k H k K k \mathcal{K}_k K k A A A
(留作课外自习, 直接利用等式 (7.1))
Lanczos过程 如果 A A A H m H_{m} H m T m T_{m} T m
T m = [ α 1 β 1 β 1 ⋱ ⋱ ⋱ ⋱ β m − 1 β m − 1 α m ] . (7.3) T _ {m} = \left[ \begin{array}{c c c c} \alpha_ {1} & \beta_ {1} & & \\ \beta_ {1} & \ddots & \ddots & \\ & \ddots & \ddots & \beta_ {m - 1} \\ & & \beta_ {m - 1} & \alpha_ {m} \end{array} \right]. \tag {7.3} T m = α 1 β 1 β 1 ⋱ ⋱ ⋱ ⋱ β m − 1 β m − 1 α m . ( 7.3 ) 与Arnoldi过程类似, 我们有下面的性质
A V m = V m T m + β m v m + 1 e m ⊤ , (7.4) A V _ {m} = V _ {m} T _ {m} + \beta_ {m} v _ {m + 1} e _ {m} ^ {\top}, \tag {7.4} A V m = V m T m + β m v m + 1 e m ⊤ , ( 7.4 ) V m T A V m = T m . (7.5) V _ {m} ^ {\mathsf {T}} A V _ {m} = T _ {m}. \tag {7.5} V m T A V m = T m . ( 7.5 ) 考察关系式 (7.4) 两边的第 j j j
β j v j + 1 = A v j − α j v j − β j − 1 v j − 1 , j = 1 , 2 , … , \beta_ {j} v _ {j + 1} = A v _ {j} - \alpha_ {j} v _ {j} - \beta_ {j - 1} v _ {j - 1}, \quad j = 1, 2, \dots , β j v j + 1 = A v j − α j v j − β j − 1 v j − 1 , j = 1 , 2 , … , 这里我们令 v 0 = 0 v_{0} = 0 v 0 = 0 β 0 = 0 \beta_0 = 0 β 0 = 0
算法7.2.Lanczos过程 1: Set v 0 = 0 v_{0} = 0 v 0 = 0 β 0 = 0 \beta_0 = 0 β 0 = 0 v 1 = r / ∥ r ∥ 2 v_{1} = r / \| r \|_{2} v 1 = r /∥ r ∥ 2 j = 1 j = 1 j = 1 m − 1 m - 1 m − 1 z = A v j z = A v_{j} z = A v j α j = ( v j , z ) \alpha_{j} = (v_{j},z) α j = ( v j , z ) z = z − α j v j − β j − 1 v j − 1 z = z - \alpha_{j}v_{j} - \beta_{j - 1}v_{j - 1} z = z − α j v j − β j − 1 v j − 1 β j = ∥ z ∥ 2 \beta_{j} = \| z\|_{2} β j = ∥ z ∥ 2 β j = 0 \beta_{j} = 0 β j = 0
10: end ifv j + 1 = z / β j v_{j + 1} = z / \beta_j v j + 1 = z / β j
可以证明, 由Lanczos过程得到的向量组 { v 1 , v 2 , … , v m } \{v_{1}, v_{2}, \ldots, v_{m}\} { v 1 , v 2 , … , v m }
定理7.3 设 { v 1 , v 2 , … , v m } \{v_{1}, v_{2}, \ldots, v_{m}\} { v 1 , v 2 , … , v m }
( v i , v j ) = δ i j = { 1 , i = j , 0 , i ≠ j , i , j = 1 , 2 , … , m . (v _ {i}, v _ {j}) = \delta_ {i j} = \left\{ \begin{array}{l l} 1, & i = j, \\ 0, & i \neq j, \end{array} \right. i, j = 1, 2, \ldots , m. ( v i , v j ) = δ ij = { 1 , 0 , i = j , i = j , i , j = 1 , 2 , … , m . (留作练习)
7.1.2 Krylov子空间方法一般格式 Krylov子空间迭代法的一般过程 (1) 令 m = 1 m = 1 m = 1 K m \mathcal{K}_m K m x ( 0 ) + K m x^{(0)} + \mathcal{K}_m x ( 0 ) + K m m ← m + 1 m \gets m + 1 m ← m + 1
在实际计算中, 为了充分利用迭代初值中所包含的有用信息, 我们通常是在仿射空间 x ( 0 ) + K m x^{(0)} + \mathcal{K}_m x ( 0 ) + K m
算法7.3.Krylov子空间迭代算法 1: 选取初始向量 x ( 0 ) x^{(0)} x ( 0 ) r 0 = b − A x ( 0 ) r_0 = b - Ax^{(0)} r 0 = b − A x ( 0 ) v 1 = r 0 / ∥ r 0 ∥ 2 v_{1} = r_{0} / \| r_{0}\|_{2} v 1 = r 0 /∥ r 0 ∥ 2 x ( 1 ) ∈ x ( 0 ) + K 1 = x ( 0 ) + span { v 1 } x^{(1)} \in x^{(0)} + \mathcal{K}_1 = x^{(0)} + \operatorname{span}\{v_1\} x ( 1 ) ∈ x ( 0 ) + K 1 = x ( 0 ) + span { v 1 } x ( 1 ) x^{(1)} x ( 1 ) m = 2 m = 2 m = 2 n n n v m v_{m} v m x ( m ) ∈ x ( 0 ) + K m = x ( 0 ) + span { v 1 , … , v m } x^{(m)} \in x^{(0)} + \mathcal{K}_m = x^{(0)} + \operatorname{span}\{v_1, \ldots, v_m\} x ( m ) ∈ x ( 0 ) + K m = x ( 0 ) + span { v 1 , … , v m } x ( m ) x^{(m)} x ( m )
11: 终止迭代
子空间的选取和更新问题可以通过Krylov子空间来解决. 下面需要考虑如何寻找方程组在仿射空间 x ( 0 ) + K m x^{(0)} + \mathcal{K}_m x ( 0 ) + K m x ( m ) x^{(m)} x ( m ) x ( m ) x^{(m)} x ( m )
我们很自然地想到用近似解与真解之间的距离来衡量“最佳”，即使得 x ( m ) − x ∗ x^{(m)} - x_{*} x ( m ) − x ∗ ∥ x ( m ) − x ∗ ∥ 2 \| x^{(m)} - x_{*}\|_{2} ∥ x ( m ) − x ∗ ∥ 2 x ∗ x_{*} x ∗
实用的“最佳”判别方式有：
(1) 使得 ∥ r m ∥ 2 = ∥ b − A x ( m ) ∥ 2 \| r_{m}\|_{2} = \| b - Ax^{(m)}\|_{2} ∥ r m ∥ 2 = ∥ b − A x ( m ) ∥ 2 r m = b − A x ( m ) r_m = b - Ax^{(m)} r m = b − A x ( m ) A A A A A A A A A ∥ r m ∥ A − 1 \| r_{m}\|_{A^{-1}} ∥ r m ∥ A − 1
∥ r m ∥ A − 1 ≜ ( r m ⊤ A − 1 r m ) 1 2 = ( ( b − A x ( m ) ) ⊤ A − 1 ( b − A x ( m ) ) ) 1 2 = ( ( A − 1 b − x ( m ) ) T A ( A − 1 b − x ( m ) ) ) 1 2 = ( ( x ∗ − x ( m ) ) T A ( x ∗ − x ( m ) ) ) 1 2 = ∥ x ∗ − x ( m ) ∥ A . \begin{array}{l} \left\| r _ {m} \right\| _ {A ^ {- 1}} \triangleq \left(r _ {m} ^ {\top} A ^ {- 1} r _ {m}\right) ^ {\frac {1}{2}} = \left(\left(b - A x ^ {(m)}\right) ^ {\top} A ^ {- 1} \left(b - A x ^ {(m)}\right)\right) ^ {\frac {1}{2}} \\ = \left(\left(A ^ {- 1} b - x ^ {(m)}\right) ^ {\mathsf {T}} A \left(A ^ {- 1} b - x ^ {(m)}\right)\right) ^ {\frac {1}{2}} \\ = \left(\left(x _ {*} - x ^ {(m)}\right) ^ {\mathsf {T}} A \left(x _ {*} - x ^ {(m)}\right)\right) ^ {\frac {1}{2}} \\ = \left\| x _ {*} - x ^ {(m)} \right\| _ {A}. \\ \end{array} ∥ r m ∥ A − 1 ≜ ( r m ⊤ A − 1 r m ) 2 1 = ( ( b − A x ( m ) ) ⊤ A − 1 ( b − A x ( m ) ) ) 2 1 = ( ( A − 1 b − x ( m ) ) T A ( A − 1 b − x ( m ) ) ) 2 1 = ( ( x ∗ − x ( m ) ) T A ( x ∗ − x ( m ) ) ) 2 1 = x ∗ − x ( m ) A . 因此CG方法也可以看作是极小化能量范数 ∥ x ∗ − x ( m ) ∥ A \| x_{*} - x^{(m)}\|_{A} ∥ x ∗ − x ( m ) ∥ A