7.2 最小平方擬合#

概念說明#

拉格朗日插值要求曲線恰好穿過每個資料點，但實測資料總有雜訊——硬要穿過，等於把雜訊也學進去了。最小平方擬合（Least Squares Fitting）改問一個更健康的問題：在選定的函數族（如直線、二次式）中，找一條使殘差平方和最小的曲線：

minimize E(c) = Σᵢ (yᵢ − Σⱼ cⱼφⱼ(xᵢ))²

其中 φⱼ 是基底函數（直線擬合時 φ₀ = 1, φ₁ = x），c 是待定係數。取平方而不是絕對值，是因為它可微、好解，且對應高斯雜訊假設下的最大概似估計（見第 9 章）。

把 E 對每個 cₖ 偏微分並令為零，得到正規方程（normal equations）：

(AᵀA)c = Aᵀy

其中 A 的第 i 列是 (φ₀(xᵢ), …, φₘ(xᵢ))。理論上解這個 m+1 階線性系統就好（5.1 節的技術全用上），但數值上 AᵀA 會把條件數平方（κ(AᵀA) = κ(A)²），雜訊大或基底病態時是災難。實務首選是 QR 分解：A = QR，正規方程變成 Rc = Qᵀy，條件數不加倍。最穩的是 SVD，連秩不足都能處理。

最小平方也是統計學「線性迴歸」的計算核心，與第 9 章的假設檢定直接相連。

演算法步驟#

選定模型（直線、多項式、或其他基底組合）。
建立設計矩陣 A（每列是 [φ₀(xᵢ), …, φₘ(xᵢ)]）。
解 min ‖Ac − y‖₂：
教學版：解正規方程 (AᵀA)c = Aᵀy。
穩健版：QR 分解 A = QR，解 Rc = Qᵀy。
檢驗殘差：畫殘差圖、算 R² = 1 − SSE/SST，判斷模型是否恰當。

Python 程式碼範例#

import numpy as np

rng = np.random.default_rng(1)

# --- 資料：帶雜訊的直線 ---
n = 50
x = np.linspace(0, 10, n)
y = 2.5 * x - 1.0 + rng.normal(0, 1.5, n)

# 方法一：正規方程（教學用）
A = np.column_stack([np.ones(n), x])
c_normal = np.linalg.solve(A.T @ A, A.T @ y)
print("正規方程解  [截距, 斜率] =", c_normal)

# 方法二：QR 分解（穩健）
Q, R = np.linalg.qr(A)
c_qr = np.linalg.solve(R, Q.T @ y)
print("QR 分解解   [截距, 斜率] =", c_qr)

# 方法三：內建 lstsq（SVD）
c_lstsq, res, rank, sv = np.linalg.lstsq(A, y, rcond=None)
print("lstsq 解    [截距, 斜率] =", c_lstsq)

# 殘差與 R²
yhat = A @ c_qr
sse = np.sum((y - yhat) ** 2)
sst = np.sum((y - y.mean()) ** 2)
print(f"R² = {1 - sse/sst:.4f}，殘差標準差 = {np.sqrt(sse/(n-2)):.3f}")

# --- 二次多項式擬合與條件數陷阱 ---
y2 = 0.5 * x**2 - 3 * x + 2 + rng.normal(0, 2, n)
A2_raw = np.column_stack([np.ones(n), x, x**2])
A2_fit = np.column_stack([np.ones(n), x, x**2]) # 未置換基底
c2 = np.linalg.lstsq(A2_fit, y2, rcond=None)[0]
print("二次擬合係數 =", c2)
print("κ(A)        =", f"{np.linalg.cond(A2_raw):.1f}")

# 多項式擬合用置換基底更穩（x 中心化）
xbar = x - x.mean()
A2_well = np.column_stack([np.ones(n), xbar, xbar**2])
print("κ(置換後)   =", f"{np.linalg.cond(A2_well):.1f}")
c2w = np.linalg.lstsq(A2_well, y2, rcond=None)[0]
print("置換基底係數 =", c2w, "（展開後應與 c2 一致）")

誤差分析#

最小平方有完整的幾何與統計理論。幾何上，擬合值 ŷ = Ac 是 y 在 A 的行空間上的正交投影，殘差 y − ŷ 與投影方向垂直——QR 分解正是這個幾何事實的代數化身。統計上，若雜訊 εᵢ 獨立同分布、方差 σ²，則係數估計的協方差矩陣為 σ²(AᵀA)⁻¹，標準誤 SE(cₖ) = σ√[(AᵀA)⁻¹]ₖₖ，這是第 9 章信賴區間與 t 檢定的基礎。R² 越接近 1 表示模型解釋力越強，但高 R² 不代表模型正確——殘差檢驗（有無系統性彎曲）不可省。數值層面：正規方程的相對誤差 ≈ ε·κ(A)²，QR 版 ≈ ε·κ(A)，SVD 版可達 ε·κ(A) 且在秩不足時給出最小範數解。模型選擇的經驗法則：基底數目 m 應遠小於資料筆數 n（建議 m ≤ n/10），否則開始「背答案」，這也是線性版的過擬合問題，正規化（ridge：加 λ‖c‖²）是標準藥方。

練習#

產生 y = sin(x) + 雜訊的 100 點資料，分別用 3、10、60 次多項式做最小平方擬合，畫出擬合曲線與殘差圖，指出哪一組開始過擬合。
實作 ridge 迴歸：解 (AᵀA + λI)c = Aᵀy，在上一題的 60 次擬合上加 λ = 0.1, 1, 10，觀察係數範數與測試誤差的變化。
用 QR 分解手推並驗證：對同一組資料，np.linalg.lstsq 與「正規方程 + Cholesky」的係數差異，並計算兩者的條件數，解釋差異來源。