計算方法3 線性方程組求解
前置
這一節主要是玩矩陣，為了偷懶只討論實線性空間
內積
二元函式 \(\left<,\right>\) 被稱為內積，則其滿足：
- \(\left<x,y\right>=\left<y,x\right>\)
- \(\left<ax+by,z\right>=a\left<x,z\right>+b\left<y,z\right>\)
- \(\left<x,x\right>\geq 0\)，等號僅在 \(x=\bold{0}\) 時取到
向量範數
一元函式 \(\norm{\cdot}\) 被稱為範數，則其滿足：
-
\(\norm{x}\geq 0\)，等號僅在 \(x=\bold0\) 時取到
- \(\norm{kx}=\abs{k}\norm{x}\)
- \(\norm{x+y}\leq \norm{x}+\norm{y}\)
在內積空間上有一個天然的範數 \(\norm{x}=\sqrt{\left<x,x\right>}\)，容易驗證滿足上述三條要求。
矩陣範數
運算元範數
矩陣可以看成線性變換，因此可以由向量範數來衡量矩陣作為線性變換（運算元）的“長度”。定義為
\[\norm{A}=\max_{\norm{x}=1} \norm{Ax} \]注意這裡的 \(Ax,x\) 都是向量，因此 RHS 出現的全都是向量範數，LHS 則是定義出來的運算元範數。
有了這個定義，我們就可以寫出這樣的不等式
\[\norm{Ax}\leq \norm{A}\norm{x} \]同時會引入新的定義，稱滿足下述要求的運算元範數具有相容性：
\[\norm{AB}\le \norm{A}\norm{B} \]有了相容性同樣可以做一些界的估計
矩陣範數
矩陣同樣也可以看成線性空間中的元素，因此可以單獨賦予範數的定義，當然這就和向量範數沒什麼關係了。
一個例子是這樣的，容易驗證其滿足三條要求
\[\norm{A}=\sum_{i,j}\abs{A_{i,j}} \]高斯消元
對於給定的線性方程組 \(Ax=b\)，可以用簡單 \(O(n^3)\)
解的穩定性
不妨假設 \(A=xb\) 中 \(\abs{A}\neq0\)，則有 \(x=A^{-1}b\)
通常 \(A\) 是給定的，而 \(b\) 是若干計算和觀察的結果，因此解的誤差主要來源於 \(b\) 引入的誤差，可以寫成
\[\hat x=A^{-1}\hat b \]考慮相對誤差的計算，即為
\[\max_{\hat b,b\neq0} {\frac{\norm{A^{-1}\hat b}}{\norm{A^{-1}b}}}/{\frac{\norm{\hat b}}{\norm{b}}} \]即後向相對誤差與前向相對誤差的比值，稱這個值為方程組 \(A\)（矩陣 \(A\)）的條件數 \(\text{cond}(A)\)
\[\begin{aligned} \text{cond}(A)&=\max_{\hat b,b\neq 0} {\frac{\norm{A^{-1}\hat b}}{\norm{A^{-1}b}}}/{\frac{\norm{\hat b}}{\norm{b}}} \\ &=\max_{\hat b\neq 0} \frac{\norm{A^{-1}\hat b}}{\norm{\hat b}}\max_{b\neq 0}\frac{\norm{b}}{\norm{A^{-1}b}} \end{aligned} \]注意到 \(A^{-1}b=x\)，且 \(b=Ax\)，帶入即得
\[\begin{aligned} \text{cond}(A)&=\max_{\hat b\neq 0}\frac{\norm{A^{-1}\hat b}}{\norm{\hat b}}\max_{x\neq 0} \frac{\norm{Ax}}{\norm{x}} \\ &=\norm{A^{-1}}\norm{A} \end{aligned} \]矩陣的條件數反映了矩陣所對應線性方程組的不穩定程度。條件數越大則不穩定程度越大，誤差的傳遞放大也就越嚴重。
迭代演算法
\(O(n^3)\) 太昂貴，考慮迭代演算法。一般而言迭代的次數是人為規定的，不少演算法能夠保證在 \(\Omega(n)\) 次迭代之後必然得到精確解。
Jacobi 迭代
對於矩陣 \(A\)，將其分解為 \(A=L+D+U\)，其中 \(L,U\) 分別為上三角矩陣和下三角矩陣，\(D\) 為對角陣。
那麼有
\[Ax=(L+U+D)x=b \\ x_{k+1}=D^{-1}(b-(L+U)x_k) \]收斂性
給出一個充分條件：若 \(A\) 是主對角線佔優矩陣，則迭代必然收斂。
只需聯立如下方程
\[\left\{ \begin{aligned} x^*&=D^{-1}(b-(L+U)x^*) \\ x_{k+1}&=D^{-1}(b-(L+U)x_{k}) \end{aligned} \right. \]兩式相減即得
\[x_{k+1}-x^*=-D^{-1}(L+U)(x_k-x^*) \]記 \(W=D^{-1}(L+U)\)，由對角佔優可知 \(Wx\) 的每一項絕對值嚴格小於 \(x\)，因此迭代收斂。
正確性
假設收斂，則有不動點 \(x\)，簡單替換即得不動點 \(x\) 是原方程的一個解。
結合收斂性的充分條件即得：若 \(A\) 為主對角線佔優矩陣，則解唯一（\(A\) 可逆），且迭代收斂至唯一解。
看起來還是比較好的。
Gauss-Seidel 迭代
用到了一個觀察，即在計算解向量 \(x_{k+1}\) 的第 \(r\) 項時，它的前 \(r-1\) 項都已經算出來了（廢話）
因此可以寫成
\[x_{k+1}=D^{-1}(b-Ux_k-Lx_{k+1}) \]證明和上面是類似的，結論也是類似的。
寫程式碼可以發現這兩個方法沒有絕對的好壞之分，迭代速度也沒有一般性的結論（至少俺不知道）。
譜半徑
為了分析一類迭代演算法的收斂性，引入譜半徑的概念
定義 \(A\) 的譜半徑為其絕對值最大的特徵值 \(\abs{\lambda_\max}\)，記為 \(\rho(A)\)
對於這樣的迭代演算法
\[x_{k+1}=Ax_k+b \]取不動點做差得
\[x_k-x^*=A^k(x_0-x^*) \]如果能說明 \(\lim\limits_{k\to\infty} A^k=O\)，那麼就能說明迭代是收斂的，且收斂到方程組的解。
有定理：\(\lim\limits_{k\to\infty}A^k=O\) 當且僅當 \(\rho(A)<1\)
分析 Jacobi 和 Gauss-Seidel 迭代矩陣的譜半徑可以知道，它們以任意初始向量開始迭代都會收斂。