逆冪法結合原點平移加速求解特征值的精確化計算

周玉娟 秦菲菲 江山

【摘要】針對非奇異方陣，使用逆冪法求其接近于某給定值的誤差最小近似特征值，結合原點平移方案，選擇合理的平移值p，從數學理論、數值算例兩方面驗證方法的有效性，故而達到加速收斂、精確計算特征值的目的。

【關鍵詞】冪法 逆冪法 原點平移 加速收斂 特征值精確化

【基金項目】國家自然科學基金資助項目（11301462）；江蘇省高校青藍工程優秀青年骨干教師資助項目；“南通大學教學改革研究課題”。

【中圖分類號】G64 【文獻標識碼】A 【文章編號】2095-3089（2018）37-0124-02

特征值是代數學的一個重要概念，在數學、物理、化學、生物、計算機、設計優化等領域有廣泛應用和悠久歷史，近來又在擾動分析、微分方程、并行計算等方面展現新的研究價值。若非奇異矩陣A的階數為n，則有n個特征值與其對應的特征向量。當矩陣為高階時，要精確地求出所有特征值十分困難、或代價太大而變得沒有必要。在實際工程計算中，往往采用數值方法精確化求出其中具有代表性的部分特征值，比如用冪法求最大特征值、用逆冪法求最小特征值，用對分法或QR分解求特殊矩陣的全部特征值等。而這些數值法有時收斂速度過慢、迭代次數過多，在本文中我們考慮結合原點平移方案，合理選擇平移值，從而實現特征值的加速收斂與計算精確。

一、方法的數學理論

矩陣的特征值問題Ax=？姿x，若直接求解行列式方程det（A-？姿I）=0計算量會非常巨大，且對高于五次的含？姿多項式不存在通用的求根公式。因而對于高階矩陣A，人們往往追求其滿足精度要求、足夠精確的近似特征值及其對應的特征向量。

冪法的過程是假設特征值滿足|？姿1|>|？姿2|≥|？姿3|≥…|？姿n|，取x0=？琢1 1+？琢2 2+…+？琢n n，使？琢1≠0，由迭代格式xk=Axk-1=？姿1k？琢1 1+？姿2k？琢2 2+…+？姿nk？琢n n=？姿1k？琢1 1+ ？琢2 2+…+ ？琢n n。則當迭代次數k→∞且 <1（j=2，3，…，n）時有xk→？姿1k？琢1 1，故Axk≈？姿1xk。這意味著？姿1是冪法求出的按模最大特征值，與此同時還求出了對應的特征向量xk。

逆冪法是冪法的逆思路，由代數學可知原矩陣A與逆矩陣A-1的特征值互為倒數、且特征向量相同。于是，逆冪法的本質是用冪法求A-1的按模最大特征值，然后再求其倒數就得到A的按模最小特征值。為防止計算機運算溢出，引入中間向量yk，將逆冪法迭代公式改進為yk= xk+1=A-1ykk=0，1…，這樣迭代一定次數后能求出滿足精度要求的最小特征值？姿n及其對應的特征向量xk。

然而，上述方法的收斂速度都取決于因子 的大小，當其小于1但接近于1時速度會很慢。我們希望加快收斂速度、減少迭代次數，可以考慮結合原點平移加速的方案。設已知矩陣A和另一矩陣B=A-pI（p為待定常數），可知兩者的特征值有如下關系：若？姿i是A的特征值，則？滋i=？姿i-p就是B的特征值，且對應的特征向量相同。

關于選取p的原則，類似地要滿足|？姿1-p|>|？姿j-p|（j=2，3，…，n）且使得 << 越小越好。這樣，新矩陣B的按模最大特征值？姿1-p：xk=（A-pI）xk-1=（？姿1-p）k？琢1 1+ k？琢j j，只要保證 << 成立，那么收斂過程就會得到加速，這個技巧稱為結合原點平移的加速。上述方案的不足之處在于p的選取并非事先預知，通常需要在對特征值分布有一定了解的基礎上，才能大致估算出p值并合理選擇，并通過計算不斷地進行精確化。

二、算例與編程實現

在這部分，我們給出算例來驗證方法的有效性。在算例中我們采用逆冪法并不是直接求最小的特征值，而是間接求接近于某給定值的、誤差最小的特征值，并最終通過數值編程來檢驗正確與效率。

例.求矩陣A=2 1 01 2 10 1 2最接近于1的特征值及對應的特征向量。

解：因要求最接近于1的特征值，取p=1，令B=A-pI=1 1 01 1 10 1 1，有B-1=0 1 -11 -1 1-1 1 2。

由結合原點平移的逆冪法公式yk= xk+1=（A-pI）-1ykk=0，1…，取初值x0=（1，0，0）T，列表計算：

執行10次迭代后，知A-pI的最大特征值？滋1=-2.4142，最小的為倒數 =-0.4142，故最接近于1的特征值為 +p=0.5858，對應的特征向量為x10=（1.7076，-2.4142，1.7066）T。

接下來，我們再給出Matlab編程來驗證上述結果的正確與效率。

A=[2 1 0； 1 2 1； 0 1 2] %直接求A的特征值

A_inv=inv（A）

x0=[1 0 0]'

k=0；

while k<16 %迭代次數

x0_max=max（abs（x0））；

x0=x0/x0_max；

x0=A_inv*x0

k=k+1；

end

k

mu_1=max（abs（x0）） %正負號判斷

mu_n=1/mu_1

lambda_n=mu_n

p=1； %選擇p值，可以改變，但需符合要求

B=A-p*eye（3） %間接求B的特征值

B_inv=inv（B）

x0=[1 0 0]'

k=0；

while k<10 %迭代次數

x0_max=max（abs（x0））；

x0=x0/x0_max；

x0=B_inv*x0

k=k+1；

end

k

mu_1=-max（abs（x0）） %正負號判斷

mu_n=1/mu_1

lambda_n=mu_n+p

運行結果為：

k =

16

mu_1 =

1.7071067811912

mu_n =

0.58578643762531

lambda_n =

0.58578643762531

k =

10

mu_1 =

-2.41421319796954

mu_n =

-0.41421362489487

lambda_n =

0.58578637510513

因為A的三個特征值為0.585786437626905，2，3.4142135 6237309，Matlab運行結果再次驗證了前者直接求需要16次迭代，而后者結合原點平移技巧達到同樣的精度只需10次迭代，充分體現了該方法精確化求解特征值及其對應特征向量的有效性。

參考文獻：

[1]楊淑娥， 陳立萍. 求實矩陣全部特征值的投影冪法[J]. 北京工業大學學報， 1994（2）： 77-82.

[2]李磊. 快速并行乘冪法及反冪法[J]. 計算數學， 1995（3）： 253-259.

[3]張青， 茍國楷， 呂崇德. 乘冪法的改進算法[J]. 應用數學與計算數學學報， 1997（1）： 51-55.

[4]侯風波， 汪永高. 求n階實方陣的任意個模最大的特征值及其相應特征向量的規范化冪法[J]. 工科數學， 2001（6）： 41-45.

[5]楊長青， 侯建花. 一般矩陣小擾動的特征值近似計算的改進[J]. 科學技術與工程， 2006（20）： 3337-3338.

[6]陳靜，顧晨宇，錢椿林.一類微分方程廣義特征值的算法[J]. 蘇州科技學院學報， 2006（2）： 9-13.

[7]曹芳芳，呂全義，聶玉峰.求實對稱矩陣部分特征值的并行算法[J]. 計算機工程與設計， 2010（22）： 4820-4823.

作者簡介：

江山（1980-），男，湖南湘潭人，副教授，博士，主要研究微分方程數值解及其應用。