Krylov 子空間法

本文的閱讀等級:中級

A 為一 n\times n 階複矩陣,\mathbf{b} 為一 n 維非零向量。1931年,俄國應用數學家、海軍工程師克雷洛夫 (Aleksey Krylov) 提出一個創新的想法[1]:運用向量序列 \{\mathbf{b},A\mathbf{b},\ldots,A^{r-1}\mathbf{b}\},稱為 Krylov 序列,計算 A 的特徵多項式。Krylov 序列的擴張稱為 Krylov 子空間 (或循環子空間),記為

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

明顯地,\mathcal{K}_r(A,\mathbf{b})\mathbb{C}^n 的一個子空間,故必存在最小正整數 k 使得 A^k\mathbf{b} 可表示為 \mathbf{b},A\mathbf{b},\ldots,A^{k-1}\mathbf{b} 的線性組合。如果

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

定義 k 次多項式

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

因為 f(A)\mathbf{b}=\mathbf{0},我們說 f(t)\mathbf{b} 相對於 A 的消滅多項式 (annihilating polynomial)。運用類似最小多項式 (minimal polynomial) 的論證方式可證明 (見“最小多項式 (上)”):給定任何矩陣—向量對 (A,\mathbf{b}),次數最小的首一 (領先係數為 1) 消滅多項式唯一存在,記為 m(t),就有 m(A)\mathbf{x}=\mathbf{0}。既然 k 是最小正整數使得 \mathbf{b},A\mathbf{b},\ldots,A^k\mathbf{b} 線性相關,可知 m(t)=f(t)

 
下面的定理說明矩陣 A 的最小多項式 m_A(t) 與向量 \mathbf{b} 相對於 A 的最小多項式 m(t) 之間的關係 (證明見“最小多項式的計算方法”)。

 
最小公倍式定理:令 \{\mathbf{b}_1,\ldots,\mathbf{b}_n\}\mathbb{C}^n 的一組基底。若 m_i(t)\mathbf{b}_i 相對於 n\times n 階矩陣 A 的最小多項式,則 A 的最小多項式 m_A(t)m_1(t),\ldots,m_n(t) 的最小公倍式。

除了直接應用於計算 A 的最小多項式 m_A(t),這個定理也可以用來推導相伴矩陣 (companion matrix) 的最小多項式和特徵多項式 (另一個證法見“多項式的相伴矩陣”)。

 
相伴矩陣

給定首一多項式

\displaystyle p(t)=t^n+a_{n-1}t^{n-1}+\cdots+a_1t+a_0

我們定義 p(t) 的相伴矩陣為下列 n\times n 階矩陣:

\displaystyle  C=\begin{bmatrix}  0&0&\cdots&0&-a_0\\  1&0&\cdots&0&-a_1\\  \vdots&\ddots&\ddots&\vdots&\vdots\\  0&\cdots&1&0&-a_{n-2}\\  0&0&\cdots&1&-a_{n-1}  \end{bmatrix}

下面證明首一多項式 p(t) 為其相伴矩陣 C 的特徵多項式和最小多項式[2]。欲證明 \det (tI-C)=p(t),寫出 C=N-\mathbf{a}\mathbf{e}_n^T,其中

\displaystyle  N=\begin{bmatrix}  0&&&\\  1&\ddots&&\\  &\ddots&\ddots&\\  &&1&0  \end{bmatrix},~~\mathbf{a}=\begin{bmatrix}  a_0\\  a_1\\  \vdots\\  a_{n-1}  \end{bmatrix},~~\mathbf{e}_n=\begin{bmatrix}  0\\  \vdots\\  0\\  1  \end{bmatrix}

因為 N^n=0

\displaystyle\begin{aligned}  (tI-N)^{-1}&=t^{-1}\left(I-\frac{N}{t}\right)^{-1}=t^{-1}\left(I+\frac{N}{t}+\cdots+\frac{N^{n-1}}{t^{n-1}}\right)\\  &=\frac{I}{t}+\frac{N}{t^2}+\cdots+\frac{N^{n-1}}{t^n}.  \end{aligned}

使用矩陣和的行列式公式 (見“矩陣和之行列式 (上)”)

\displaystyle  \det(A+\mathbf{u}\mathbf{v}^T)=(\det A)(1+\mathbf{v}^TA^{-1}\mathbf{u})

推得

\displaystyle\begin{aligned}  \det(tI-C)&=\det(tI-N+\mathbf{a}\mathbf{e}_n^T)\\  &=\det(tI-N)\left(1+\mathbf{e}_n^T(tI-N)^{-1}\mathbf{a}\right)\\  &=t^n\left(1+\mathbf{e}_n^T\left(\frac{I}{t}+\frac{N}{t^2}+\cdots+\frac{N^{n-1}}{t^n}\right)\mathbf{a}\right)\\  &=t^n+\mathbf{e}_n^TI\mathbf{a}t^{n-1}+\mathbf{e}_n^TN\mathbf{a}t^{n-2}+\cdots+\mathbf{e}_n^TN^{n-1}\mathbf{a}\\  &=t^n+a_{n-1}t^{n-1}+a_{n-2}t^{n-2}+\cdots+a_0=p(t),  \end{aligned}

故證明伴隨矩陣 C 的特徵多項式為 p_C(t)=p(t)。考慮 \mathbb{C}^n 的標準基底 \boldsymbol{\beta}=\{\mathbf{e}_1,\ldots,\mathbf{e}_n\}。令 m_i(t)\mathbf{e}_i 相對於 C 的最小多項式,1\le i\le n。因為 C\mathbf{e}_i=\mathbf{e}_{i+1}1\le i\le n-1

\displaystyle  \{\mathbf{e}_1,C\mathbf{e}_1,\ldots,C^{n-1}\mathbf{e}_1\}=\{\mathbf{e}_1,\mathbf{e}_2,\ldots,\mathbf{e}_n\}

構成一線性獨立集,即有

\displaystyle  C^n\mathbf{e}_1=CC^{n-1}\mathbf{e}_1=C\mathbf{e}_n=-\mathbf{a}=-\sum_{i=0}^{n-1}a_i\mathbf{e}_{i+1}=-\sum_{i=0}^{n-1}a_iC^i\mathbf{e}_1

根據定義,m_1(t)=p(t)。我們知道 m_1(t) 整除 m_1(t),\ldots,m_n(t) 的最小公倍式,由最小公倍式定理可知 m_1(t) 整除 m_C(t)。另一方面,最小多項式 m_C(t) 必整除特徵多項式 p_C(t),即證得 m_C(t)=p(t)

 
Krylov 子空間法

最小公倍式定理表明任一 n 維非零向量 \mathbf{b} 相對於 n\times n 階矩陣 A 的最小多項式 m(t) 整除 A 的最小多項式 m_A(t),因此也整除 A 的特徵多項式 p_A(t)。自然地,我們設想 p_A(t) 可否表示為某些 \mathbf{b}_i 相對於 A 的最小多項式 m_i(t) 的乘積?事實上,克雷洛夫正是根據這個思路發展出一個特徵多項式算法,以下稱為 Krylov 子空間法。

 
首先,隨意挑選一個非零向量 \mathbf{b} 作為 Krylov 序列的「種子」。令 m_1(t)=t^k-\sum_{i=0}^{k-1}c_it^i\mathbf{b} 相對於 A 的最小多項式。將線性獨立的 Krylov 序列 \{\mathbf{b},A\mathbf{b},\ldots,A^{k-1}\mathbf{b}\} 合併為一 n\times k 階矩陣 K_1=\begin{bmatrix}  \mathbf{b}&A\mathbf{b}&\cdots&A^{k-1}\mathbf{b}  \end{bmatrix},稱為 Krylov 矩陣,\text{rank}K_1=k。令 C_1 為對應多項式 m_1(t)k\times k 階相伴矩陣。直接計算可證明

\displaystyle\begin{aligned}  K_1C_1&=\begin{bmatrix}  \mathbf{b}&A\mathbf{b}&\cdots&A^{k-2}\mathbf{b}&A^{k-1}\mathbf{b}  \end{bmatrix}\begin{bmatrix}  0&0&\cdots&0&c_0\\  1&0&\cdots&0&c_1\\  \vdots&\ddots&\ddots&\vdots&\vdots\\  0&\cdots&1&0&c_{k-2}\\  0&0&\cdots&1&c_{k-1}  \end{bmatrix}\\  &=\begin{bmatrix}  A\mathbf{b}&A^2\mathbf{b}&\cdots&A^{k-1}\mathbf{b}&\sum_{i=0}^{k-1}c_iA^i\mathbf{b}  \end{bmatrix}\\  &=\begin{bmatrix}  A\mathbf{b}&A^2\mathbf{b}&\cdots&A^{k-1}\mathbf{b}&A^k\mathbf{b}  \end{bmatrix}\\  &=AK_1.  \end{aligned}

k=n,則 K_1^{-1}AK_1=C_1,也就是說,A 相似於 C_1。因為相似矩陣有相同的特徵多項式,且 m_1(t)C_1 的特徵多項式,立知 p_A(t)=m_1(t)。若 k<n,將 K_1 擴充為一 n\times n 階可逆矩陣 \tilde{K}_1=\begin{bmatrix}  K_1&K'_1  \end{bmatrix},其中 n\times(n-k) 階分塊 K'_1 可由 Steinitz 替換原則求得 (見“基底與維數常見問答集”)。使用上面導出的關係式,

\displaystyle  A\tilde{K}_1=\begin{bmatrix}  AK_1&AK'_1  \end{bmatrix}=\begin{bmatrix}  K_1C_1&A{K}'_1  \end{bmatrix}=\begin{bmatrix}  K_1&{K}'_1  \end{bmatrix}\begin{bmatrix}  C_1&X\\  0&A_2  \end{bmatrix}=\tilde{K}_1\begin{bmatrix}  C_1&X\\  0&A_2  \end{bmatrix}

其中 \displaystyle  \begin{bmatrix}  X\\  A_2  \end{bmatrix}=\tilde{K}_1^{-1}A{K}'_1,故得

\displaystyle  \tilde{K}_1^{-1}A\tilde{K}_1=\begin{bmatrix}  C_1&X\\  0&A_2  \end{bmatrix}

換句話說,A 相似於 \begin{bmatrix}  C_1&X\\  0&A_2  \end{bmatrix},即知

\displaystyle  p_A(t)=m_1(t)\det(tI_{n-k}-A_2)

繼續對 (n-k)\times(n-k) 階矩陣 A_2 重複上述步驟,可得最小多項式 m_2(t)。若 A_2 的 Krylov 序列產生可逆的 Krylov 矩陣 K_2,則 K_2^{-1}A_2K_2=C_2p_A(t)=m_1(t)m_2(t),其中 C_2 是對應 m_2(t) 的相伴矩陣。使用以上結果,寫出

\displaystyle  \begin{bmatrix}  C_1&X\\  0&A_2  \end{bmatrix}=\begin{bmatrix}  C_1&X\\  0&K_2C_2K_2^{-1}  \end{bmatrix}=\begin{bmatrix}  I&0\\  0&K_2  \end{bmatrix}\begin{bmatrix}  C_1&XK_2\\  0&C_2  \end{bmatrix}\begin{bmatrix}  I&0\\  0&K_2^{-1}  \end{bmatrix}

也就有

\displaystyle  K^{-1}AK=\begin{bmatrix}  C_1&XK_2\\  0&C_2  \end{bmatrix}

其中 K=\tilde{K}_1\begin{bmatrix}  I&0\\  0&K_2  \end{bmatrix}。若 A_2 的 Krylov 序列產生不可逆 Krylov 矩陣 K_2,則 p_A(t)=m_1(t)m_2(t)\det(tI-A_3),其中 A_3 的階數必小於 n-k。執行前述演算程序,直到產生可逆 Krylov 矩陣即告結束。假設經過 q 個步驟,最後 A 相似於一分塊上三角矩陣,稱為 Frobenius 矩陣 (見“不使用行列式的特徵值和特徵向量算法 (中)”):

\displaystyle  K^{-1}AK=\begin{bmatrix}  C_1&\cdots&\ast\\  &\ddots&\vdots\\  &&C_q  \end{bmatrix}=H

其中 C_i 是相伴矩陣,對應的特徵多項式為 m_i(t)。矩陣 A 的特徵多項式為 p_A(t)=m_1(t)\cdots m_q(t) (另一個推導方式請見“利用循環子空間計算特徵多項式”)。

 
下面舉一例展示 Krylov 子空間算法:

\displaystyle  A=\left[\!\!\begin{array}{rrcr}  1&-3&5&11\\  2&-1&3&1\\  -1&1&0&0\\  1&-1&2&2  \end{array}\!\!\right]

(1) 選擇 \mathbf{b}=(1,0,0,0)^T,計算得到

\displaystyle  A\mathbf{b}=\left[\!\!\begin{array}{r}  1\\  2\\  -1\\  1  \end{array}\!\!\right],~~A^2\mathbf{b}=\left[\!\!\begin{array}{r}  1\\  -2\\  1\\  -1  \end{array}\!\!\right]

因為 A^2\mathbf{b}+A\mathbf{b}-2\mathbf{b}=\mathbf{0},故 \mathbf{b} 相對於 A 的最小多項式為 m_1(t)=t^2+t-2,相伴矩陣為 C_1=\left[\!\!\begin{array}{cr}  0&2\\  1&-1  \end{array}\!\!\right]

(2) 寫出

\displaystyle  \tilde{K}_1=\begin{bmatrix}  K_1&K'_1  \end{bmatrix}=\left[\!\!\begin{array}{crccc}  1&1&\vline&0&0\\  0&2&\vline&1&0\\  0&-1&\vline&0&1\\  0&1&\vline&0&0  \end{array}\!\!\right]

並計算

\displaystyle  \begin{bmatrix}  X\\  A_2  \end{bmatrix}=\tilde{K}_1^{-1}AK'_1=\left[\!\!\begin{array}{rr}  -2&3\\  -1&2\\\hline  1&-1\\  0&2  \end{array}\!\!\right]

(3) 選擇 \mathbf{b}_2=(0,1)^T,計算得到

\displaystyle  A_2\mathbf{b}_2=\left[\!\!\begin{array}{r}  -1\\  2  \end{array}\!\!\right],~~A_2^2\mathbf{b}_2=\left[\!\!\begin{array}{r}  -3\\  4  \end{array}\!\!\right]

因為 A_2^2\mathbf{b}_2-3A_2\mathbf{b}_2+2\mathbf{b}_2=\mathbf{0},故 \mathbf{b}_2 相對於 A_2 的最小多項式為 m_2(t)=t^2-3t+2,相伴矩陣為 C_2=\left[\!\!\begin{array}{cr}  0&-2\\  1&3  \end{array}\!\!\right]。因為 K_2=\left[\!\!\begin{array}{cr}  0&-1\\  1&2  \end{array}\!\!\right] 是可逆矩陣,停止計算,所求的特徵多項式為 p_A(t)=m_1(t)m_2(t)=(t^2+t-2)(t^2-3t+2)

 
觀察發現所有的相伴矩陣 C_i 都是上 Hessenberg 矩陣 (見“特殊矩陣 (19):Hessenberg 矩陣”),以 5\times 5 階矩陣為例,形式如下:

\displaystyle  \begin{bmatrix}  \ast&\ast&\ast&\ast&\ast\\  \ast&\ast&\ast&\ast&\ast\\  0&\ast&\ast&\ast&\ast\\  0&0&\ast&\ast&\ast\\  0&0&0&\ast&\ast  \end{bmatrix}

顯然,包含主對角分塊 C_i 的分塊上三角矩陣 H 也是一上 Hessenberg 矩陣。所以說,Krylov 子空間法的作為即在透過數個 Krylov 序列,建立一相似變換使得 A 相似於一上 Hessenberg 矩陣。就理論的角度而言,Krylov 子空間法確實頗具吸引力,但實際執行時卻會遭遇困難。Krylov 序列的生成機制與 Power 迭代法相同,A^k\mathbf{b} 終將收斂至 A 的一個特徵向量 (對應絕對值最大的特徵值,見“特徵向量是甚麼物,恁麼來?”)。問題出在這裡,當 k 增大時,Krylov 序列生成的向量方向逐漸趨於一致,Krylov 矩陣 K 跟著傾向病態 (ill-conditioned,見“病態系統”),嚴重危害相似變換的數值穩定性。此外,在許多應用場合,A 是一個稀疏矩陣,然而 Krylov 矩陣的多數元通常不為零,這意味未加修飾的 Krylov 子空間法難免耗費多餘的計算。

 
1950年代,QR 分解的出現使得 Krylov 子空間法出現了轉機,通過對 n\times k 階 Krylov 矩陣 K_1=\begin{bmatrix}  \mathbf{b}&A\mathbf{b}&\cdots&A^{k-1}\mathbf{b}  \end{bmatrix} 執行 QR 分解 K_1=QR 可克服病態問題,其中 Q 是一 n\times k 階矩陣滿足 Q^\ast Q=I_kR 是一 k\times k 階可逆上三角矩陣 (見“線代膠囊──QR 分解”)。將 QR 分解代入 AK_1=K_1C_1,可得 AQR=QRC_1,則

\displaystyle  Q^\ast AQ=RC_1R^{-1}=H_1

其中 H_1 雖不再是相伴矩陣,但仍為上 Hessenberg 矩陣 (上 Hessenberg 矩陣與上三角矩陣的乘積仍是上 Hessenberg 矩陣)。引入 QR 分解的價值在於 Q 有單範正交 (orthonormal) 的行向量,因此 Q 是一個完美的良置 (well-conditioned) 矩陣。當 k<nQ 並非方陣,Q^\ast AQ=H_1 不是相似變換,也就是說,AH_1 的特徵值不同。不過這點遺憾並未阻饒 Krylov 子空間法的發展,日後我們將陸續介紹 Krylov 子空間法衍生的特徵值和線性方程算法。特別一提的是,從1970年代起,Krylov 子空間法已逐漸成為線性方程迭代算法的主流[3]

 
參考來源:
[1] 維基百科:Krylov subpsace
[2] Carl Meyer, Matrix Analysis and Applied Linear Algebra, SIAM, 2000, pp 648.
[3] The Best of the 20th Century: Editors Name Top 10 Algorithms

相關閱讀:
廣告
本篇發表於 特徵分析, 線性代數專欄 並標籤為 , , , , , 。將永久鏈結加入書籤。

9 Responses to Krylov 子空間法

  1. 張浩軒 說道:

    感謝這篇對Krylov子空間法的解釋 收穫很大 希望未來老師能寫一些關於Arnoldi和GMRES的方法

    • ccjou 說道:

      好的,沒想到這麼冷門的主題也有讀者關注。

      • 張浩軒 說道:

        感謝周老師,因為最近在讀Trefethen的Numerical Linear Algebra中的Iterative chapter,一直無法有觀念貫通的感覺。我主要是無法理解用Arnoldi是如何loacte eigenvalue(Trefethen p 261)。 還有在Trefethen p255:Arnoldi在第n階時,找到Hessenberg matrix,那此Hessenberg的eigenvalue要如何找?而此Hessenberg的eigenvalue和欲求A矩陣的eigenvalue有多近似的關係?實際運作上如何知道n要多大,才足夠逼近欲求A矩陣的eigenvalue?
        我已看過老師的"利用循環子空間計算特徵多項式"、"最小多項式"、"Krylov子空間法" 但還是不能理解我上述的疑問。希望老師有時間的話可以指點迷津,感謝

  2. Xin Liu 說道:

    “日後我們將陸續介紹 Krylov 子空間法衍生的特徵值和線性方程算法”,我是大陆的学生,看到老师的文章实在是受益匪浅,非常希望老师继续给出关于Krylov子空间法求解线性方程的算法,十分渴望,忘老师不吝赐教!

    • ccjou 說道:

      你說的應該是GMRES算法,以及Conjugate Gradient算法(應用於正定矩陣)。等我回覆了先前另一位讀者的問題後再討論這個主題。

  3. Vahi Chen 說道:

    周老师,如果可能的话,希望能听您讲讲线性方程组的迭代解法(FOM、GMRES、Lanczos、CG等),很是期待!

  4. sowter 說道:

    老師您好

    請問這方法 跟 固定點跌代法做比較

    這種方法好處在哪

發表迴響

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

WordPress.com Logo

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

Twitter picture

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

Facebook照片

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

Google+ photo

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

連結到 %s