矩阵特征值求解.docx
《矩阵特征值求解.docx》由会员分享，可在线阅读，更多相关《矩阵特征值求解.docx（17页珍藏版）》请在冰豆网上搜索。
矩阵特征值求解
矩阵特征值求解的分值算法
12组
1.1矩阵计算的基本问题
（1）求解线性方程组的问题•即给定一个n阶非奇异矩阵A和n维向量b，求一个n维向量X,使得
Ax=b（1.1.1）
（2）线性最小二乘问题，即给定一个mxn阶矩阵A和m维向量b,求一个n维向量
X,使得
|AX-b|=min{|Ay-比严Rn}（1.1.2）
（3）矩阵的特征问题，即给定一个n阶实（复）矩阵A,求它的部分或全部特征值以及对应的特征向量，也就是求解方程
（1.1.3）
一对解（4X）,其中R（C）,x-Rn（Cn）,即A为矩阵A的特征值，X为矩阵
Ax=Zx
A的属于特征值A的特征向量。
在工程上，矩阵的特征值具有广泛的应用，如大型桥梁或建筑物的振动问题：
机械和机件的振动问题；飞机机翼的颤振问题；无线电电子学及光学系统的电磁
振动问题;调节系统的自振问题以及声学和超声学系统的振动问题•又如天文、地
震、信息系统、经济学中的一些问题都与矩阵的特征值问题密切相关。
在科学上，计算流体力学、统计计算、量子力学、化学工程和网络排队的马尔可夫链模拟等实际问题，最后也都要归结为矩阵的特征值问题.由于特征值问题在许多科学和工程领域中具有广泛的应用，因此对矩阵的特征值问题的求解理论研究算法的开发软件的制作等是当今计算数学和科学与工程计算研究领域的重大课题，国际上这方面的研究工作十分活跃。
1.2矩阵的特征值问题研究现状及算法概述
对一个nxn阶实（复）矩阵A,它的特征值问题，即求方程（1.1.3）式的非平凡解，是数值线性代数的一个中心问题•这一问题的内在非线性给计算特征值带来许多计算问题•为了求（1.1.3）式中的A,—个简单的想法就是显式地求解特征方程
det（A一几I）二0（121）
除非对于个别的特殊矩阵，由于特征方程的系数不能够用稳定的数值方法由行列式的计算来求得，既使能精确计算出特征方程的系数，在有限精度下，其特征多项式f〃）二det（A-ZJ）的根可能对多项式的系数非常敏感能•因此，这个方法只
在理论上是有意义的，实际计算中对一般矩阵是不可行的数_._.人
较大，则行列式det（A-几I）的计算量将非常大；其次，根据•首先，右矩即AfbJ阳数大于四的多项式求根不存在一种通用的方法，基于上述原Galois理论对于次因，人们只能寻求其它途径•因此，如何有效地！
精确地求解’
矩阵特征值问题，就成为数值线性代数领域的一个中心问题.
目前，求解矩阵特征值问题的方法有两大类：
一类称为变换方法，另一类称为向
量迭代方法•变换方法是直接对原矩阵进行处理，通过一系列相似变换，使之变换成一个易于求解特征值的形式，如Jacobi算法,Givens算法,QR算法等。
变换方法由于要存储矩阵元素，因而它只适用于求解中小型矩阵，它一般和向量迭代方法结合起
女口Lanczos
来使用•向量迭代方法是通过一系列矩阵向量乘积而求得特征值和特征向量的.由于向量迭代方法可采用压缩存储技术，因而它适合于求大规模矩阵的特征值问题，尤其是大型稀疏矩阵的部分特征值和特征向量问题方法,Amoldi方法Qavidson方法等，现在这类问题仍是比较热的研究课题。
2分治方法的基础及理论研究
2.1分治方法的概述
考虑对称三对角矩阵Tn的特征值问题
（2.1.1）
TnX二ZX
其中
（2.1.2）
1981年Cuppen提出一种求上述对称三对角矩阵几所有特征值和特征向量
的分而治之方法（divide一*and—conquermethod）.其基本思想是先将对称三对角矩阵&分割为两个分别为kXk阶和（n-kp<（n-k）阶低阶对称三对角子矩阵
T（9）和T
（1）.T⑹和T（l））可以用同样的方法也分别分割为两个更低阶的子矩阵,递
归的采用这种分割技术可以把矩阵分割为一些能直接求出特征值的足够小的子矩阵（比如2阶或1阶矩阵），或者按照某种标准分割到适当阶数（如小于等于25阶）后，结合其它求矩阵特征值的方法，如QR算法，求出其特征值。
在求出低阶矩阵特征值的基础上，开始胶合过程。
在胶合阶段，分割前的矩阵『的特征值的求出（所谓的“治之”）是建立在其两个子矩阵T⑹和T⑴的特征值的基础上的，其中严和T⑴•是在分割阶段由『分割出的低阶子矩阵•随后的数值分析表明，Cuppen的方法存在着数值不稳定的危险，特别是当存在特征值束时，计算出的特征向量可能不正交。
Gu和Eisenstat对Cuppen的方法作了改进，极大地降低了数值不稳定的危险性。
Cuppen的方法在计算Tn的特征值的同时也需要计算对应的特征向量，并且
•根据文中，当用残量
Cuppen的方法比二
是在T⑼和T（l）的特征值和特征向量的基础上进行计算的
TnX-AX|和正交性xTX-In作为检验准确性的标准时,
分法或多分法精确的多。
但文［3中，，如果用扎-kTn作为衡量特征值准确性
的标准时，二分法或多分法精确些•此外Cuppen的分而治之方法要求矩阵乘积，
存储量为0（£）,而二分法或多分法的存储量仅为0（n）,比前者少•应此，当只需计算特征值时，通常选用后者.1987年Dongarra和Sorensen把分治思想应用到求对称三对角矩阵特征值的并行计算，取得了不错的效果，再次引起了人们对分治方法的极大关注。
分割胶合方法（split-mergemethod）是后来提出用分治方式求对称三对
角矩阵特征值Tn的方法，不同于分而治之方法在分割过程中采用矩阵的秩1扰动，
它采用矩阵的秩2扰动•与二分法和多分法相似，在分割胶合方法中，特征值的计算独立于特征向量的计算•如果需要计算特征向量，可采用反幕法。
由于分割胶合方法计算特征值时采用具有三次收敛的Laguerre迭代，数值试验表明，其计算
速度和精度都明显优于二分法和多分法•并且文中给出的用于计算Laguerre迭
代的非线性三项递归式可以避免上溢和下溢问题。
2.1.1分割策略
分治方法的第一步就是把原来高阶的对称三对角矩阵特征值问题转化为两个低阶对称三对角矩阵特征值问题示如，即所谓的分割阶段•设对称三对角矩阵Tn表下
不妨假设pjHO（j二1,2,…，n-l）,即称Tn为不可约矩阵。
否则，若存在某些PJ•二0,
则Tn就可以约化为若干个低阶主子矩阵特征值问题，Tn的特征值就由这若干个
低阶主子矩阵的特征值构成•当Tn为不可约对称三对角矩阵时，对其作如下分割
fr（0））
TndT⑴丿
（2.1.1.2）
记Tn二Tn+E，其中T°和T⑴）分别为k沢k阶和（n一k），zn~k）或
（n-k-1）x（n-kT）阶对称三对角矩阵，通常W=[n/2]o
PP0
pepe2
PvvT时，其中V
（o；
T,…，0）且
pl
0
Qk卄用，
Pk卅
Pl
+
〃\
T）=
PkH4
ak*一
■
+
+Pk.
++
J
J
JJ
Pn_
PkA叫-P」
J丿
称此为秩扰动矩阵•此时有
T
（o）
⑵当
PekeJ
（0
称此为秩2
扰动矩阵•此时有
（3）
Pk
时，称此为秩3扰动矩阵•此时有
务
1
pi
g七Pkpk七七
宀
洛1P1
Pl+
■J-
T（l）=
J
■
+
J
■
id
rd
Pk4叫丿
Pn4J丿
fo
Pk
nJ
人们关注的问题是：
对于上述三种分治策略几与Tn的特征值之间的关系.
利用Hoffman-Wielandt定理的结论，我们可以得到如下定理：
定理2.1设A和B是两个n阶Hermite矩阵，它们的特征值分别是7Z'＞沁才,心
和已”2A〃＞巴，贝IJ
n
[J仏j-Pj）]~〈A_BF（2.1.1.3）
其中Hf为Frobenius范数。
根据上述定理，设几与Tn的特征值分别为必舒口
山心”工扎；则有下面关系成立
（2.1.1.4）
因此，只要（2.1.1.4）式右端越小，kj就越接近Aj，即Tn与的特征值就越接
近。
对于秩1和秩2扰动,人与Tn的特征值之间关系更详细的描述，将在后文给出。
2.1.2胶合
在矩阵Tn分割后，先求出矩阵Tn的特征值，即求出T（9）和T⑴的特征值。
剩下的工作就是如何由Tn的特征值出发求出Tn的特征值，这就是胶合阶段主要任务。
在这个阶段可以采用不同的迭代方法，如割线法!
Newton迭代、Laguerre迭代以及路径跟踪等，以的某个特征值或的特征值构成的区间内的点为初始点经过若千步迭代，最后收敛到几的某个特征值•本文中采用Laguerre迭代从特征区间提取特征值。
2.2Laguerre迭代
2.2.1Laguerre迭代及其计算
根据文中，对于不可约对称三对角矩阵Tn,它的特征值或特征多项式
fa）二detTn-Mn）的根全为互不相同的实数•因此，适合求多项式的根都是实
（2.2.1.1）
单根的具有三次收敛的Laguerre迭代
Ljx）=X+
-「罟卯V需j需）]
非常适合用来从特征区间提取出Tn的特征值•这个方法早在13世纪就已经提出，为了避免上述（2.2.1.3）式、（2.2.1.4）式和（221.5）式在计算中可能出现的上溢或下溢问题，文中给出如下的改进的非线性三项递归式：
近年来结合分治策略又重新焕发出生机•文中首次对用Laguerre迭代求不可约
对称三对角矩阵的特征值的实用性进行了研究，而后又被用于求矩阵的奇异值问
题。
为了利用Laguerre迭代求Tn的特征值,我们需要计算Tn的特征多项式
f仏）二det（Tn-讥）以及其一阶导数f，仏）二阶导数f〃仏）•由文中知，计算这些值的有效的方法是三项递归式•众所周知，当矩阵阶数比较大时，三项递归式可能出现上溢或下溢问题•文中给出了一种改进的非线性的三项递归式，数值试验表明，除极其个别的情形，改进的三项递归式可以避免上溢或下溢问题。
对称三对角矩阵Tn的特征多项式及其一、二阶倒数的计算公式
设对称三对角矩阵Tn如（221.1）式所示,Tk表示其k阶顺序主子式，表
示如下
P1
匕PkjPk4«k>
（2.2.1.2）
fk（Q二det（Tk-）
・10为Tk的特征多项式，则特征多项式及其导数计算公式如
下（三项递归式）:
rfoG）二1,fl（Q=%
tfi仏）=（%-Qfi4仏）—Pi!
仁（4）,i=2,3,…，n
（221.3
rfo©）二1,时
（二）=T
（fi（Q=（%-Qfi：
仏）-巧仁（Q,i=2,3,…，n
（221.4
'fOS）二l,flG）二0
（221.5fi（Q=（W—几）fi:
⑷一2fi：
仏）—巧fi：
（4）,i=2,3,…，n
f〃）=fnWf⑷=fnO）f⑺二f；G）,
（2・2.1・3）式两边都除以得f厂GJ得
（221.6）
匕1（扌L）=a1_A
几）二ctj_几一i二2,3,•••,n
=4仏）
亠⑷一尚和匚©）—15“
（2・2.1.4）式和（22.1.5）式两边除以f「仏）并注意下述关系
可得下面的递归关系
SOW]
（1）
1
3（4）二帀”叽⑷〃
占pi/仏）]J=2,3,…,n
（2.2.1.7）
j£仏）二0,勺仏）二0
J彳垄
Fi（Q二y*[（%-Q九⑷+1-jriTI
（2.2.1.8）
分别对
式和（2.2.1.8）式进
仏）],i=2,3,…,nIi（九）
行简单的推导，我们可以得到下面两个关
系式-需S需干⑴
我们称（2.2.1.6）式为改进的三项递归
进的计算特征多项式一阶导数和二阶导数的计算公式式・<221.7）式和（2.2.1.8）式分别为改
这些是我们在调用
Laguerre迭代时所需要计算的•根据文中知，除了个别极其特殊的情形，它们可
以避免计算过程中出现的上溢或下溢现象•但是，值得注意的是，如果对于某个
i（l勻Vn）使得耳⑷二0时，上述改进的各式就会出现中断现象•为了保证排除中断的发生，在计算中作如下处理
如果©
（1）=0,则令匕（扎）二Pig・i二1,2,…，n—1如果=0,则令En（）J=
这里E为机器精度•上述处理相当于对8或an引入一个小扰动IPiS或I注2.1：
我们由三项递归式（2.2.1.3）式经过改进得到（221.6piTsl。
式,在计算改进的非线性三项递归式（2.2.1.6）式时，同时可以得到在某个给定值A处，不大于该固定值A的特征值的个数，即q（IXQ（对所有1兰i兰n）的个数.
222Laguerre迭代的性质
对于不可约对称三对角矩阵Tn（如（2.1.1.1）式），其特征多项式
f仏）=detQn-几In）,（2.221）
的所有根为互不相等的实根。
设f（心的根为
并令扎。
二二和扎卄二。
对X亡仏，AS,Laguerre迭代L+（x）和L_（x）分别定
义如下
L』x）二x+ia（2.222）
/f（x）\,Lm/小/F（x）2zf〃（x）f（x）V
f（x）f（x）
L-（x）二X+4[zoooI
（q）—J（nT）[（n—1）（—3）"—n（d）]f（x）Yf（x）（222.3）
f（x）
并且有下面的命题成立.
命题2.1对仏m,kmSm=0,1,…，n,则有下述结果成立：
（D几mVL_（X）CXCL+（X）CAm+；
（2）存在常数C+和c_使得
当x趋于几m+时，有L+（X）-km+兰cjx-Am4t3；
当X趋于丸血时，有L_（X）m兰C」X-珀1。
根据命题2.1,如果令
Jk
=L；（x）=L+（L+（…L+（x）…））
k
—I
L（x）・〃））,
I
xT二Ll（x）二L_（L_（・
则有
（1）
（2）
（2.224）
（2）
（1）
一yVY（2.22.4）式说明当X忘仏m,）uJ时,L_（x）从X的左侧趋向于Am,L+（x）从X的
右侧趋向于文中指出只要X非常接近丸"（或非常接近Am半），一般调用
（2..2.2.3）式（或（2..2.2.2）式）两！
三次便可收敛于〈（或亦丿。
2.2.3病态接近特征值及改进的Laguerre迭代
根据文献，当多项式有重根时Laguerre迭代只是线性收敛•虽然不可约对称
三对角矩阵的特征值为两两互不相同的实数，但是,它们中的某些特征值也可能非常接近以致于在数值上无法区分，例如Wiikinson矩阵。
当矩阵有一个由r个非常接近的特征值构成的束时，在计算中被看作为一个r重特征值，在这种情形
下，Laguerre迭代的收敛率仅为线性的。
为了在上述情形一卜'加速迭代Laguerre收敛，对Laguerre迭代（2..2.22）式和
（2..2.23）式作重新定义，对正整数r（0Lr+（x）二X+d°（2.2.3.1）
（223.2）
Lr_（x）二X+4
（f\x））（n一ri）（f\x））2n（_匚厂7）_』[（n_Dp）
f（x）Vrf（x）f（x）
我们称（2.2.3.1）式和（2.2.3.2）式为改进的Laguerre迭代•我们需要注意的是:
当r=1时,（2.2.3.1）式和（2.2.3.2）式就分别是（2.2.2.2）式和（2.2.2.3）式•因此，我们在调用Laguerre迭代时，通常是先令r=1,然后计算r来调用（2.2。
3.1）式或（2.2.3.2）式。
由上面的两个式子，可以得到两个序列
X（k；=LkJX）=LrMLr4CLr"X）・・・））
—北—
X当X”为f仏）的r重根时，如果X并且f（兀）在区间（X,X”）内无根，则序列
X黑，k=1,2,-以三次收敛速度递增地收敛于xj类似地有，如果Yx并且f仏）在区间（x：
x）内无根，则序列x（1,k二1,2,…，以三次收敛速度递减地收敛于x*。
数值试验表明，即使在病态根束的情形下，（223.1）式和（223.2）式也能
达到三次收敛•对上述改进的Laguerre迭代1/（x）和1/（x）,有下面定理：
定理2.2设禹V扎2成・・・X-C~m,/~mG和任一正整数rcn,我们有
⑴如果-则AmCXCL+（X）〈L2+（X）〈•…（J如果-沏，则5<\）「V〃丄（心宀，为了方便起见，记％二亠和几卄二g。
注2.2:
设仏口如十）为Laguerre迭代的初始点并用它逼近特征值几m十，在某
些情形下，在离扎m+右边非常近的地方可能存在一组特征值，即入m十V扎m七<兀叶它们相对于x离的距离，非常靠近几皿卅。
在这种情形下，在达到Am寺的能够具有三次收敛速度的邻域前，实际计算需要花费许多次调用Laguerre!
代•这时，我们可以选择其它初始点来重新开始
Laguerre迭代或者调用改进的Laguerre迭代Lrdx）来加速收敛过程。
3求反对称三对角矩阵特征值的方法及分而治之算法
3.1反对称三对角矩阵的算法介绍
反对称矩阵特征值问题在许多实际问题中有广泛应用，如防潮陀螺仪系统的模型分析［231等，对这类具有特殊结构的矩阵，已有不少求特征值的算法［21'251它们基本上是把反对称矩阵特征值问题经过正交相似变换，如HousheuodI变换或
Gviens变换，转化为反对称三对角矩阵特征值问题•众所周知，在相似变换下不改变原矩阵的谱集•对一般矩阵，在采用QL算法之前,通常也是先把原矩阵化为上
Hessnehegr型，当矩阵为反对称矩阵时，根据其反对称性，其上Hessnehegr型矩阵
也就是三对角矩阵•因此，我们在下面讨论的反对称矩阵总是已三对角化了的反对称三对角矩阵，形如
-P1
0
-P„.1.1）
韵。
丿
不失一般性假设PjHO,j=1,2,…，n-1,式称之为不可约反对称三对角矩阵，简记为LATS,且近一步假设pj>0,j二1,2,…，nT,若存在pj=0,03.1.1求IAST特征值的QL算法
考虑IAST（3.1.1）,如果矩阵阶数n为奇数，由于IAST的特征值为互为共扼的纯虚数或零，因此可以通过位移为零的QL迭代把Jn化为
（3.1.1.1）
其中JnA为偶数阶矩阵，其特征值为互为共辘的纯虚数。
不失一般性，我们以后不妨假设矩阵的阶数n为偶数。
令Jn二J,计算j的特征值的
QL算法描述如下:
⑴计算J2+肾1的第n列
a—（Js+pl:
en=（0,0,…，Pn4,Pn20,Pl:
-P
⑵计算Househoulder矩阵po使得
PoQ.—|a6n,
这里Po=I—2wWT，其中w=（0,0,・・;，4pnjh_2^2,0,叮-p（-T丁/h
2
广邛nrpsz+（障_P②
h2邛2邛2+（障邛边7二茁-fnT-P2nJsign（i）二-signer一P_nj）
⑶计算J=PoTJPo,即
0-Pi
0*
0*0*
*0*0
0*0*
然后利用一系列旋转变换把J化为新的LAST
具体地，注意J只是比原矩阵J在（n,n-3）和（n-3,n）位置上多两个非零元素•用旋转矩阵
（In-4
-s
把P2RJP1P2的（n,n~4）和5-4,n）位置上多两个非零元素,这时在（n,n~5）和（n-5,n）位置上多两个非零元素。
类似地，通过旋转矩阵R,匕,・・・，Pn-3,把J化为
J二Pn-…PJPJJPlP2・・・Pn,此时#为一个新的LASTo此时对哒行压缩。
在对啲压缩后的矩阵执行双位
移QL算法，经过若干步上面的循环后，最终可以收敛于J的标准型冏
其中
‘01、
12二
.一10丿.
3.1.2由LAST的叉积的压缩矩阵求其特征值的方法
1971年，M.H.C.Paardekooper曾用于对称矩阵的Jacobi思想推广到反对称矩阵情形，给出了一种求反对称矩阵特征值得算法。
但他的算法仅考虑了反对称
性，未利用其纯虚数特征值共辘成对的性质，廉庆荣教授等利用这一性质，构造
出矩阵阶数为n=2N的jjT（其中J为LAST）的压缩矩阵（N阶不可约对称三对角矩阵），简记CIST（condenseirreduciblesymmetrietridiagonalmatrix）o卜面的定
理可以说明：
CIST的特征值的开方是对LAST的特征值的虚部。
有关
算法是在构造出JjT的CIST后，采用对称QL算法（带W订kinson位移）求岀CIST的特征值，然后在对其特征值开方求得LAST特征值的虚部。
不失一般性，仍考虑形如（3.1.1）式的n阶不可约反对称三对角矩阵J（不妨设n=2N为偶数,n为奇数时，N二h/2］）。
令二胪二妬），J必为特殊的五对角对称矩阵：
a*=肾,Ctnn=P：
」,ajj=Pj2」+Pj2,''=2-3n—
P2j,j-P2J-1.N出亠一\\,Pj二0（其
匕），
为不可约对称三对角矩阵，T二&ij）迂R'•环
何如2"jj\jH•pjHj,j=2,3,…，n-1,a*=0,其他.
定理3.1设J—JF如（3.121）所示，
（3.1.2.1）
若选择排列阵P二（Pij）€才满足
则P打”P二T©〒，其中T和T均
且aj=a：
i,2j,T=（%j空R'邓且
Gj=a2iA2j」0—j—
）•
注3.1当n=2N+1时，T为N+1阶，定理仍成立。
本文仍采用文［22］中的叫法，称T和T为r的压缩矩阵，其元素可以从r中直接得出，以T为例：
fa-f5nN—Gjj邛：
丄+時，j二1,2,…，N-1,L2N?
=-p2jJP2jj=l,2；・•,N-l,
心j卅,jRj‘j十二
（3.1.2.2）
严j二0,其他.
F面的定理给出定理3.1
升级会员 