矩陣的 QR 分解在電腦運算中可能造成誤差，本文探討一下一種改進版本的 Gram-Schmidt 正交化方法。

經典的 Gram-Schmidt 方法可能造成數值不穩定性。在電腦中，舍入誤差可能會累積，造成得到的正交基並不具有正交性。

經典的 Gram-Schmidt

我們先來回顧一下經典的 Gram-Schmidt 方法：

1
2
3
4
5
6

for j = 1 : n
    v_j = x_j
    for k = 1 : j - 1
        v_j = v_j - ( (v_k^T x_j) / (v_k^T v_k) ) * v_k
    endfor
endfor

最終我們將 $v_j$ 歸一化得到標準正交基 $q_j = (v_j)/(| v_j |)$ 。

注意：CGS 始終使用原始向量 $x_j$ 與之前的基向量計算投影

數學證明

簡單進行數學證明

定義第 $j$ 步生成的向量 $v_j$

$$v_j = x_j - \sum_{k=1}^{j-1} ( (v_k^T x_j) / (v_k^T v_k) ) v_k$$

選取任意一個之前的基向量 $v_m$（其中 $1 \le m < j$），計算它與 $v_j$ 的內積 $v_m^T v_j$：

$$v_m^T v_j = v_m^T ( x_j - \sum_{k=1}^{j-1} ((v_k^T x_j) / (v_k^T v_k)) v_k )$$

利用內積的線性性質，將 $v_m^T$ 乘進去：

$$v_m^T v_j = v_m^T x_j - \sum_{k=1}^{j-1} ((v_k^T x_j) / (v_k^T v_k)) (v_m^T v_k)$$

由於前提假設 $v_1, \dots, v_{j-1}$ 是兩兩正交的，所以在求和符號 $\sum$ 中：

• 當 $k \neq m$ 時，$v_m^T v_k = 0$ 。
• 當且僅當 $k = m$ 時，$v_m^T v_k = v_m^T v_m \neq 0$ 。

因此，求和項中只剩下 $k=m$ 這一項：

$$v_m^T v_j = v_m^T x_j - ((v_m^T x_j) / (v_m^T v_m)) (v_m^T v_m)$$

分子分母中的標量 $v_m^T v_m$ 相互抵消：

$$v_m^T v_j = v_m^T x_j - v_m^T x_j$$

$$v_m^T v_j = 0$$

證畢

誤差分析

如果在 CGS（經典格拉姆-施密特正交化）中計算 $q_2$ 時發生誤差，導致 $q_1^T q_2 = \delta$ 是一個很小但非零的數值，這個誤差將不會在隨後的任何計算中被修正：

$$v_3 = x_3 - (q_1^T x_3)q_1 - (q_2^T x_3)q_2$$

$$q_2^T v_3 = q_2^T x_3 - (q_1^T x_3)\delta - (q_2^T x_3) = -(q_1^T x_3)\delta$$

$$q_1^T v_3 = q_1^T x_3 - (q_1^T x_3) - (q_2^T x_3)\delta = -(q_2^T x_3)\delta$$

我們可以看出 $v_3$ 與 $q_1$ 或 $q_2$ 都不正交。

改進的 Gram-Schmidt (MGS)

1
2
3
4
5
6
7
8
9
10

for j = 1 : n
    v_j = x_j
endfor

for j = 1 : n
    q_j = v_j / ||v_j||_2
    for k = j + 1 : n
        v_k = v_k - (q_j^T v_k) q_j
    endfor
endfor

最終我們將 $v_j$ 歸一化得到標準正交基 $q_j = (v_j)/(| v_j |)$ 。

區別在於，一旦算出來一個向量 $q$，立即用它去更新後面所有的向量（減去向量在 $q$ 上的投影），使得後面所有的向量與這個向量 $q$ 正交。

誤差分析

假設 MGS 出現：$q_1^T q_2 = \delta$ 很小但非零。

對於第三個向量 $v_3$ ：

• 初始狀態：$v_3^{(0)} = x_3$
• $j = 1: v_3^{(1)} = v_3^{(0)} - (q_1^T v_3^{(0)})q_1$
• $j = 2: v_3 = v_3^{(1)} - (q_2^T v_3^{(1)})q_2$

此時 $v_3$ 是第三個向量在歸一化之前的最終形式。

讓我們檢查正交性，假設不再產生其他誤差：

$q_2^T v_3 = q_2^T v_3^{(1)} - (q_2^T v_3^{(1)}) = 0$

所以，我們保留了對 $q_2$ 的正交性。

接下來看 $q_1$：

$$q_1^T v_3 = q_1^T v_3^{(1)} - (q_2^T v_3^{(1)})\delta$$

$$q_2^T v_3^{(1)} = q_2^T v_3^{(0)} - (q_1^T v_3^{(0)})\delta$$

$$q_1^T v_3^{(1)} = q_1^T v_3^{(0)} - (q_1^T v_3^{(0)}) = 0$$

由此可得：

$$q_1^T v_3 = -(q_2^T v_3^{(0)} - (q_1^T v_3^{(0)})\delta)\delta= -q_2^T v_3^{(0)}\delta + q_1^T v_3^{(0)}\delta^2$$

由於 $\delta$ 非常小，$\delta^2$ 就更小了。因此 $q_1^T v_3$ 中的誤差和 CGS 差不多，但我們消除了 $q_2^T v_3$ 中的誤差，使得誤差更小。

總結

在電腦中，由於浮點運算的特殊性，有許多數學方法需要重新考慮，儘量提高數值穩定性。這使得我們有必要重新審視我們學習過的數學知識，針對電腦的特殊性進行重新設計與優化。