矩阵特征值求解的分值算法12组1. 1矩阵计算的基本问题(1) 求解线性方程组的问题.即给定一个n 阶非奇异矩阵A 和n 维向量b ,求 一个n 维向量X,使得Ax =b (1.1.1 )(2) 线性最小二乘问题,即给定一个mx n 阶矩阵A 和m 维向量b ,求一个n 维向量使得 |A X -b | =min{ |Ay -比严 R n }(3) 矩阵的特征问题,即给定一个n 阶实(复)矩阵A ,求它的部分或全部特 征值以及对应的特征向量,也就是求解方程Ax = Z xA 的属于特征值A 的特征向量。
在工程上,矩阵的特征值具有广泛的应用,如大型桥梁或建筑物的振动问题: 机械和机件的振动问题;飞机机翼的颤振问题 ;无线电电子学及光学系统的电磁 振动问题;调节系统的自振问题以及声学和超声学系统的振动问题 .又如天文、地 震、信息系统、经济学中的一些问题都与矩阵的特征值问题密切相关。
在科学上,计算流体力学、统计计算、量子力学、化学工程和网络排队的马 尔可夫链模拟等实际问题,最后也都要归结为矩阵的特征值问题.由于特征值问 题在许多科学和工程领域中具有广泛的应用,因此对矩阵的特征值问题的求解理 论研究算法的开发软件的制作等是当今计算数学和科学与工程计算研究领域的 重大课题,国际上这方面的研究工作十分活跃。
1.2矩阵的特征值问题研究现状及算法概述对一个nxn 阶实(复)矩阵A,它的特征值问题,即求方程(I.1.3)式的非平凡 解,是数值线性代数的一个中心问题.这一问题的内在非线性给计算特征值带来 许多计算问题.为了求(1.1.3)式中的A , —个简单的想法就是显式地求解特征方 程det(A —几I) = 0除非对于个别的特殊矩阵,由于特征方程的系数不能够用稳定的数值方法由 行列式的计算来求得,既使能精确计算出特征方程的系数,在有限精度下,其特征 多项式f ") =det(A-ZJ)的根可能对多项式的系数非常敏感 能在理论上是有意义的,实际计算中对一般矩阵是不可行的 数较大,则行列式det(A -几I)的计算量将非常大;其次,根据 数大于四的多项式求根不存在一种通用的方法 ,基于上述原因,人们只能寻求其 它途径.因此,如何有效地!精确地求解矩阵特征值问题,就成为数值线性代数领 域的一个中心问题.目前,求解矩阵特征值问题的方法有两大类:一类称为变换方法,另一类称为 向X,(1.1.2 )(1.1.3 ) 一对解(4 X),其中R(C),x- R n (C n ),即A 为矩阵A 的特征值,X 为矩阵(121 ).因此,这个方法只 .首先,若矩阵A 的阶 Galois 理论,对于次量迭代方法.变换方法是直接对原矩阵进行处理,通过一系列相似变换,使之变换成一个易于求解特征值的形式,如Jacobi算法,Givens算法,QR算法等。
变换方法由于要存储矩阵元素,因而它只适用于求解中小型矩阵,它一般和向量迭代方法结合起来使用.向量迭代方法是通过一系列矩阵向量乘积而求得特征值和特征向量的.由于向量迭代方法可采用压缩存储技术,因而它适合于求大规模矩阵的特征值问题,尤其是大型稀疏矩阵的部分特征值和特征向量问题,如Lanczos 方法,Amoldi方法Qavidson方法等,现在这类问题仍是比较热的研究课题。
2分治方法的基础及理论研究2.1分治方法的概述考虑对称三对角矩阵T n的特征值问题TnX = Z X(2.1.1 ) 其中(2.1.2 )P n.1981年Cuppen提出一种求上述对称三对角矩阵T n所有特征值和特征向量的分而治之方法(divide 一and —conquer method).其基本思想是先将对称三对角矩阵T n分割为两个分别为k X k阶和(n - kp<(n - k)阶低阶对称三对角子矩阵T⑼和T⑴.T (0)和T⑴)可以用同样的方法也分别分割为两个更低阶的子矩阵,递归的采用这种分割技术可以把矩阵分割为一些能直接求出特征值的足够小的子矩阵(比如2阶或1阶矩阵),或者按照某种标准分割到适当阶数(如小于等于25 阶)后,结合其它求矩阵特征值的方法,如QR算法,求出其特征值。
在求出低阶矩阵特征值的基础上,开始胶合过程。
在胶合阶段,分割前的矩阵T1的特征值的求出(所谓的“治之”)是建立在其两个子矩阵T(0)和T(1)的特征值的基础上的,其中T(0)和T⑴.是在分割阶段由T1分割出的低阶子矩阵.随后的数值分析表明,Cuppen的方法存在着数值不稳定的危险,特别是当存在特征值束时,计算出的特征向量可能不正交。
Gu和Eisenstat对Cuppen的方法作了改进,极大地降低了数值不稳定的危险性。
Cuppen的方法在计算T n的特征值的同时也需要计算对应的特征向量,并且是在T ⑼和T ⑴的特征值和特征向量的基础上进行计算的 .根据文中,当用残量|TnX -A X |和正交性|x T X -ln |作为检验准确性的标准时,的标准时,二分法或多分法精确些.此外Cuppen 的分而治之方法要求矩阵乘积, 存储量为0(n 2),而二分法或多分法的存储量仅为 0(n),比前者少.应此,当只需 计算特征值时,通常选用后者.1987年Dongarra 和Sorensen 把分治思想应用到 求对称三对角矩阵特征值的并行计算,取得了不错的效果,再次引起了人们对分 治方法的极大关注。
分割胶合方法(split-merge method )是后来提出用分治方式求对称三对 角矩阵特征值T n 的方法,不同于分而治之方法在分割过程中采用矩阵的秩 1扰动, 它采用矩阵的秩2扰动.与二分法和多分法相似,在分割胶合方法中,特征值的计 算独立于特征向量的计算.如果需要计算特征向量,可采用反幕法。
由于分割胶 合方法计算特征值时采用具有三次收敛的Laguerre 迭代,数值试验表明,其计算速度和精度都明显优于二分法和多分法.并且文中给出的用于计算 Laguerre 迭代的非线性三项递归式可以避免上溢和下溢问题。
2.1.1分割策略分治方法的第一步就是把原来高阶的对称三对角矩阵特征值问题转化为两 个低阶对称三对角矩阵特征值问题 示如下(2.1.1.1 )不妨假设P j HO(j =1,2,…,n-1),即称T n 为不可约矩阵。
否则,若存在某些P j =0, 则T n 就可以约化为若干个低阶主子矩阵特征值问题,T n 的特征值就由这若干个低阶主子矩阵的特征值构成.当T n 为不可约对称三对角矩阵时,对其作如下分割Cuppen 的方法比二分法或多分法精确的多。
但文[3中,,如果用扎-kT n 作为衡量特征值准确性,即所谓的分割阶段.设对称三对角矩阵T n 表P n Jan丿fr (0))Tnd T ⑴丿(2.1.1.2)记T n =T n+E ,其中T (0)和T ⑴)分别为k 沢k 阶和(n —k )" n-k )或(n - k -1)x( n - k -1)阶对称三对角矩阵,通常k =[ n/2]。
(0扰动矩阵.此时有名1 p 1、、宀P 1+ J ■+J ■,T ⑴=+ J ■+ J ■Lk』 Ln』P k 4 叫丿P n 4 J 丿f oP knJan/p 1Q k卄用2P k卅、P 1+ ■T ")=P k H 4a k * -+ J+ J P k .+ +JJP n」P k A叫-p 」J 丿1扰动矩阵.此时有 T (0)P peP0 pe 2Pvv T时,其中 v =(o ;•-,T,…,0)且称此为秩(2)当Pe k eJ^称此为秩2(3)P k时,称此为秩3扰动矩阵.此时有务1 P 1g 七 P k 七P k七 叭七人们关注的问题是:对于上述三种分治策略T n 与T n 的特征值之间的关系.利用Hoffman-Wielandt 定理的结论,我们可以得到如下定理:定理2.1设A 和B 是两个n 阶Hermite 矩阵,它们的特征值分别是几’> 沁才,心 和已>^2 A" >巴,则其中H F 为Frobenius 范数。
根据上述定理,设T n 与T n 的特征值分别为Uh "和A 4 > A 2 A"工扎n,则有下面关系成立(2.1.1.4 )因此,只要(2.1.1.4)式右端越小,k j 就越接近A j ,即T n与T n 的特征值就越接近。
对于秩1和秩2扰动,T n 与T n 的特征值之间关系更详细的描述,将在后文 给出。
2.1.2胶合在矩阵T n 分割后,先求出矩阵T n 的特征值,即求出T ⑼和T ⑴的特征值。
剩下的工作就是如何由T n 的特征值出发求出T n的特征值,这就是胶合阶段主要任务。
在这个阶段可以采用不同的迭代方法,如割线法!Newton 迭代、Laguerre 迭 代以及路径跟踪等,以T n 的某个特征值或T n 的特征值构成的区间内的点为初 始点经过若千步迭代,最后收敛到T n 的某个特征值.本文中采用Laguerre 迭代从 特征区间提取特征值。
2.2 Laguerre 迭代2.2.1 Laguerre 迭代及其计算根据文中,对于不可约对称三对角矩阵 T n ,它的特征值或特征多项式f a ) =detT n -M n)的根全为互不相同的实数.因此,适合求多项式的根都是实n[J仏j -P j) ]2 < A -B |F (2.1.1.3)单根的具有三次收敛的Laguerre迭代L』x) = X +---------- 「…(2.2.1.1)-「罟卯V需j需)]非常适合用来从特征区间提取出T n 的特征值.这个方法早在13世纪就已经提出, 为了避免上述(2.2.1.3)式、(2.2.1.4)式和(221.5)式在计算中可能出现的上 溢或下溢问题,文中给出如下的改进的非线性三项递归式:令2治/2…,n,(2 . 2. 1.3)式两边都除以得fi^GJ 得近年来结合分治策略又重新焕发出生机 .文中首次对用Laguerre 迭代求不可约对称三对角矩阵的特征值的实用性进行了研究 题。
,而后又被用于求矩阵的奇异值问为了利用Laguerre 迭代求T n 的特征值,我们需要计算T n 的特征多项式f 仏)=det (T n -讥)以及其一阶导数f ‘仏)、 二阶导数f "仏).由文中知,计算这 些值的有效的方法是三项递归式.众所周知,当矩阵阶数比较大时,三项递归式 可能出现上溢或下溢问题.文中给出了一种改进的非线性的三项递归式 ,数值试 验表明,除极其个别的情形,改进的三项递归式可以避免上溢或下溢问题。
对称三对角矩阵T n 的特征多项式及其一、二阶倒数的计算公式设对称三对角矩阵T n 如 (221.1 )式所示,T k 表示其k 阶顺序主子式,表 示如下P 1匕P k j P k 4 «k >(2.2.1.2)f k (Q = det (T k -)・l k )为T k 的特征多项式,则特征多项式及其导数计算公式如下(三项递归式):rfo G ) =1,f 1(Q=% 八t f i 仏)=(% -Q f i 4仏)—P i !仁⑷,i =2,3,…,n(221.3r f 。