逆幂法结合原点平移加速求解特征值的精确化计算
2018-11-07周玉娟秦菲菲江山
周玉娟 秦菲菲 江山
【摘要】针对非奇异方阵，使用逆幂法求其接近于某给定值的误差最小近似特征值，结合原点平移方案，选择合理的平移值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-），男，湖南湘潭人，副教授，博士，主要研究微分方程数值解及其应用。
