從投影到一個方陣系統
上一篇留給我們一幅乾淨的圖：在 最小平方法問題 裡，我們無法把高瘦的系統 A x = b 解得恰好成立，因為 b 幾乎從不落在 A 的行空間（column space）裡。我們能做到的最好程度，就是落在那個行空間中最接近的點上，而那個最接近的點，正是 b 到行空間的 正交投影。這個投影的關鍵特徵是：剩下的殘差 r = b - A x，會筆直地指出行空間之外——它垂直於 A 的每一行（每一個 column）。
現在把這句幾何的話翻成代數。「殘差垂直於 A 的每一行」意思是：每一行與 r 做內積都得零。把這些內積疊起來，恰好就是矩陣乘積 A^T r = 0。代入 r = b - A x，就得到 A^T (b - A x) = 0，重新整理成著名的 法方程：A^T A x = A^T b。我們從一個有 m 個方程式、n 個未知數、且無解的系統出發；最後落到一個有 n 個方程式、n 個未知數、而且有解的系統。
為什麼這道食譜看起來無法抗拒
在紙上，法方程簡直是一場美夢。矩陣 A^T A 是方陣、對稱，而且（當 A 的各行線性獨立時）正定——正好是最乾淨、最快的直接解法所鍾愛的那一類矩陣。你甚至不需要通用方法：一個 Cholesky 分解——它只用下三角把 A^T A 寫成 A^T A = L L^T——解這個系統所花的工作量，大約只是標準消去法的一半。對於參數個數 n 很小的情形，成本微不足道，程式碼三行就寫完。
- 組出小小的方陣 M = A^T A 與向量 c = A^T b——對 n 個參數而言，這只是 n×n 與長度 n 而已。
- 對它做 Cholesky 分解 M = L L^T，只用到對稱正定矩陣 M 的下三角部分。
- 先以前向代入解 L y = c，再以回代解 L^T x = y；得到的 x 即最小化 ||A x - b||_2。
所以法方程是每一門課最先教的東西，而就理解幾何而言，它們完美無瑕。麻煩純粹出在數值上，而且整個藏在那看似無辜的第一行裡：M = A^T A。組出這個乘積的動作——還沒等任何求解器跑起來——就已經悄悄地把問題搞壞了。要看清這是怎麼回事，我們得談談條件數。
隱藏的危險：把條件數平方
回想條件數那一階梯的主規則：準確度 = 條件數 × 穩定度。即便是一個無瑕、後向穩定 的求解器，也只能交付問題的條件數所留下來、未被吞掉的那些位數；一個接近 1e8 的條件數，就已經要花掉你約 16 位 雙精度 數字裡的大約 8 位。一個最小平方擬合的條件數，由 A 的條件數（記作 cond(A)）所主宰。而殘酷的事實是：你真正交給求解器的那個矩陣 A^T A，其條件數為 cond(A^T A) = cond(A)^2。
把那個指數慢慢讀，因為它就是整篇的全部。把一個條件數平方，會讓它吃掉的位數加倍。如果 A 輕度病態、cond(A) = 1e4——對一個真實的擬合來說，這是再普通不過的值——那麼你原本的問題損失約 4 位，這沒問題，你還留著 12 位。但 A^T A 的條件數是 1e8，所以走法方程這條路會損失約 8 位。你白白丟掉了 4 位好的準確度——不是因為資料、不是因為求解器，而純粹是因為「組出 A^T A」這個代數步驟。在 Cholesky 還沒碰到任何一個數之前，傷害就已經造成了。
一個你幾乎能手算的小例子
這裡有個經典的迷你例子，把危險變得具體——Läuchli 矩陣。讓 A 的各行幾乎平行：取 A = [[1, 1], [eps, 0], [0, eps]]，其中 eps 很小，比如 eps = 1e-8。各行線性獨立，所以原則上這個擬合是良置（well posed）的，而 cond(A) 大約是 1/eps = 1e8——大，但在雙精度（約帶 16 位數字）下還撐得住。
A = [ 1 1 ] eps = 1e-8
[ eps 0 ]
[ 0 eps]
A^T A = [ 1 + eps^2 1 ]
[ 1 1 + eps^2 ]
In double precision, 1 + eps^2 = 1 + 1e-16 -> ROUNDS to exactly 1.
so the stored matrix becomes [ 1 1 ] <- SINGULAR (rank 1!)
[ 1 1 ]
cond(A) ~ 1e8 (the honest problem: solvable)
cond(A^T A) ~ 1e16 (what you actually solve: numerically singular)看看食譜的第一行做了什麼。A^T A 的對角線是 1 + eps^2 = 1 + 1e-16，而在雙精度下，這個加法會直接捨回成 1——回想浮點數無法表示 1 + 1e-16，而且加法並不精確。非對角線則恰好是 1。於是你存下來的矩陣是 [[1,1],[1,1]]，它是奇異的：兩列一模一樣。Cholesky 會直接失敗、或回傳垃圾，然而底層的擬合本來是可解的。區分這兩行的資訊原本住在 eps^2 那些項裡，而組出 A^T A 把它捨進了虛無。
這就是危險的一幀全貌。你想解的問題條件數是 1e8，舒舒服服地落在雙精度的能力範圍內。你選來解的那個矩陣條件數是 1e16，正卡在全面崩潰的邊緣。資料本身一點問題都沒有；這份損失完全是演算法自己造成的。在一個不穩定的表述上跑一個穩定的求解器，你拿到的依然是壞答案——因為準確度等於條件數乘以穩定度，而你親手把條件數加倍了。
什麼時候沒問題，以及該改走哪條路
這一切並不代表法方程被禁用。當 A 良置時——cond(A) 大約落在 1 到 1e3 上下——把它平方仍會給你留下充足的位數，這時法方程快速、簡單，而且完全體面。當你有數百萬筆資料、卻只有寥寥幾個參數時，它的精簡是實打實的優點，因為 A^T A 很小，可以一列一列地累加，完全不必把整個 A 存下來。規則不是「永遠別用」，而是「在用之前，先弄清楚條件數」。
但這個危險絕非邊角案例，因為最常見的擬合任務之一就會一頭撞上它。多項式迴歸——用一條 d 次曲線去穿過資料——是用 1、x、x^2、x^3 等等這些行來組出它的矩陣 A。這些冪次行隨著次數升高，會愈長愈像，於是這個 范德蒙矩陣 會變得驚人地病態；即便只是等距點上的 10 次，就能把 cond(A) 推過 1e10。把它丟進法方程，你會把它平方成 1e20，遠遠超出任何雙精度所能挽救的範圍。下一篇就專門講這個解藥。
逃生路線是：根本就別去組 A^T A。穩定的方法不先把問題平方再去解，而是直接在 A 上動工：對 A 做 QR 分解，或在最棘手的情形下用 奇異值分解 與 偽逆。這些方法以 cond(A)、而非 cond(A)^2 的代價換來最小平方解——它們保住問題真正允許的全部位數。那正是接下來兩篇要鋪的路，所以請把這個危險記在心裡，當作它們存在的理由。