計算數學入門系列 - 線性代數問題（一）

上文提到如何在程式內建構一個矩陣，在這篇文章我們討論一下如何解決一些矩陣的數值問題。最簡單的數值問題，是去找尋N 元一次方程組的答案，數學上來說我們希望找出Ax=b 的根。方程式裏面A 代表一個矩陣，b 代表着方程組的一個已知向量。對於純數學家來說（在本科生課程裏面的線性代數），我們只會問這個方程組有答案嗎？有答案的話，他是唯一的嗎？要真的找出這個方程組的答案，課程裏面，我們也只會用高斯消去法（Gaussian Elimination）一步一步的把答案找出來。對於一些不是很大的問題（就是說N 並不會太大），這個高斯消去法當然不會有太大的問題。方法 一般都可以把確切解（Exact Solution）找出。

可是這裏有兩個問題。第一個，確切解有多exact 其實是可以被質疑的。主要原因是因為前面文章有提過，在電腦裏面我們根本沒有辦法把所有數字確切的標示出來。在那麼多的加減乘除後，我們真的有辦法知道計算出來的數字到底有多準確？當然我們可以說，如果答案準確至小數點後十多個位，誤差一般都不會很大。在一般應用上來說，這個問題應該不會有太大影響。可是這個解釋，對應用數學家來說並沒有太大的吸引力。我們需要的是一個更嚴謹的證明去說明每一步驟微小的不準確，是如何影響着答案的準確度。在數學上來說，這是一個關於算法（Algorithm）本身的穩定性（Stability）的問題。一般在本科生的課程裏面並不會提到。同學們如果有興趣可以參考一些數值線性代數的參考書（譬如說Golub 和Van Loan 的Matrix Computations）。這裏就不再仔細討論。可是比較大的問題，是這個算法的計算量實在太大。簡單來說，要找出一個N 元一次方程組的答案，高斯消去法就差不多需要N的3次方那麼多的計算量。對於應用數學上的問題，這個方法就變得不切實際。

在程式內，很多同學都會第一時間使用MATLAB 的內置函數inv 把矩陣A 的逆

（Inverse）找出來

>> x=inv(A)*b;

這裏沒有錯， x的解（如果我們可以把A的逆找得到）就真的等於這個表達式。 可是，要找出這個矩陣的inverse，可能就真的需要用一些像高斯消去法的方法。所以當矩陣很大時，這個方法就不可行。而且，儘管A本身可能是非常稀疏，他的inverse 一般來說都會是一個全矩陣。這個表達式就先需要把這個全矩陣儲存起來，然後再乘 以方程式右手邊的向量。所以不單是計算量巨大，而且對電腦的儲存也有一定要求。 這個方法也有一個好處，就是如果A並不可逆，MATLAB 會告訴我們

>> A=ones(2,2);

>> b=[1;2];

>> inv(A)*b

Warning: Matrix is singular to working precision. 

 

ans =

   Inf

   Inf