數(shù)值分析原理第八章

1、第八章 矩陣特征值和特征向量的計算對于nn階的實矩陣，線性代數(shù)理論中是通過求解特征多項式的零點而得到特征值l，然后通過求得齊次線性方程組的非零向量而得到矩陣的相應于特征值l的特征向量當矩陣階數(shù)較高時，這種方法計算量很大，故常用數(shù)值方法近似求解特征值與特征向量目前常用的數(shù)值方法有迭代法(冪法)和變換法(Jacobi方法等)兩類8.1 乘冪法與反冪法一、乘冪法乘冪法是求矩陣按模最大的特征值(主特征值)和相應的特征向量的一種迭代法設，初始向量，令(8.1)生成迭代向量序列由遞推公式(8.1)，知(8.2)這表明等于用矩陣的k次冪左乘，故稱此方法為乘冪法下面分析當k時，向量序列的變化規(guī)律設，為矩陣的n

2、個特征值，且滿足(8.3)相應于特征值，的n個線性無關(guān)的特征向量構(gòu)成向量空間上的一組基任取非零的初始向量，則可由這組特征向量線性表出(8.4)其中為線性組合系數(shù)將式(8.4)代入式(8.2)，得(8.5)由和式(8.5)，得(8.6)當時，由式(8.3)知，特征值下面針對進行討論由式(8.6)有由于，故若，當k充分大時，此時有(8.7)上式表明，與只近似相差一個常數(shù)因子，故可取作為矩陣的相應于主特征值的近似特征向量當k充分大時，若，則有(8.8)這表明主特征值可由式(8.8)近似求得如果矩陣的特征值滿足則根據(jù)式(8.6)有(8.9)則當k充分大時，由于，故有 (8.10)由于都是矩陣的特征值對

3、應的特征向量，故也是矩陣的特征值對應的特征向量由式(8.10)知，k較大時，就是與主特征值的對應的近似特征向量類似于式(8.8)，可求得主特征值的近似由于此時的特征向量子空間不是一維的，故由式(8.10)得到的近似特征向量只是該子空間的一個特征向量，對于不同的初始向量可能得到的特征向量空間中線性無關(guān)的近似特征向量對于矩陣的其它主特征值情形，如，等，同樣可以用乘冪法求解，具體過程可參閱文獻6以上討論說明了乘冪法的基本原理通過上述對乘冪法過程的分析可知，乘冪法是一種迭代法，公式計算簡單，便于上機實踐，可以方便地用于近似求矩陣按模最大的一個(或幾個)特征值及相應的特征向量需要注意的是：(1) 從理論

4、上講，對于任意給定的初始向量，有可能使式(8.4)中的，但因舍入誤差的存在，隨著迭代過程的進行，等效于從的出發(fā)進行迭代(2) 在用乘冪法(8.1)進行迭代計算時，迭代向量的分量的絕對值可能會出現(xiàn)非常大(當)或者非常小(當)的現(xiàn)象，甚至出現(xiàn)溢出為此，實用中每進行m步就需要對迭代向量進行一次規(guī)范化，即用(其中表示向量的按模最大的分量)代替繼續(xù)迭代由于特征向量允許相差一個常數(shù)因子，故按前面乘冪法的理論依然得到正確的近似特征向量在每次規(guī)范化后，用乘冪計算前后兩個向量的分量的比值作為主特征值的近似，這種規(guī)范化并不影響主特征值的近似計算。規(guī)范化的乘冪法避免了溢出的可能性，至于m取多少，取決于實際情形，可以

5、取m=5或m=1等等下面給出乘冪法的具體算法算法8.1 (乘冪法)輸入 矩陣、誤差限、任意初始向量、m；輸出 主特征值的近似值及相應的近似特征向量；Step 1 利用式(8.1)求，并由式(8.8)求，置k=1；Step 2 利用式(8.1)求，并利用式(8.8)求；Step 3 判斷是否滿足？若滿足，則取，結(jié)束；否則，轉(zhuǎn)向Step4；Step 4 置k=k+1，判斷mod(k1,m)=0？若滿足，取；轉(zhuǎn)向Step2例8.1 用乘冪法計算矩陣的主特征值及相應特征向量，要求解 取初始向量，用乘冪法公式進行計算，且取，每迭代5步進行一次規(guī)范化，計算結(jié)果見表8.1表8.1 乘冪法的計算結(jié)果k0123

6、453637(1.0000000 1.0000000 1.0000000)(-6.0000000 2.0000000 8.0000000)(102.0000000 -32.0000000 34.0000000)(-1218.0000000 206.0000000 608.0000000)(17058.0000000 -4664.0000000 190.0000000)(-218118.0000000 46130.0000000 61832.0000000)(-13.2201800 3.1081369 2.2688627)(174.7731588 -41.0901285 -28.9947753)

7、-6.0000000-17.0000000-11.9411765-14.0049261-12.7868449-13.2201800-13.2201800取-13.2201800作為主特征值的近似值，相應于的特征向量的近似值取為(174.7731588, -41.0901285, -28.9947753)T二、原點平移法乘冪法的收斂速度取決于或的程度，r1時收斂速度較快，r1時收斂速度較慢對收斂較慢的乘冪法，可以采用原點平移法加快其收斂速度設矩陣，這里p是可以選取的參數(shù)當?shù)奶卣髦禐闀r，的特征值，且與有相同的特征向量若的主特征值為，則要選擇適當?shù)膮?shù)p，使其滿足是的主特征值，即 ()，且對矩陣應用

8、乘冪法，可以使得計算的主特征值的過程得到加速這種方法通常稱為原點平移法對于參數(shù)p的選擇依賴于對矩陣的特征值分布的大致了解通?？梢杂肎erschgorin圓盤定理得到矩陣的特征值分布情況定理8.1(Gerschgorin圓盤定理)19 設為n階實矩陣，則(1) 的每一個特征值必定屬于下述n個閉圓盤(稱為Gerschgorin圓) (8.11)的并集；(2) 由矩陣的所有Gerschgorin圓組成的連通部分中任取一個，如果它是由k個Gerschgorin圓構(gòu)成，則在這個連通部分中有且僅有的k個特征值(Gerschgorin圓相重時重復計數(shù)，特征值相同時也重復計算)求得矩陣的主特征值后，可得矩陣的

9、主特征值，同時得到對應的特征向量的近似例8.2 對例8.1的矩陣，取p=4.6，用原點平移法求其主特征值及相應的特征向量解 對應用乘冪法算法8.1，每迭代1步進行一次規(guī)范化，計算結(jié)果見表8.2由表8.2可知，的主特征值= -17.8201809+4.6= -13.2201809，相應的特征向量的近似值取為(1.0000000, -0.2351055, -0.1716216)T表8.2 原點平移法的計算結(jié)果k01231112(1.0000000 1.0000000 1.0000000)(-10.6000000 -2.6000000 3.4000000)(-16.8264160 2.7584906

10、 1.7396227)(-17.4019737 3.7969499 3.0797486)(-17.8201809 4.1896219 3.0583200)(-17.8201809 4.1896216 3.0583203)(1.0000000 1.0000000 1.0000000)(1.0000000 0.2452830 -0.3207547)(1.0000000 -0.1639281 -0.1033864)(1.0000000 -0.2181908 -0.1769770)(1.0000000 -0.2351055 -0.1716212)(1.0000000 -0.2351055 -0.171

11、6216)-10.6000000-16.8264160-17.4019737-17.8201809-17.8201809事實上，如果對于矩陣的特征值能夠分離的很清楚，可以利用原點平移法求得矩陣的所有特征值及其相應的特征向量但需要說明的是，雖然常常能夠選擇合適的p值使乘冪法得到加速，但設計一個自動選擇適當參數(shù)的過程非常困難原點平移法的價值不在于直接使用它使迭代過程加速，而在于把原點平移法與其它方法(如反冪法等)結(jié)合使用以獲得更好的效果三、反冪法反冪法又稱逆冪法，它是近似求矩陣按模最小的特征值及其相應特征向量的一種方法設，且無零特征值，則若l為矩陣的特征值，則必為矩陣的特征值，且與具有相同的特征向

12、量如果的n個特征值滿足則的n個特征值滿足因此，若用做乘冪矩陣，由乘冪迭代格式 (8.12)便可求出的按模最大特征值，進而得到矩陣的按模最小的特征值因此，對任取初始向量，稱式(8.12)為求矩陣按模最小特征值的反冪法在應用式(8.12)計算時，高效的算法并不是先計算的逆矩陣，而是通過求解如下線性方程組求得 (8.13)由于在迭代過程中求解的線性方程組(8.13)的系數(shù)矩陣不變，故可以利用矩陣的三角分解，則每次迭代只需解兩個三角形方程組 (8.14)下面給出反冪法的具體算法算法8.2 (反冪法)輸入 輸入矩陣、誤差限、任意初始向量；輸出 按模最小的特征值的近似值及相應的近似特征向量；Step 1

13、利用式(8.14)求，并由式(8.8)求，置k=1；Step 2 利用式(8.14)求，并利用式(8.8)求；Step 3 判斷是否滿足？若滿足，則取，結(jié)束；否則，置k=k+1，轉(zhuǎn)向Step2根據(jù)反冪法與乘冪法的上述關(guān)系，若p是的相對分離較好的近似值()，可以結(jié)合原點平移法與反冪法來求得的更精確的特征值與相應的特征向量的近似此時只需利用反冪法算法8.2求解的按摸最小的特征值及其相應的特征向量即可8.2 Jacobi方法Jacobi方法是近似求實對稱矩陣全部特征值及相應特征向量的一種變換方法其基本思想是將對稱矩陣經(jīng)一系列正交相似變換約化為一個近似對角陣，從而該對角陣的對角元就是的近似特征值，由各

14、個正交變換陣的乘積可得對應的特征向量在本節(jié)中，將用到下列線性代數(shù)知識：(1) 若矩陣與相似，即存在可逆陣使，且與具有完全相同的特征值(2) 若矩陣滿足，則稱為正交矩陣顯然，且當，是正交矩陣時，其乘積仍為正交矩陣(3) 實對稱矩陣的特征值均為實數(shù)(4) 對任何實對稱矩陣，總存在正交陣，使得=，其中為的特征值，的各列為相應的特征向量(5) 設為實對稱矩陣，為正交陣，則(6) 稱矩陣為n階的Givens旋轉(zhuǎn)矩陣，它是在單位陣的p行、q行和p列、q列的四個交叉位置上分別置上cosq，sinq，-sinq和cosq而成的容易驗證旋轉(zhuǎn)矩陣是正交矩陣，即，用它作相似變換時十分方便Jacobi方法就是用這種旋

15、轉(zhuǎn)矩陣對實對稱陣作一系列的旋轉(zhuǎn)相似變換，從而將約化為近似對角陣一、Jacobi方法首先考慮的二階實對稱矩陣的對角化取二階矩陣對矩陣做旋轉(zhuǎn)變換，即其中 (8.15)由式(8.15)的最后一式知，要使的相似矩陣成為對角陣，只須適當選取q，使即時取 (8.16)當時，取確定q后，旋轉(zhuǎn)矩陣則隨之確定此時得到矩陣的特征值為，它們對應的特征向量分別是的各列，即對應于的特征向量分別是可以看出，對二階實對稱矩陣，用適當?shù)恼幌嗨谱儞Q，一次即可把化為對角陣當為n階實對稱矩陣時，需使用n階Givens旋轉(zhuǎn)矩陣容易驗證，對于n階Givens旋轉(zhuǎn)矩陣來說，只改變的第p行與第q行元素，只改變的第p列與第q列元素，只改變

16、的第p行、第q行、第p列與第q列元素Jacobi方法是通過一系列旋轉(zhuǎn)相似變換逐漸將實對稱矩陣化為對角陣的過程，即 (8.17)恰當?shù)剡x取每個旋轉(zhuǎn)矩陣，就可以使趨于對角化設，其中p和q分別指向矩陣的嚴格上三角陣中絕對值最大的元素的行號和列號變換過程中產(chǎn)生的也是實對稱陣與的差別僅在于第p、q行與第p、q列的元素由矩陣乘法可得 (8.18) (8.19) (8.20)由式(8.18)(8.20)知，再結(jié)合正交矩陣的性質(zhì)，可得若，選取q使，只需q滿足 (8.21)或 (8.22)此時則有 (8.23)引入記號 (8.24)式中表示矩陣的對角元的平方和，表示矩陣的非對角元的平方和由于，故有 (8.25)

17、這說明，只要，則按上述方法構(gòu)造的Givens旋轉(zhuǎn)矩陣對變換后就會使對角線元素的平方和增加，而非對角線元素的平方和減少由于的對稱性，實際上只要計算的上三角元素，而下三角元素由對稱性獲得，這樣即節(jié)省了計算量，又避免了舍入誤差對對稱性的影響需要指出的是，若在某一步已有，則可能又變?yōu)榉橇阍匾虼耍⒉荒鼙ＷC通過有限次旋轉(zhuǎn)相似變換將矩陣化為對角陣Jacobi方法中每一步都有兩個主要步驟(1) Givens旋轉(zhuǎn)矩陣的確定選取Givens旋轉(zhuǎn)矩陣，使得，則需使q滿足式(8.21)或式(8.22)考慮到舍入誤差的影響，當時，取 (8.26)當時，記 (8.27)由式(8.21)及式(8.27)有 (8.28)

18、 (8.29)為避免相近數(shù)相減(當非常大時)，取 (8.30)進而有 (8.31)(2) 特征向量的計算由式(8.17)知，若記，則有下述遞推關(guān)系 (8.32)如果趨于對角陣，則的各列就是近似特征向量的計算可由式(8.32)得到，具體地其元素間的關(guān)系為 (8.33)下面給出Jacobi方法求矩陣特征問題的具體算法算法8.3 (Jacobi算法)輸入 輸入矩陣、誤差限；輸出 矩陣的所有特征值及其相應的特征向量；Step 1 置，k=0；Step 2 選主元，即確定，使；Step 3 按式(8.27)或式(8.31)計算；Step 4 按式(8.33)計算新正交陣的元素；Step 5 按式(8.1

19、8)(8.20)計算Ak+1的元素；Step 6 計算；若，則停止計算，否則，置k=k+1，轉(zhuǎn)向Step2下面給出Jacobi方法收斂性的結(jié)論定理8.2 設為實對稱矩陣，則由出發(fā)，Jacobi方法所產(chǎn)生的矩陣序列收斂于對角陣，且就是實對稱矩陣的全部特征值證明 由式(8.26)可知由知從而有進而有遞推得到當時，顯然有，即非對角元的平方和趨于零，趨于對角陣，Jacobi方法收斂例8.3 用Jacobi方法求例8.1中實對稱矩陣的全部特征值與特征向量，取解 按Jacobi算法8.3編程上機計算，結(jié)果見表8.3從表8-3可以看出，的特征值分別為，對應的特征向量分別為表8.3 Jacobi方法的計算結(jié)果k0AI1267，二、Jacobi法的變形使用Jacobi方法進行計算時，每次先要在非對角元中循環(huán)掃描以挑選主元，這要花費較多的機時實用中，常常對Jacobi方法進行修正在開始掃描之前，首先確定一個閥值，在實對稱矩陣的嚴格