QR方法計算一般矩陣特征值和特征向量的算法研究

趙 洋 劉曉宙,2

(1.南京大學(xué)物理學(xué)院 江蘇 南京 210093 2.近代聲學(xué)教育部重點實驗室，南京大學(xué)聲學(xué)研究所 江蘇 南京 210093)

引言

計算矩陣的特征值和特征向量是計算方法中一個常見的問題.一般常用的方法有乘冪法、反冪法，雅克比法和QR方法[1-2].乘冪法用于求矩陣的模最大的特征值和對應(yīng)的特征向量[3]，反冪法用于求矩陣的模最小的特征值和對應(yīng)的特征向量[4]，雅克比法計算實對稱矩陣的特征值和特征向量[5].QR法通常僅用來求一般矩陣的特征值，獲得特征值后采用反冪法來獲得特征向量[6]，一般來說并不直接通過QR法獲得特征向量.但是反冪法存在著其缺陷，它對輸入的試探向量的值有一定要求，在只進行一次嘗試的情況下得到的不一定就是想要的特征向量[7][8],需要對冪法做一些改進和補充才能得到想要的結(jié)果[9].為了避免迭代法中初值的影響，本文探索直接采用QR算法獲得特征值和特征向量的方法，討論了用QR方法在不同情況下特征向量的一些算法，并對算法進行了驗證.

1 QR方法計算特征值和特征向量的一般方法

1.1 QR分解

QR分解的基本步驟是：令A(yù)0=A,對k=1,2,…

分解Ak-1=Qk-1Rk-1

(1)

其中Qk-1為正交矩陣，Rk-1為上三角矩陣.

令A(yù)k=Rk-1Qk-1

(2)

(3)

從而Ak～Ak-1～…～A0

(4)

假設(shè)對于矩陣Ak，有特征值λk,特征向量Xk，

即 AkXk=λkXk

(5)

代入(3)式，得：

(6)

即有：

Ak-1(Qk-1Xk)=λk(Qk-1Xk)

(7)

由(4)式可知，Ak-1和Ak具有相同特征值λk.Ak-1的特征向量為Xk-1=Qk-1Xk,

遞推可得：Ai(i=1,2,…,k)具有相同特征值λk，A0的特征向量為X0=Q0Q1…Qk-2Qk-1Xk.

因此，如果Ak的特征值和特征向量有辦法求出，那么就能獲得A0的特征值和特征向量.下面分別對特征值和特征向量兩部分的求解做介紹.

1.2 矩陣特征值的導(dǎo)出

先介紹幾個定理[9]：

定理1 設(shè)A=(aij)∈Rn×n

1°如果A的特征值滿足：|λ1|>|λ2|>…>|λn|>0

2° A有標準型A=XDX-1其中D=diag(λ1,λ2,…,λn)，且設(shè)X-1有三角分解X-1=LU(L為單位下三角陣，U為上三角陣)，則由QR算法產(chǎn)生的{Ak}本質(zhì)上收斂于上三角陣.

定理2 設(shè)A=(aij)∈Rn×n，如果A的等模特征值中只有實重特征值或者多重復(fù)的共軛特征值，則QR算法產(chǎn)生的{Ak}本質(zhì)收斂于分塊上三角陣(對角塊為一階和二階子塊)且對角塊每一個2×2子塊給出A的一對共軛復(fù)特征值，每一個對角子塊給出A的實特征值.

本文討論的矩陣均滿足上述定理，可以通過QR算法收斂于分塊上三角陣，求得特征值.

1.3 矩陣特征向量的導(dǎo)出

在不借助冪法的情況下，采用直接由特征值特征向量的定義計算.

特征值特征向量的定義為：存在某個向量ξ，使得Aξ=λξ，那么就稱λ為矩陣的一個特征值，ξ為其對應(yīng)的特征向量.所以在特征值已知的情況下，只要求解：

(A-λE)ξ=θ

(8)

即可求出特征值λ對應(yīng)的特征向量ξ，即求新矩陣A-λE的零空間.

下面說明零空間的求法：

任意一個矩陣均可以通過行變換化為行最簡形矩陣.行最簡形矩陣的定義為非零行的第一個非零元素全是1，且非零行的第一個元素1所在列的其余元素全為零.例如：

(9)

接下來就將借助行最簡形矩陣求出矩陣的零空間.

對于一個矩陣A和它的行最簡形矩陣U，Ax=θ和Ux=θ應(yīng)當為同解方程組，求出U的零空間就求出了A的零空間.

對U做適當?shù)牧薪粨Q，可以把U化成如下U'的形式

(10)

I是單位矩陣，K是一個一般的非零矩陣.

由于列交換會將解x中的變量順序交換，所以求出U′x′=θ的解x'后，將其中的變量依照列變換順序反變換回去，就得到了我們想要的解x.

2 不同類型矩陣特征值和特征向量的具體計算

前文已經(jīng)得到了特征值和特征向量求解的主要算法.下面對不同情況的矩陣進行進一步的討論，主要分類為實對稱矩陣和一般實系數(shù)矩陣.在一般實系數(shù)矩陣的討論中又分為可對角化矩陣和不可完成對角化的矩陣.

2.1 實對稱矩陣

當A0為實對稱陣時，特征值均為實數(shù)，且經(jīng)過正交相似變換后得到的仍然是一個實對稱陣，所以經(jīng)過QR算法最終一定收斂到一個對角陣A∞.

(11)

很容易知道單位陣I中任意一個列向量均能成為該對角陣的特征向量，所以可以認為它的特征值組成的矩陣就是單位陣I.由此A0的特征向量組成的矩陣可表示為:

Q0Q1……Qk-2Qk-1I=Q0Q1…Qk-2Qk-1=P

(12)

取P的列向量即得到A的特征向量.

2.2 一般實系數(shù)矩陣

對于非實對稱的矩陣，它和實對稱矩陣主要有如下區(qū)別：

(1)矩陣不一定可以對角化.即在特征值出現(xiàn)重根的情況下，n階方陣不一定有n個線性無關(guān)的特征向量.

(2)特征值不一定都是實數(shù).所以經(jīng)過QR法計算出來最終收斂的矩陣可能是一個分塊上三角陣，此時的特征值會出現(xiàn)共軛復(fù)根.

(3)(擬)上三角化后得到的矩陣A∞的特征向量不再具有單位陣這樣的簡單形式，需要通過求矩陣零空間的方法將特征向量求出.

所以對一般實系數(shù)矩陣求特征值的步驟為：使用QR法讓其收斂到一個(分塊)上三角陣，將正交相似變換矩陣Q存下來，求出該(分塊)上三角陣的特征值和特征向量，求出的特征值就是對應(yīng)原矩陣的特征值，特征向量經(jīng)過矩陣Q的作用得到原矩陣的特征向量.

由于零空間的算法具有普適性，所以不論矩陣是否可以對角化，總能求出所有線性無關(guān)的解，可以根據(jù)最終求得的特征向量的個數(shù)來判斷是否可以對角化.

3 計算實例

3.1 實對稱陣

先選取了一個特征值均為單重根的實對稱矩陣的情況：

各符號代表的含義是：λ0為matlab內(nèi)置函數(shù)計算出來的特征值，λ為QR算法求出來的上三角陣的對角線元素，Q為算法中的正交矩陣不斷累乘得到的變換矩陣，在此情況下，Q的每一列就對應(yīng)著A的一個特征向量.后續(xù)經(jīng)過特征向量的定義驗證該結(jié)果正確.

選取特征值含重根的實對稱矩陣驗證，計算結(jié)果也是正確的.

3.2 一般實系數(shù)的特征值和特征向量

特征值均為單根且可以對角化的例子：

其中A為待求矩陣，λ0為matlab內(nèi)置函數(shù)求出的特征值，λ為QR方法計算后的分塊上三角陣的特征值.P是分塊上三角陣的線性無關(guān)的特征向量，D為正交相似變換矩陣，P1=DP就是A的特征向量.

考慮特征值出現(xiàn)重根且矩陣不可對角化的例子：

符號含義與上例相同，和上例的區(qū)別在于求出來的P只有三個線性無關(guān)的列向量，矩陣無法完成對角化.使用matlab內(nèi)置函數(shù)驗證時，特征向量有兩列是線性相關(guān)的，這也驗證了該矩陣無法完成對角化.

4 結(jié)論

采用QR方法計算一般實系數(shù)矩陣的特征值和特征向量，此方法的優(yōu)點在于不需要輸入試探向量進行迭代，避免了冪法中最終結(jié)果的正確性與輸入的試探向量直接相關(guān)的問題.而且該算法具有普適性，不論矩陣是否可以對角化，只要能夠滿足定理條件，用QR算法最終收斂到一個(擬)上三角陣，就可以求得所有的特征值和對應(yīng)的特征向量.利用該(擬)上三角陣求解零空間時比較方便，因為擬三角陣近似一個行梯形，在變成行最簡型矩陣時更為方便.整個算法的缺點是需要在QR分解的時候?qū)⒚看蔚恼幌嗨谱儞Q陣Q累乘起來，加大了一定的計算量.