QR 分解的數值計算方法比較

本文的閱讀等級:中級

QR 分解是數值線性代數中具備多種用途的計算工具,主要應用於線性方程、最小平方法和特徵值問題。常見的 QR 分解的計算方法包括 Householder 變換、Givens 旋轉以及 Gram-Schmidt 正交法。本文先回顧 QR 分解的主要性質,介紹 QR 分解於計算最小平方解的應用,並討論上述三種算法的運算量與數值穩定性。

 
A 是一個 m\times n 實矩陣,以下假設 m\ge n,因為 QR 分解的多數應用均符合此情況。矩陣 A 的 QR 分解如下:

A=QR=\begin{bmatrix}    Q_1&Q_2    \end{bmatrix}\begin{bmatrix}    R_1\\    0    \end{bmatrix}=Q_1R_1

其中 Q=\begin{bmatrix}    Q_1&Q_2    \end{bmatrix}m\times m 階正交矩陣,Q^T=Q^{-1}Q_1m\times n 階分塊,Q_2m\times(m-n) 階分塊,R_1=[r_{ij}]n\times n 階上三角矩陣。為便於區分,有人稱完整型 A=QR 為「胖」QR 分解,經濟型 A=Q_1R_1 為「瘦」QR 分解。

 
正交矩陣 Q 給出 A 的行空間 (column space) C(A) 的一組單範正交 (orthonormal) 基底與正交補餘 C(A)^{\perp} (即左零空間 N(A^T)) 的一組單範正交基底。將 AQ 以行向量表示,A=\begin{bmatrix}    \mathbf{a}_1&\cdots&\mathbf{a}_n    \end{bmatrix}Q=\begin{bmatrix}    \mathbf{q}_1&\cdots&\mathbf{q}_n&\mathbf{q}_{n+1}&\cdots\mathbf{q}_m    \end{bmatrix}。若 \mathrm{rank}A=n,亦即 A 有線性獨立的行向量,則

\mathrm{span}\{\mathbf{a}_1,\ldots,\mathbf{a}_j\}=\mathrm{span}\{\mathbf{q}_1,\ldots,\mathbf{q}_j\},~~j=1,\ldots,n

也就推得

C(A)=C(Q_1),~~C(A)^{\perp}=C(Q_2)

證明於下。直接乘開 A=QR 可得

\displaystyle  \mathbf{a}_j=\sum_{i=1}^jr_{ij}\mathbf{q}_i,~~j=1,\ldots,n

故知 \mathrm{span}\{\mathbf{a}_1,\ldots,\mathbf{a}_j\}\subseteq\mathrm{span}\{\mathbf{q}_1,\ldots,\mathbf{q}_j\}。因為 AQ 有線性獨立行向量,\mathrm{span}\{\mathbf{a}_1,\ldots,\mathbf{a}_j\}\mathrm{span}\{\mathbf{q}_1,\ldots,\mathbf{q}_j\} 的維數都等於 j,由此推論 \mathrm{span}\{\mathbf{a}_1,\ldots,\mathbf{a}_j\}=\mathrm{span}\{\mathbf{q}_1,\ldots,\mathbf{q}_j\}。當 j=n,上式等同於 C(A)=C(Q_1)。因為 Q 是正交矩陣,\mathrm{span}\{\mathbf{q}_{1},\ldots,\mathbf{q}_n\}^{\perp}=\mathrm{span}\{\mathbf{q}_{n+1},\ldots,\mathbf{q}_m\},也就證明了 C(A)^{\perp}=C(Q_2)

 
\mathrm{rank}A=nA 有唯一的「瘦」QR 分解 A=Q_1R_1。利用 Cholesky 分解即可證明。利用 Q_1^TQ_1=I_n,就有

A^TA=(Q_1R_1)^T(Q_1R_1)=R_1^TQ_1^TQ_1R_1=R_1^TR_1

得知 G=R_1^TA^TA 的 Cholesky 分解矩陣,即 A^TA=GG^TG 是下三角矩陣且主對角元都是正數。因為 A 滿行秩,R_1 是可逆的,A^TA 是正定矩陣,故 A^TA 必定有唯一的 Cholesky 分解 (見“Cholesky 分解”),因此證明 Q_1=AR_1^{-1} 唯一存在。

 
計算最小平方解是 QR 分解的一個典型應用。考慮 A\mathbf{x}=\mathbf{b},最小平方解 \hat{\mathbf{x}} 滿足正規方程 (normal equation) A^TA\hat{\mathbf{x}}=A^T\mathbf{b}。若 \mathrm{rank}A=n (若此條件不成立,則存在無窮多組解),則 \hat{\mathbf{x}}=(A^TA)^{-1}A^T\mathbf{b}。然而,A^TA 很容易引入誤差。為方便說明,假設 A 是一個方陣,採用2-範數計算條件數 (後面詳述):

\begin{aligned}  \kappa(A^TA)&=\Vert A^TA\Vert_2\cdot\Vert (A^TA)^{-1}\Vert_2=\frac{\sigma_{\max}(A^TA)}{\sigma_{\min}(A^TA)}\\  &=\frac{(\sigma_{\max}(A))^2}{(\sigma_{\min}(A))^2}=\left(\Vert A\Vert_2\cdot\Vert A^{-1}\Vert_2\right)^2=(\kappa(A))^2,  \end{aligned}

其中 \sigma_{\max}(X)\sigma_{\min}(X) 分別表示矩陣 X 最大與最小奇異值 (計算詳見“每週問題 April 13, 2015”)。交互乘積 A^TA 的條件數是 A 的條件數平方,換句話說,A^TA 很可能使原本為良置的系統變成病態。因此,我們應該直接處理矩陣 A 而非 A^TA,作法如下:將 QR 分解 A=Q_1R_1 代入正規方程式 A^TA\hat{\mathbf{x}}=A^T\mathbf{b},就有

R^T_1R_1\hat{\mathbf{x}}=R_1^TQ_1^T\mathbf{b}

左乘 (R_1^T)^{-1},立得

R_1\hat{\mathbf{x}}=Q_1^T\mathbf{b}

其中 R_1 是上三角矩陣,使用反向迭代即可解出 \hat{\mathbf{x}}。在求最小平方解的過程中,我們其實並不需要明確地知道 Q_1 矩陣,此性質對減少運算量有很大的助益。

 
一般通用的 QR 分解計算方法有 Householder 變換 (見“Householder 變換於 QR 分解的應用”)、Givens 旋轉 (見“Givens 旋轉於 QR 分解的應用”) 和 Gram-Schmidt 正交化 (見“Gram-Schmidt 正交化改良版”)。這些算法的基本原理介紹請參閱相關文章,在此不再贅述。欲進一步了解算法的詳細編程,則請參閱文獻[1],下面列出 m\times n 階矩陣 A 的 QR 分解運算量,包含加法和乘法運算。

  • Householder 變換:計算 R_1 使用 2mn^2-2n^3/3 次浮點運算,若要得到 Q_1,還必須再增加 2mn^2-2n^3/3 次浮點運算。
  • Givens 旋轉:計算 R_1 使用 3mn^2-n^3 次浮點運算,略大於 Householder 變換。
  • 改良 Gram-Schmidt 正交化 (MGS):此法一併得到 Q_1R_1,共需 2mn^2 次浮點運算。所以,如欲得到 Q_1R_1,MGS 約耗費 Householder 變換的一半運算量。

 
QR 分解算法的精確度及穩定性分析相當複雜,有興趣深入了解的讀者可以參考文獻[2],以下僅說明主要結果。MGS 計算出的 Q_1,記為 \hat{Q}_1,滿足

\hat{Q}_1^T\hat{Q}_1=I+E_{MGS}

其中 \Vert E_{MGS}\Vert\approx u\cdot\kappa(A),這裡 u 代表單位捨去 (chopped) 誤差或捨入 (rounded) 誤差,\kappa(A)A 的條件數 (見“條件數”)。對於一般矩陣 A,包括非方陣,以2-範數定義的條件數為 \kappa(A)=\sigma_{\max}/\sigma_{min}\sigma_{\max}\sigma_{\min} 分別是 A 的最大和最小奇異值。另一方面,Householder 變換得到的 \hat{Q}_1 滿足

\hat{Q}_1^T\hat{Q}_1=I+E_{H}

其中 \Vert E_{H}\Vert\approx u。這指出 Householder 變換 (Givens 旋轉亦同) 是無條件穩定算法,因為 E_H 完全由浮點運算的精確度決定。MGS 則是有條件穩定算法,穩定與否取決於條件數 \kappa(A),只有當 \kappa(A) 不很大,也就是 A 的行向量「相當獨立」時,MGS 才適用於計算 Q_1

 
最後我們想知道為何 Householder 變換和 Givens 旋轉能夠擁有無條件數值穩定性?我提供一個簡單的解釋[3]。假設 A 是我們要處理的矩陣,E 代表因捨去誤差而產生的微小擾動矩陣,即 \Vert E\Vert_F 遠比 \Vert A\Vert_F 來得小,這裡 Frobenius 範數 \Vert A\Vert_F 的算式為 (見“矩陣範數”)

\Vert A\Vert_F=\sqrt{\mathrm{trace}(A^TA)}

P_i 為任意正交矩陣,考慮

P_k\cdots P_1(A+E)=P_k\cdots P_1A+P_k\cdots P_1E=PA+PE

其中正交矩陣乘積 P=P_k\cdots P_1 亦為正交矩陣,故

\Vert PE\Vert_F=\sqrt{\mathrm{trace}(E^TP^TPE)}=\sqrt{\mathrm{trace}(E^TE)}=\Vert E\Vert_F

因此證明執行一連串正交矩陣乘法的 Householder 變換和 Givens 旋轉並未放大浮點運算誤差。

 
本文參考:
[1] Gene H. Golub and Charles F. Van Loan, Matrix Computations, 2nd ed., The John Hopkins University Press, 1989, pp 211-219.
[2] Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, Society of Industrial and Applied Mathematics, 1996, pp 363-386.
[3] Carl Meyer, Matrix Analysis and Applied Linear Algebra, Society for Industrial and Applied Mathematics, 2000, pp 345-350.

Advertisements
本篇發表於 線性代數專欄, 數值線性代數 並標籤為 , , , , 。將永久鏈結加入書籤。

8 則回應給 QR 分解的數值計算方法比較

  1. Beck 說道:

    感謝這篇分析,目前我用GS就遇到穩定性的問題

  2. 李孟豪 說道:

    老師你好,我在你的教學網頁上常看到運算量,但不清楚是怎麼算出來的,我對這個很感興趣,想請問您運算量要怎麼去計算呢?

    • ccjou 說道:

      我不記得是怎麼得到這些數字,很可能抄錄自文獻[1]。

    • Meiyue Shao 說道:

      运算量可以直接对算法里的循环进行计数得到,前提是你知道算法应当如何实现。

      • 李孟豪 說道:

        感謝Meiyue Shao的回答!請問有沒有參考文獻或關鍵字可以讓我去學習這部分的相關內容呢?

        • Meiyue Shao 說道:

          可以去找本《数值线性代数》或者《矩阵计算》的书看看。当然,直接找算法的书看也行,但那里只有复杂度分析的介绍,不会专门讨论数值算法。

發表迴響

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

WordPress.com Logo

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

Twitter picture

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

Facebook照片

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

Google+ photo

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

連結到 %s