逆冪法結(jié)合原點(diǎn)平移加速求解特征值的精確化計(jì)算

周玉娟 秦菲菲 江山

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

【關(guān)鍵詞】冪法 逆冪法 原點(diǎn)平移 加速收斂 特征值精確化

【基金項(xiàng)目】國家自然科學(xué)基金資助項(xiàng)目（11301462）；江蘇省高校青藍(lán)工程優(yōu)秀青年骨干教師資助項(xiàng)目；“南通大學(xué)教學(xué)改革研究課題”。

【中圖分類號】G64 【文獻(xiàn)標(biāo)識碼】A 【文章編號】2095-3089（2018）37-0124-02

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

一、方法的數(shù)學(xué)理論

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

冪法的過程是假設(shè)特征值滿足|？姿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。則當(dāng)?shù)螖?shù)k→∞且 <1（j=2，3，…，n）時(shí)有xk→？姿1k？琢1 1，故Axk≈？姿1xk。這意味著？姿1是冪法求出的按模最大特征值，與此同時(shí)還求出了對應(yīng)的特征向量xk。

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

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

關(guān)于選取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，只要保證 << 成立，那么收斂過程就會(huì)得到加速，這個(gè)技巧稱為結(jié)合原點(diǎn)平移的加速。上述方案的不足之處在于p的選取并非事先預(yù)知，通常需要在對特征值分布有一定了解的基礎(chǔ)上，才能大致估算出p值并合理選擇，并通過計(jì)算不斷地進(jìn)行精確化。

二、算例與編程實(shí)現(xiàn)

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

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

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

由結(jié)合原點(diǎn)平移的逆冪法公式y(tǒng)k= xk+1=（A-pI）-1ykk=0，1…，取初值x0=（1，0，0）T，列表計(jì)算：

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

接下來，我們再給出Matlab編程來驗(yàn)證上述結(jié)果的正確與效率。

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

A_inv=inv（A）

x0=[1 0 0]'

k=0；

while k<16 %迭代次數(shù)

x0_max=max（abs（x0））；

x0=x0/x0_max；

x0=A_inv*x0

k=k+1；

end

k

mu_1=max（abs（x0）） %正負(fù)號判斷

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 %迭代次數(shù)

x0_max=max（abs（x0））；

x0=x0/x0_max；

x0=B_inv*x0

k=k+1；

end

k

mu_1=-max（abs（x0）） %正負(fù)號判斷

mu_n=1/mu_1

lambda_n=mu_n+p

運(yùn)行結(jié)果為：

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

因?yàn)锳的三個(gè)特征值為0.585786437626905，2，3.4142135 6237309，Matlab運(yùn)行結(jié)果再次驗(yàn)證了前者直接求需要16次迭代，而后者結(jié)合原點(diǎn)平移技巧達(dá)到同樣的精度只需10次迭代，充分體現(xiàn)了該方法精確化求解特征值及其對應(yīng)特征向量的有效性。

參考文獻(xiàn)：

[1]楊淑娥， 陳立萍. 求實(shí)矩陣全部特征值的投影冪法[J]. 北京工業(yè)大學(xué)學(xué)報(bào)， 1994（2）： 77-82.

[2]李磊. 快速并行乘冪法及反冪法[J]. 計(jì)算數(shù)學(xué)， 1995（3）： 253-259.

[3]張青， 茍國楷， 呂崇德. 乘冪法的改進(jìn)算法[J]. 應(yīng)用數(shù)學(xué)與計(jì)算數(shù)學(xué)學(xué)報(bào)， 1997（1）： 51-55.

[4]侯風(fēng)波， 汪永高. 求n階實(shí)方陣的任意個(gè)模最大的特征值及其相應(yīng)特征向量的規(guī)范化冪法[J]. 工科數(shù)學(xué)， 2001（6）： 41-45.

[5]楊長青， 侯建花. 一般矩陣小擾動(dòng)的特征值近似計(jì)算的改進(jìn)[J]. 科學(xué)技術(shù)與工程， 2006（20）： 3337-3338.

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

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

作者簡介：

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