9.2 條件數與數值穩定性#

本節目標#

建立「條件數」的直覺：輸入誤差被放大的倍數
學會用 np.linalg.cond 量測矩陣的病態程度
用 Hilbert 矩陣親眼看到條件數隨維度爆炸、解的誤差跟著失控
學會四個實務對策：正規化、用 lstsq、先檢查 \(\kappa\) 再相信答案

問題：同一個方程組，答案為什麼這麼脆弱？#

考慮方程組 \(A\mathbf{x} = \mathbf{b}\) 。理論上 \(\mathbf{x} = A^{-1}\mathbf{b}\) 。但實務上有兩個誤差來源：

\(\mathbf{b}\) 本身來自測量或前序計算，天生帶誤差（上節學過浮點誤差無所不在）
求解過程自身又有捨入誤差

關鍵問題：輸入的小誤差，會在輸出被放大成多大的誤差？ 答案由矩陣 \(A\) 本身決定，這個放大倍數就是條件數：

\[ \kappa(A) = \lVert A \rVert \, \lVert A^{-1} \rVert \]

直覺解讀：

\[ \frac{\lVert \delta \mathbf{x} \rVert}{\lVert \mathbf{x} \rVert} \;\lesssim\; \kappa(A) \cdot \frac{\lVert \delta \mathbf{b} \rVert}{\lVert \mathbf{b} \rVert} \]

也就是說：** \(\mathbf{b}\) 有 \(10^{-16}\) 的相對誤差， \(\mathbf{x}\) 的相對誤差最壞可達 \(\kappa(A) \times 10^{-16}\) **。 \(\kappa = 10^3\) 無妨， \(\kappa = 10^{16}\) 答案完全不可信。

矩陣自身還有一個事實： \(\kappa(A) \geq 1\) 恆成立，且 \(A\) 奇異（不可逆）時 \(\kappa(A) = \infty\) 。條件數因此可以看作「距離奇異有多遠」的量尺。

用 NumPy 量測條件數#

import numpy as np

rng = np.random.default_rng(0)

I = np.eye(3)
print(np.linalg.cond(I))          # 1.0 —— 完美狀態

D = np.diag([1.0, 1.0, 1e-8])
print(np.linalg.cond(D))          # 1e+08 —— 對角線差 8 個數量級

singular = np.array([[1.0, 2.0],
                     [2.0, 4.0]])  # 第二列 = 2 倍第一列
print(np.linalg.cond(singular))   # inf（伴隨 warning）

A = rng.standard_normal((10, 10))
print(np.linalg.cond(A))          # 隨機矩陣通常在 10~100 之間，算健康

np.linalg.cond 預設用 2-範數，內部透過 SVD 計算： \(\kappa_2(A) = \sigma_{\max} / \sigma_{\min}\) （最大奇異值除以最小奇異值）。最小奇異值越接近 0，矩陣越接近奇異，條件數越大。

病態矩陣明星：Hilbert 矩陣#

Hilbert 矩陣定義為 \(H_{ij} = 1/(i+j-1)\) ，元素全是「無害」的簡單分數，卻是出了名的病態。我們來看 \(\kappa\) 隨維度爆炸，以及解的實際誤差：

import numpy as np
from scipy.linalg import hilbert   # 若沒有 scipy，下面也附了手刻版

def hilbert_manual(n):
    i, j = np.indices((n, n))
    return 1.0 / (i + j + 1)       # H[i,j] = 1/(i+j+1)

for n in [3, 5, 8, 10, 12]:
    H = hilbert(n)
    kappa = np.linalg.cond(H)
    print(f"n={n:2d}  kappa={kappa:.2e}")

典型輸出：

n= 3  kappa=5.24e+02
n= 5  kappa=4.77e+05
n= 8  kappa=1.53e+10
n=10  kappa=1.60e+13
n=12  kappa=1.75e+16

\(12 \times 12\) 的小矩陣，條件數已經達到 float64 的極限（約 \(10^{16}\) ）。接著看這對「解方程」的實際殺傷力：

x_true = np.ones(10)              # 我們心裡先想好答案
H = hilbert(10)
b = H @ x_true                    # 製造出對應的 b（過程帶入一點捨入誤差）

x = np.linalg.solve(H, b)
print(np.linalg.norm(x - x_true) / np.linalg.norm(x_true))
# 典型輸出：1e-3 左右 —— 明明 b 幾乎完美，答案卻錯了 3 位有效數字！

# 理論預測：誤差上限約 kappa * eps = 1.6e13 * 2.2e-16 ≈ 3.5e-3
print(np.linalg.cond(H) * np.finfo(float).eps)

觀察重點：誤差實測值和 \(\kappa \cdot \varepsilon_{\text{mach}}\) 的理論預測同一量級。這不是巧合，是數值線性代數最核心的洞見： 答案可不可信，決定權在矩陣的條件數，不在電腦貴不貴。

實驗：給 b 一點擾動，看解差多少#

直接驗證「輸入誤差被放大 \(\kappa\) 倍」：

import numpy as np

rng = np.random.default_rng(1)

for n in [4, 8, 12]:
    H = hilbert(n)
    x_true = rng.standard_normal(n)
    b = H @ x_true

    delta = 1e-10 * rng.standard_normal(n)   # 極小的 b 擾動
    x1 = np.linalg.solve(H, b)
    x2 = np.linalg.solve(H, b + delta)

    rel_b = np.linalg.norm(delta) / np.linalg.norm(b)
    rel_x = np.linalg.norm(x2 - x1) / np.linalg.norm(x1)
    kappa = np.linalg.cond(H)
    print(f"n={n:2d} kappa={kappa:.1e}  "
          f"b 相對誤差={rel_b:.1e}  x 相對誤差={rel_x:.1e}  "
          f"放大倍數={rel_x/rel_b:.1e}")

結果會顯示：n=4 時放大倍數數百倍，n=12 時可以放大到 \(10^{8}\) 倍以上。同一行程式碼、同一個擾動大小，矩陣一病態，一切走樣。

順帶一提：不只需要解 \(\mathbf{x}\) 的時候要小心。 算逆矩陣本身就是病態操作—— \(\kappa(A)\) 大時， \(A^{-1}\) 的元素可能大到荒謬，這就是下一節實務準則「別用 inv」的伏筆。

奇異矩陣與 κ = ∞#

A = np.array([[1.0, 2.0],
              [2.0, 4.0]])
print(np.linalg.det(A))       # 0.0（或極小的浮點殘值）
print(np.linalg.cond(A))      # inf

# NumPy 不會沉默地給爛答案，它會警告你：
try:
    np.linalg.solve(A, np.array([3.0, 6.0]))
except np.linalg.LinAlgError as e:
    print("捕捉到：", e)      # Singular matrix

# 但「接近奇異」比「完全奇異」更危險 —— 它不報錯，只是給你錯的答案
A_close = A + 1e-15 * np.eye(2)
print(np.linalg.cond(A_close))    # ~1e+15，超大但不無限，solve 不會報錯！
print(np.linalg.solve(A_close, np.array([3.0, 6.0])))

完全奇異會被偵測到；接近奇異不會。這正是必須主動檢查條件數的原因。

實務對策#

對策一：先正規化#

條件數是「尺度」的敏感指標。把身高用公尺還是公里表示， \(\kappa\) 可以差好幾個數量級：

# 特徵尺度差異巨大的設計矩陣：先標準化各欄
rng = np.random.default_rng(2)
X = np.column_stack([rng.standard_normal(100) * 1000,      # 大尺度特徵
                     rng.standard_normal(100) * 0.001])    # 小尺度特徵
print(np.linalg.cond(X))          # 很大

Xs = (X - X.mean(0)) / X.std(0)
print(np.linalg.cond(Xs))         # 顯著下降

機器學習裡「特徵標準化」不只是收斂技巧，本質上就是在壓低條件數。

對策二：用 lstsq，不要用 inv#

x = np.linalg.lstsq(H, b, rcond=None)[0]   # 用 SVD 求解，天生數值穩定
# np.linalg.inv(H) @ b 則是兩次病態操作的疊加，誤差更大

lstsq 基於 SVD，在接近奇異時的行為也優雅（最小奇異值方向自然被壓制），還能順便解「行滿秩不足」的最小平方問題。

對策三：檢查 κ，再相信答案#

把條件數檢查變成習慣：

def solve_report(A, b):
    kappa = np.linalg.cond(A)
    print(f"kappa = {kappa:.2e}, 可信有效位數 ≈ {max(0, 16 - np.log10(kappa)):.1f}")
    return np.linalg.solve(A, b)

一個實用的經驗法則：**可信的位數 ≈ 16 − \(\log_{10} \kappa\) **。 \(\kappa = 10^8\) 的系統，答案大約只有 8 位可信。

對策四：改問題，不硬解#

\(\kappa\) 實在太大時，換方法：正則化（Ridge 的 \(\lambda I\) 就是抬高 \(\sigma_{\min}\) 、壓低 \(\kappa\) ）、降維、或重新設計模型。硬解病態系統是在跟數學定律對賭。

小結與練習#

條件數 \(\kappa(A) = \lVert A\rVert\lVert A^{-1}\rVert\) ：輸入相對誤差被放大的最壞倍數
\(\kappa \geq 1\) ；奇異時 \(\kappa = \infty\) ；「接近奇異」不報錯但答案不可信，最危險
Hilbert 矩陣： \(n=12\) 時 \(\kappa \sim 10^{16}\) ，解的有效數字全毀
誤差預測：輸出相對誤差 \(\lesssim \kappa \times\) 輸入相對誤差，實測吻合
對策：正規化、用 lstsq 而非 inv、先查 \(\kappa\) 再相信答案、必要時改問題

練習：

對 \(n = 2\) 到 \(13\) 的 Hilbert 矩陣畫出 \(\log_{10} \kappa\) 對 \(n\) 的曲線，驗證它幾乎是直線（成長速度約每年 1.1 個數量級）。
構造一個 \(\kappa = 10^{12}\) 的 \(2 \times 2\) 矩陣（提示：對角矩陣），用實驗量測 \(\mathbf{b}\) 擾動 \(10^{-8}\) 時解的相對誤差，與理論放大倍數比較。
用 np.linalg.solve 和 np.linalg.inv(H) @ b 分別解 \(n=10\) 的 Hilbert 系統，比較兩者與真解的誤差，驗證「別用 inv」。
生成隨機 \(100 \times 100\) 矩陣 1000 個，統計條件數的分佈（中位數、95 百分位），建立「隨機矩陣通常多健康」的直覺。