Gram-Schmidt 正交化改良版

本文的閱讀等級:中級

A 為一個 m\times n 階實矩陣,\mathrm{rank}A=n,QR 分解 A=QR 滿足以下兩個條件:Qm\times n 階矩陣,其行向量 (column vector) 組成單範正交 (orthonormal) 向量集,Rn\times n 階上三角矩陣。Gram-Schmidt 正交化是最常見於一般線性代數教科書的 QR 分解演算法 (見“Gram-Schmidt 正交化與 QR 分解”),以下稱之為古典 (classical) Gram-Schmidt 正交化,簡稱 CGS。不幸的是,當數值計算引入捨入 (roundoff) 誤差時,CGS 最終產生的 Q 矩陣的行向量正交性可能變的很糟,就此觀點而言,CGS 是數值不穩定的。本文介紹 CGS 的一個修改版本,稱為改良 (modified) Gram-Schmidt 正交化,簡稱 MGS[1,2]

 
我們從 QR 分解推導 Gram-Schmidt 正交化過程。將 AQ 以行向量表示為 A=\begin{bmatrix}  \mathbf{a}_1&\cdots&\mathbf{a}_n  \end{bmatrix}Q=\begin{bmatrix}  \mathbf{q}_1&\cdots&\mathbf{q}_n  \end{bmatrix}。因為 \mathbf{a}_j,\mathbf{q}_j\in\mathbb{R}^mR=[r_{ij}] 是上三角矩陣,乘開 A=QR,如下:

\displaystyle \begin{bmatrix} \mathbf{a}_1&\mathbf{a}_2&\cdots&\mathbf{a}_n \end{bmatrix}=\begin{bmatrix} \mathbf{q}_1&\mathbf{q}_2&\cdots&\mathbf{q}_n \end{bmatrix}\begin{bmatrix} r_{11}&r_{12}&\cdots&r_{1n}\\ 0&r_{22}&\cdots&r_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&r_{nn} \end{bmatrix}

比較等號兩邊即得

\mathbf{a}_j=r_{1j}\mathbf{q}_1+r_{2j}\mathbf{q}_2+\cdots+r_{jj}\mathbf{q}_j,~~j=1,\ldots,n

\mathrm{rank}A=n 可知 \mathrm{rank}R=n,故 R 不含零主對角元,就有

\mathbf{q}_j=\displaystyle\frac{1}{r_{jj}}\left(\mathbf{a}_j-\sum_{i=1}^{j-1}r_{ij}\mathbf{q}_i\right)

\mathbf{q}_j 當成與下列向量同方向的單位向量:

\mathbf{u}_j=\mathbf{a}_j-\displaystyle\sum_{i=1}^{j-1}r_{ij}\mathbf{q}_i

故得 r_{jj}=\Vert\mathbf{u}_j\Vert。若 \{\mathbf{q}_1,\ldots,\mathbf{q}_{j-1}\} 為一個單範正交集,為了保證 \mathbf{u}_j\in\mathrm{span}\{\mathbf{q}_1,\ldots,\mathbf{q}_{j-1}\}^{\perp},令

r_{ij}=\mathbf{q}_i^T\mathbf{a}_j,~~i=1,\ldots,j-1

對於 k=1,\ldots,j-1,就有

\begin{aligned} \mathbf{u}_j^T\mathbf{q}_k&=\displaystyle\left(\mathbf{a}_j-\sum_{i=1}^{j-1}r_{ij}\mathbf{q}_i \right)^T\mathbf{q}_k\\ &=\mathbf{a}_j^T\mathbf{q}_k-\sum_{i=1}^{j-1}r_{ij}\mathbf{q}_i^T\mathbf{q}_k\\ &=r_{kj}-r_{kj}=0\end{aligned}

將以上結果整理成演算法,如下:

古典 Gram-Schmidt 正交化


  1. 對於 j=1,\ldots,n
  2.     \mathbf{u}_j\leftarrow\mathbf{a}_j
  3.     對於 i=1,\ldots,j-1
  4.         r_{ij}\leftarrow\mathbf{q}_i^T\mathbf{a}_j
  5.         \mathbf{u}_j\leftarrow\mathbf{u}_j-r_{ij}\mathbf{q}_i
  6.         r_{jj}\leftarrow\Vert\mathbf{u}_j\Vert
  7.         \mathbf{q}_j\leftarrow\mathbf{u}_j/r_{jj}

 
接下來我用一個例子說明 CGS 的數值不穩定現象[3]。考慮下列 4\times 3 矩陣

A=\begin{bmatrix}  1&1&1\\  \epsilon&0&0\\  0&\epsilon&0\\  0&0&\epsilon  \end{bmatrix}

\epsilon 為極小正數,數值計算誤差使得 1+\epsilon^2\approx 1,此處 \approx 代表浮點運算的捨入過程,CGS 計算步驟如下:

j=1

\begin{aligned} \mathbf{u}_1&\leftarrow\mathbf{a}_1=\begin{bmatrix}  1\\  \epsilon\\  0\\  0  \end{bmatrix}\\  r_{11}&\leftarrow\Vert\mathbf{u}_1\Vert=\sqrt{1+\epsilon^2}\approx 1,~~\mathbf{q}_1\leftarrow\mathbf{u}_1/r_{11}=\begin{bmatrix}  1\\  \epsilon\\  0\\  0  \end{bmatrix}\end{aligned}

j=2

\begin{aligned} \mathbf{u}_2&\leftarrow\mathbf{a}_2=\begin{bmatrix}  1\\  0\\  \epsilon\\  0  \end{bmatrix}\\  r_{12}&\leftarrow\mathbf{q}_1^T\mathbf{a}_2=1,~~\mathbf{u}_2\leftarrow\mathbf{u}_2-r_{12}\mathbf{q}_1=\begin{bmatrix}  0\\  -\epsilon\\  \epsilon\\  0  \end{bmatrix}\\  r_{22}&\leftarrow\Vert\mathbf{u}_2\Vert=\sqrt{2}\epsilon,~~\mathbf{q}_2\leftarrow\mathbf{u}_2/r_{22}=\begin{bmatrix}  0\\  -1/\sqrt{2}\\  1/\sqrt{2}\\  0  \end{bmatrix}\end{aligned}

j=3

\begin{aligned} \mathbf{u}_3&\leftarrow\mathbf{a}_3=\begin{bmatrix}  1\\  0\\  0\\  \epsilon  \end{bmatrix}\\  r_{13}&\leftarrow\mathbf{q}_1^T\mathbf{a}_3=1,~~\mathbf{u}_3\leftarrow\mathbf{u}_3-r_{13}\mathbf{q}_1=\begin{bmatrix}  0\\  -\epsilon\\  0\\  \epsilon  \end{bmatrix}\\  r_{23}&\leftarrow\mathbf{q}_2^T\mathbf{a}_3=0,~~\mathbf{u}_3\leftarrow\mathbf{u}_3-r_{23}\mathbf{q}_2=\begin{bmatrix}  0\\  -\epsilon\\  0\\  \epsilon  \end{bmatrix}\\  r_{33}&\leftarrow\Vert\mathbf{u}_3\Vert=\sqrt{2}\epsilon,~~\mathbf{q}_3\leftarrow\mathbf{u}_3/r_{33}=\begin{bmatrix}  0\\  -1/\sqrt{2}\\  0\\  1/\sqrt{2}  \end{bmatrix}\end{aligned}

最後得到

Q=\begin{bmatrix}  1&0&0\\  \epsilon&-1/\sqrt{2}&-1/\sqrt{2}\\  0&1/\sqrt{2}&0\\  0&0&1/\sqrt{2}  \end{bmatrix},~ R=\begin{bmatrix}  1&1&1\\  0&\sqrt{2}\epsilon&0\\  0&0&\sqrt{2}\epsilon  \end{bmatrix}

明顯地,\mathbf{q}_2\mathbf{q}_3 偏離正交。

 
上例顯示 CGS 產生新正交基底 (即 \mathbf{u}_3) 的計算出了問題,合併

r_{ij}=\mathbf{q}_i^T\mathbf{a}_j,~~i=1,\ldots,j-1

\mathbf{u}_j=\mathbf{a}_j-\displaystyle\sum_{i=1}^{j-1}\mathbf{q}_ir_{ij}

可得

\mathbf{u}_j=\displaystyle\mathbf{a}_j-\sum_{i=1}^{j-1}\mathbf{q}_i\mathbf{q}_i^T\mathbf{a}_j=(I-\mathbf{q}_1\mathbf{q}_1^T-\cdots-\mathbf{q}_{j-1}\mathbf{q}_{j-1}^T)\mathbf{a}_j

這個算式之所以發生錯誤的原因在於 \mathbf{q}_1,\ldots,\mathbf{q}_{j-1} 並非完全正交,如上例 \mathbf{q}_1^T\mathbf{q}_2=-\epsilon/\sqrt{2},若 \mathbf{u}_j 的各元數值都很小,經正規化而得的 \mathbf{q}_j 很容易遠離 \mathrm{span}\{\mathbf{q}_1,\ldots,\mathbf{q}_{j-1}\}^{\perp}。為了克服此問題,我們應設法降低不全然正交的基底 \mathbf{q}_1,\ldots,\mathbf{q}_{j-1} 造成的破壞。

 
P_i=I-\mathbf{q}_{i}\mathbf{q}_i^Ti\ge 1P_0\equiv I,若 \{\mathbf{q}_i\} 為一個單範正交集 (實際上,此條件並不成立),對於 j\neq k,就有

P_jP_k=(I-\mathbf{q}_j\mathbf{q}_j^T)(I-\mathbf{q}_k\mathbf{q}_k^T)=I-\mathbf{q}_j\mathbf{q}_j^T-\mathbf{q}_k\mathbf{q}_k^T

所以

P_{j-1}\cdots P_1=I-\mathbf{q}_1\mathbf{q}_1^T-\cdots-\mathbf{q}_{j-1}\mathbf{q}_{j-1}^T

Gram-Schmidt 正交化的主要公式可表示為

\mathbf{u}_j=P_{j-1}\cdots P_1\mathbf{a}_j=(I-\mathbf{q}_{j-1}\mathbf{q}_{j-1}^T)\cdots(I-\mathbf{q}_1\mathbf{q}_1^T)\mathbf{a}_j

j=3 為例說明,

\begin{aligned} \mathbf{u}_3&=(I-\mathbf{q}_2\mathbf{q}_2^T)(I-\mathbf{q}_1\mathbf{q}_1^T)\mathbf{a}_3\\ &=(I+\mathbf{q}_2(\mathbf{q}_2^T\mathbf{q}_1)\mathbf{q}_1^T-\mathbf{q}_2\mathbf{q}_2^T-\mathbf{q}_1\mathbf{q}_1^T)\mathbf{a}_3\end{aligned}

上式的意思是先將 \mathbf{a}_3 投影至 \mathbf{q}_1 的分量扣除,再將殘餘量 (I-\mathbf{q}_1\mathbf{q}_1^T)\mathbf{a}_3 投影至 \mathbf{q}_2 的分量扣除,這個步驟可以「投影消滅」因 \mathbf{q}_2^T\mathbf{q}_1\neq 0 而產生的誤差。下面是根據 \mathbf{u}_j 新算式整理而成的改良法 MGS,比較 CGS 和 MGS,兩者唯一的差別在於 r_{ij} 的計算方式不同。

改良 Gram-Schmidt 正交化


  1. 對於 j=1,\ldots,n
  2.     \mathbf{u}_j\leftarrow\mathbf{a}_j
  3.     對於 i=1,\ldots,j-1
  4.         r_{ij}\leftarrow\mathbf{q}_i^T\mathbf{u}_j
  5.         \mathbf{u}_j\leftarrow\mathbf{u}_j-r_{ij}\mathbf{q}_i
  6.         r_{jj}\leftarrow\Vert\mathbf{u}_j\Vert
  7.         \mathbf{q}_j\leftarrow\mathbf{u}_j/r_{jj}

 
將 MGS 應用於計算上例 A 矩陣的 QR 分解。當 j=1j=2,MGS 和 CGS 有相同結果,以下只列出 j=3 的情況:

\begin{aligned} \mathbf{u}_3&\leftarrow\mathbf{a}_3=\begin{bmatrix}  1\\  0\\  0\\  \epsilon  \end{bmatrix}\\  r_{13}&\leftarrow\mathbf{q}_1^T\mathbf{u}_3=1,~~\mathbf{u}_3\leftarrow\mathbf{u}_3-r_{13}\mathbf{q}_1=\begin{bmatrix}  0\\  -\epsilon\\  0\\  \epsilon  \end{bmatrix}\\  r_{23}&\leftarrow\mathbf{q}_2^T\mathbf{u}_3=\displaystyle\frac{\epsilon}{\sqrt{2}},~~\mathbf{u}_3\leftarrow\mathbf{u}_3-r_{23}\mathbf{q}_2=\begin{bmatrix}  0\\  -\epsilon/{2}\\  -\epsilon/{2}\\  \epsilon  \end{bmatrix}\\  r_{33}&\leftarrow\Vert\mathbf{u}_3\Vert=\sqrt{3/2}\epsilon,~~\mathbf{q}_3\leftarrow\mathbf{u}_3/r_{33}=\begin{bmatrix}  0\\  -1/\sqrt{6}\\  -1/\sqrt{6}\\  \sqrt{2/3}  \end{bmatrix}\end{aligned}

這個 \mathbf{q}_3 向量和 CGS 的結果完全不同,MGS 算出的 \mathbf{q}_1, \mathbf{q}_2, \mathbf{q}_3 幾乎正交,QR 分解矩陣為

Q=\begin{bmatrix}  1&0&0\\  \epsilon&-1/\sqrt{2}&-1/\sqrt{6}\\  0&1/\sqrt{2}&-1/\sqrt{6}\\  0&0&\sqrt{2/3}  \end{bmatrix},~ R=\begin{bmatrix}  1&1&1\\  0&\sqrt{2}\epsilon&\epsilon/\sqrt{2}\\  0&0&\sqrt{3/2}\epsilon  \end{bmatrix}

 
MGS 雖然可以降低數值計算造成的誤差,但仍不能保證最終得到的 \{\mathbf{q}_i\} 具備良好的正交性,也就是說,不論 CGS 或 MGS 都不適用於 QR 分解。現實應用中,我們經常採用 Householder 變換或 Givens 旋轉計算 QR 分解 (見“Householder 變換於 QR 分解的應用”和“Givens 旋轉於 QR 分解的應用”),他日將另文解說各方法的數值穩定性和運算量。

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

延伸閱讀:
This entry was posted in 線性代數專欄, 數值線性代數 and tagged , , , . Bookmark the permalink.

10 則回應給 Gram-Schmidt 正交化改良版

  1. levinc 說:

    妙哉! 期待續集!

  2. ccjou 說:

    為了凸顯 A 不必是方陣,我將原先的例子更換為 4\times 3 階矩陣。

  3. Liang Dai 說:

    这种矩阵分解的方法从原始GS正交化到MGS的讲解方法让我感觉很新鲜,比单纯从算法流程的角度讲解深刻。谢谢周老师!

  4. ccjou 說:

    其實我寫的有些心虛,總覺得從CGS到MGS不夠慎密,後來翻閱一些材料也沒有具體收穫,暫時將就先放著,有更好的主意時再修改好了。

  5. vivian 說:

    請問 如果A的維度m*n 且m>n
    則Q的最後m-n個行向量要怎麼算呢??

  6. vivian 說:

    謝謝你的回覆 另外如果我想計算QR分解的複雜度 (需要做幾次乘法) (要包含Q最後的m-n個行向量)
    這查書查的到嗎??
    真是個頭疼的問題啊..

  7. Meiyue Shao 說:

    其实历史上Laplace解最小二乘问题时就已经用MGS了,后来Schmidt使用的反而被称为CGS。

發表迴響

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

WordPress.com Logo

你正使用 WordPress.com 帳號留言。 登出 / 變更 )

Twitter picture

你正使用 Twitter 帳號留言。 登出 / 變更 )

Facebook照片

你正使用 Facebook 帳號留言。 登出 / 變更 )

Google+ photo

你正使用 Google+ 帳號留言。 登出 / 變更 )

連結到 %s