Krylov 子空間法──線性方程的數值解法 (一):Arnoldi 與 Lanczos 算法

本文的閱讀等級:高級

A 為一 n\times n 階複矩陣,\mathbf{b}\in\mathbb{C}^n 為一非零向量。向量序列 \{\mathbf{b},A\mathbf{b},\ldots,A^{j-1}\mathbf{b}\} 稱為 Krylov 序列,此序列所生成的子空間稱為 Krylov 子空間 (見“Krylov 子空間法”),記為

\displaystyle  \mathcal{K}_j(A,\mathbf{b})=\text{span}\{\mathbf{b},A\mathbf{b},\ldots,A^{j-1}\mathbf{b}\}

因為 \mathcal{K}_j(A,\mathbf{b}) 是有限維空間 \mathbb{C}^n 的一個子空間,當 j 不斷增大時,Krylov 序列 \{\mathbf{b},A\mathbf{b},\ldots,A^{j-1}\mathbf{b}\} 最終會是一個線性相關集。設 k 為最小的正整數使得 \mathcal{K}_{k+1}(A,\mathbf{b})=\mathcal{K}_k(A,\mathbf{b}),也就是說 \{\mathbf{b},A\mathbf{b},\ldots,A^{k-1}\mathbf{b}\} 為一個線性獨立集且 A^k\mathbf{b}\in \mathcal{K}_k(A,\mathbf{b})。因此,存在唯一數組 c_0,c_1,\ldots,c_{k-1} 滿足

A^k\mathbf{b}=c_0\mathbf{b}+c_1A\mathbf{b}+\cdots+c_{k-1}A^{k-1}\mathbf{b}

定義 k 次多項式

\displaystyle  m(t)=t^k-c_{k-1}t^{k-1}-\cdots-c_1t-c_0

因為 m(A)\mathbf{b}=\mathbf{0},我們稱 m(t)\mathbf{b} 相對於 A 的最小 (消滅) 多項式 (minimal polynomial)。以下考慮 A 為可逆矩陣。我們可以斷定 c_0\neq 0,否則有 A^{k-1}\mathbf{b}=c_1\mathbf{b}+\cdots+c_{k-1}A^{k-2}\mathbf{b},即存在次數小於 k 的消滅多項式。改寫上面的線性組合式,

\displaystyle  A\left(\frac{A^{k-1}\mathbf{b}-c_{k-1}A^{k-2}\mathbf{b}-\cdots-c_1\mathbf{b}}{c_0}\right)=\mathbf{b}

故線性方程 A\mathbf{x}=\mathbf{b} 的精確解為

\displaystyle  \mathbf{x}_\ast=\frac{1}{c_0}\left(A^{k-1}\mathbf{b}-c_{k-1}A^{k-2}\mathbf{b}-\cdots-c_1\mathbf{b}\right)

理論上,一旦算出 Krylov 序列 \{\mathbf{b},A\mathbf{b},\ldots,A^{k-1}\mathbf{b}\} 以及 \mathbf{b} 相對於 A 的最小多項式 m(t) 便得到線性方程 A\mathbf{x}=\mathbf{b} 的解。不過,實現這個想法必須克服一個困難。當 j 增大時,Krylov 序列 \{\mathbf{b},A\mathbf{b},\ldots,A^{r-1}\mathbf{b}\} 生成的向量方向逐漸趨於一致 (對應 A 的絕對值最大的特徵值的特徵向量,見“Power 迭代法”),也就是說,Krylov 序列幾乎線性相關 (稱為病態),將嚴重危害數值計算的穩定性。

 
解決 Krylov 子空間基底蘊含的病態問題的最直接方法是正交化 Krylov 向量序列。給定 n\times n 階矩陣 A 與種子向量 \mathbf{v}\mathbf{v}\neq\mathbf{0},假設 \{\mathbf{q}_1,\ldots,\mathbf{q}_j\} 是 Krylov 子空間 \mathcal{K}_j(A,\mathbf{v})=\hbox{span}\{\mathbf{v},A\mathbf{v},\ldots,A^{j-1}\mathbf{v}\} 的一組單範正交基底 (orthonormal basis)。如果 A^j\mathbf{v}\notin\mathcal{K}_j(A,\mathbf{v}),欲建構 \mathbf{q}_{j+1},使用 Gram-Schmidt 正交化程序:扣除 A^j\mathbf{v} 至基底向量 \mathbf{q}_1,\ldots,\mathbf{q}_j 的投影量 (見“Gram-Schmidt 正交化與 QR 分解”),

\displaystyle  \mathbf{y}_j=A^j\mathbf{v}-\sum_{i=1}^j(\mathbf{q}_i^\ast A^j\mathbf{v})\mathbf{q}_i

再予以歸一化,\mathbf{q}_{j+1}=\mathbf{y}_j/\Vert\mathbf{y}_j\Vert。因此,\{\mathbf{q}_1,\ldots,\mathbf{q}_j,\mathbf{q}_{j+1}\}\mathcal{K}_{j+1}(A,\mathbf{v}) 的一組單範正交基底。上面的正交化算式可以再化簡。設 \mathbf{q}_1=\mathbf{v}/\Vert\mathbf{v}\VertA\mathbf{q}_1=d_1\mathbf{q}_1+d_2\mathbf{q}_2,可得

\displaystyle  \begin{aligned}  \mathcal{K}_{j+1}(A,\mathbf{v})&=\hbox{span}\{\mathbf{v},A\mathbf{v},A^2\mathbf{v},\ldots,A^j\mathbf{v}\}\\  &=\hbox{span}\{\mathbf{q}_1,A\mathbf{q}_1,A^2\mathbf{q}_1,\ldots,A^j\mathbf{q}_1\}\\  &=\hbox{span}\{\mathbf{q}_1,d_1\mathbf{q}_1+d_2\mathbf{q}_2,A(d_1\mathbf{q}_1+d_2\mathbf{q}_2),\ldots,A^{j-1}(d_1\mathbf{q}_1+d_2\mathbf{q}_2)\}\\  &=\hbox{span}\{\mathbf{q}_1,\mathbf{q}_2,A\mathbf{q}_2,\ldots,A^{j-1}\mathbf{q}_2\}\\  &=\cdots\\  &=\hbox{span}\{\mathbf{q}_1,\mathbf{q}_2,\ldots,\mathbf{q}_j,A\mathbf{q}_j\}.  \end{aligned}

在 Gram-Schmidt 正交化過程中,我們可用 A\mathbf{q}_j 取代 A^j\mathbf{v},算式如下:

\displaystyle  \mathbf{u}_j=A\mathbf{q}_j-\sum_{i=1}^j(\mathbf{q}_i^\ast A\mathbf{q}_j)\mathbf{q}_i

\mathbf{u}_j=\mathbf{0},停止計算,表示 \mathcal{K}_{j+1}(A,\mathbf{v})=\mathcal{K}_j(A,\mathbf{v});否則令 \mathbf{q}_{j+1}=\mathbf{u}_j/\Vert\mathbf{u}_j\Vert,並重複上述步驟繼續擴充 Krylov 子空間。因為 \mathbf{q}_{j+1}\perp\hbox{span}\{\mathbf{q}_1,\ldots,\mathbf{q}_j\}

\displaystyle  \mathbf{q}_{j+1}^\ast\mathbf{u}_j=\mathbf{q}_{j+1}^\ast\left(A\mathbf{q}_j-\sum_{i=1}^j(\mathbf{q}_i^\ast A\mathbf{q}_j)\mathbf{q}_i\right)=\mathbf{q}_{j+1}^\ast A\mathbf{q}_j

另一方面,\mathbf{q}_{j+1}^\ast\mathbf{u}_j=\mathbf{u}_j^\ast\mathbf{u}_j/\Vert\mathbf{u}_j\Vert=\Vert\mathbf{u}_j\Vert,故 \Vert\mathbf{u}_j\Vert=\mathbf{q}_{j+1}^\ast A\mathbf{q}_j。合併以上結果,單範正交基底向量的關係式如下:

\displaystyle\begin{aligned}  A\mathbf{q}_j&=\mathbf{u}_j+\sum_{i=1}^j(\mathbf{q}_i^\ast A\mathbf{q}_j)\mathbf{q}_i\\  &=\Vert\mathbf{u}_j\Vert\mathbf{q}_{j+1}+\sum_{i=1}^j(\mathbf{q}_i^\ast A\mathbf{q}_j)\mathbf{q}_i\\  &=\sum_{i=1}^{j+1}(\mathbf{q}_i^\ast A\mathbf{q}_j)\mathbf{q}_i\\  &=\sum_{i=1}^{j+1}h_{ij}\mathbf{q}_i,  \end{aligned}

其中 h_{ij}=\mathbf{q}_i^\ast A\mathbf{q}_j1\le i\le j+1。特別的是,h_{j+1,j}=\mathbf{q}_{j+1}^\ast A\mathbf{q}_j=\Vert\mathbf{u}_j\Vert

 
為了提升數值計算的穩定性,我們以改良 Gram-Schmidt 正交化 (簡稱 MGS) 取代古典 Gram-Schmidt 正交化,兩者的差異在於 h_{ij} 的計算方式 (見“Gram-Schmidt 正交化改良版”)。底下的演算法稱為 Arnoldi-MGS 算法:對於一 n\times n 階矩陣 A 和非零向量 \mathbf{v},找出 Krylov 子空間 \mathcal{K}_k(A,\mathbf{v}) 的一組單範正交基底 \{\mathbf{q}_1,\ldots,\mathbf{q}_k\}

Arnoldi-MGS 算法


  1. \mathbf{q}_1\leftarrow\mathbf{v}/\Vert\mathbf{v}\Vert
  2. 對於 j=1,2,\ldots,n
  3.     \mathbf{u}\leftarrow A\mathbf{q}_j
  4.     對於 i=1,\ldots,j
  5.         h_{ij}\leftarrow\mathbf{q}_i^\ast\mathbf{u}
  6.         \mathbf{u}\leftarrow\mathbf{u}-h_{ij}\mathbf{q}_i
  7.     h_{j+1,j}\leftarrow\Vert\mathbf{u}\Vert
  8.     若 h_{j+1,j}=0,停止計算。
  9.     \mathbf{q}_{j+1}\leftarrow\mathbf{u}/h_{j+1,j}

當 Arnoldi-MGS 算法於 j=k 停止時,種子向量 \mathbf{v} 相對於 A 的最小多項式的次數為 k\mathbf{v} 生成的最大線性獨立 Krylov 序列所生成的子空間為 \mathcal{K}_k(A,\mathbf{v})=\hbox{span}\{\mathbf{q}_1,\ldots,\mathbf{q}_k\}

 
Arnoldi-MGS 算法可用簡明的矩陣運算表示。令 n\times j 階矩陣 Q_j=\begin{bmatrix}  \mathbf{q}_1&\cdots&\mathbf{q}_j  \end{bmatrix} 以及 j\times j 階上 Hessenberg 矩陣 H_j=[h_{il}],即 h_{il}=0i>l+1 (見“特殊矩陣 (19):Hessenberg 矩陣”)。寫出單範正交基底的關鍵式

\displaystyle  A\mathbf{q}_l=\begin{bmatrix}  \mathbf{q}_1&\cdots&\mathbf{q}_l  \end{bmatrix}\begin{bmatrix}  h_{1l}\\  \vdots\\  h_{ll}  \end{bmatrix}+\mathbf{q}_{l+1}h_{l+1,l},~~l=1,\ldots,j

合併可得

\displaystyle\begin{aligned}  AQ_j&=\begin{bmatrix}  A\mathbf{q}_1&A\mathbf{q}_2&\cdots &A\mathbf{q}_j  \end{bmatrix}\\  &=\begin{bmatrix}  \mathbf{q}_1&\mathbf{q}_2&\cdots&\mathbf{q}_j  \end{bmatrix}\begin{bmatrix}  h_{11}&h_{12}& h_{13}&\cdots&h_{1j}\\  h_{21}&h_{22}&h_{23}&\cdots&h_{2j}\\  0&h_{32}&h_{33}&\cdots&h_{3j}\\  \vdots&\ddots&\ddots&\ddots&\vdots\\  0&\cdots&0&h_{j,j-1}&h_{jj}  \end{bmatrix}+\begin{bmatrix}  \mathbf{0}&\cdots&\mathbf{0}&\mathbf{q}_{j+1}h_{j+1,j}  \end{bmatrix}\\  &=Q_jH_j+\mathbf{q}_{j+1}h_{j+1,j}\mathbf{e}_j^T,  \end{aligned}

稱為 Arnoldi 關係式,其中 \mathbf{e}^T_j=(0,\ldots,0,1)j 維單位向量。Arnoldi 關係式也可以表示為

AQ_j=Q_{j+1}\tilde{H_j}

其中 Q_{j+1}=\begin{bmatrix}  \mathbf{q}_1&\cdots&\mathbf{q}_{j+1}  \end{bmatrix}\tilde{H}_j(j+1)\times j 階矩陣,

\displaystyle  \tilde{H}_j=\begin{bmatrix}  H_j\\  h_{j+1,j}\mathbf{e}_j^T  \end{bmatrix}=\begin{bmatrix}  h_{11}&h_{12}& h_{13}&\cdots&h_{1j}\\  h_{21}&h_{22}&h_{23}&\cdots&h_{2j}\\  0&h_{32}&h_{33}&\cdots&h_{3j}\\  \vdots&\ddots&\ddots&\ddots&\vdots\\  0&\cdots&0&h_{j,j-1}&h_{jj}\\  0&\cdots&\cdots&0&h_{j+1,j}  \end{bmatrix}

因為 Q_j^\ast Q_j=I_jQ_j^\ast\mathbf{q}_{j+1}=\mathbf{0},上式左乘 Q_j^\ast,可得

Q_j^\ast AQ_j=Q_j^\ast Q_{j+1}\tilde{H_j}=H_j

上式可以解釋為 A 至 Krylov 子空間 \mathcal{K}_j(A,\mathbf{v}) 的正交投影,Q_j 的行 (column) 向量為 \mathcal{K}_j(A,\mathbf{v}) 的基底。Arnoldi-MGS 算法的作用在產生 Krylov 子空間基底和正交投影 H_j。當計算程序於 j=k 停止時,h_{k+1,k}=0,即有

AQ_k=Q_kH_k

 
A 是 Hermitian 矩陣,即 A^\ast=A,則

H_j^\ast=(Q_j^\ast AQ_j)^\ast=Q_j^\ast AQ_j=H_j

上 Hessenberg 矩陣 H_j 是 Hermitian,表示 H_j 是三對角矩陣 (見“特殊矩陣 (11):三對角矩陣”)。為容易識別,我們以 T_j 替換 H_j,如下:

\displaystyle  T_j=\begin{bmatrix}  \alpha_1&\beta_1& & & \\  \beta_1&\alpha_2&\beta_2&&\\   & \beta_2&\alpha_3&\ddots&\\  && \ddots&\ddots&\beta_{j-1}\\  &&&\beta_{j-1}&\alpha_j  \end{bmatrix}

使用 \mathbf{q}_j^\ast A\mathbf{q}_j=\alpha_j\mathbf{q}_{j-1}^\ast A\mathbf{q}_j=\beta_{j-1}\mathbf{q}_i^\ast A\mathbf{q}_j=0i<j-1,Gram-Schmidt 正交化計算可以大幅簡化:

\displaystyle  \mathbf{u}_j=A\mathbf{q}_j-\sum_{i=1}^j(\mathbf{q}_i^\ast A\mathbf{q}_j)\mathbf{q}_i=A\mathbf{q}_j-\alpha_j\mathbf{q}_j-\beta_{j-1}\mathbf{q}_{j-1}

類似前面的推導步驟,

\begin{aligned}  \Vert\mathbf{u}_j\Vert&=\mathbf{q}_{j+1}^\ast\mathbf{u}_j=\mathbf{q}_{j+1}^\ast\left(A\mathbf{q}_j-\alpha_j\mathbf{q}_j-\beta_{j-1}\mathbf{q}_{j-1}\right)\\  &=\mathbf{q}_{j+1}^\ast A\mathbf{q}_j=(\mathbf{q}_j^\ast A\mathbf{q}_{j+1})^\ast=\bar\beta_j,\end{aligned}

可知 \beta_j 是非負實數。因此,\beta_j=\Vert\mathbf{u}_j\Vert\mathbf{u}_j=\beta_j\mathbf{q}_{j+1}。合併以上結果,

A\mathbf{q}_j=\beta_{j-1}\mathbf{q}_{j-1}+\alpha_j\mathbf{q}_j+\beta_j\mathbf{q}_{j+1}

我們將推導結果整理成一個算法,稱為對稱 Lanczos 算法 (或簡稱 Lanczos 算法):對於一 n\times n 階 Hermitian 矩陣 A 和非零向量 \mathbf{v},找出 Krylov 子空間 \mathcal{K}_k(A,\mathbf{v}) 的一組單範正交基底 \{\mathbf{q}_1,\ldots,\mathbf{q}_k\}

對稱 Lanczos 算法


  1. \mathbf{q}_1\leftarrow\mathbf{v}/\Vert\mathbf{v}\Vert\mathbf{q}_0\leftarrow\mathbf{0}\beta_0\leftarrow 0
  2. 對於 j=1,2,\ldots,n
  3.     \mathbf{u}\leftarrow A\mathbf{q}_j
  4.     \alpha_j\leftarrow \mathbf{q}_j^\ast\mathbf{u}
  5.     \mathbf{u}\leftarrow\mathbf{u}-\alpha_j\mathbf{q}_j-\beta_{j-1}\mathbf{q}_{j-1}
  6.     \beta_j\leftarrow\Vert\mathbf{u}\Vert
  7.     若 \beta_j=0,停止計算。
  8.     \mathbf{q}_{j+1}\leftarrow\mathbf{u}/\beta_j

 
當 Lanczos 算法於 j=k 停止時,種子向量 \mathbf{v} 相對於 A 的最小多項式的次數為 k\mathbf{v} 生成的最大線性獨立 Krylov 序列所生成的子空間為 \mathcal{K}_k(A,\mathbf{v})=\hbox{span}\{\mathbf{q}_1,\ldots,\mathbf{q}_k\}。類似 Arnoldi 關係式,Lanczos 關係式為

AQ_j=Q_jT_j+\mathbf{q}_{j+1}\beta_j\mathbf{e}_j^T

AQ_j=Q_{j+1}\tilde{T}_j

其中 \tilde{T}_j=\begin{bmatrix}  T_j\\  \beta_j\mathbf{e}_j^T  \end{bmatrix}。Lanczos 算法停止時,\beta_k=0,可得

AQ_k=Q_kT_k

 
A 為大尺寸矩陣,我們可以採用迭代法求解。定常迭代法是較古老且相對簡單的線性方程迭代算法,但僅保證有限的 A 矩陣具備收斂性,如對角佔優矩陣與正定矩陣 (見“定常迭代法──線性方程迭代解法”)。即便 Krylov 序列未直接提供線性方程 A\mathbf{x}=\mathbf{b} 的精確解 \mathbf{x}_\ast=A^{-1}\mathbf{b},但至少我們知道精確解 \mathbf{x}_\ast 屬於 Krylov 子空間 \mathcal{K}_k(A,\mathbf{b}),其中 k\mathbf{b} 相對 A 的最小多項式的次數。下文將介紹通過 Krylov 序列的線性方程的迭代解法,統稱為 Krylov 子空間法,適用於任何型態的可逆矩陣,是目前數值線性代數中最成功的一類算法。

繼續閱讀:
廣告
本篇發表於 線性代數專欄, 數值線性代數 並標籤為 , , , , , , , , 。將永久鏈結加入書籤。

2 Responses to Krylov 子空間法──線性方程的數值解法 (一):Arnoldi 與 Lanczos 算法

  1. Ran 說道:

    请问老师 subspace iteration 和 Lanczos method 有什么异同呢?

    • ccjou 說道:

      我還沒寫完,後續會介紹。簡單講,Krylov subspace methods由兩個子空間(U_jV_j)定義:
      1. x_j-x_0\in U_j,
      2. b-Ax_j\perp V_j.

      兩類方法的最大差異在於V_j

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s