第8章矩陣特征值計算

1、上頁上頁下頁下頁第第8章章 矩陣特征值計算矩陣特征值計算 8.1 特征值性質(zhì)和估計特征值性質(zhì)和估計 8.2 冪法及反冪法冪法及反冪法 8.3 正交變化與矩陣分解正交變化與矩陣分解 8.4 QR方法方法上頁上頁下頁下頁8.1 特征值性質(zhì)和估計特征值性質(zhì)和估計 工程技術(shù)中有多種振動問題，如橋梁或建筑物的工程技術(shù)中有多種振動問題，如橋梁或建筑物的振動，機(jī)械零件、飛機(jī)機(jī)翼的振動，及一些穩(wěn)定性分振動，機(jī)械零件、飛機(jī)機(jī)翼的振動，及一些穩(wěn)定性分析和相關(guān)分析在數(shù)學(xué)上都可轉(zhuǎn)化為求矩陣特征值與特析和相關(guān)分析在數(shù)學(xué)上都可轉(zhuǎn)化為求矩陣特征值與特征向量的問題征向量的問題. 下面先復(fù)習(xí)一些矩陣的特征值和特征向量的基礎(chǔ)下面

2、先復(fù)習(xí)一些矩陣的特征值和特征向量的基礎(chǔ)知識及性質(zhì)知識及性質(zhì).上頁上頁下頁下頁8.1.1 特征值問題及性質(zhì)特征值問題及性質(zhì)求求A的特征值問題的特征值問題(1.1)等價于求等價于求A的特征方程的特征方程設(shè)矩陣設(shè)矩陣A Rnn，特征值問題是求，特征值問題是求 C和非零向和非零向量量x Rn，使，使( )det()0(1.2)pIA的根的根. 在在5.1.3節(jié)已經(jīng)給出特征值的一些性質(zhì)，下面節(jié)已經(jīng)給出特征值的一些性質(zhì)，下面再補(bǔ)充一些基本性質(zhì)再補(bǔ)充一些基本性質(zhì). Ax= x (1.1)其中其中 是矩陣是矩陣A的的特征值特征值，x是矩陣是矩陣A屬于特征值屬于特征值 的的特征向量特征向量. 本章討論計算矩陣特

3、征值的數(shù)值方法本章討論計算矩陣特征值的數(shù)值方法, 在在科學(xué)和工程技術(shù)中很多問題在數(shù)學(xué)上都?xì)w結(jié)為求矩科學(xué)和工程技術(shù)中很多問題在數(shù)學(xué)上都?xì)w結(jié)為求矩陣的特征值問題陣的特征值問題.上頁上頁下頁下頁 定理定理1 設(shè)設(shè) 為為ARnn的特征值的特征值, 且且Ax= x, x 0，則有則有 (2) - - 為為A- - I的特征值，即的特征值，即(A- - I)x=( - - )x ； (1) c 為的為的cA特征值特征值(c 0為常數(shù)為常數(shù))； (3) k為為Ak的特征值，即的特征值，即Akx= kx ； (4) 設(shè)設(shè)A為非奇異矩陣，那么為非奇異矩陣，那么 0 , 且且 - -1為為A- -1的特的特征值，

4、即征值，即A- -1x= - -1x .上頁上頁下頁下頁 定理定理2 (1) 設(shè)矩陣設(shè)矩陣ARnn可對角化，即存在非可對角化，即存在非奇異矩陣奇異矩陣P使使的充分必要條件是的充分必要條件是A具有具有n個線性無關(guān)的特征向量個線性無關(guān)的特征向量. (2) 如果矩陣如果矩陣ARnn有有m個個(mn)不同的特征不同的特征值值 1, 2, m，則對應(yīng)的特征向量，則對應(yīng)的特征向量 x1, x2, , xm 線線性無關(guān)性無關(guān).121,nPAP 上頁上頁下頁下頁 定理定理3(對稱矩陣的正交約化對稱矩陣的正交約化) 設(shè)設(shè)ARnn為對稱為對稱矩陣，則矩陣，則 (3) 存在一個正交矩陣存在一個正交矩陣P使的使的且且

5、 1, 2, n為為A的特征值，而的特征值，而P(u1,u2,un)的列的列向量向量uj為為A的對應(yīng)于的對應(yīng)于 j 的單位特征向量的單位特征向量.12,TnP AP (1) A的特征值均為實數(shù)；的特征值均為實數(shù)； (2) A有有n個線性無關(guān)的特征向量；個線性無關(guān)的特征向量；上頁上頁下頁下頁 定義定義 設(shè)設(shè)ARnn為對稱矩陣，對于任一非零向為對稱矩陣，對于任一非零向量量x 0，稱，稱,),(),()(xxxAxxR 為對應(yīng)于向量為對應(yīng)于向量x的的瑞利瑞利(Rayleigh)商商. 定理定理4 設(shè)設(shè)ARnn為對稱矩陣為對稱矩陣(其特征值依次記其特征值依次記為為 1 2 n)，則，則 (1). (對

6、任何非零對任何非零xRn)；1(, )( , )nAx xx x (2). ； (1.3)10(, )max( , )nx RxAx xx x 0(, )min.( , )nnx RxAx xx x 上頁上頁下頁下頁 證明證明 只證只證(1)，關(guān)于，關(guān)于(2)自己作練習(xí)自己作練習(xí). 由于由于A為實對稱矩陣，可將為實對稱矩陣，可將 1, 2, n 對應(yīng)的特對應(yīng)的特征向量征向量 x1,x2,xn 正交規(guī)范化，則有正交規(guī)范化，則有 (xi, xj)=ij，設(shè)，設(shè)x 0為為Rn中任一向量，則有中任一向量，則有. 0,1221 niiniiixxx 于是于是.),(),(11212 niiniiinxx

7、xAx從而從而(1)成立成立. 結(jié)論結(jié)論(1)說明說明瑞利商瑞利商必位于必位于 n和和 1之間之間.上頁上頁下頁下頁8.1.2 特征值估計與擾動特征值估計與擾動 定義定義1 設(shè)設(shè)n階矩陣階矩陣A=(aij)，令，令).,(niarnijjiji211 (1) ； (2) 集合集合 稱為復(fù)平面上以稱為復(fù)平面上以aii為圓心，以為圓心，以ri為半徑的為半徑的n階矩陣階矩陣A的的n個個格什戈林格什戈林(Gerschgorin)圓盤圓盤. |,(1,2, )iiiiDzzarzCin 上頁上頁下頁下頁 定理定理5 (Gerschgorin圓盤定理圓盤定理) 特別地，如果特別地，如果A的一個圓盤的一個圓

8、盤Di是與其它圓盤分離是與其它圓盤分離(即即孤立圓盤孤立圓盤)，則，則Di中精確地包含中精確地包含A的一個特征值的一個特征值.1(1,2,. ).(1.4)niiiijjj iarain (1) 設(shè)設(shè)n階矩陣階矩陣A(aij)，則，則A的每一個特征值必屬的每一個特征值必屬于下面某個圓盤之中于下面某個圓盤之中 (2) 如果如果A有有m個個圓盤圓盤組成一個連通的并集組成一個連通的并集S，且，且S與余下與余下n- -m個個圓盤圓盤是分離的，則是分離的，則S內(nèi)恰包含內(nèi)恰包含A的的m個個特征值特征值.或者說或者說 A的特征值都在的特征值都在n個圓盤的個圓盤的并集并集中中.上頁上頁下頁下頁 證明證明 只就

9、只就(1)給出證明給出證明. 設(shè)設(shè) 為為A的特征值，即的特征值，即.1knkjjkjkkraa Ax= x，其中，其中 x=(x1,x2, xn)T 0.或或 記記 ，考慮，考慮Ax= x的第的第k個個方程，即方程，即0max1 xxxinik,1knjjkjxxa ,)( kjjkjkkkxaxa 于是于是即即, kjkjkkjjkjkkkaxxaxa 上頁上頁下頁下頁這說明，這說明，A的每一個特征值必位于的每一個特征值必位于A的一個圓盤的一個圓盤中，并且相應(yīng)的特征值中，并且相應(yīng)的特征值 一定位于第一定位于第k個圓盤中個圓盤中(其中其中k是對應(yīng)特征向量是對應(yīng)特征向量x絕對值最大的分量的下標(biāo)絕

10、對值最大的分量的下標(biāo)).利用相似矩陣性質(zhì)，有時可以獲得利用相似矩陣性質(zhì)，有時可以獲得A的特征值進(jìn)的特征值進(jìn)一步的估計，即適當(dāng)選取非奇異對角陣一步的估計，即適當(dāng)選取非奇異對角陣,112111 nD 并做相似變換并做相似變換 .適當(dāng)選取適當(dāng)選取 可使某些圓盤半徑及連通性發(fā)生變化可使某些圓盤半徑及連通性發(fā)生變化.nnijijaADD 1), 2 , 1(nii 上頁上頁下頁下頁.411101014 A 例例1 估計矩陣估計矩陣A的特征值范圍，其中的特征值范圍，其中 解解 矩陣矩陣A的的3個圓盤為個圓盤為. 24:, 2:, 14:321 DDD 由定理由定理8，可知，可知A的的3個特征值位于個特征值

11、位于3個圓盤的并個圓盤的并集中，由于集中，由于D1是孤立圓盤，所以是孤立圓盤，所以D1內(nèi)恰好包含內(nèi)恰好包含A的的一個特征值一個特征值 1(為實特征值為實特征值)，即，即. 531 A的其它兩個特征值的其它兩個特征值 2, 3包含在包含在D2, D3的并集中的并集中.上頁上頁下頁下頁現(xiàn)在取對角陣現(xiàn)在取對角陣,9 . 0000100011 D做相似變換做相似變換.49 . 09 . 00101491011 ADDAA 矩陣矩陣A1的的3個圓盤為個圓盤為. 8 . 14:,919:, 14:321 EEE上頁上頁下頁下頁顯然，顯然，3個圓盤都是孤立圓盤，所以，每一個圓個圓盤都是孤立圓盤，所以，每一個

12、圓盤都包含盤都包含A的一個特征值的一個特征值(為實特征值為實特征值)，且有估計，且有估計 . 2 . 28 . 5,919919, 53321 上頁上頁下頁下頁 定理定理6(Bauer-Fike定理定理) 設(shè)設(shè) 是是A+IRnn的一個的一個特征值，且特征值，且P- -1AP=D=diag( 1, 2, n)，則有，則有1()min,(1.5)pppAPPI 其中其中|p為矩陣的為矩陣的p范數(shù)，范數(shù)，p=1,2, . 證明證明 由于由于 (A)時顯然成立時顯然成立,故只考慮故只考慮 (A).這這時時D- - I非奇異，設(shè)非奇異，設(shè)x是是A+I對應(yīng)于對應(yīng)于 的特征向量，由的特征向量，由(A+I-

13、- I)x=0左乘左乘P- -1可得可得111()()()(),DIPxPIPPx 1111() ()(),PxDIPIPPx 而對角矩陣而對角矩陣(D- - I)- -1的范數(shù)為的范數(shù)為1()1(),min,pADImm 所以有所以有1.pppPIPm 這就得到這就得到(1.5)式式. 這時總有這時總有中的一個中的一個 取到取到m值值.上頁上頁下頁下頁 由定理由定理6可知可知|P- -1| |P|=cond(P)是特征值擾動的是特征值擾動的放大系數(shù)，但將放大系數(shù)，但將A對角化的相似變化矩陣不是唯一的，對角化的相似變化矩陣不是唯一的，所以取所以取cond(P)的下確界的下確界 112()inf

14、()(,) ,(1.6)nAcond P PAPdiag 稱為特征值問題的條件數(shù)稱為特征值問題的條件數(shù). 只要只要v(A)不很大，矩陣微不很大，矩陣微小擾動只帶來特征值的微小擾動小擾動只帶來特征值的微小擾動. 但是但是v(A)難以計算難以計算,有時只對一個有時只對一個P，用，用cond(P)代替代替v(A).特征值問題的條件數(shù)和解線性方程組時的矩陣是特征值問題的條件數(shù)和解線性方程組時的矩陣是兩個不同的概念，對于一個矩陣兩個不同的概念，對于一個矩陣A，兩者可能一大一，兩者可能一大一小，例如二階矩陣小，例如二階矩陣A=diag(1, 10- -10)，有，有v(A) =1，但解，但解線性方程組的矩

15、陣條件數(shù)線性方程組的矩陣條件數(shù)cond(A)=1010.上頁上頁下頁下頁 關(guān)于計算矩陣關(guān)于計算矩陣A的特征值問題，當(dāng)?shù)奶卣髦祮栴}，當(dāng)n=2,3時，我們時，我們還可按行列式展開的辦法求特征方程還可按行列式展開的辦法求特征方程p( )=0的根的根. 但但當(dāng)當(dāng)n較大時，如果按展開行列式的辦法較大時，如果按展開行列式的辦法, 首先求出首先求出p( )的系數(shù)，再求的系數(shù)，再求p( )的根，工作量就非常大，用這種辦的根，工作量就非常大，用這種辦法求矩陣特征值是不切實際的，由此需要研究法求矩陣特征值是不切實際的，由此需要研究A的特的特征值及特征向量的數(shù)值方法征值及特征向量的數(shù)值方法. 本章將介紹一些計算機(jī)上

16、常用的兩類方法，一類本章將介紹一些計算機(jī)上常用的兩類方法，一類是冪法及反冪法是冪法及反冪法(迭代法迭代法)，另一類是正交相似變化的，另一類是正交相似變化的方法方法(變化法變化法).上頁上頁下頁下頁8.2 冪法及反冪法冪法及反冪法 冪法與反冪法都是求冪法與反冪法都是求實矩陣實矩陣的特征值和特征向的特征值和特征向量的量的向量迭代法向量迭代法，所不同的是冪法是計算矩陣的，所不同的是冪法是計算矩陣的主主特征值特征值(矩陣矩陣按模最大的特征值按模最大的特征值稱為稱為主特征值主特征值，其模，其模就是該矩陣的就是該矩陣的譜半徑譜半徑)和和相應(yīng)特征向量相應(yīng)特征向量的一種向量迭的一種向量迭代法，特別適用于大型稀

17、疏矩陣代法，特別適用于大型稀疏矩陣. 而反冪法則是計而反冪法則是計算非奇異算非奇異(可逆可逆)矩陣矩陣按模最小的特征值按模最小的特征值和和相應(yīng)特征相應(yīng)特征向量向量的一種向量迭代法，特別是計算海森伯格陣或的一種向量迭代法，特別是計算海森伯格陣或三對角陣的對應(yīng)一個給定近似特征值的特征向量的三對角陣的對應(yīng)一個給定近似特征值的特征向量的有效方法之一有效方法之一. 下面分別介紹冪法與反冪法下面分別介紹冪法與反冪法. 上頁上頁下頁下頁8.2.1 冪法冪法(又稱又稱乘冪法乘冪法)現(xiàn)討論求現(xiàn)討論求 1及及x1的方法的方法.), 2 , 1(nixAxiii 設(shè)實矩陣設(shè)實矩陣A=(aij)有一個有一個完全的特征

18、向量組完全的特征向量組，即，即A有有n個線性無關(guān)的特征向量個線性無關(guān)的特征向量，設(shè)矩陣，設(shè)矩陣A的特征值為的特征值為 1, 2, n, 相應(yīng)的特征向量為相應(yīng)的特征向量為x1,x2,xn. 已知已知A的的主特征值主特征值 1是實根，且滿足條件是實根，且滿足條件)1 . 2(|,|21n 上頁上頁下頁下頁)3 . 2(),0(122110 axaxaxavnn設(shè)設(shè) 冪法的冪法的基本思想基本思想是是: 任取任取非零非零的初始向量的初始向量v0 , 由矩由矩陣陣A構(gòu)造一向量序列構(gòu)造一向量序列vk稱為迭代向量，由假設(shè)，稱為迭代向量，由假設(shè)，v0可唯一表示為可唯一表示為)2 . 2(.,.,0110212

19、01 vAAvvvAAvvAvvkkk上頁上頁下頁下頁)3 . 2(),0(122110 axaxaxavnn設(shè)設(shè)于是于是).()/(1112111122211101kkniikiiknknnkkkkkxaxaxaxaxaxavAAvv 其中其中.)/(21 niikiikxa 由假設(shè)由假設(shè) 故故 從而從而), 3 , 2( 1/1nii , 0lim kk )4 . 2(.lim111xavkkk 為為 1的特征向量的特征向量.上頁上頁下頁下頁所以當(dāng)所以當(dāng)k充分大時，有充分大時，有)5 . 2(,111xavkk 即為矩陣即為矩陣A的的對應(yīng)特征值對應(yīng)特征值 1 的一個近似特征向量的一個近似特

20、征向量.用用(vk)i 表示表示vk的第的第i個分量，則當(dāng)個分量，則當(dāng)k充分大時，有充分大時，有 )7 . 2(.11 ikikvv即為即為A的的主特征值主特征值 1的近似值的近似值.)6 . 2(,111111kkkkvxaAvv 由于由于 這種由已知非零向量這種由已知非零向量v0及矩陣及矩陣A的乘冪的乘冪Ak構(gòu)造向構(gòu)造向量序列量序列vk以計算以計算A的的主特征值主特征值 1(利用利用(2.7)式式)及相應(yīng)及相應(yīng)特征向量特征向量(利用利用(2.5)式式)的方法就稱為的方法就稱為(乘乘)冪法冪法.上頁上頁下頁下頁 迭代公式實質(zhì)上是由矩陣迭代公式實質(zhì)上是由矩陣A的乘冪的乘冪 Ak與非零向與非零向

21、量量v0相乘來構(gòu)造向量序列相乘來構(gòu)造向量序列vk=Akv0，從而計算主特，從而計算主特征值征值 1及其對應(yīng)的特征向量，這就是及其對應(yīng)的特征向量，這就是冪法冪法的思想的思想. ).(11 kvvikik 的收斂速度由比值的收斂速度由比值,12 r來確定，來確定，r 越小收斂越快，但當(dāng)越小收斂越快，但當(dāng)r 1 1時收斂可能很慢時收斂可能很慢.總結(jié)上述討論，有下面的定理總結(jié)上述討論，有下面的定理.上頁上頁下頁下頁 定理定理7 設(shè)設(shè)ARnn有有n個線性無關(guān)的特征向量，個線性無關(guān)的特征向量，主特征值主特征值 1滿足滿足 | 1| 2| | n|，則對任何非零向量則對任何非零向量v0(a1 0)，(2.4

22、)式和式和(2.7)式成立式成立.又設(shè)又設(shè)A有有n個線性無關(guān)的特征向量，個線性無關(guān)的特征向量， 1對應(yīng)的對應(yīng)的r個線性個線性無關(guān)的特征向量為無關(guān)的特征向量為x1,x2,xr，則由，則由(2.2)式有式有 如果如果A的主特征值為實的重根的主特征值為實的重根, 即即 1= 2= r, 且且 | r| r+1| | n|，,)/(11110 nriikiiriiikkkxaxavAv 上頁上頁下頁下頁).0(lim111 riiiriiikkkxaxav設(shè)設(shè) 為為 1的特征向量，這說明當(dāng)?shù)奶卣飨蛄?，這說明當(dāng)A的主特征值是實的重根的主特征值是實的重根時，定理時，定理7的結(jié)論還是正確的的結(jié)論還是正確的.

23、 應(yīng)用應(yīng)用冪法冪法計算計算A的主特征值的主特征值 1及其對應(yīng)的特征向及其對應(yīng)的特征向量時，如果量時，如果| 1|1(或或| 1| 2| | n|，則對任意非零初始，則對任意非零初始向量向量v0=u0(a1 0)，有冪法計算公式為，有冪法計算公式為(uk,vk)則有則有 (1),)max(lim11xxukk (2).lim1 kk上頁上頁下頁下頁例例2 用冪法計算矩陣用冪法計算矩陣 0 . 225. 05 . 025. 00 . 10 . 15 . 00 . 10 . 1A的主特征值和相應(yīng)的特征向量的主特征值和相應(yīng)的特征向量.解解取取 v0=u0=(1,1,1)T, 則則 TTvuv1 ,81

24、82. 0,9091. 01,75. 2,75. 2 ,25. 2 , 5 . 211111 ), 2 , 1(max1 k/vuvAuvkkkkkkk 上頁上頁下頁下頁計算結(jié)果如下表計算結(jié)果如下表k0(1, 1, 1) 1(0.9091, 0.8182, 1) 2.75000005(0.7651, 0.6674, 1)2.558791810(0.7494, 0.6508, 1)2.538002915(0.7483, 0.6497, 1)2.536625616(0.7483, 0.6497, 1)2.536584017(0.7482, 0.6497, 1)2.536559818(0.7482,

25、 0.6497, 1)2.536545619(0.7482, 0.6497, 1)2.536537420(0.7482, 0.6497, 1)2.5365323)(規(guī)規(guī)范范化化向向量量Tkumaxkkv 上頁上頁下頁下頁,5365323. 21 這個結(jié)果是用這個結(jié)果是用8位浮點數(shù)字進(jìn)行運算得到的位浮點數(shù)字進(jìn)行運算得到的, uk的分量值是舍入值的分量值是舍入值. 于是得到于是得到及相應(yīng)的特征向量及相應(yīng)的特征向量(0.7482, 0.6497, 1)T. 1和相應(yīng)的特和相應(yīng)的特征向量的真值征向量的真值(8位數(shù)字位數(shù)字)為為,5365258. 21 .)1,64966116. 0,74822116.

26、 0(1Tx 上頁上頁下頁下頁例例 用冪法計算矩陣用冪法計算矩陣 1634310232A的主特征值與其對應(yīng)的特征向量的主特征值與其對應(yīng)的特征向量.解解取取 v0=u0=(0,0,1)T , 則則 TTvuv25. 0 , 1, 5 . 01, 4,1 , 4 , 211111 ), 2 , 1(max1 k/vuvAuvkkkkkkk 上頁上頁下頁下頁直到直到k=8 時的計算結(jié)果見下表時的計算結(jié)果見下表k1 2, 4, 1, 4 0.5, 1, 0.252 4.5, 9, 7.75 90.5, 1, 0.86113 5.7222, 11.4444, 8.36111.44440.5, 1, 0.

27、73604 5.4621, 10.9223, 8.2306 10.92230.5, 1, 0.75365 5.5075, 11.0142, 8.2576 11.01420.5, 1, 0.74946 5.4987, 10.9974, 8.2494 10.99740.5, 1, 0.75017 5.5002, 11.0005, 8.2501 11.00050.5, 1, 0.75008 5.5000, 11.0000, 8.2500 11.00000.5, 1, 0.7500TkvTku Tx7500. 0,0 . 1,5 . 0,0000.1111 從而從而k 上頁上頁下頁下頁8.2.2 加速

28、方法加速方法1、原點平移法、原點平移法 由前面討論知道，應(yīng)用冪法計算由前面討論知道，應(yīng)用冪法計算A的主特征值的的主特征值的收斂速度主要由比值收斂速度主要由比值 r=| 2/ 1|來決定，但當(dāng)來決定，但當(dāng)r 接近于接近于1時，收斂可能很慢時，收斂可能很慢. 這時，一個補(bǔ)救辦法是采用加速這時，一個補(bǔ)救辦法是采用加速收斂的方法收斂的方法.其中其中p為參數(shù)，設(shè)為參數(shù)，設(shè)A的特征值為的特征值為 i，則對矩陣，則對矩陣B的特征的特征值為值為 i- -p ，而且，而且A, B的特征向量相同的特征向量相同. 引進(jìn)矩陣引進(jìn)矩陣 B=A- -pI .上頁上頁下頁下頁 如果要計算如果要計算A的主特征值的主特征值 1

29、, 只要只要選擇合適的數(shù)選擇合適的數(shù)p，使使 1- -p為矩陣為矩陣B=A- -pI 的主特征值，且的主特征值，且 1212max ppini那么，對矩陣那么，對矩陣B=A- -pI應(yīng)用冪法求其主特征值應(yīng)用冪法求其主特征值 1- -p, 收收斂速度將會加快斂速度將會加快. 這種通過求這種通過求B=A- -pI的主特征值和特的主特征值和特征向量，而得到征向量，而得到A的主特征值和特征向量的方法叫的主特征值和特征向量的方法叫原原點平移法點平移法. 對于對于A的特征值的某種分布，它是十分有的特征值的某種分布，它是十分有效的效的.上頁上頁下頁下頁例例3 設(shè)設(shè)AR44有特征值有特征值),4 , 3 ,

30、2 , 1(15 jji 比值比值r=| 2/ 1| 0.9. 做變換做變換B=A- -12I (p=12),則則B的特征值為的特征值為. 1, 0, 1, 24321 應(yīng)用冪法計算應(yīng)用冪法計算B的主特征值的主特征值1的收斂速度的比值為的收斂速度的比值為. 9 . 021121212 pp 雖然常常能夠選擇有利的雖然常常能夠選擇有利的p值值, 使冪法得到加速使冪法得到加速, 但設(shè)計一個自動選擇適當(dāng)參數(shù)但設(shè)計一個自動選擇適當(dāng)參數(shù)p的過程是困難的的過程是困難的. 1212max142 ii上頁上頁下頁下頁 下面考慮當(dāng)下面考慮當(dāng)A的特征值是實數(shù)時，怎樣選擇的特征值是實數(shù)時，怎樣選擇p使采使采用冪法計

31、算用冪法計算 1得到加速得到加速.ppn 1且使且使收斂速度的比值收斂速度的比值.min,max112 ppppn 設(shè)設(shè)A的特征值都是實數(shù)，且滿足的特征值都是實數(shù)，且滿足)10. 2(,121nn 則對實數(shù)則對實數(shù)p，使矩陣，使矩陣A- -pI的主特征值為的主特征值為 1- -p或或 n- -p時，時，當(dāng)當(dāng)我們計算我們計算 1及及x1時，首先應(yīng)選取時，首先應(yīng)選取p使使上頁上頁下頁下頁顯然，當(dāng)顯然，當(dāng) 2- -p=- -( n- -p )時，即時，即P=( 2+ n)/2P*時時 為為最小值，這時最小值，這時收斂速度的比值收斂速度的比值為為.2212112nnnpppp 當(dāng)希望計算當(dāng)希望計算 n

32、時，應(yīng)選取時，應(yīng)選取 p=( 1+ n-1)/2P* 使得應(yīng)用冪法計算使得應(yīng)用冪法計算 n得到加速得到加速. 當(dāng)當(dāng)A的特征值都是實數(shù)，滿足的特征值都是實數(shù)，滿足(2.10)式式且且 2, n能初步估計出來，我們就能確定能初步估計出來，我們就能確定P*的近似值的近似值.nn 121上頁上頁下頁下頁 例例4 用原點平移加速法求用原點平移加速法求例例2中矩陣中矩陣A的主特征值的主特征值與其對應(yīng)的特征向量與其對應(yīng)的特征向量.25. 125. 05 . 025. 025. 00 . 15 . 00 . 125. 0 B,0 . 225. 05 . 025. 00 . 10 . 15 . 00 . 10

33、. 1 A對對B應(yīng)用冪法，仍應(yīng)用冪法，仍取取 v0=(1,1,1)T , 則則 .1,75. 0,875. 01, 2,2, 5 . 1,75. 111111TTvuv ), 2 , 1(max1 k/vuvBuvkkkkkkk 解解 取取p=0.75, 做平移變換做平移變換B=A- -pI，則，則上頁上頁下頁下頁k0(1, 1, 1) 1(0.8750, 0.7500, 1) 2.00000005(0.7516, 0.6522, 1)1.79140116(0.7491, 0.6511, 1)1.78884437(0.7488, 0.6501, 1)1.78733008(0.7484, 0.6

34、499, 1)1.78691529(0.7483, 0.6497, 1)1.786658710(0.7482, 0.6497, 1)1.7865914)(規(guī)規(guī)范范化化向向量量Tkumaxkkv 計算結(jié)果如下表計算結(jié)果如下表由此得由此得B的主特征值為的主特征值為 1 1.7865914，A的主特征值的主特征值 1為為 1 = 1 +0.75 2.5365914.上頁上頁下頁下頁 原點位移的加速方法，是一個矩陣變換方法原點位移的加速方法，是一個矩陣變換方法. 這這種變換容易計算，又不破壞矩陣種變換容易計算，又不破壞矩陣A的稀疏性，但的稀疏性，但p的的選擇選擇依賴依賴對對A的的特征值分布的大致了特征

35、值分布的大致了解解.及相應(yīng)的特征向量為及相應(yīng)的特征向量為 x1 (0.7482, 0.6497, 1)T. 與例與例2結(jié)果比較，上述結(jié)果比例結(jié)果比較，上述結(jié)果比例2迭代迭代15次還好次還好. 若若迭代迭代15次次 1 1.7865258(相應(yīng)的相應(yīng)的 1 2.5365258)就是真就是真值值.上頁上頁下頁下頁 設(shè)設(shè)ARnn為為對稱矩陣對稱矩陣，稱，稱.),(),()(xxxAxxR 為向量為向量x的的瑞利商瑞利商，其中，其中(x, x)=xTx為內(nèi)積為內(nèi)積. 由定理由定理4知道，實對稱矩陣知道，實對稱矩陣A的特征值的特征值 1及及 n可用可用瑞利商瑞利商的的極限值表示極限值表示. 下面我們將下

36、面我們將瑞利商瑞利商應(yīng)用到用冪法計算應(yīng)用到用冪法計算實對稱矩陣實對稱矩陣A的主特征值的加速上來的主特征值的加速上來.2、瑞利商、瑞利商(Rayleigh)加速加速上頁上頁下頁下頁 定理定理9 設(shè)設(shè)ARnn為為對稱矩陣對稱矩陣，特征值滿足，特征值滿足對應(yīng)的特征向量對應(yīng)的特征向量xi滿足滿足(xi, xj)= ij (單位正交向量單位正交向量), 應(yīng)應(yīng)用冪法公式用冪法公式(2.9)計算計算A的主特征值的主特征值 1，則規(guī)范化向，則規(guī)范化向量量uk的的瑞利商瑞利商給出給出 1的較好的近似值為的較好的近似值為,121nn kkkkkkuuuAuuR2121, 由此可見，由此可見，R(uk)比比 k更快

37、的收斂于更快的收斂于 1.上頁上頁下頁下頁 證明證明 由由(2.8)式及式及,)max(,)max(00100uAuAAuvuAuAukkkkkkk )11. 2(.),(),(),(),()(2121122112200001 knjkjjnjkjjkkkkkkkkkaauAuAuAuAuuuAuuR 得得上頁上頁下頁下頁 冪法的冪法的瑞利商加速迭代公式瑞利商加速迭代公式可以寫為可以寫為 kkkkkkkkkkvukuuuvAuv /), 2 , 1(),(),(1111其中其中A為為n階實對稱矩陣階實對稱矩陣.,11kkux 對給定的誤差限對給定的誤差限 ，當(dāng)，當(dāng)| k k- -1| 時，取近

38、似值時，取近似值上頁上頁下頁下頁8.2.3 反冪法反冪法 反冪法是用于求非奇異矩陣反冪法是用于求非奇異矩陣A的的按模最小的特征按模最小的特征值和對應(yīng)特征向量值和對應(yīng)特征向量的方法的方法. 而結(jié)合原點平移法的反冪而結(jié)合原點平移法的反冪法則可以求矩陣法則可以求矩陣A的任何一個的任何一個給定近似特征值對應(yīng)的給定近似特征值對應(yīng)的特征向量特征向量.設(shè)矩陣設(shè)矩陣A非奇異非奇異,其特征值其特征值 i (i=1,2,n) ,滿足滿足0121 nn 其相應(yīng)的特征向量其相應(yīng)的特征向量x1,x2,xn線性無關(guān)，則線性無關(guān)，則 A- -1 的特征的特征值為值為1/ i ，對應(yīng)的特征向量仍為，對應(yīng)的特征向量仍為xi (

39、i=1,2,n).iiiiiixxAxAx11 上頁上頁下頁下頁此時此時, A- -1的特征值滿足的特征值滿足11111 nn因此因此, 對對A- -1應(yīng)用冪法應(yīng)用冪法,可求出其可求出其主特征值主特征值 1/ n k 和和特征向量特征向量 xn uk .從而求得從而求得A的的按模最小特征值按模最小特征值 n 1/ k 和對應(yīng)的和對應(yīng)的特征向量特征向量 xn uk ，這種求這種求A- -1的方法就稱為的方法就稱為反冪法反冪法.上頁上頁下頁下頁為了避免求為了避免求A- -1, 可通過解線性方程組可通過解線性方程組Avk=uk- -1得到得到vk，采用采用LU分解法，即先對分解法，即先對A進(jìn)行進(jìn)行L

40、U分解分解A=LU, 此時此時反冪反冪法的迭代公式法的迭代公式為為 , 2 , 1/max,1 kvuvvzUvzuLzkkkkkkkkkkk 求求出出解解求求出出解解), 2 , 1()max(11 k/vuvuAvkkkkkkk knknux ,1 反冪法的迭代公式反冪法的迭代公式為為:任取任取v0=u0 0 0，構(gòu)造，構(gòu)造上頁上頁下頁下頁對給定的誤差對給定的誤差 ，當(dāng)，當(dāng)|kk- -1| n|0，則對任意非零初始向量則對任意非零初始向量u0(an 0) ，由反冪法計算公，由反冪法計算公式構(gòu)造的向量序列式構(gòu)造的向量序列vk，uk滿足滿足 (1)lim,max()nkknxux (2).1)

41、max(limnkkv 收斂速度收斂速度的比值的比值r=| n/ n- -1|. 上頁上頁下頁下頁 在反冪法中也可以用在反冪法中也可以用原點平移法原點平移法加速迭代過程加速迭代過程，或或求其它特征值與其對應(yīng)的特征向量求其它特征值與其對應(yīng)的特征向量. 如果矩陣如果矩陣(A- -pI)- -1存在，顯然其特征值為存在，顯然其特征值為,1,1,121pppn 對應(yīng)的特征向量仍然是對應(yīng)的特征向量仍然是x1,x2,xn，現(xiàn)對矩陣，現(xiàn)對矩陣(A- -pI)- -1應(yīng)用冪法，得到反冪法的迭代公式應(yīng)用冪法，得到反冪法的迭代公式00110,(),(1,2,.).(2.12)/ max()kkkkkuvvApIu

42、kuvv 初初始始向向量量上頁上頁下頁下頁 如果如果p是是A的特征值的特征值 j的一個近似值，且設(shè)的一個近似值，且設(shè) j與其它與其它特征值是分離的，即特征值是分離的，即),(jippij 就是說就是說1/( j- -p)是矩陣是矩陣 (A- -pI)- -1的主特征值，可用反冪的主特征值，可用反冪法法(2.12)計算特征值計算特征值 j及特征向量及特征向量xj. 設(shè)設(shè)ARnn有有n個線性無關(guān)的特征向量個線性無關(guān)的特征向量x1,x2, xn,則則),0(10 jniiiaxau,)max()(0)1(0upIAupIAvkkk 上頁上頁下頁下頁,)max()(00upIAupIAukkk 其中其

43、中.)()(10 niikiikxpaupIA 同理可得下面的定理同理可得下面的定理.上頁上頁下頁下頁 定理定理11 設(shè)設(shè)ARnn有有n個線性無關(guān)的特征向量，個線性無關(guān)的特征向量， 矩陣矩陣A的特征值及對應(yīng)的特征向量分別記為的特征值及對應(yīng)的特征向量分別記為 i 及及xi (i=1,2,n)，而，而p為為 j的近似值，的近似值，(A- -pI)- -1存在，且存在，且 (1),)max(limjjkkxxu (2)jkjkkvppv )max(1,1)max(lim即即則對任意非零初始向量則對任意非零初始向量u0(aj 0) ，由反冪法計算公式，由反冪法計算公式(2.12)構(gòu)造的向量序列構(gòu)造的向

44、量序列vk，uk滿足滿足. )( |jippij . |min/ |pprijij 且收斂速度為且收斂速度為上頁上頁下頁下頁 由該定理知，對由該定理知，對A- -pI(其中其中p j)應(yīng)用反冪法，可應(yīng)用反冪法，可用來計算特征向量用來計算特征向量xj，只要選擇，只要選擇p是是 j的一個較好的近的一個較好的近似且特征值分離情況較好，一般似且特征值分離情況較好，一般r很小，常常只要迭很小，常常只要迭代一二次就可完成特征向量的計算代一二次就可完成特征向量的計算.反冪法迭代公式中的反冪法迭代公式中的vk是通過解方程組是通過解方程組.)(1 kkuvpIA求得的求得的, 為了節(jié)省工作量為了節(jié)省工作量, 可

45、以先將可以先將A- -pI進(jìn)行三角分解進(jìn)行三角分解.)(LUpIAP 于是求于是求vk相對于解兩個三角形方程組相對于解兩個三角形方程組.,1kkkkyUvPuLy 上頁上頁下頁下頁實驗表明實驗表明, 按下述方法選擇按下述方法選擇u0是較好的是較好的: 選選u0使使)13. 2()1 , 1 , 1(011 PuLUv用回代求解三角形方程組用回代求解三角形方程組(2.13)即得即得v1，然后再按公，然后再按公式式(2.12)進(jìn)行迭代進(jìn)行迭代.上頁上頁下頁下頁反冪法計算公式：反冪法計算公式：1.分解計算分解計算P(A- -pI)=LU,且保留且保留L, U及及P信息信息2.反冪迭代法反冪迭代法(1

46、) 解解Uv1=(1,1,1)T求求v1.),max(kkkkk/vuv (2) 解解k=2,3, 解解Lyk=Puk- -1求求yk 解解Uvk=yk求求vk k=max(vk) 計算計算uk=vk/ k上頁上頁下頁下頁例例5 用反冪法求用反冪法求 410131012A的對應(yīng)于計算特征值的對應(yīng)于計算特征值 =1.2679(精確特征值為精確特征值為 )的特征向量的特征向量(用用5位浮點數(shù)進(jìn)行計算位浮點數(shù)進(jìn)行計算).333 解解 用部分選主元的三角分解將用部分選主元的三角分解將A- -pI(其中其中p=1.2679)分解為分解為 P(A- -pI)=LU,其中其中,126807. 07321.

47、0010001 L上頁上頁下頁下頁,1029405. 0007321. 21017321. 113 U.001100010 P由由Uv1=(1,1,1)T,得得,)8 .3400, 3 .9290,12692(1Tv ,)26796. 0,73206. 0, 1(1Tu 由由LUv2=Pu1 ,得得,)4 .5467,14937,20404(2Tv ,)26796. 0,73206. 0, 1(2Tu 3對應(yīng)的特征向量是對應(yīng)的特征向量是,)26795. 0,73205. 0, 1()32,31, 1(3TTx 由此看出由此看出u2是是x3的相當(dāng)好的近似的相當(dāng)好的近似.特征值特征值 3的真值為的

48、真值為,26794901. 1/12679. 123 .26794912. 1333 上頁上頁下頁下頁8.3 正交變化與矩陣分解正交變化與矩陣分解正交變化是計算矩陣特征值的有力工具，本節(jié)介正交變化是計算矩陣特征值的有力工具，本節(jié)介紹豪斯紹豪斯霍爾德霍爾德(Householder)變換變換和和吉文斯吉文斯(Givens)變變換換，并利用它們討論，并利用它們討論矩陣分解矩陣分解. 主要討論主要討論實矩陣實矩陣和和實實向量向量.上頁上頁下頁下頁8.3.1 豪斯霍爾德變換豪斯霍爾德變換定義定義2 設(shè)向量設(shè)向量w Rn，且，且wTw=1，稱矩陣，稱矩陣)1 . 3(.212222122221)(2212

49、221212121 nnnnnwwwwwwwwwwwwwwwwHTwwIwH2)( 為為初等反射矩陣初等反射矩陣，這個矩陣也稱為，這個矩陣也稱為豪斯霍爾德變換豪斯霍爾德變換. 如果記如果記w=(w1,w2, ,wn)，則有，則有上頁上頁下頁下頁定理定理12 設(shè)有初等反射矩陣，其中設(shè)有初等反射矩陣，其中H=I- -2wwT，其，其中中wTw=1，則，則(1) H是對稱矩陣，即是對稱矩陣，即HT=H.(2) H是正交矩陣，即是正交矩陣，即H- -1=H.(3) 設(shè)設(shè)A為對稱矩陣，那么為對稱矩陣，那么 A1=H- -1AH=HAH 即亦即亦是對稱矩陣是對稱矩陣.證明證明 只證只證H的正交性，其中顯然

50、的正交性，其中顯然. HTH=H2=(I- -2wwT)(I- -2wwT)=I- -4wwT+4w(wTw)wT=I.設(shè)向量設(shè)向量u 0，則顯然，則顯然是一個初等反射矩陣是一個初等反射矩陣.222.TuuHIu上頁上頁下頁下頁vwSxv 圖圖 8-18-1Oy下面考察下面考察初等反射矩陣初等反射矩陣的的幾何意義幾何意義. 參見參見圖圖8-1，考，考慮以為慮以為w法向量且過原點法向量且過原點O的超平面的超平面S:wTx=0. 設(shè)任意向設(shè)任意向量量v Rn, 則則v=x+y, 其中其中x S, y ST. 于是于是 Hx=(I- -2wwT)x=x- -2wwTx=x.對于對于y ST，易知，易

51、知Hy=- -y, 從而對從而對任意向量任意向量v Rn，總有，總有 Hv=x- -y=v .其中其中v 為為v關(guān)于平面關(guān)于平面S的鏡面反射的鏡面反射(見圖見圖8-1).初等反射矩陣在計算上的意義是它能用來約化矩初等反射矩陣在計算上的意義是它能用來約化矩陣，例如設(shè)向量陣，例如設(shè)向量Hx 0，可選擇一初等反射矩陣，可選擇一初等反射矩陣H，使，使Hx=y.上頁上頁下頁下頁定理定理13 設(shè)設(shè)x,y為兩個不相等的為兩個不相等的n維向量維向量, |x|2=|y|2,則存在一個初等反射矩陣則存在一個初等反射矩陣H，使，使Hx=y.2222().TTTxyHIwwIxyxy 2xywxy 證明證明 令令 ,

52、 則得到一個初等反射矩陣則得到一個初等反射矩陣而且而且2222()()2()2.TTTTxyxyx xy xHxxxyxxxyxy 因為因為22() ()2().TTTxyxyxyx xy x 所以所以().Hxxxyy 容易說明，容易說明，w是使是使Hx=y成立的唯一長度等于成立的唯一長度等于1的的向量向量(不計符號不計符號).上頁上頁下頁下頁定理定理14(約化定理約化定理) 設(shè)設(shè)x=(x1,x2,xn)T 0, 則存在則存在初等反射矩陣初等反射矩陣H，使，使Hx=- -e1，其中，其中1121212,sgn(),(3.2),1().2THIuuxxuxeux 證明證明 記記y=- -e1，

53、設(shè)，設(shè)x y，取，取=|x|2，則有，則有|x|2=|y|2，于是由定理，于是由定理13存在存在H變換變換 H=I- -2wwT ，其中其中 ，使，使Hx=y=- -e1. .112xewxe 上頁上頁下頁下頁記記u=x+ +e1=(u1,u2,un)T，于是，于是,(0).xdxxdd 1222.TTuuHIIuuu 其中記其中記u=(x1+ +,x2,xn)T， . 顯然顯然2212u 如果如果和和x1異號，那么計算異號，那么計算x1+ +時有效數(shù)字可能時有效數(shù)字可能損失，我們?nèi)p失，我們?nèi)『秃蛒1有相同的符號，即取有相同的符號，即取1/221121sgn()sgn().niixxxx 在

54、計算在計算時，可能上溢或下溢，為了避免溢出，時，可能上溢或下溢，為了避免溢出，將將x規(guī)范化規(guī)范化上頁上頁下頁下頁則有則有H 使使 H x = e1，其中，其中221221(9, 5,1,1) ,108,54,2Tuxeuu 12(),/,/,/,.THIu ud uu ddHH 例例6 設(shè)設(shè)x=(3, 5, 1, 1)T，則，則|x|2=6. 取取=6,12745994529551,549553195153THIuu 可直接驗證可直接驗證Hx=- -e1=(- -6, 0, 0, 0)T.上頁上頁下頁下頁8.3.2 吉文斯變換吉文斯變換設(shè)向量設(shè)向量x,y R2，則變換，則變換1122cossi

55、n,sincosyxyx .yPx 或或是平面上向量的一個旋轉(zhuǎn)變換，其中是平面上向量的一個旋轉(zhuǎn)變換，其中cossin( ),sincosP 為正交變換為正交變換.Rn中變換中變換.yPx 其中其中x=(x1,x2,xn)T, y=(y1,y2,yn)T，而，而上頁上頁下頁下頁11cossin1( , , )(3.3)1sincos11iPP i jj 第第 行行第第 行行稱為稱為Rn中平面中平面xi, xj的的旋轉(zhuǎn)變換旋轉(zhuǎn)變換，也稱為，也稱為吉文斯變換吉文斯變換. P=P(i, j, )=P(i, j)稱為平面旋轉(zhuǎn)矩陣稱為平面旋轉(zhuǎn)矩陣.上頁上頁下頁下頁顯然，顯然，P=P(i, j, )具有性質(zhì)

56、：具有性質(zhì)：(1) P與單位陣與單位陣I只是在只是在(i, i), (i, j), (j, i), (j, j)位置元位置元素不一樣，其它相同素不一樣，其它相同.(2) P為正交矩陣為正交矩陣(P- -1=PT).(3) P(i, j)A(左乘左乘)只需計算第只需計算第i行與第行與第j行元素，即行元素，即對對A=(aij)mn有有cs,1,2, .scililjljlaalnaa 其中其中c=cos , s=sin .(4) AP(i, j)(右乘右乘)只需計算第只需計算第i列與第列與第j列元素列元素cs(,)(,),1,2,.scliljliljaaaalm 利用平面旋轉(zhuǎn)變換，可得向量利用平

57、面旋轉(zhuǎn)變換，可得向量x中的指定元素變?yōu)榱阒械闹付ㄔ刈優(yōu)榱?上頁上頁下頁下頁定理定理15(約化定理約化定理) 設(shè)設(shè)x=(x1,xi,xj,xn)T, 其其中不全為零，則可選擇平面旋轉(zhuǎn)陣中不全為零，則可選擇平面旋轉(zhuǎn)陣P=P(i, j, )，使，使1(,0,) ,TinPxxxxij 其中其中22,arctan(/).iijjixxxxx 于是，由于是，由c, s的取法得的取法得證明證明 取取c=cos =xi/x i , s=sin =xj/x i ，由，由P(i, j, )x =x =(x 1,x i,x j,x n)T, 利用矩陣乘法利用矩陣乘法, 顯然有顯然有, .iijjijkkxcxs

58、xxsxcxxxki j 22,0.iijjxxxx上頁上頁下頁下頁8.3.3 矩陣的矩陣的QR分解與舒爾分解分解與舒爾分解 定理定理16 設(shè)設(shè)ARnn非奇異，則存在正交矩陣非奇異，則存在正交矩陣P，使使PA=R，其中，其中R為上三角矩陣為上三角矩陣. 證明證明 我們先用吉文斯變換給出構(gòu)造我們先用吉文斯變換給出構(gòu)造P的方法的方法 (1) 第第1步約化，由設(shè)有步約化，由設(shè)有j(j=1,2, ,n)使使aj1 0，則，則可選擇吉文斯變換可選擇吉文斯變換P(1, j)，將，將aj1處的元素化為零處的元素化為零. 即即若若aj1 0 (j=2,3, ,n)，則存在，則存在P(1, j)使得使得可簡記為

59、可簡記為P1A=A(2)，其中，其中P1=P(1,n)P(1,2).)2 , 1(), 1()2()2()2(2)2(2)2(2211211AaaaarrrAPnPnnnnn 上頁上頁下頁下頁 (2) 第第k步約化，設(shè)上述過程已完成第步約化，設(shè)上述過程已完成第1步至第步至第k-1步，于是有步，于是有 由設(shè)有由設(shè)有j(n j k)使使ajk(k) 0 (j=k+1,n)，則可選擇，則可選擇吉文斯變換吉文斯變換P(k, j) (j=k+1,n)使使1112112222( )221( )( )( )( ).knknkkkkkkknkknknnrrrrrrrPP P AAaaaa PkA(k)=P(k

60、, n)P(k, k+1)A(k)=PkPk- -1P1A=A(k+1),其中其中 Pk=P(k, n)P(k, k+1) .上頁上頁下頁下頁 (3) 繼續(xù)上述約化過程，最后則有繼續(xù)上述約化過程，最后則有Pn- -1P2P1A=A(n)=R (上三角形矩陣上三角形矩陣). 令令P=Pn- -1P2P1, 它是一個正交矩陣它是一個正交矩陣, 有有PA=R. 也可以用豪斯霍爾德變換構(gòu)造出正交矩陣也可以用豪斯霍爾德變換構(gòu)造出正交矩陣P，記，記A(0)=A，它的第一列記為，它的第一列記為u1，不妨設(shè)，不妨設(shè)a1(0) 0，可按公，可按公式式(3.2)找到矩陣找到矩陣 使使111111,n nTHRHI