第8章矩陣特征值問題的數(shù)值解法

1、上頁上頁下頁下頁第第8章章 矩陣特征值問題的數(shù)值解法矩陣特征值問題的數(shù)值解法 8.1 引言引言 8.2 冪法及反冪法冪法及反冪法 8.3 雅可比方法雅可比方法 8.4 QR算法算法上頁上頁下頁下頁8.1 引引 言言 工程技術中有多種振動問題，如橋梁或建筑物的工程技術中有多種振動問題，如橋梁或建筑物的振動，機械零件、飛機機翼的振動，及一些穩(wěn)定性分振動，機械零件、飛機機翼的振動，及一些穩(wěn)定性分析和相關分析在數(shù)學上都可轉化為求矩陣特征值與特析和相關分析在數(shù)學上都可轉化為求矩陣特征值與特征向量的問題征向量的問題. 下面先復習一些矩陣的特征值和特征向量的基礎下面先復習一些矩陣的特征值和特征向量的基礎知識

2、知識.上頁上頁下頁下頁 定義定義1 已知已知n階矩陣階矩陣A=(aij)，則，則)2()(det)det()(12211212222111211的的項項次次數(shù)數(shù) naaaaaaaaaaaaAInnnnnnnnnn 稱為稱為A的的特征多項式特征多項式. 一般有一般有n個根個根(實的或復的，復根按重數(shù)計算實的或復的，復根按重數(shù)計算)稱為稱為A的的特征值特征值. 用用(A)表示表示A的所有特征值的集合的所有特征值的集合. A的特征方程的特征方程)1 . 1(0)det()( AI 上頁上頁下頁下頁 設設為為A的特征值，相應的齊次方程組的特征值，相應的齊次方程組 注：注：當當A為實矩陣時，為實矩陣時，

3、 ()=0為實系數(shù)為實系數(shù)n次代數(shù)次代數(shù)方程，其復根是共軛成對出現(xiàn)方程，其復根是共軛成對出現(xiàn).)2 . 1(0)( xAI 的的非零解非零解x稱為矩陣稱為矩陣A的對應于的對應于的的特征向量特征向量.上頁上頁下頁下頁 定理定理1 設設為為ARnn的特征值的特征值, 且且Ax=x (x 0)，則有則有 - -p為為A- -pI的特征值，即的特征值，即(A- -pI)x=(- -p)x ； c為的為的cA特征值特征值(c0為常數(shù)為常數(shù))； 下面敘述有關特征值的一些下面敘述有關特征值的一些結論結論： k為為Ak的特征值，即的特征值，即Akx=kx ； 設設A為非奇異矩陣，那么為非奇異矩陣，那么0 ,

4、且且- -1為為A- -1的特的特征值，即征值，即A- -1x=- -1x .上頁上頁下頁下頁 定理定理2 設設i(i=1,2,n)為為n階矩陣階矩陣A=(aij)的特征值，的特征值，則有則有)(Atraniiinii 11 稱為稱為A的的跡跡； .nA21 定理定理3 設設ARnn，則有，則有. )()(AAT 上頁上頁下頁下頁 定義定義2 設設n階矩陣階矩陣A=(aij)，令，令 下面討論下面討論矩陣特征值界的估計矩陣特征值界的估計.).,(niarnijjiji211 ； 集合集合 稱稱為復平面上以為復平面上以aii為圓心，以為圓心，以ri為半徑的為半徑的n階矩陣階矩陣A的的n個個Ger

5、schgorin圓盤圓盤. ),(,|niCzrazzDiiii21 上頁上頁下頁下頁 定理定理4 (Gerschgorin圓盤定理圓盤定理) 特別地，如果特別地，如果A的一個圓盤的一個圓盤Di是與其它圓盤分離是與其它圓盤分離(即即孤立圓盤孤立圓盤)，則，則Di中精確地包含中精確地包含A的一個特征值的一個特征值.).,(niaranijjijiii211 設設n階矩陣階矩陣A(aij)，則，則A的每一個特征值必屬的每一個特征值必屬于下面某個圓盤之中于下面某個圓盤之中 如果如果A有有m個個圓盤圓盤組成一個連通的并集組成一個連通的并集S，且，且S與余下與余下n- -m個個圓盤圓盤是分離的，則是分離

6、的，則S內(nèi)恰包含內(nèi)恰包含A的的m個個特征值特征值.或者說或者說 A的特征值都在的特征值都在n個圓盤的個圓盤的并集并集中中.上頁上頁下頁下頁 證明證明 只就給出證明只就給出證明. 設設為為A的特征值，即的特征值，即.1knkjjkjkkraa Ax=x，其中，其中x=(x1,x2, xn)T 0.或或 記記 ，考慮，考慮Ax=x的第的第k個個方程，即方程，即0max1 xxxinik,1knjjkjxxa ,)( kjjkjkkkxaxa 于是于是即即, kjkjkkjjkjkkkaxxaxa 上頁上頁下頁下頁這說明，這說明，A的每一個特征值必位于的每一個特征值必位于A的一個圓盤的一個圓盤中，并

7、且相應的特征值中，并且相應的特征值一定位于第一定位于第k個圓盤中個圓盤中(其中其中k是對應特征向量是對應特征向量x絕對值最大的分量的下標絕對值最大的分量的下標).上頁上頁下頁下頁.411101014 A 例例2 估計矩陣估計矩陣A的特征值范圍，其中的特征值范圍，其中 解解 矩陣矩陣A的的3個圓盤為個圓盤為. 24:, 2:, 14:321 DDD 由定理由定理8，可知，可知A的的3個特征值位于個特征值位于3個圓盤的并個圓盤的并集中，由于集中，由于D1是孤立圓盤，所以是孤立圓盤，所以D1內(nèi)恰好包含內(nèi)恰好包含A的的一個特征值一個特征值1(為實特征值為實特征值)，即，即. 531 A的其它兩個特征值

8、的其它兩個特征值2, 3包含在包含在D2, D3的并集中的并集中.上頁上頁下頁下頁 定義定義3 設設ARnn為對稱矩陣，對于任一非零為對稱矩陣，對于任一非零向量向量x，稱，稱,),(),()(xxxAxxR 為對應于向量為對應于向量x的的瑞利瑞利(Rayleigh)商商. 定理定理5 設設ARnn為對稱矩陣為對稱矩陣(其特征值次序記其特征值次序記為為12n)，則，則 1. (對任何非零對任何非零xRn)；1 ),(),(xxxAxn 2. ；),(),(maxxxxAxxRxn01 3. .),(),(minxxxAxxRxnn0 上頁上頁下頁下頁 證明證明 只證只證1，關于，關于2, 3自己

9、作練習自己作練習. 由于由于A為實對稱矩陣，可將為實對稱矩陣，可將 1,2,n 對應的特對應的特征向量征向量 x1,x2,xn 正交規(guī)范化，則有正交規(guī)范化，則有 (xi, xj)=ij，設，設x 0為為Rn中任一向量，則有中任一向量，則有. 0,1221 niiniiixxx 于是于是.),(),(11212 niiniiinxxxAx從而從而1成立成立. 結論結論1說明說明瑞利商瑞利商必位于必位于n和和1之間之間.上頁上頁下頁下頁 關于計算矩陣關于計算矩陣A的特征值問題，當?shù)奶卣髦祮栴}，當n2,3時，我時，我們還可按行列式展開的辦法求們還可按行列式展開的辦法求 ()=0的根的根. 但當?shù)攏

10、較較大時，如果按展開行列式的辦法，首先求出大時，如果按展開行列式的辦法，首先求出 ()的系的系數(shù)，再求數(shù)，再求 ()的根，工作量就非常大，用這種辦法求的根，工作量就非常大，用這種辦法求矩陣的特征值是不切實際的，由此需要研究矩陣的特征值是不切實際的，由此需要研究求求A的特的特征值及特征向量的數(shù)值解法征值及特征向量的數(shù)值解法. 本章將介紹一些計算機上常用的本章將介紹一些計算機上常用的兩類方法兩類方法，一，一類是類是冪法及反冪法冪法及反冪法（迭代法），另一類是（迭代法），另一類是正交相似正交相似變換的方法變換的方法（變換法）（變換法）.上頁上頁下頁下頁 冪法與反冪法都是求冪法與反冪法都是求實矩陣實矩

11、陣的特征值和特征的特征值和特征向量的向量的向量迭代法向量迭代法，所不同的是冪法是計算矩陣，所不同的是冪法是計算矩陣的的主特征值主特征值(矩陣矩陣按模最大的特征值按模最大的特征值稱為稱為主特征值主特征值，其模就是該矩陣的其模就是該矩陣的譜半徑譜半徑)和和相應特征向量相應特征向量的一種的一種向量迭代法，而反冪法則是計算非奇異向量迭代法，而反冪法則是計算非奇異(可逆可逆)矩陣矩陣按模最小的特征值按模最小的特征值和和相應特征向量相應特征向量的一種向量迭的一種向量迭代法代法. 下面分別介紹冪法與反冪法下面分別介紹冪法與反冪法.8.2 冪法及反冪法冪法及反冪法上頁上頁下頁下頁現(xiàn)討論求現(xiàn)討論求1及及x1的方

12、法的方法.), 2 , 1(nixAxiii 設實矩陣設實矩陣A=(aij)有一個有一個完全的特征向量組完全的特征向量組，即，即A有有n個線性無關的特征向量個線性無關的特征向量，設矩陣，設矩陣A的特征值為的特征值為1,2,n, 相應的特征向量為相應的特征向量為x1,x2,xn. 已知已知A的主的主特征值特征值1是實根，且滿足條件是實根，且滿足條件)1 . 2(|,|21n 8.2.1 冪法（又稱冪法（又稱乘冪法乘冪法）上頁上頁下頁下頁)3 . 2(),0(122110 axaxaxavnn設設 冪法的冪法的基本思想基本思想是是: 任取任取非零非零的初始向量的初始向量v0 , 由矩由矩陣陣A構造

13、一向量序列構造一向量序列vk稱為迭代向量，由假設，稱為迭代向量，由假設，v0可唯一表示為可唯一表示為)2 . 2(.,.,011021201 vAAvvvAAvvAvvkkk上頁上頁下頁下頁)3 . 2(),0(122110 axaxaxavnn設設于是于是).()/(1112111122211101kkniikiiknknnkkkkkxaxaxaxaxaxavAAvv 其中其中.)/(21 niikiikxa 由假設由假設 故故 從而從而), 3 , 2( 1/1nii , 0lim kk )4 . 2(.lim111xavkkk 為為1的特征向量的特征向量.上頁上頁下頁下頁所以當所以當k充

14、分大時，有充分大時，有)5 . 2(,111xavkk 即為矩陣即為矩陣A的的對應特征值對應特征值 1 的一個近似特征向量的一個近似特征向量.用用(vk)i 表示表示vk的第的第i個分量，則當個分量，則當k充分大時，有充分大時，有 )7 . 2(.11 ikikvv即為即為A的的主特征值主特征值 1的近似值的近似值.)6 . 2(,111111kkkkvxaAvv 由于由于上頁上頁下頁下頁 迭代公式實質上是由矩陣迭代公式實質上是由矩陣A的乘冪的乘冪 Ak與非零向與非零向量量v0相乘來構造向量序列相乘來構造向量序列vk=Akv0，從而計算主特，從而計算主特征值征值1及其對應的特征向量，這就是及其

15、對應的特征向量，這就是冪法冪法的思想的思想. ).(11 kvvikik 的收斂速度由比值的收斂速度由比值,12 r來確定，來確定，r越小收斂越快，但當越小收斂越快，但當r1 1時收斂可能很慢時收斂可能很慢.上頁上頁下頁下頁 定理定理6 設設ARnn有有n個線性無關的特征向量，個線性無關的特征向量，主特征值主特征值1滿足條件滿足條件 |1|2|n|，則對任何非零向量則對任何非零向量v0(a1 0)，冪法的算式成立，冪法的算式成立.又設又設A有有n個線性無關的特征向量，個線性無關的特征向量，1對應的對應的r個線性個線性無關的特征向量為無關的特征向量為x1,x2,xr，則由，則由(2.2)式有式有

16、 如果如果A的主特征值為實的重根的主特征值為實的重根, 即即1=2=r, 且且 |r|r+1|n|，,)/(11110 nriikiiriiikkkxaxavAv 上頁上頁下頁下頁).0(lim111 riiiriiikkkxaxav設設 為為A的特征向量，這說明當?shù)奶卣飨蛄?，這說明當A的主特征值是實的重根的主特征值是實的重根時，定理時，定理6的結論還是正確的的結論還是正確的. 應用應用冪法冪法計算計算A的主特征值的主特征值1及其對應的特征向及其對應的特征向量時，如果量時，如果|1|1(或或|1|2|n|，則對任意非零初始，則對任意非零初始向量向量v0=u0(a1 0)，有冪法計算公式為，有冪

17、法計算公式為則有則有 ,)max(lim11xxukk .lim1 kk上頁上頁下頁下頁例例1 用冪法計算矩陣用冪法計算矩陣 1634310232A的主特征值與其對應的特征向量的主特征值與其對應的特征向量.解解取取 v0=u0=(0,0,1)T , 則則 TTvuv25. 0 , 1, 5 . 01, 4,1 , 4 , 211111 ), 2 , 1(max1 k/vuvAuvkkkkkkk 上頁上頁下頁下頁直到直到k=8 時的計算結果見下表時的計算結果見下表k1 2, 4, 1, 4 0.5, 1, 0.252 4.5, 9, 7.75 90.5, 1, 0.86113 5.7222, 1

18、1.4444, 8.36111.44440.5, 1, 0.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,00

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

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

21、上頁下頁下頁例例4 設設AR44有特征值有特征值),4 , 3 , 2 , 1(15 jji 比值比值r=|2/1|0.9. 做變換做變換B=A- -12I (p=12),則則B的特征值為的特征值為. 1, 0, 1, 24321 應用冪法計算應用冪法計算B的主特征值的主特征值1的收斂速度的比值為的收斂速度的比值為. 9 . 021121212 pp 雖然常常能夠選擇有利的雖然常常能夠選擇有利的p值值, 使冪法得到加速使冪法得到加速, 但設計一個自動選擇適當參數(shù)但設計一個自動選擇適當參數(shù)p的過程是困難的的過程是困難的.上頁上頁下頁下頁 定理定理 設設ARnn為為對稱矩陣對稱矩陣，特征值滿足，特

22、征值滿足對應的特征向量對應的特征向量xi滿足滿足(xi, xj)=ij (單位正交向量單位正交向量) ，應用冪法公式應用冪法公式(2.9)計算計算A的主特征值的主特征值 1，則規(guī)范化，則規(guī)范化向量向量uk的的瑞利商瑞利商給出給出 1的較好的近似值為的較好的近似值為,121nn kkkkkkuuuAuuR2121, 由此可見，由此可見，R(uk)比比k更快的收斂于更快的收斂于 1.2、瑞利商加速、瑞利商加速上頁上頁下頁下頁 證明證明 由由(2.8)式及式及,)max(,)max(00100uAuAAuvuAuAukkkkkkk )11. 2(.),(),(),(),()(212112211220

23、0001 knjkjjnjkjjkkkkkkkkkaauAuAuAuAuuuAuuR 得得上頁上頁下頁下頁 冪法的冪法的瑞利商加速迭代公式瑞利商加速迭代公式可以寫為可以寫為 kkkkkkkkkkvukuuuvAuv /), 2 , 1(),(),(1111其中其中A為為n階實對稱矩陣階實對稱矩陣.,11kkux 對給定的誤差限對給定的誤差限 ，當，當| kk- -1| 時，取近似值時，取近似值上頁上頁下頁下頁反冪法反冪法 反冪法是用于求非奇異矩陣反冪法是用于求非奇異矩陣A的的按模最小的特征按模最小的特征值和對應特征向量值和對應特征向量的方法的方法. 而結合原點平移法的反冪而結合原點平移法的反冪

24、法則可以求矩陣法則可以求矩陣A的任何一個的任何一個具有先了解的特征值和具有先了解的特征值和對應的特征向量對應的特征向量。設矩陣設矩陣A非奇異非奇異,其特征值其特征值 i (i=1,2,n) ,滿足滿足0121 nn 其相應的特征向量其相應的特征向量x1,x2,xn線性無關，則線性無關，則 A- -1 的特征的特征值為值為1/ i ，對應的特征向量仍為，對應的特征向量仍為xi (i=1,2,n).iiiiiixxAxAx11 上頁上頁下頁下頁此時此時, A- -1的特征值滿足的特征值滿足11111 nn因此因此, 對對A- -1應用冪法應用冪法,可求出其可求出其主特征值主特征值1/ n k 和和

25、特征向量特征向量 xn uk .從而求得從而求得A的的按模最小按模最小特征值特征值 n 1/k 和對應的和對應的特征向量特征向量 xn uk ，這種求這種求A- -1的方法就稱的方法就稱為為反冪法反冪法.上頁上頁下頁下頁為了避免求為了避免求A- -1, 可通過解線性方程組可通過解線性方程組Avk=uk- -1得到得到vk，采用采用LU分解法，即先對分解法，即先對A進行進行LU分解分解A=LU, 此時此時反冪反冪法的迭代公式法的迭代公式為為 , 2 , 1/max,1 kvuvvzUvzuLzkkkkkkkkkkk 求求出出解解求求出出解解 ), 2 , 1(max11 k/vuvuAvkkkk

26、kkk knknux ,1 反冪法的迭代公式反冪法的迭代公式為為上頁上頁下頁下頁對給定的誤差對給定的誤差 ，當，當|kk- -1| n|0，則對任意非零初始向量則對任意非零初始向量u0(an 0) ，由反冪法計算公，由反冪法計算公式構造的向量序列式構造的向量序列vk，uk滿足滿足 ,)max(limnnkkxxu .1)max(limnkkv 上頁上頁下頁下頁 在反冪法中也可以用在反冪法中也可以用原點平移法原點平移法加速迭代過程加速迭代過程，或或求其它特征值與其對應的特征向量求其它特征值與其對應的特征向量. 如果矩陣如果矩陣(A- -pI)- -1存在，顯然其特征值為存在，顯然其特征值為,1,

27、1,121pppn 對應的特征向量仍然是對應的特征向量仍然是x1,x2,xn. 如果如果p是是A的特征值的特征值 j的一個近似值，且設的一個近似值，且設 j與其它與其它特征值是分離的，即特征值是分離的，即),(jippij 就是說就是說1/( j- -p)是矩陣是矩陣 (A- -pI)- -1的主特征值，可用反冪的主特征值，可用反冪法應用于矩陣法應用于矩陣(A- -pI)- -1計算計算 j及對應的特征向量及對應的特征向量.上頁上頁下頁下頁現(xiàn)對矩陣現(xiàn)對矩陣(A- -pI)- -1應用反冪法，得到如下迭代公式應用反冪法，得到如下迭代公式 00110,(),(1,2,).(2.12)max/ ma

28、x()kkkkkkkuvvApIukvuv 初初始始向向量量上頁上頁下頁下頁 定理定理9 設設ARnn有有n個線性無關的特征向量，個線性無關的特征向量， 矩陣矩陣A的特征值及對應的特征向量分別記為的特征值及對應的特征向量分別記為 i 及及xi (i=1,2,n)，而，而p為為 j的近似值，的近似值，(A- -pI)- -1存在，且存在，且 ,)max(limjjkkxxu jkjkkvppv )max(1,1)max(lim即即則對任意非零初始向量則對任意非零初始向量u0(aj 0) ，由反冪法計算公式，由反冪法計算公式(2.12)構造的向量序列構造的向量序列vk，uk滿足滿足. )( |ji

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

30、行三角分解.)(LUpIAP 于是求于是求vk相對于解兩個三角形方程組相對于解兩個三角形方程組.,1kkkkyUvPuLy 上頁上頁下頁下頁實驗表明實驗表明, 按下述方法選擇按下述方法選擇u0是較好的是較好的: 選選u0使使)13. 2()1 , 1 , 1(011 PuLUv用回代求解用回代求解(2.13)即得即得v1,然后再按公式然后再按公式(2.12)進行迭代進行迭代.上頁上頁下頁下頁例例5 求矩陣求矩陣A最接近于最接近于1.9的特征值和相應的特征向量的特征值和相應的特征向量 3 1- 2-1- 4 3 2- 3 7 A取取Ty) 1 , 1 , 1 (0作迭代，結果如表：作迭代，結果如

31、表：上頁上頁下頁下頁上頁上頁下頁下頁1222,1,(1) , diag(,) (1,2, ),2(),(), Tniijn nnTijn nijiji ji jAUU AUDinAUAaUBU AUbab 任任意意實實對對稱稱矩矩陣陣 可可通通過過正正交交相相似似變變換換化化成成對對角角陣陣 即即存存在在 正正交交矩矩陣陣使使得得其其中中是是 的的特特征征值值中中各各列列即即為為相相應應的的特特征征向向量量。（ ）在在正正交交相相似似變變換換下下，矩矩陣陣元元素素的的平平方方和和不不變變。設設為為正正交交矩矩陣陣，記記則則理理論論基基礎礎：1n Jacobi方方法法用用來來求求實實對對稱稱矩矩

32、陣陣的的全全部部特特征征值值及及相相應應特特征征向向量量。8.3 Jacobi方法方法上頁上頁下頁下頁,A通通過過一一次次正正交交變變換換 將將 中中一一對對非非零零的的非非對對角角元元素素化化成成零零 并并且且使使得得非非對對角角元元素素的的平平方方和和減減少少。反反復復進進行行上上述述過過程程，使使變變換換后后的的矩矩陣陣的的非非對對角角元元素素的的平平方方和和趨趨于于零零，從從而而使使該該矩矩陣陣近近似似為為對對角角矩矩陣陣，得得到到全全部部特特征征值值和和特特征征向向量量。Jacobi 方法的基本思路方法的基本思路上頁上頁下頁下頁矩陣的旋轉變換矩陣的旋轉變換cos ,1sin ,pqp

33、pqqpqqppqUuuunuU 其其中中，的的主主對對角角線線元元素素中中其其余余為為 ；而而其其非非主主對對角角線線元元素素中中維維空空間間中中的的二二維維坐坐其其余余為為標標旋旋0 0。稱稱為為轉轉矩矩陣陣。( )pqU 坐坐標標旋旋轉轉矩矩陣陣是是正正交交矩矩陣陣. .11cossin1( )1sincos11pqU 上頁上頁下頁下頁(1)(1)(1)22(1)22(1)(1)0, ()cossinsin2 sincossin2 cossin (, ) TpqqppqpqijppppqqpqqqppqqpqpiippiqiAaaAU AUaaaaaaaaaaaaaip qa 設設 為為

34、實實對對稱稱矩矩陣陣，且且若若記記則則有有(1)(1)(1)(1)(1)(1)sincos ( , )1()sin2cos2 2()cot2 2qiiqpiqiijjiijpqqpqqpppqppqqpqaaaaaai jp qaaaaaaaa 如如果果取取 使使得得(1)(1)(/4)0,.pqqppqqpaaAaa 則則有有得得到到一一個個使使 中中非非零零的的非非對對角角元元素素變變成成零零的的正正交交相相似似變變換換上頁上頁下頁下頁 (1)(2)( )(1)2(1)(1)2(1)(1)(1)(1)(1)(1), ( ),()()( , ),cossin ,kijijijijijjiij

35、piippiqiqiiqAAAAAE AaE Aaaaai jp q aaaaaa 對對重重復復上上述述過過程程得得到到一一個個矩矩陣陣序序列列?？煽勺C證，雖雖然然這這種種變變換換不不一一定定能能使使矩矩陣陣中中非非對對角角元元中中零零元元素素的的個個數(shù)數(shù)單單調調增增加加，但但可可保保證證非非對對角角元元的的平平方方和和遞遞減減。以以 與與為為例例：設設，則則由由(1)(1)(1)(1)2(1)2(1)2(1)2(1)2,2222,sincos ,(, )=0()()()2()() +2() 2()( )2( ) piqipqqpijijpiqipqijijip qi jp qijpiqipq

36、ijip qi jp qaaip qaaE AaaaaaaaaE AaE A ，上上式式表表明明，在在上上述述旋旋轉轉變變換換下下，非非對對角角元元的的平平方方和和嚴嚴格格單單調調遞遞減減，因因而而對對角角元元的的平平方方和和單單調調遞遞增增。上頁上頁下頁下頁 ( )(0)2( )( )Jacobi, lim ()=lim0,Jacobikkkijkki jAnAAAAE Aa 設設 為為 階階實實對對稱稱陣陣，對對 用用法法得得到到序序列列其其中中則則即即法法收收斂斂。Jacobi 迭代法收斂定理迭代法收斂定理 上頁上頁下頁下頁說說明明o1 Jacobi;定定理理表表明明，法法是是收收斂斂的

37、的Jacobi法法是是求求中中小小型型稠稠密密實實對對稱稱矩矩陣陣的的全全部部特特征征值值與與特特征征向向量量的的較較好好方方法法。o2 JacobiAnAn當當 的的階階數(shù)數(shù) 不不太太高高時時，算算法法的的收收斂斂速速度度很很快快；但但當當 的的 階階數(shù)數(shù) 變變得得較較大大時時，其其收收斂斂速速度度將將會會變變慢慢，即即法法 為為適適合合計計算算中中等等規(guī)規(guī)模模的的實實對對稱稱矩矩陣陣的的特特征征值值問問題題；o3 對對中中等等規(guī)規(guī)模模問問題題，具具有有較較好好的的數(shù)數(shù)值值穩(wěn)穩(wěn)定定性性；求求得得的的結結果果 的的精精度度也也很很高高，得得到到的的特特征征向向量量正正交交性性很很好好；o4 不

38、不足足之之處處：運運算算量量大大，不不能能保保持持矩矩陣陣的的特特殊殊形形狀狀 （如如稀稀疏疏性性）。上頁上頁下頁下頁 00(0)12T(1)(0)(0)(0)0:1,2,cot20,4cos0.7071068,sin0.70710680.70710680.70710680 ( )0.70710680.70710680 ,001100.7071068030.70710680.70710680.70710682kpqUUAUA U 210121012AA 用用J Ja ac co ob bi i方方法法求求 的的特特征征值值。解：解： 例例6： 上頁上頁下頁下頁11(1)230.70710678

39、0.47765830.88807380.45970081 0 00 0.888071:2,3,cot2,cos,sin 38 -0.45970080 0.4597008 ( ) 0.88807kpqUU T(2)(1)(1)(1)38 3.0000000 0.3250576 0.6279630 0.3250576 0.6339746 -0.0000000 0.6279630 0 2.3,U660254UAA 上頁上頁下頁下頁22(2)230.50478660.55166340.85165390.52410460.8516539 0 -0.5241046 0 1 2:1,3,cot2,cos,s

40、in 00.5241046 0 0.8516( )53 9kpqUU T(3)(2)(2)(2) 3.3864461 0.2768366 -0.0000000 0.2768366 0.6339746 -0.1703642-0.0000000 -0.1703642 1.9795,UU793AA 上頁上頁下頁下頁66(6)236:1,2,cot2-693.88568,-0.00072057929,cos0.9999997,sin0.0007205792 0.9999997 0.0007205792 0 ( ) -0.0007205792 0.9999997 0 0 kpqUU T(7)(6)(6)

41、(6), 0 1 3.4142136 0 0.0000000UU -0.0000000 0.5857864 -0.0000242 0.0000000 -0.0000242 1.9999999(AAE (7)1.1668492e-009A 上頁上頁下頁下頁123123(0)(1)(6)3.4142136,0.5857864,1.9999999;223.4142136,220.5857864,20.0000001 0.5000000 0.5000121 -0.7070982-0.7071068 0.7071068 0.0000121 0.5000000 QUUU所所以以，準準確確特特征征值值最最大

42、大誤誤差差不不超超過過；又又 0.4999879 0.7071153 0.50000000.5000121 -0.7070982-0.70710680.7071068 0.0000121 0.50000000.4999879 0.7071153 所所以以，對對應應的的特特征征向向量量分分別別為為，上頁上頁下頁下頁8.4 QR算法算法8.4.1. 化矩陣為化矩陣為Hessenberg形形8.4.2 QR 算法及其收斂性算法及其收斂性（1）鏡面反射矩陣）鏡面反射矩陣 （2）化一般實矩陣為上）化一般實矩陣為上Hessenberg矩陣矩陣 （3）化對稱矩陣為三對角矩陣）化對稱矩陣為三對角矩陣上頁上頁下

43、頁下頁 當當A為實矩陣，如果限制用正交相似變換，由于為實矩陣，如果限制用正交相似變換，由于A有復的特征值，有復的特征值， A不能用正交相似變換約化為上三不能用正交相似變換約化為上三角陣角陣. 用正交相似變換能約化到什么程度呢？用正交相似變換能約化到什么程度呢？(Schur定理定理) 設設ARnn，則存在酉矩陣，則存在酉矩陣U使使其中其中rii(i=1,2,n)為為 A的特征值的特征值. 下面給出理論上有關通過酉相似變換及正交變下面給出理論上有關通過酉相似變換及正交變換可以約化一般矩陣換可以約化一般矩陣A到什么程度的問題到什么程度的問題.),(22211211上三角陣上三角陣RrrrrrrAUU

44、nnnnH 上頁上頁下頁下頁其中其中Rii(i=1,2,m)為一階或二階方陣，且每個一階為一階或二階方陣，且每個一階Rii是是A的實特征值，每個二階對角塊的實特征值，每個二階對角塊Rii的兩個特征值的兩個特征值是是 A的兩個共軛復特征值的兩個共軛復特征值. (實實Schur分解分解) 設設ARnn，則存在正交矩陣，則存在正交矩陣Q使使,22211211 mmmmTRRRRRRAQQ上頁上頁下頁下頁63（1）鏡面反射矩陣）鏡面反射矩陣?；蚍Q為或定義：矩陣r r變變換換H Ho ou us se eh ho ol ld de e鏡鏡面面反反射射變變換換,RuuuuuEHuuuuEHnTTT),2(

45、222。，使，變換存在、為對稱正交變換；即、：TnTexeHxH,xRxHHHH)0 , 0 , 1 (, rHouseholde0,2 ,1 1211性性質質8.4.1. 化矩陣為化矩陣為Hessenberg形形鏡面反射矩陣的意義是鏡面反射矩陣的意義是“成批成批”消去向量的非零元素消去向量的非零元素。上頁上頁下頁下頁.,11221yxeyxT滿足變成求變換矩陣將, 1121623. 5,1623. 3)(, 1121121TTexuxxsignx9387.00613.01225.03162.00613.09387.01225.03162.01225.01225.0755.06325.0316

46、2.03162.06325.06325.02uuuuEHTT例：例：解解:THxy0001623. 3#上頁上頁下頁下頁 (1)用用初等反射矩陣作正交相似變換初等反射矩陣作正交相似變換約化一般約化一般實矩陣實矩陣A為為上上Hessenberg矩陣矩陣.化矩陣為上化矩陣為上Hessenberg矩陣矩陣討論討論兩個問題兩個問題 (2)用用初等反射矩陣作正交相似變換初等反射矩陣作正交相似變換約化對稱約化對稱矩陣矩陣A為為對稱三對角矩陣對稱三對角矩陣. 于是，于是，求原矩陣特征值問題求原矩陣特征值問題，就，就轉化為轉化為求上求上Hessenberg矩陣矩陣或或對稱三對角矩陣的特征值對稱三對角矩陣的特征

47、值問題問題.上頁上頁下頁下頁（2）化一般實矩陣為上）化一般實矩陣為上Hessenberg矩陣矩陣 設設ARnn，下面來說明，可選擇初等反射矩，下面來說明，可選擇初等反射矩陣陣U1,U2,Un- -2使使A經(jīng)正交相似變換約化為一個上經(jīng)正交相似變換約化為一個上Hessenberg陣陣.(1) 設設,)1(221)1(1211212222111211 AcAaaaaaaaaaaAnnnnnn上頁上頁下頁下頁21/212112111 111121sgn()(),(3.1)().niiaaucea 其中其中c1=(a21,an1)TRn- -1 ,不妨設不妨設c10，否則這一步不，否則這一步不需要約化需

48、要約化. 于是于是, 可選擇初等反射陣可選擇初等反射陣11111TRIu u 1 11 1R ce 令令,111 RU使使 其中其中上頁上頁下頁下頁則則 1)1(221111)1(12111112RARcRRAaUAUA(2)(2)(2)1112131(2)(2)(2)122232(2)(2)1112(2)(2)(2)32333(2)122(2)(2)(2)23,000nnnnnnnaaaaaaaAAaaacAaaa 其中其中.,),()2()2()2(222)2(2)2(322 nnnTnRARaac上頁上頁下頁下頁(2) 第第k步約化：重復上述過程，設對步約化：重復上述過程，設對A已完成第

49、已完成第1步步,第第k- -1步正交相似變換，即有步正交相似變換，即有,111 kkkkUAUA或或,11111 kkkUUAUUA且且 )()(1,)()(, 1)(1, 1)(, 1)()(1,)(1)(2)(1,2)(2)1(1,2)2(221)(1)(1, 1)(1)1(1, 1)2(12)1(11knnkknknkknkkkkkkkkknkkkkkkkknkkkkkkknkkkkkkkaaaaaaaaaaaaaaaaaaaaA 上頁上頁下頁下頁,0)(22)(12)(11knkAcAAknkkkkk 其中其中 為為k階上階上Hessenberg陣，陣，設設ck0, 于是可選擇初等反射

50、陣于是可選擇初等反射陣Rk使使 其中，其中，Rk計算公式為計算公式為1,kkkR ce )(11)()(, 1,),(kknTknkkkkkARaac .)()()(22knknkRA 上頁上頁下頁下頁( )( )2 1/21,11( )1,1sgn()() ),(3.2)(),.nkkkkkiki kkkkkkkkkkTkkkkaauceaRIu u 令令, kkRIU則則)3 . 3(00)1(221)1(12)1(11)(22)(12)1(111 kkkkkkkkkkkkkkkkAcAARARcRRAAUAUA上頁上頁下頁下頁221122 nnUUAUUUU其中其中 為為k+1階上階上H

51、essenberg陣，第陣，第k步約化只需步約化只需計算計算 及及 (當當A為對稱矩陣時，只需要為對稱矩陣時，只需要計算計算 ).)1(11 kAkkRA)(12kkkRAR)(22kkkRAR)(22(3) 重復上述過程，則有重復上述過程，則有.1)(1)1(1, 12)3(332)2(22111 nnnnnnnnnAaaaaa 上頁上頁下頁下頁 定理定理11 (Household約化矩陣為上約化矩陣為上Hessenberg陣陣) 設設ARnn則存在初等反射矩陣則存在初等反射矩陣U1,U2,Un- -2 使使.00221122HAUUUUAUUUUTnn 為為上上Hessenberg矩陣矩陣

52、.總結上述結論，有總結上述結論，有上頁上頁下頁下頁例例7 用用Household方法方法將矩陣將矩陣 7242327341AA約化為約化為上上Hessenberg陣陣. 解解 選取初等反射陣選取初等反射陣R1使使 , 其中其中c1=(2,4)T.1111ecR (1) 計算計算R1:)() 1 , 5 . 0(, 4) 4 , 2max(11規(guī)規(guī)范范化化Tcc .,472136. 4,809017. 1)5 . 0(,)1,618034. 1(,118034. 125. 11111111111TTuuIRecu 上頁上頁下頁下頁.1111ecR 則有則有(2) 約化計算約化計算:,111 RU

53、則得到則得到上上Hessenberg陣陣為為.200000. 2399999. 00400000. 0799999. 7472136. 4447214. 0602631. 74112HAUUA 上頁上頁下頁下頁（3）化對稱矩陣為對稱三對角矩陣）化對稱矩陣為對稱三對角矩陣推論推論 (Household約化對稱矩陣為對稱三對角矩陣約化對稱矩陣為對稱三對角矩陣) 設設ARnn為對稱矩陣，則存在初等反射矩陣為對稱矩陣，則存在初等反射矩陣U1,U2,Un- -2使使.111222111221122CcbbcbbcbbcUUAUUUUnnnnnn 為為對稱三對角矩陣對稱三對角矩陣.上頁上頁下頁下頁8.4.

54、2 QR 算法及其收斂性算法及其收斂性QR算法算法是一種變換方法，是計算一般矩陣是一種變換方法，是計算一般矩陣(中小中小型矩陣型矩陣)全部特征值全部特征值問題的問題的最有效方法之一最有效方法之一. (1) 上上Hessenberg矩陣矩陣的的全部特征值全部特征值問題；問題； (2) 計算計算對稱三對角矩陣對稱三對角矩陣的的全部特征值全部特征值問題，問題， 目前目前QR算法算法主要主要用來計算：用來計算：且且QR算法算法具有收斂快，算法穩(wěn)定等特點具有收斂快，算法穩(wěn)定等特點.上頁上頁下頁下頁 對于一般矩陣對于一般矩陣ARnn (或對稱矩陣或對稱矩陣)，則首先用，則首先用Household方法將方法

55、將A化為上化為上Hessenberg陣陣B(或對稱三或對稱三對角陣對角陣)，然后再用，然后再用QR算法算法計算計算B的全部特征值的全部特征值. 設設ARnn，且，且A對進行對進行QR分解，即分解，即A=QR，其中其中R為上三角陣為上三角陣, Q為正交陣為正交陣, 上頁上頁下頁下頁于是可得到一新矩陣于是可得到一新矩陣B=RQ=QTAQ.顯然，顯然，B是由是由A經(jīng)過正交相似變換得到，因此經(jīng)過正交相似變換得到，因此B與與A的的特征值相同特征值相同. 再對再對B進行進行QR分解，又可得一新矩陣分解，又可得一新矩陣,重復這一過程可得到矩陣序列：重復這一過程可得到矩陣序列： 設設A=A1 將將A1進行進行

56、QR分解分解A1=Q1R1 作矩陣作矩陣A2=R1Q1=Q1TR1Q1 QR算法，就是利用矩陣的算法，就是利用矩陣的QR分解，按上述遞分解，按上述遞推法則構造矩陣序列推法則構造矩陣序列Ak的過程的過程. 只要只要A為非奇異矩為非奇異矩陣，則由陣，則由QR算法就完全確定算法就完全確定Ak.上頁上頁下頁下頁 定理定理13 (基本基本QR算法算法) 設設A=A1Rnn, 構造構造QR算法算法:)1 . 4(), 2 , 1(),(1 kQRARIQQRQAkkkkkTkkkk為為上上三三角角陣陣其其中中則則有有記記,1221RRRRQQQQkkkk ;,111kTkkkkAQQAAA 即即相相似似于

57、于）（;)()()2(1211211kTkkTkkQAQQQQAQQQA .)3(kkkkRQAA 的的分分解解式式為為上頁上頁下頁下頁.)()(11111111111111111kkkkkkTkkkkkkkkkkkkkkAAARQARQAQQRAQRRAQQRRRQQQRQ 于是于是 證明證明 (1),(2)顯然，現(xiàn)證顯然，現(xiàn)證(3). 用歸納法，顯然，當用歸納法，顯然，當k=1時有時有 , 設設Ak- -1有分解式有分解式1111RQRQA ,111 kkkRQA上頁上頁下頁下頁 定理定理14 (QR算法的收斂性算法的收斂性) 設設A=(aij)Rnn, (1) 如果如果A的特征值滿足的特征值滿足:; 021 n (2) A有標準型有標準型A=XDX- -1, 其中其中D=diag( 1, 2, n),且設且設X- -1有三角分解有三角分解X- -1=LU(L為單位下三