趙 洋 劉曉宙,2

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

引言

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

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

1.1 QR分解

QR分解的基本步驟是：令A0=A,對k=1,2,…

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

(1)

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

令Ak=Rk-1Qk-1

(2)

(3)

從而Ak～Ak-1～…～A0

(4)

假設對于矩陣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 矩陣特征值的導出

先介紹幾個定理[9]：

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

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

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

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

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

1.3 矩陣特征向量的導出

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

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

(A-λE)ξ=θ

(8)

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

下面說明零空間的求法：

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

(9)

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

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

對U做適當的列交換，可以把U化成如下U'的形式

(10)

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

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

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

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

2.1 實對稱矩陣

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

(11)

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

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

(12)

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

2.2 一般實系數矩陣

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

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

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

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

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

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

3 計算實例

3.1 實對稱陣

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

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

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

3.2 一般實系數的特征值和特征向量

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

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

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

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

4 結論

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