天天看點

《數值分析(原書第2版)》—— 2.3 誤差來源

本節書摘來自華章出版社《數值分析(原書第2版)》一 書中的第2章,第2.3節,作者:(美)timothy sauer,更多章節内容可以通路雲栖社群“華章計算機”公衆号檢視。

正如我們已經描述的,在高斯消去法中有兩個潛在的誤差來源.病态問題的概念和解對于輸入資料的敏感性相關.我們将使用第1章中的後向和前向誤差讨論條件數.在求解病态矩陣方程時,幾乎沒什麼辦法可以避免誤差,因而能夠盡可能識别并避免病态矩陣很重要.第二個誤差來源的原因是淹沒,在大部分問題中可以通過一個簡單的修正,稱為部分主元,進行避免,這部分将在2.4節中進行讨論.85

随後介紹向量和矩陣範數的概念,并用于度量誤差的大小,在目前方程組中的誤差用向量表示.我們将主重點關注一個叫做無窮範數的概念上.

在第1章中,我們發現部分方程求解問題在前向和後向誤差之間有很大差異.對于線性方程組系統也是如此.為了量化誤差,我們首先講述向量的無窮範數概念.

定義2.3 向量x=(x1,…,xn)的無窮範數或者最大範數為‖x‖∞=maxxi,i=1,…,n,即x所有元素中的最大絕對值.

後向誤差和前向誤差和定義1.8相關.後向誤差表示輸入時候的差異,或者問題的資料一方,而前向誤差表示在輸出時候的差異,在算法解的一方.

定義2.4 令xa是線性方程組ax=b的近似解.餘項是向量r=b-axa.後向誤差是餘項的範數‖b-axa‖∞,前向誤差是‖x-xa‖∞.

例2.10 找出近似解xa=[1,1]的後向誤差和前向誤差,方程組如下:11

3-4x1

x2=3

2正确的解是x=[2,1].使用無窮範數,後向誤差是‖b-axa‖∞=3

2-11

3-41

1∞=1

3∞=3前向誤差是‖x-xa‖∞=2

1-1

0∞=1在其他情況下,後向誤差和前向誤差可能具有不同的數量級.

例2.11 找出近似解[-1,3.0001]的後向誤差和前向誤差,方程組如下x1+x2=2

1.0001x1+x2=2.0001(2.17)首先,計算精确解[x1,x2].高斯消去包含步驟11|2

1.00011|2.0001→第2行減去第1行的1.0001倍→11|2

0-0.0001|-0.0001

86求解得到的方程組x1+x2=2

-0.0001x2=-0.0001得到解[x1,x2]=[1,1].

後向誤差是如下向量的無窮範數b-axa=2

2.0001-11

1.00011-1

3.0001

=2

2.0001-2.0001

2=-0.0001

0.0001誤差是0.0001.前向誤差是如下向量差的無窮範數圖2.2 例2.11對應的幾何表示.系統(2.17)表示為直線x2=2-x1和x2=2.0001-1.0001x1,它們的交點在(1,1).點(-1,3.0001)差一點就在兩條直線上,并對應一個近似解.在圖中,兩條直線之間的差異被放大,它們實際上離得非常近x-xa=1

1--1

3.0001=2

-2.0001前向誤差為2.0001.

圖2.2有助于了解為什麼小的後向誤差和大的前向誤差可以同時存在.即使“近似根”(-1,3.0001)相對遠離真實根(1,1),這個近似根也幾乎就位于兩個直線上.由于兩條直線近似平行,是以出現這種情況也是可能的.如果兩條直線不平行,前向誤差和後向誤差将在同一個數量級上.

把餘項表示為r=b-axa.系統ax=b的相對後向誤差定義為‖r‖∞‖b‖∞相對前向誤差定義為

‖x-xa‖∞‖x‖∞87條件 條件數的概念貫穿整個數值分析.在第1章讨論威爾金森多項式時,我們知道當給方程f(x)=0一個小的擾動,如何求“解根過程”的誤差放大因子.對于矩陣方程ax=b,也有相似的誤差放大因子,最大可能的放大因子是cond(a)=‖a‖‖a-1‖.

方程ax=b的誤差放大因子是二者的比率,或者誤差放大因子=相對前向誤差相對後向誤差=‖x-xa‖∞‖x‖∞‖r‖∞‖b‖∞(2.18)對于系統(2.17),相對後向誤差是0.00012.0001≈0.00005=0.005%相對前向誤差是2.00011=2.0001≈200%誤差放大因子是2.001/(0.0001/2.0001)=40004.0001.

在第1章中,我們定義條件數的概念,條件數對應在預先定義的輸入誤差範圍中最大的誤差放大倍數.“預定義範圍”依賴于上下文.現在我們對于目前線性方程組系統将更精确定義這個概念.對于一個固定的矩陣a,考慮對不同向量b求解ax=b.在這個上下文中,b是輸入,對應的解x是輸出.輸入中小的變化對應b的小變化,這對應一個誤差放大因子.我們做如下定義:

定義2.5 方陣a的條件數cond(a)為求解ax=b時,對于所有右側向量b,可能出現的最大誤差放大因子.

令人驚訝的是,對于方陣有一個關于條件數的緊緻的公式.和向量範數類似,定義n×n矩陣a的矩陣範數為‖a‖∞=每行元素絕對值之和的最大值(2.19)即每行元素絕對值求和,并把n行求和的最大值作為矩陣a的範數.

定理2.6 n×n矩陣a的條件數是cond(a)=‖a‖·‖a-1‖後面證明的定理允許我們計算例2.11的系數矩陣的條件數.根據(2.19),矩陣a=11

1.0001188的範數為‖a‖=2.0001.a的逆為a-1=-1000010000

10001-10000對應的範數‖a-1‖=20001.a的條件數是cond(a)=(2.0001)(20001)=40004.0001這正是我們在例2.11中看到的條件數,如此定義的條件數顯然對應最糟的情況.在這個系統中,對于任何其他的b,誤差放大因子将小于或者等于40004.0001.習題3要求計算其他的誤差放大因子.

條件數的重要性和第1章相同.誤差放大因子可能具有cond(a)的數量級.在浮點算術中,相對後向誤差不可能小于εmach,這是由于b的元素的存儲已經引入了和εmach差不多大的誤差.依據(2.18),在求解ax=b可能出現的相關前向誤差是εmach·cond(a).換句話講,如果cond(a)≈10k,我們在計算x時,将丢掉k位數字精度.

在例2.11中,cond(a)≈4×104,因而在雙精度中求解x時,我們應該期望得到16-4=12個正确數位.我們可以通過引入matlab中的通用方程求解器“\”來驗證.

在matlab中,反斜線指令x=a\b使用最先進的lu分解方法求解線性方程組,我們将在第2.4節中進行描述.現在,我們将它作為例子,來看看在浮點算術中我們從最好的算法可以得到什麼結果.下面的matlab指令可以得到例2.10的計算機解xa:

和正确解x=[1,1]相比,計算機求解有大約11個正确的數位,和條件數預測的結果相似.

希爾伯特矩陣h的元素是hij=1/(i+j-1),其對應的條件數非常大.

例2.12 令h表示n×n希爾伯特矩陣.使用matlab的\指令從n=6和10計算hx=b的解,其中b=h·[1,…,1]t.

選擇右側的b,使得對應的解向量中的所有n個元素都是1,這使得前向誤差的檢查十分簡單.matlab使用無窮範數找出條件數,并計算解:

《數值分析(原書第2版)》—— 2.3 誤差來源

條件數大約是107,可以預計在最壞的情況下得到16-7=9個正确的數位;在求解中有大約9位正确的數位.現在計算n=10的情況:

《數值分析(原書第2版)》—— 2.3 誤差來源

由于條件數是1013,在解中隻出現16-13=3個正确的數位.

對于比10大一點兒的n,希爾伯特矩陣的條件數比1016大,在求解xa時不能保證具有正确的數位.

對于病态問題即使最優的軟體也可能無能為力.提高精度會有所幫助;在擴充精度中εmach=2-64≈5.42×10-20,我們開始計算時具有20位而不是16位.但是,希爾伯特矩陣的條件數最終随着n的增長而快速變大,使得任何有意義的有限精度都變得無能為力.

幸運的是,大條件數的希爾伯特矩陣并不常見.n個方程n個未知變量構成的良态線性系統通常在雙精度中可以解到n=104或者更大規模.但同時,知道病态問題是否存在也很重要,而條件數有助于診斷病态問題是否出現.在程式設計問題1~4還有更多關于誤差放大和條件數的例子.

在本節中無窮範數以一種簡單方式給向量定義了一個長度.它是向量範數‖x‖的一個例子,向量範數滿足如下的屬性:

(i) ‖x‖≥0,當且僅當x=[0,…,0]時等号成立.

(ii) 對于每個标量α和向量x,‖αx‖=α·‖x‖.

(iii) 對于向量x、y,‖x+y‖≤‖x‖+‖y‖.

此外,‖a‖∞是矩陣範數的例子,矩陣範數滿足三個和上面相似的性質:

(i) ‖a‖≥0,當且僅當a=0時等号成立.

(ii) 對于每個标量α和矩陣a,‖αa‖=α·‖a‖.

(iii) 對于矩陣a、b,‖a+b‖≤‖a‖+‖b‖.

作為一個不同的例子,向量x=[x1,…,xn]的1-範數是‖x‖1=x1+…+xn.n×n矩陣a的矩陣1-範數是‖a‖1=最大絕對列和,即列向量的1-範數的最大值.關于這些範數定義的驗證,可以參見習題9和習題10.90

剛剛讨論的誤差放大因子、條件數和矩陣範數可用于定義任何向量或者矩陣範數.我們将僅關注稱為算子範數的矩陣範數,這意味着矩陣範數使用特定的向量範數進行定義‖a‖=max‖ax‖‖x‖對于所有非零向量x取最大值.然後,根據該定義,矩陣範數和相關的向量範數一緻,對于任意矩陣a和向量x滿足‖ax‖≤‖a‖·‖x‖(2.20)(2.20)中定義的無窮範數‖a‖∞不僅是矩陣範數,而且是算子範數,驗證見習題10和習題11.

該事實允許我們證明前面所講的對于cond(a)的簡單表示.該證明對于無窮範數和任何算子範數都成立.

定理2.6的證明 我們使用等式a(x-xa)=r以及ax=b.由相容性質(2.20),‖x-xa‖≤‖a-1‖·‖r‖以及1‖b‖≥1‖a‖‖x‖把兩個不等式放在一起得到‖x-xa‖‖x‖≤‖a‖‖b‖‖a-1‖·‖r‖表明‖a‖‖a-1‖是所有誤差放大因子的上界.第二,通常可以取到這個上界.選擇x滿足‖a‖=‖ax‖/‖x‖,以及r滿足‖a-1‖=‖a-1r‖/‖r‖,二者都可以根據算子矩陣範數的定義獲得.令xa=x-a-1r,因而x-xa=a-1r.對于特定選擇的x和r還需要驗證如下等式:‖x-xa‖‖x‖=‖a-1r‖‖x‖=‖a-1‖‖r‖‖a‖‖ax‖2.3.2 淹沒

經典高斯消去法的第二個主要誤差來源可以用更加簡單的方式來修正.我們使用下面的例子來展示淹沒.

例2.13 考慮方程組10-20x1+x2=1

x1+2x2=4我們将三次求解這個方程組:一次使用完全精度,第二次模拟ieee雙精度算術進行求解,第三次我們首先交換方程的順91序.

1.精确解.在表格形式中,高斯消去過程如下10-201|1

12|4→第2行減去第1行的1020倍→10-201|1

02-1020|4-1020底端的方程為(2-1020)x2=4-1020→x2=4-10202-1020頂端的方程得到10-20x1+4-10202-1020=1

x1=10201-4-10202-1020

x1=-2×10202-1020精确解是[x1,x2]=2×10201020-2,4-10202-1020≈[2,1]2.ieee雙精度.計算機版本的的高斯消去有一些不同:10-201|1

02-1020|4-1020ieee雙精度中,2-1020由于舍入等于-1020.類似地,4-1020也被儲存為-1020.現在底端的方程是-1020x2=-1020→x2=1頂端方程的機器算術版本為10-20x1+1=1因而x1=0.得到的計算解是[x1,x2]=[0,1]和真實解相比,這個解有非常大的相對誤差.

3.ieee雙精度,使用行交換.我們改變了兩個方程的順序,重複進行高斯消去法的計算機版本的計算.12|4

10-201|1→第2行減去第1行的10-20倍

→12|4

01-2×10-20|1-4×10-20在ieee雙精度中,1-2×10-20儲存為1,1-4×10-20儲存為1.方程變為x1+2x2=4

x2=1

92得到的計算解是x1=2,x2=1.很顯然,這不是精确解,但是它具有大約16位的精确數位,這是我們使用52位浮點數字所能得到的最大精度.

前面的兩次計算之間的差異很大.第3種方法給出一個可以接受的結果,但是第2種方法卻沒有.關于第2種方法為什麼出了問題,經過分析發現問題出在消去過程中的乘子1020.從底部的方程減去頂部方程的1020倍,底部方程會被抑制,或者稱為“swamp,”.盡管初始的時候有兩個獨立的方程或者源資訊,但是經過第2種計算的消去,這裡僅僅留下了頂部方程的兩個副本.由于底部方程消失,出于所有實際的目的,我們不能指望計算結果可以滿足底部的方程,實際上得到的解也不能滿足對應的方程.

而另一方面,第3種方法,在消去過程中沒有産生覆寫,因為乘子是10-20.在消去之後,原始的兩個方程仍然存在,并且變為三角形式.結果是更加精确的近似解.

例2.13的主要意思是在高斯消去的過程中保證乘子盡可能小,同時避免淹沒.幸運的是,通過對樸素的高斯消去法的簡單修正,可以使得高斯消去法中的乘子的絕對值不大于1.這種新的原則,包括在表格中明智的行交換,稱為部分主元方法,在下一節中将進行詳細的描述.

2.3節習題

1.對于下面的矩陣計算‖a‖∞.

(a) a=12

34(b) a=151

-12-3

1-70

2.找出下述矩陣的條件數(無窮範數).

34(b) a=12.01

36(c) a=63

42

3.找出例2.11中方程組的如下的近似解xa的前向和後向誤差,以及誤差放大因子(使用無窮範數):

(a) [-1,3] (b) [0,2] (c) [2,2] (d) [-2,4] (e) [-2,4.0001]

4.找出x1+2x2=1,2x1+4.01x2=2方程組在如下近似解時的前向和後向誤差,以及誤差放大因子:

(a) [-1,1] (b) [3,-1] (c) [2,-1/2]

5.對于方程組x1-2x2=3,3x1-4x2=7,在如下的近似解處,找出相對的前向和後向誤差以及誤差放大因子:(a) [-2,-4](b) [-2,-3](c) [0,-2](d) [-1,-1](e) 系數矩陣的條件數是多少?

6.對于方程組x1+2x2=3,2x1+4.01x2=6.01,在如下的近似解處,找出相對的前向和後向誤差以及誤差放大因子:(a) [-10,6],(b) [-100,52],(c) [-600,301],(d) [-599,301].(e) 系數矩陣的條件數是多少?93

7.計算5×5希爾伯特矩陣的‖h‖∞.

8.(a) 找出方程組系數矩陣的條件數關于δ>0的函數11

1+δ1x1

x2=2

2+δ(b) 找出近似解xa=[-1,3+δ]處的誤差放大因子.

9.(a) 證明:無窮範數‖x‖∞是向量範數.(b) 證明:1-範數‖x‖1是向量範數.

10.(a) 證明:無窮範數‖a‖∞是矩陣範數.(b) 證明:1-範數‖a‖1是矩陣範數.

11.證明:矩陣無窮範數是向量無窮範數的算子範數.

12.證明:矩陣1-範數是向量1-範數的算子範數.

13.對于習題1中的矩陣,找出滿足‖a‖∞=‖ax‖∞/‖x‖∞的向量.

14.對于習題1中的矩陣,找出滿足‖a‖1=‖ax‖1/‖x‖1的向量.

15.找出矩陣a的lu分解:a=10201

11.996

0501所需的最大數量級的乘子lij是多少?

2.3節程式設計問題

1.對于元素為aij=5/(i+2j-1)的n×n矩陣,令x=[1,…,1]t,b=ax.使用程式設計問題2.1.1的matlab程式或者matlab的反斜線指令計算雙精度解xc.找出ax=b的前向誤差和誤差放大因子的無窮範數,并與a的條件數進行比較:(a) n=6(b) n=10.

2.對于元素是aij=1/(i-j+1)的矩陣執行程式設計問題1.

3.令a是n×n的矩陣,元素是aij=i-j+1.定義x=[1,…,1]t,b=ax.對于n=100,200,300,400,500,使用程式設計問題2.1.1中的matlab程式或者matlab的反斜線指令計算雙精度的計算解xc.對于每個解計算前向誤差的無窮範數.找出問題ax=b的5個誤差放大因子,并和對應的條件數進行比較.

4.對于元素是aij=(i-j)2+n/10的矩陣執行程式設計問題3中的步驟.

5.n取多少,可以使得程式設計問題1中的解沒有正确的有效數字?

6.使用程式設計問題2.1.1中的matlab程式計算例2.13第2和第3種方法的雙精度實作,并和課本中的理論解進行比較.94

繼續閱讀