5.1 高斯消去與 LU 分解#

概念說明#

解線性聯立方程組 Ax = b 是所有數值計算中最基本的任務。高中教的是「消去法」，高斯消去就是它的系統化版本：用第 1 列去消掉下面各列的第 1 個未知數，再處理第 2 個未知數，直到矩陣變成上三角，最後用反向代換從最後一列一路解回第一列。

LU 分解則是高斯消去的「產品化」表述：把消去過程記錄下來，寫成 A = LU，其中 L 是單位下三角（記錄了每一步乘了什麼倍數），U 是上三角（消去後的結果）。這樣一來：

解 Ax = b 變成解 Ly = b（前向代換）再解 Ux = y（反向代換），各花 O(n²)。
若要解多組不同 b，分解只需做一次 O(n³)，之後每組 b 只要 O(n²)——這是 LU 相對「每次都重新消去」的最大賣點。

但純消去有個陷阱：若主對角線元素（主軸，pivot）接近零，除以它會放大誤差甚至除以零。部分軸選（partial pivoting）在每一步選該行絕對值最大的元素當主軸並交換列，把穩定性救回來。含軸選的分解寫成 PA = LU（P 是排列矩陣）。

演算法步驟#

對 k = 0, …, n−1：
在第 k 行、第 k 列以下找絕對值最大的元素，交換該列至第 k 列（記錄到 P）。
對 i > k：計算倍數 mᵢₖ = aᵢₖ/aₖₖ，存入 L；第 i 列減去 mᵢₖ 倍第 k 列。
前向代換解 Ly = Pb。
反向代換解 Ux = y。

Python 程式碼範例#

import numpy as np

def lu_decompose(A):
    """含部分軸選的 LU 分解，回傳 P, L, U 使 PA = LU"""
    A = A.astype(float).copy()
    n = len(A)
    L = np.eye(n)
    P = np.eye(n)
    for k in range(n):
        # 部分軸選
        p = k + np.argmax(np.abs(A[k:, k]))
        if p != k:
            A[[k, p]] = A[[p, k]]
            P[[k, p]] = P[[p, k]]
            L[[k, p], :k] = L[[p, k], :k]
        if abs(A[k, k]) < 1e-14:
            raise ValueError("矩陣奇異")
        m = A[k+1:, k] / A[k, k]
        L[k+1:, k] = m
        A[k+1:, k:] -= np.outer(m, A[k, k:])
    return P, L, A

def lu_solve(A, b):
    P, L, U = lu_decompose(A)
    b = P @ b
    y = np.zeros_like(b, dtype=float)
    for i in range(len(b)):                    # 前向代換
        y[i] = b[i] - L[i, :i] @ y[:i]
    x = np.zeros_like(b, dtype=float)
    for i in range(len(b)-1, -1, -1):          # 反向代換
        x[i] = (y[i] - U[i, i+1:] @ x[i+1:]) / U[i, i]
    return x

A = np.array([[2., 1., 1.], [4., -6., 0.], [-2., 7., 2.]])
b = np.array([5., -2., 9.])
x = lu_solve(A, b)
print("解 x =", x)
print("檢驗 Ax - b =", A @ x - b)       # 應接近 0
print("與 np.linalg.solve 比對:", np.linalg.solve(A, b))

複雜度與誤差分析#

分解階段的三層迴圈共做約 n³/3 次乘加，時間 O(n³)，記憶體可原地覆蓋（L 與 U 共用一個陣列）。兩次代換各 O(n²)。因此 n = 1000 時約 3×10⁸ 次浮點運算，現代電腦一秒內可完成；這也是為什麼 BLAS 庫會拼命最佳化這個核心。誤差方面，誤差界由條件數 κ(A) = ‖A‖·‖A⁻¹‖ 控制：相對誤差約為 ε·κ(A)（ε 為機器精度）。部分軸選保證增長因子有界（一般情況 ≤ 2ⁿ⁻¹，實務上遠小），使 LU 成為稠密線性系統的預設方法。對稱正定矩陣可改用不需軸選的 Cholesky 分解（A = LLᵀ），成本再減半；稀疏大矩陣則改用迭代法（共軛梯度、GMRES），詳見第 10 章的延伸討論。

練習#

把 lu_decompose 改成原地運算：把 L 的下三角部分存進 U 陣列被清零的位置，最後只回傳一個矩陣與排列陣列。
對 Hilbert 矩陣 H（Hᵢⱼ = 1/(i+j−1)，n = 10）解 Hx = 1，比較計算解與精確解，計算 κ(H) = np.linalg.cond(H)，體會條件數 10¹³ 的災難。
同一矩陣 A、100 組隨機 b：比較「100 次 np.linalg.solve」與「一次 LU 分解 + 100 次代換」的總時間，驗證 O(n³) 對 O(n²) 的差異。