QR 演算法 (上)

本文的閱讀等級:中級

1920至1930年代,量子力學之矩陣力學 (matrix mechanics) 首次將矩陣理論推向應用領域 (見“海森堡不確定原理的矩陣證明”)。此後,特徵值與特徵向量遂成為線性代數所處理的四個核心問題之一 (見“線性代數的核心問題”)。令 A 為一個 n\times n 階實矩陣或複矩陣,\lambda 為一特徵值,\mathbf{x} 為對應的 (非零) 特徵向量,滿足特徵方程 A\mathbf{x}=\lambda\mathbf{x}。1940年代,許多數學家和工程師相繼加入特徵值與特徵向量算法的研究行列。最明顯的特徵值算法,即線性代數教科書所述的標準方法,包含二個步驟:先計算 A 的特徵多項式的係數,

\displaystyle  p_A(t)=\det(tI-A)=t^n+a_{n-1}t^{n-1}+\cdots+a_1t+a_0

再求出 p_A(t) 的所有根。多項式的求根問題在當時曾經是一個熱門的研究題目。不幸的是,除了小尺寸矩陣 (n<10),運用有限位元的電子計算機實現上述算法最終是一場無可逃避的災難,因為多項式求根在本質上是一個病態問題,根的位置很容易受到多項式係數的微小擾動而發生劇烈改變 (Wilkinson 多項式說明了這種情況,詳見“Power 迭代法”)。然而,在多數的情況下,矩陣的特徵值對於 n^2 個元的微小擾動並不敏感。這表示將矩陣 An^2 個元替換為特徵多項式 p_A(t)n 個係數過度壓縮資料,故而危害特徵值的計算。

 
既然直接求解特徵多項式的根不可行,研究工作者於是改採其他途徑。他們所考慮的問題形式並非特徵方程 A\mathbf{x}=\lambda\mathbf{x},而是 Schur 分解定理 (見“矩陣三角化的 Schur 定理”):任一方陣 A 必可三角化為 A=UTU^\ast,其中 T 是一個上三角矩陣,U 是一個么正 (unitary) 矩陣,滿足 U^\ast=U^{-1}。因為 A 相似於 T,兩個相似矩陣擁有相同的特徵值,而上三角矩陣的主對角元為其特徵值,據此可知 T 的主對角元即為 A 的特徵值。換句話說,我們的目標要尋找一個么正矩陣 U 使得 T=U^\ast AU 為上三角矩陣。那麼該從何處著手呢?假設 A=XY,其中 X 是可逆矩陣,則 B=YX 相似於 A,因為 B=YX=X^{-1}AX (見“AB 和 BA 有何關係?”)。這個觀察結果指引了一個方向,接下來的挑戰在於找出適當的 XY 分解形式與設計算法程序。1954年,瑞士數學家、計算機科學家盧蒂斯豪賽爾 (Heinz Rutishauser) 提出一個基於 LU 分解的迭代算法,名為 LR 演算法[1]。這是特徵值算法的一個重大進步,並為後續發展豎立典範。1961年,英國計算機科學家弗朗西斯 (John G.F. Francis) 和俄國數學家卡布諾夫斯凱雅 (Vera Kublanovskaya) 獨立提出了 QR 演算法。經過五十年的演化,QR 演算法逐漸成為當今特徵值和特徵向量數值計算的骨幹,並被視為20世紀最重要的十個演算法之一[2]。另一方面,由於 LR 演算法的穩定性不如 QR 演算法,今天已鮮少被提及。

 
顧名思義,QR 演算法建立在 QR 分解之上。給定一個 n\times n 矩陣 A,QR 分解記為 A=QR,其中 Qn\times n 階么正矩陣 (實矩陣則為正交矩陣),Q^\ast=Q^{-1}Rn\times n 階上三角矩陣。任一矩陣皆存在 QR 分解,而且可逆矩陣有唯一的標準 QR 分解 (即 R 的主對角元皆為正數,見“線代膠囊──QR 分解”)。目前常用的 QR 分解算法包括 Gram-Schmidt 正交化、Householder 變換,以及 Givens 旋轉 (見“QR 分解的數值計算方法比較”)。基本 QR 演算法的主要思路是利用迭代技巧,以 QR 分解衍生 (相似) 上三角矩陣:初始化 A_0=A,先分解出 A_{k-1}=Q_kR_k,再計算 A_k=R_kQ_k,詳細如下。

基本 QR 演算法


  1. A_0\leftarrow AU_0\leftarrow I
  2. 對於 k=1,2,\ldots
  3.     A_{k-1}\rightarrow Q_kR_k (QR 分解)
  4.     A_k\leftarrow R_kQ_k
  5.     U_k\leftarrow U_{k-1}Q_k
  6. T\leftarrow A_\inftyU\leftarrow U_\infty

 
如前所述,A_k=R_kQ_k 相似於 A_{k-1}=Q_kR_k,合併二式可得

\displaystyle  A_k=R_kQ_k=Q_k^\ast A_{k-1}Q_k

繼續往回迭代,

\displaystyle\begin{aligned}  A_k&=Q_k^{\ast}A_{k-1}Q_k=Q_k^{\ast}Q_{k-1}^{\ast} A_{k-2}Q_{k-1}Q_k=\cdots\\  &=Q_k^{\ast}Q_{k-1}^\ast\cdots Q_1^{\ast} A_0 Q_1\cdots Q_{k-1}Q_k=U_k^{\ast} A_0 U_k,\end{aligned}

上面定義了 U_k=Q_1\cdots Q_k。不難驗證 U_k^\ast U_k=I。矩陣序列 \{A_k\}_0^\infty 構成一個「相似家族」,因此特徵值維持不變。然而,\{A_k\} 果真收斂至一個上三角矩陣嗎?

 
基本 QR 演算法的收斂條件並不嚴苛。假設 A 的特徵值 \lambda_1,\ldots,\lambda_n 滿足

\displaystyle  \vert\lambda_1\vert>\vert\lambda_2\vert > \cdots > \vert\lambda_n\vert>0

底下證明 \{Q_k\} 收斂至一個對角矩陣,因此證明 \{A_k\} 收斂至一個上三角矩陣。我們將推導 A^k 的兩個 QR 分解表達式[3]。令上三角矩陣 P_k=R_kR_{k-1}\cdots R_1。使用 U_kA_k=A_0U_k,即 Q_1\cdots Q_k A_k=A_0Q_1\cdots Q_k,計算得到

\displaystyle\begin{aligned}  U_kP_k&=Q_1\cdots Q_{k-1}Q_kR_kR_{k-1}\cdots R_1\\  &=Q_1\cdots Q_{k-1}A_{k-1}R_{k-1}\cdots R_1\\  &=A_0Q_1\cdots Q_{k-1}R_{k-1}\cdots R_1\\  &=A_0U_{k-1}P_{k-1}=\cdots\\  &=A_0^{k-1}U_1P_1=A_0^{k-1}Q_1R_1=A_0^{k}.\end{aligned}

上式說明 A^k=A_0^k 的 QR 分解為 A^k=U_kP_k。這個式子留待稍後使用。

 
因為 A 有相異的特徵值,故可對角化為 A=X\Lambda X^{-1}=X\Lambda Y,其中 Y=X^{-1}\Lambda=\hbox{diag}(\lambda_1,\ldots,\lambda_n),也就有 A^k=X\Lambda^k X^{-1}=X\Lambda^kY。令 X 的 QR 分解與 Y 的 LU 分解分別為

X=QR,~~Y=LW

其中 RW 是上三角矩陣,L=[l_{ij}] 是下三角矩陣且主對角元為 1Q 是么正矩陣。注意,R 的可逆性來自 X 的可逆性。可逆矩陣 Y 存在 LU 分解的條件為 Y 的所有 k\times k 階領先主子陣 (leading principal submatrix) 都是可逆矩陣 (見“LU 分解”)。注意,Q, R, L, W 不因迭代程序而改變,也就是與 k 無關。將上面兩式代入 A^k=X\Lambda^k Y,可得

\displaystyle  A^k=QR\Lambda^kLW=QR\left(\Lambda^kL\Lambda^{-k}\right)\Lambda^kW

其中 \Lambda^{-k}=(\Lambda^k)^{-1}\Lambda^kL\Lambda^{-k} 是單位下三角矩陣 (即主對角元為 1),且

\displaystyle  \left(\Lambda^kL\Lambda^{-k}\right)_{ij}=l_{ij}\left(\frac{\lambda_i}{\lambda_j}\right)^k,~~i > j

\Lambda^kL\Lambda^{-k}=I+E_k。當 i>j,已知條件表明 \left|\lambda_i/\lambda_j\right|<1,故 \lim_{k\to\infty}E_k=0。合併以上結果,

\displaystyle\begin{aligned}  A^k&=QR(I+E_k)\Lambda^kW\\  &=Q(I+RE_kR^{-1})R\Lambda^kW\\  &=Q(I+F_k)R\Lambda^kW,  \end{aligned}

其中 F_k=RE_kR^{-1}。明顯地,\lim_{k\to\infty}F_k=\lim_{k\to\infty}RE_kR^{-1}=0。設 I+F_k 的 QR 分解為 I+F_k=\tilde{Q}_k\tilde{R}_k。當 k\to\infty,我們有 \tilde{Q}_k\tilde{R}_k\to I,即知 \tilde{Q}_k\to I\tilde{R}_k\to I。因此最後得到

\displaystyle  A^k=(Q\tilde{Q}_k)(\tilde{R}_kR\Lambda^kW)

其中 Q\tilde{Q}_k 是么正矩陣,\tilde{R}_kR\Lambda^kW 是上三角矩陣。上式很容易轉換為標準 QR 分解,

\displaystyle  A^k=(Q\tilde{Q}_k D_k)(D_k^{-1}\tilde{R}_kR\Lambda^kW)

其中 D_k=\hbox{diag}(e^{j\theta_{k1}},\ldots,e^{j\theta_{kn}}) 是么正矩陣,j=\sqrt{-1},選擇旋轉角度 \theta_{k1},\ldots,\theta_{kn} 使得 D_k^{-1}\tilde{R}_kR\Lambda^kW 的主對角元為正數。

 
比較冪矩陣 A^k 的兩種 QR 分解表達式:

\displaystyle  A^k=U_kP_k,~~  A^k=(Q\tilde{Q}_k D_k)(D_k^{-1}\tilde{R}_kR\Lambda^kW)

可得 U_k=Q\tilde{Q}_k{D}_k。當 k\to\infty\tilde{Q}_k\to I,因此 U_k\to Q{D}_k。寫出

\displaystyle  Q_k=(Q_1\cdots Q_{k-1})^\ast(Q_1\cdots Q_{k-1})Q_k=U^\ast_{k-1}U_k=D_{k-1}^\ast \tilde{Q}_{k-1}^\ast \tilde{Q}_kD_k

即知 \{Q_k\} 收斂至 {D}_{k-1}^\ast{D}_k。因為 \tilde{R}_kR 的主對角元已經是正數,但 \Lambda^kW 則否,故可推論 {D}_{k-1}^\ast{D}_k=D_{k-1}^{-1}D_k=\hbox{diag}(e^{j\phi_1},\ldots,e^{j\phi_n}),其中 \phi_i=\measuredangle\lambda_i1\le i\le n。所以,

\displaystyle T=\lim_{k\to\infty}A_k=\lim_{k\to\infty}R_kQ_k=\lim_{k\to\infty}R_k\hbox{diag}(e^{j\phi_1},\ldots,e^{j\phi_n})

為一個上三角矩陣。

 
下文將探討 QR 演算法的理論與實現,我們要回答兩個問題:第一,QR 演算法的基本原理與幾何意義是甚麼?第二,針對一般矩陣,基本 QR 演算法的收斂速度並不理想,如何修改使之成為收斂速度佳且穩定性高的實用算法?

 
註解:
[1] 維基百科:QR algorithm
[2] 弗朗西斯 (John G.F. Francis),英國人,生於1934年 (維基百科)。1954年,任職於隸屬英國政府的國家研究發展機構 (National Research Development Corporation,NRDC)。1955-56年,就讀劍橋大學,未獲得學位,隨後返回 NRDC 工作。1961-62年,先後發表兩篇 QR 演算法的理論與實務期刊論文。1961年,弗朗西斯離開 NRDC,從此脫離數值計算的研究領域。2007年,美國史丹佛教授格魯布 (Gene H. Golub) 聯繫上弗朗西斯,他全然不知自己早年的研究成果受到廣泛地引用,並對 QR 演算法被譽為20世紀最重要的十個演算法之一深表驚訝。弗朗西斯已經退休,2007年他在空中大學攻讀學位。(Gene H. Golub and Frank Uhlig, The QR algorithm: 50 years later–its genesis by John Francis and Vera Kublanovskaya, and subsequent developments, IMA Journal of Numerical Analysis, 2009, 29 (3): 467–485.)
[3] 取自 James H. Wilkinson 的 The Algebraic Eigenvalue Problem,Clarendon Press,1965,頁517-519。

This entry was posted in 線性代數專欄, 數值線性代數 and tagged , , , . Bookmark the permalink.

3 則回應給 QR 演算法 (上)

  1. kaerfuka 說:

    周老师,看完意犹未尽,为什么只有(上),没找到(下)?

    • ccjou 說:

      你沒找到(下)是因為我還沒寫。當初寫完(上)之後,我發現自己並未掌握QR算法的全貌故而停筆,不意日子久了便忘記這件事。如果你想深入了解,請參考下文:
      Understanding the QR algorithm
      https://www.jstor.org/stable/2029539?seq=1#page_scan_tab_contents

      我目前工作繁忙,估計短時間內不會有(下)。

      • Meiyue Shao 說:

        我觉得Watkins的文章不太适合理解QR算法。我的建议是先找本矩阵计算的教材看看,把算法全貌了解一下,等到QR算法以及子空间迭代法都熟悉了再去看Watkins的文章。而且即便如此仍然很难完全领悟QR算法的精妙。

發表迴響

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

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / 變更 )

Twitter picture

You are commenting using your Twitter account. Log Out / 變更 )

Facebook照片

You are commenting using your Facebook account. Log Out / 變更 )

Google+ photo

You are commenting using your Google+ account. Log Out / 變更 )

連結到 %s