Gram-Schmidt 正交化改良版

本文的閱讀等級：中級

設 $A$ 為一個 $m\times n$ 階實矩陣， $\mathrm{rank}A=n$ ，QR 分解 $A=QR$ 滿足以下兩個條件： $Q$ 是 $m\times n$ 階矩陣，其行向量 (column vector) 組成單範正交 (orthonormal) 向量集， $R$ 是 $n\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 正交化過程。將 $A$ 和 $Q$ 以行向量表示為 $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}^m$ ， $R=[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 正交化

對於 $j=1,\ldots,n$
$\mathbf{u}_j\leftarrow\mathbf{a}_j$
對於 $i=1,\ldots,j-1$
$r_{ij}\leftarrow\mathbf{q}_i^T\mathbf{a}_j$
$\mathbf{u}_j\leftarrow\mathbf{u}_j-r_{ij}\mathbf{q}_i$
$r_{jj}\leftarrow\Vert\mathbf{u}_j\Vert$
$\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^T$ ， $i\ge 1$ ， $P_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 正交化

對於 $j=1,\ldots,n$
$\mathbf{u}_j\leftarrow\mathbf{a}_j$
對於 $i=1,\ldots,j-1$
$r_{ij}\leftarrow\mathbf{q}_i^T\mathbf{u}_j$
$\mathbf{u}_j\leftarrow\mathbf{u}_j-r_{ij}\mathbf{q}_i$
$r_{jj}\leftarrow\Vert\mathbf{u}_j\Vert$
$\mathbf{q}_j\leftarrow\mathbf{u}_j/r_{jj}$

將 MGS 應用於計算上例 $A$ 矩陣的 QR 分解。當 $j=1$ 和 $j=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, 2^nd 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.

10 Responses to Gram-Schmidt 正交化改良版

levinc says:

06/01/2011 at 10:26 pm

妙哉! 期待續集!

ccjou says:

06/02/2011 at 12:29 pm

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

Liang Dai says:

06/02/2011 at 4:50 pm

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

ccjou says:

06/02/2011 at 7:23 pm

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

vivian says:

01/17/2014 at 12:01 am

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

- ccjou says:
  
  01/17/2014 at 7:29 am
  
  你可以使用Steinitz 替換原則將mxn階矩陣擴充為mxm階，然後再執行正交化。Steinitz 替換原則的細節請見下文：
  https://ccjou.wordpress.com/2010/06/15/%E5%9F%BA%E5%BA%95%E8%88%87%E7%B6%AD%E5%BA%A6-%E5%B8%B8%E8%A6%8B%E5%95%8F%E7%AD%94%E9%9B%86/
  
vivian says:

01/17/2014 at 8:43 pm

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

- ccjou says:
  
  01/17/2014 at 10:19 pm
  
  請見下文：
  https://ccjou.wordpress.com/2011/06/03/qr-%E5%88%86%E8%A7%A3%E7%9A%84%E6%95%B8%E5%80%BC%E8%A8%88%E7%AE%97%E6%96%B9%E6%B3%95%E6%AF%94%E8%BC%83/
  這些數據是從Gene H. Golub and Charles F. Van Loan, Matrix Computations找到的。
  
Meiyue Shao says:

01/13/2016 at 10:25 am

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

- ccjou says:
  
  01/13/2016 at 11:18 am
  
  謝謝分享，有興趣深入了解的讀者請參閱下文：
  
  Click to access 0907.4695v1.pdf

	陳倍恩 on 線性代數的第一堂課──矩陣乘法的定義
	輕鬆談如何教學二項式定理？… on 牛頓的二項式定理 (上)
	madhouse on 高斯消去法
	WishMobile on 翻轉 LU 分解
	周子傑 on Karush-Kuhn-Tucker (KKT) 條件
	Cloud Huang on 線性泛函與伴隨