1、線性代數(shù)機算的速度和精度線性代數(shù)機算的速度和精度機算中精度和速度的重要性 在用筆算時，通常都用整數(shù)矩陣來演示和做題，幾乎不涉及誤差，也就沒有精度問題。至于矩陣計算速度的低下，即使在二、三階的低價矩陣中，就已暴露無遺，但那從來不是純粹進行理論探討的數(shù)學家所關心的內容，甚至很怕讀者觸及這個敏感問題，這是用筆算解矩陣方程的根本弱點。承認這個弱點必然導致計算機的引入和課程的改造。 線性代數(shù)進入應用領域，必定要使用計算機，速度和精度兩個現(xiàn)實問題就不可避免地擺在我們面前。在這里，我們將注重精度問題，并只著重于MATLAB命令中包含的、或會對計算精度作出提示的問題做一介紹。1雙精度浮點數(shù)的精度雙精度浮點數(shù)的

2、精度 按照IEEE標準，表達一個數(shù)需要8個字節(jié)，也就是64個二進制位。雙精度浮點數(shù)用下式表示 其中M是一個小于一大于1/2的二進制分數(shù)，稱為尾數(shù)，占用52個二進制位表示； 而指數(shù)E是一個帶符號的二進制整數(shù)。占11個二進制位，總共可表示2048個整數(shù)，即可以表示從-1023到+1024的數(shù)集，它決定了數(shù)的動態(tài)范圍。數(shù)的正負號反映在S上，它只占一位。總計1+52+11=64位。ESM 2) 1(12/1 MMATLAB中數(shù)的精度 浮點數(shù)的量化步長可以代表它的相對精度。它是由M的位數(shù)決定的。52位二進制數(shù)的量化步長是2 -52=2.220410 16。該數(shù)的動態(tài)范圍取決于指數(shù)部分。因為2 10231

3、0 307及2102410308。所以MATLAB中數(shù)的動態(tài)范圍為2.225110 -3081.797710+308。在MATLAB命令窗中，鍵入eps, realmin, realmax，系統(tǒng)就會給出上面的幾個數(shù)。MATLAB中運算的精度 在做浮點加法、乘法、除法運算時，相對誤差可以基本不變。由于多次運算造成誤差的積累，MATLAB中的大部分運算結果的誤差比eps 要大一些。因此正常情況下（即系統(tǒng)不給出警告時），計算結果至少可以有12位十進制有效數(shù)字的可信度。但是減法有時會有問題，如果遇到兩個很接近的大數(shù)相減，把有效數(shù)位中的前N位都消掉了，那么數(shù)字的精度就減少了N位。比如12345-1234

4、4=1，這兩個數(shù)具有5位有效數(shù)位，但它們的差只有一位有效。 2高斯消元法中的精度問題高斯消元法中的精度問題主主元交換法（元交換法（partial Pivoting） 設方程組如下，寫成矩陣形式： 用消元法化簡增廣矩陣C=A,B 由此得知 x=10000/10001=0.9999; y=10000/10001=0.9999，這是手算的精確解。0.000110.0001 1x1*0-110 xyxyy A X = B0.0001 111 10000100001 1000010000-110-1100 10001100001 10000100001 0(1-1/10001)0110000/10001

5、0 110000/10001C =A,B計算機舍入誤差造成消元誤差 假如我們采用的是一個很粗糙的計算機，只能保持三位有效數(shù)字，那么它遇到10001時，只會四舍五入解讀為10000。則上面的消元過程將成為： 這個結果將解讀為 x=0, y=1，x的誤差達到了100%，完全不能采用。如果計算機的精度是N位，其結果雖然不至于完全無用，但也將使有效數(shù)位減少4位。 0.0001 111 10000100001 1000010000-110-1100 10000100001 10000100001 000110 11C =A,B主元太小是造成誤差的根源 從計算過程看，其實問題就出在回代過程中出現(xiàn)了相接近的

6、兩個大數(shù)相減。兩個大數(shù)出現(xiàn)于把微小的主元0.0001化為1的消元過程中，主元太小是造成誤差的根源。如果在做消元法的時候，檢查主元行中各個元素，把最大的主元通過列交換放到主元的位置，那樣這個問題就可以解決，例如上題中先交換第一行中的兩列，可得： 0.0001111 0.000111 0.00011-1101-100 -1.0001-11 0.000111 0 1 1/1000101 1/1.00010 11/1.0001C =A,B在rref函數(shù)中的換主元措施 由此得知 x=10000/10001=0.9999; y=0.9999，這和精確解完全一致，說明只要進行主元最大的列交換。計算精度是可以

7、保證的。 在rref命令中，語句p,k = max(abs(A(i:m,j);就是為了在第i行中，找到絕對值最大的元素A(i,k)，以后再進行列交換。 在MATLAB命令窗中鍵入： A=0.0001,1;-1,1,B=1;0,X=AB， 得到 X= 0.9999；0.9999鍵入type rref得出的主要程序代碼while (i = m) & (j = n) % Find value and index of largest element in the remainder of column j. p,k = max(abs(A(i:m,j); k = k+i-1; if (p n

8、,即方程數(shù)多于變量數(shù)的情況，嚴格地滿足Ax=b的解x是不存在的。要找到近似的 ，就要引入誤差e， Ax-b=e x誤差線性方程組的建立 引入誤差向量e。 e Ax b 寫出其完全的矩陣形式如下 問題是，找到解x，使e的長度或范數(shù)為最小。 1111112221nmmnmnmexbaaexbaaexbeAx -b從向量空間的視點分析 研究例6.1的超定方程組（d）： 改寫成簡寫為 選擇不同的x1和x2將得到不同的合成向量A*xx1v1x2v2q ，q必定處于v1和v2張成的平面之內。而方程中的b則一般不會在這個平面內，所以精確解通常不可能存在。 12111113123xx 12*xx12Axvvb

9、本例的向量空間圖 這時最近似的解就應該是該平面上與b點最近的點所對應的坐標A*xhat。它應該是b點向v1和v2張成的平面的投影。所以A*xhat和b的連線應該和v1和v2張成的平面垂直，也就是說必須分別與v1和v2正交。如圖8.6所示。最小二乘解的公式推導 A*xhat和b的連線向量應該是這兩個向量之差，即，它與v1和v2正交的要求可以分別表示為： 和 綜合在一起可以寫成：最后得到公式Ax - b1Tv (Ax -b)= 0T2v (Ax -b)= 0T1TT2vA (Ax - b)=(Ax - b)= 0vT-1Tx = (A A) A b最小二乘解的數(shù)字例 例8. 求題6.1(d)方程組

10、的最小二乘解。 解：MATLAB程序ag808如下：A=1,1;1,1; 1,2, b=1;3;3xhat=inv(A*A)*A*be=A*xhatb, norm(e) 運行此程序，得到 1.00001.0000,3.0000 ,1.00002.0000 xe 1 11 1 1 ,3 ,1 23 AbMATLAB中超定方程的解 在MATLAB中，把運算(ATA) -1AT單獨編成一個子程序，稱為pinv函數(shù)。求最小二乘解的公式可以寫成 xpinv(A)*b，與適定方程的解x inv(A)* b非常相似，只是pinv函數(shù)并不要求A是方陣。 最小二乘解也可用運算符表示，這就把欠定方程、適定方程和超

11、定方程用統(tǒng)一的運算格式： x A b MATLAB會自動根據(jù)系數(shù)矩陣A的行數(shù)m和列數(shù)n，來判斷采用哪個方法和程序。不過對于欠定方程，這個式子只給出了一個特解，沒給通解。 二。二。 機算速度問題機算速度問題 運算時間通常用完成某一運算所需做的乘加法次數(shù)來判斷。消元法正向消元所需的乘法次數(shù)為(n3-n)/3，回代所需的次數(shù)少，為n2/2，在n很大時，近似按n3/3估計。N=10時，為333次；n=25時為5208次。對現(xiàn)有的最新計算機（2008年測試的最高紀錄是每秒1015次，測試所用的標準就是用LINPACK軟件包解方程AX=b），這都可以在幾十個微微秒中完成。 求行列式所需的時間與此相同。但如果按行列式的數(shù)序定義來求，25階的行列式，就要做3.7*1026次乘法，用最新的計算機，也要100萬年。所以消元法是一個巨大的成就，雖然以高斯命名，其實是中國人在1000多年前的“九章算法”中最先提出的方法， 矩陣形式對計算速度影響 矩陣的構造形式對計算速度影響也很大。最顯然的是對角矩陣，它完全可以按標量計算，算一個行列式只需要做n-1次乘法，解一個方程組也是這樣。介乎滿元素矩陣和對角矩陣之間的是對角帶型矩陣，其特點是所有的元素聚集在對角元素的一側或兩側，如下圖所示：幾種函數(shù)所需的運算次數(shù) 對這種單邊帶寬為w的矩陣進行LU分解（即前向消元