跳轉到

6.2 QR 分解#

本節目標#

  • 理解 \(A = QR\) 的意義： \(Q\) 是 orthonormal 基、 \(R\) 是上三角
  • 知道為什麼 QR 在數值上特別穩定（Householder 概念一句話版）
  • 用 np.linalg.qr 實作，並用 QR 解方程組 \(R\mathbf{x} = Q^T\mathbf{b}\)
  • 認識 QR 迭代求特徵值的概念（Power iteration 預告）
  • 實驗：驗證隨機矩陣的 \(Q^T Q = I\)

\(A = QR\) 的意義：換一組好基#

延續上一節：Gram-Schmidt 把矩陣 \(A\) 的行向量扳正成 orthonormal 基 \(Q\) 。 但扳正的過程「剪掉了東西」， \(R\) 就是把剪掉的量記帳下來的上三角矩陣：

\[ A = QR, \qquad A = \begin{bmatrix} \mathbf{a}_1 & \mathbf{a}_2 & \cdots \end{bmatrix}, \quad Q = \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 & \cdots \end{bmatrix} \]

幾何上這條式子說的是：

  • \(Q\) 的行 \(\mathbf{q}_1, \mathbf{q}_2, \dots\) 是 \(A\) 的行空間的一組好基（單位長、互相垂直）
  • \(R\) 的第 \(j\) 行記錄：原來的行 \(\mathbf{a}_j\) 在這組好基下的座標，其中對角線元素 \(r_{jj}\) 是 \(\mathbf{a}_j\) 在新基上的長度貢獻，非對角線是投影係數

舉 \(2 \times 2\) 看： \(\mathbf{a}_1 = r_{11}\mathbf{q}_1\) ， \(\mathbf{a}_2 = r_{12}\mathbf{q}_1 + r_{22}\mathbf{q}_2\) 。 用 \(\mathbf{q}_1\) 去乘第一條式，因為 \(\mathbf{q}_1 \cdot \mathbf{q}_2 = 0\) ，立刻得到 \(r_{11} = \mathbf{a}_1 \cdot \mathbf{q}_1\) —— 正交性讓係數一乘內積就現形，這正是上三角的由來。

import numpy as np

A = np.array([[1., 1., 0.],
              [1., 0., 1.],
              [0., 1., 1.]])

Q, R = np.linalg.qr(A)
print(Q)          # 行是 orthonormal 基
print(R)          # 上三角
print(Q @ R)      # 等於 A

為什麼 QR 數值上特別穩定？#

一個（幾乎只要一句話的）概念：

np.linalg.qr 實際上不是用課本的 Gram-Schmidt，而是用 Householder 反射—— 一系列「把某個向量鏡射到座標軸上」的正交變換，逐步把 \(A\) 的下方元素消成零。

為什麼這更穩？Gram-Schmidt 是「相減」，數字相近時相減會把有效位數吃掉（cancellation）。 Householder 則全程只做旋轉與鏡射這類保長度、保夾角的變換，誤差不易放大。 細節在數值分析課程再學，此刻你只需記住結論：

需要穩定的正交化時，一律用 np.linalg.qr，不要自己手刻 Gram-Schmidt。

用 QR 解方程組：先換基、再回代#

解 \(A\mathbf{x} = \mathbf{b}\) ，代入 \(A = QR\) ：

\[ QR\mathbf{x} = \mathbf{b} \;\Longrightarrow\; R\mathbf{x} = Q^T\mathbf{b} \]

關鍵一步是左乘 \(Q^T\) 。因為 \(Q\) 是正交的， \(Q^{-1} = Q^T\) ——反矩陣不用算！ 於是流程變成兩段：

  1. 換基： \(\mathbf{y} = Q^T \mathbf{b}\) （就是 \(\mathbf{b}\) 在 orthonormal 基下的座標，一連串內積）
  2. 回代：解上三角系統 \(R\mathbf{x} = \mathbf{y}\) （從最後一列往上逐一解出）

上三角系統超好解：最後一列只有一個未知數，倒數第二列兩個，一路往上收。

import numpy as np

A = np.array([[2., 1.],
              [1., 3.]])
b = np.array([5., 10.])

Q, R = np.linalg.qr(A)

y = Q.T @ b                # 第一步：換基
x = np.linalg.solve(R, y)  # 第二步：解上三角（R 是上三角，solve 會很快）
print(x)

print(np.linalg.solve(A, b))   # 直接解，答案一致

好處：如果同一個 \(A\) 要配很多個不同的 \(\mathbf{b}\) （例如每天更新的資料）， QR 只做一次，之後每個 \(\mathbf{b}\) 都只需便宜的第二步。

QR 迭代求特徵值：Power iteration 預告#

QR 的另一個大招是求特徵值。概念如下：

  1. 對 \(A\) 做 QR 分解： \(A = Q_1 R_1\)
  2. 反過來乘： \(A_1 = R_1 Q_1\)
  3. 重複： \(A_1 = Q_2 R_2\) ， \(A_2 = R_2 Q_2\) ，……

每一步 \(R Q = Q^{-1} A Q = Q^T A Q\) ，是正交相似變換——特徵值完全不變（第 8 章會證）。 奇妙的是，反覆迭代之後 \(A_k\) 會越來越接近上三角，而上三角矩陣的特徵值就明明白白躺在對角線上！

import numpy as np

A = np.array([[4., 1.],
              [2., 3.]])

Ak = A.copy()
for _ in range(50):
    Q, R = np.linalg.qr(Ak)
    Ak = R @ Q

print(Ak)                  # 越來越接近上三角
print(np.diag(Ak))         # 對角線 ≈ 特徵值
print(np.linalg.eigvals(A))   # NumPy 驗證

這就是為什麼 NumPy 的 eig 底層其實是 QR 迭代。 它還有一個更樸素的親戚 Power iteration（反覆乘同一個矩陣，抓住最大特徵值方向）， Google PageRank 的原始精神正是它——第 8 章我們會親手實作。

實驗：隨機矩陣的 \(Q^T Q = I\)#

import numpy as np
from linalg_viz import plot_vectors

np.random.seed(42)

# 100x100 的隨機矩陣，看看正交性有多乾淨
A = np.random.randn(100, 100)
Q, R = np.linalg.qr(A)

err = np.abs(Q.T @ Q - np.eye(100)).max()
print(f"max |Q^T Q - I| = {err:.2e}")   # 通常 < 1e-14

# 上三角檢查：R 的左下角應為零
print(f"max |下三角部分| = {np.abs(np.tril(R, -1)).max():.2e}")

# 小矩陣畫圖看：Q 的兩行互相垂直且等長
A2 = np.random.randn(2, 2)
Q2, _ = np.linalg.qr(A2)
plot_vectors([A2[:, 0], A2[:, 1]], labels=['a1', 'a2'])
plot_vectors([Q2[:, 0], Q2[:, 1]], labels=['q1', 'q2'])

畫圖會看到： \(\mathbf{a}_1, \mathbf{a}_2\) 大概率歪歪斜斜、長短不一， 而 \(\mathbf{q}_1, \mathbf{q}_2\) 永遠是兩支等長的直角箭頭——這就是 QR 送你的禮物。

小結與練習#

本節重點 1. \(A = QR\) ： \(Q\) 把 \(A\) 的行空間換成 orthonormal 基， \(R\) 記下座標（上三角）。 2. 數值穩定來自 Householder 反射：全程只用保長度的正交變換。 3. 解 \(A\mathbf{x} = \mathbf{b}\) 兩步走： \(\mathbf{y} = Q^T\mathbf{b}\) ，再回代 \(R\mathbf{x} = \mathbf{y}\) 。 4. QR 迭代 \(A \to RQ \to \cdots\) 收斂到上三角，對角線即特徵值——eig 的底層引擎。

練習 1. 用 \(A = \begin{bmatrix} 2 & 1 \\ 1 & 3 \end{bmatrix}\) 手驗 \(Q^T Q = I\) 與 \(QR = A\) 。 2. 寫一個 qr_solve(A, b) 函式封裝「換基 + 回代」，與 np.linalg.solve 比較答案與速度。 3. 把 QR 迭代跑在 \(A = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}\) 上，觀察收斂失敗（特徵值 \(\pm 1\) 絕對值相同）——想想這呼應 Power iteration 的什麼限制。