M2221-L02 — 線性方程組、殘差、條件性與最小二乘

預備知識：消元法、矩陣範數與正交性。MATLAB 示例使用標準數值運算。

1. 求解不等同於顯式構造逆矩陣

對非奇異方陣 $A$，數學上的解為 $x=A^{-1}b$。數值計算通常通過分解後的方程組求解，而非顯式求出 $A^{-1}$ 的每一欄。MATLAB 的 A\b 根據矩陣結構選擇求解器；其文檔提醒注意尺度不良、接近奇異及秩虧的輸入。求解器執行成功或沒有警告，並不是精度保證。參見官方 mldivide 文檔。

一個標準實驗是

A = [3 1; 1 2]; b = [5;5];
x = A\b;
r = b-A*x;
relres = norm(r,inf)/norm(b,inf);
assert(norm(x-[1;2],inf) < 1e-12);

精確解為 $(1,2)^{\mathsf T}$。測試同時核對已知答案及殘差；對解未知的問題，殘差小本身並不保證前向誤差小。

2. 消元、主元選擇與可檢驗的分解

消元法從下方各行減去主元行的倍數。若主元為零，就與可用的非零行交換；數值部分主元法從當前欄的可選行中取絕對值最大的元素。這使當前每個消元乘子的絕對值不超過一，但不能保證對所有矩陣，每個中間元素都小。

通常將結果寫為

$$PA=LU,$$

其中 $P$ 置換行，$L$ 是對角元為一的下三角矩陣，$U$ 是上三角矩陣。先解 $Ly=Pb$，再解 $Ux=y$。省略 $P$ 會改變方程。在 MATLAB 中使用明確的三個輸出，並驗證：

[L,U,P] = lu(A);
assert(norm(P*A-L*U,inf) < 1e-12);
y = L\(P*b);
xlu = U\y;
assert(norm(A*xlu-b,inf) < 1e-12);

這裡 $A$ 必須如前例那樣定義，或在同一完整腳本中定義。一般輸入的容差應隨尺度變化；固定容差僅適用於這個小型演示。

例題：為甚麼換行重要。 對 $A=\begin{pmatrix}0&1\\2&3\end{pmatrix}$，交換行得到 $P=\begin{pmatrix}0&1\\1&0\end{pmatrix}$ 及 $PA=\begin{pmatrix}2&3\\0&1\end{pmatrix}$；後者已是上三角矩陣，所以 $L=I$。取 $b=(1,5)^{\mathsf T}$，解 $Ux=Pb=(5,1)^{\mathsf T}$ 得 $x=(1,1)^{\mathsf T}$。若用未置換的 $b$ 作右側，就解了另一個方程組。

主元選擇是算法保障，不能治癒數學問題本身的病態。矩陣在精確算術中可能非奇異，卻在可用精度下數值上接近秩虧。行列式大小依賴尺度，單獨使用並非一般的條件性診斷：$10^{-100}I$ 的行列式很小，但在任何誘導範數下條件數都是一。

3. 不應混淆的三種誤差

對近似解 $\widehat x$，殘差為 $r=b-A\widehat x$；前向誤差為 $\widehat x-x$，其中 $x$ 是真解。後向誤差考察要改變多少輸入數據才能使 $\widehat x$ 成為精確解。固定 $A$ 時，它恰好解 $A\widehat x=b-r$，所以當 $b\ne0$ 時，$\|r\|/\|b\|$ 是其相對右側後向誤差。

對誘導矩陣範數定義

$$\kappa(A)=\|A\|\|A^{-1}\|.$$

由 $x-\widehat x=A^{-1}r$ 和 $\|b\|\le\|A\|\|x\|$，當 $b\ne0$ 時，

$$\frac{\|x-\widehat x\|}{\|x\|}\le\kappa(A)\frac{\|r\|}{\|b\|}.$$

這是一個界，不是等式，也不自動保證估計緊。若 $b=0$，歸一化表達式須改用絕對誤差。常用的尺度敏感殘差指標是在分母非零時取 $\|r\|/(\|A\|\|\widehat x\|+\|b\|)$；它的精確後向誤差解釋取決於允許的擾動模型。

例題：敏感性實驗。 設

$$A_\delta=\begin{pmatrix}1&1\\1&1+\delta\end{pmatrix},\quad b=\binom2{2+\delta},\quad \delta>0.$$

精確解為 $(1,1)^{\mathsf T}$。將 $b$ 擾動 $(0,\varepsilon)^{\mathsf T}$，兩式相減得到解的變化

$$\Delta x=\binom{-\varepsilon/\delta}{\varepsilon/\delta}.$$

當 $\delta$ 很小時，微小輸入擾動可能引起較大解變化。無窮範數下，

$$\|A_\delta\|_\infty=2+\delta,\quad \|A_\delta^{-1}\|_\infty=(2+\delta)/\delta,\quad \kappa_\infty(A_\delta)=(2+\delta)^2/\delta.$$

這種敏感性在精確算術中已經存在；沒有額外建模信息時，準確的求解器也不能消除它。

4. 最小二乘：投影而非精確相容

對 $A\in\mathbb R^{m\times n}$ 或 $\mathbb C^{m\times n}$，其中 $m\ge n$，最小化 $\|Ax-b\|_2$。殘差必須與欄空間正交，得到正規方程

$$A^*(Ax-b)=0.$$

不借助非正式求導也可證明必要性：若一個最小殘差在欄空間上的投影非零，減去該投影將嚴格減小平方範數。反之，正交性給出

$$\|A(x+h)-b\|_2^2=\|Ax-b\|_2^2+\|Ah\|_2^2,$$

所以正規方程的每個解都是最小解。若各欄線性無關，只有 $h=0$ 才有 $Ah=0$，因此最小解唯一。

若 $A=QR$ 是薄 QR 分解，且 $Q^*Q=I$、$R$ 是非奇異方形上三角矩陣，則求解

$$Rx=Q^*b.$$

這樣避免顯式構造 $A^*A$；滿欄秩時，後者的 2-範數條件數為 $\kappa_2(A)^2$，因為其特徵值是奇異值的平方。方法代數等價不代表數值穩定性相同。秩虧的最小二乘需要考慮秩的方法，例如 SVD 或帶主元的 QR，並說明在可能存在的多個最小解中選擇哪一個。

例題：擬合直線。 用 $y=\alpha+\beta t$ 擬合 $(t,y)=(0,1),(1,2),(2,2)$。此時

$$A=\begin{pmatrix}1&0\\1&1\\1&2\end{pmatrix},\quad b=\begin{pmatrix}1\\2\\2\end{pmatrix},$$

正規方程為 $\begin{pmatrix}3&3\\3&5\end{pmatrix}(\alpha,\beta)^{\mathsf T}=(5,6)^{\mathsf T}$，解得 $\alpha=7/6$、$\beta=1/2$。殘差 $b-Ax=(-1/6,1/3,-1/6)^{\mathsf T}$ 與兩欄均正交，平方範數為 $1/6$。

5. 實驗步驟

對相同輸入矩陣分別使用直接求解、顯式帶主元 LU 求解，以及適用時的 QR 最小二乘求解。分開記錄分解殘差、解的殘差及相對已知答案的誤差。讓 $\delta$ 構成遞減正數列並固定擾動比率，重複 $A_\delta$ 實驗。將觀察誤差與解析公式 $\varepsilon/\delta$ 比較，注意何時有限精度已無法區分 $1+\delta$ 與一。

擬合直線時畫出數據和擬合線，同時輸出正交性檢驗 $A^*r$ 與殘差平方和。圖形看似合理並不證明最小二乘最優性。保留單位，不要把擬合解釋為因果證據。

6. 練習

P01 — 基礎。 求解 $\begin{pmatrix}3&1\\1&2\end{pmatrix}x=(5,5)^{\mathsf T}$ 並驗證殘差。

P02 — 核心。 對第 2 節的換行例子，計算在右側省略 $P$ 時得到的錯誤結果。

P03 — 核心。 在敏感性例中取 $\delta=10^{-4}$、$\varepsilon=10^{-6}$，計算精確的解變化及其無窮範數。

P04 — 核心。 若 $\kappa(A)\le100$ 且相對殘差至多為 $10^{-8}$，可推出甚麼前向相對誤差界？列出所需的非零條件。

P05 — 診斷。 為甚麼行列式小不足以說明方程組病態？給出一個矩陣族。

P06 — 證明。 證明 $A^*A$ 正定當且僅當 $A$ 的欄線性無關。

P07 — 核心。 計算第 4 節的最小二乘擬合直線，驗證各欄與殘差正交。

P08 — 擴展。 若擬合矩陣的兩欄相等，解釋系數向量為甚麼不唯一，擬合值為何仍可能唯一。

7. 完整解答

P01. 方程 $3x_1+x_2=5$、$x_1+2x_2=5$ 給出 $x=(1,2)^{\mathsf T}$。直接相乘得到 $b$，因此精確殘差為零。

P02. 解 $\begin{pmatrix}2&3\\0&1\end{pmatrix}x=(1,5)^{\mathsf T}$ 得 $x_2=5$、$x_1=-7$。原矩陣將此向量映為 $(5,1)^{\mathsf T}$，而不是 $(1,5)^{\mathsf T}$，因此不能丟掉置換。

P03. $\varepsilon/\delta=10^{-2}$，所以解變化為 $(-0.01,0.01)^{\mathsf T}$，其無窮範數為 $0.01$；而右側擾動的同一範數僅為 $10^{-6}$。

P04. 界為 $100\cdot10^{-8}=10^{-6}$。矩陣必須非奇異，且 $b\ne0$，這保證真解非零、相對表達式有定義。殘差須在與條件數界相兼容的同一範數下計算。

P05. 對 $c\ne0$，$A=cI_n$ 的行列式 $c^n\to0$（當 $c\to0$），但在誘導範數下 $\kappa(A)=\|cI\|\|(cI)^{-1}\|=1$。縮放改變行列式大小，卻不改變條件性。

P06. 對每個 $x$，$x^*A^*Ax=\|Ax\|^2$。它對所有非零 $x$ 嚴格為正，當且僅當 $\ker A=0$，亦即各欄線性無關。

P07. 解所示兩式得 $\beta=1/2$、$\alpha=7/6$。殘差各項為 $-1/6,1/3,-1/6$，其和為零；以 $t=0,1,2$ 加權的和為 $0+1/3-2/6=0$，故 $A^{\mathsf T}r=0$。平方範數為 $1/36+1/9+1/36=1/6$。

P08. 兩欄相等意味著存在一個非零核向量，在相應位置的分量為 $1,-1$。給最小解加上它的任意倍數不會改變 $Ax$，故系數不唯一。欄空間中離 $b$ 最近的向量是其到該子空間的唯一正交投影，所以即使系數表示不唯一，擬合值仍唯一。

8. 總結

分解殘差、數據後向誤差和解的前向誤差回答不同問題。主元選擇保護消元步驟，卻不改變問題條件性。最小二乘是正交投影問題；QR 避免顯式正規方程造成的條件數平方。實驗報告應說明每項測試實際保證了甚麼。