大型实对称矩阵特征值的数值解法

第18卷第3期2002年9月

北京建筑工程学院学报

JOURNALOFBEIJINGINSTITUTEOFCIVILENGINEERINGANDARCHTECTURE

Vol.18No.3Jun.2002

文章编号:1004-6011(2002)04-0058-03

大型实对称矩阵特征值的数值解法

刘长河 寿玉亭 马龙友 代西武 刘世祥

(基础部,北京,100044)

摘 要:本文介绍计算稀疏大型实对称矩阵特征值的方法 Davidson方法。并把它与矩阵的拟上三角化方法结合起来,得到一种求一般大型实对称矩阵特征值的方法。关键词:特征值;三对角矩阵;Davidson方法.中图分类号:O241 6 文献标识码:A

0 引言

早在19世纪中叶,Jacobi就给出了求实对称矩阵的特征值的数值解法 经典Jacobi方法。此后,人们研究出关于矩阵的特征值及特征向量的许多新的数值解法。如幂方法,Krylov方法,Lanczos方法,Frame方法,QR方法等。特别是近几十近来,随着现代科学的发展,不断地提出一些大型的矩阵计算问题。同时,计算机技术的飞跃、计算能力的增强,使解决这些问题成为现实。Davidson方法便是目前应用较广的计算稀疏大型实对称矩阵特征值的方法之一。本文利用矩阵的拟上三角化方法,对David son方法进行改进,使其能够用来计算一般的大型实对称矩阵的特征值。

r)

(1)若a(ir(i=r+2,r+3,!n)全为零,则令

A

(r+1)

=A

(r)

,转(5);否则转(2)。

(2)计算

dr=

(r)2

(air)=i+1

r)(r)

cr=-sign(a(r+1,r)dr,(若ar+1,r=0,则取cr=

dr)

(r)hr=c2c-crar+1

r)(r)(r)Tn

(3)令ur=(0,!,0,a(r+1,r-cr,ar+2,r,!,anr) R

(4)计算

pr=A

(r)T

ur/hr

qr=A(r)Tur/hrtr=pTrur/hr r=qr-trurA

(r+1)

(r)

T

T

1 矩阵的拟上三角化

定义:实矩阵Vk称为拟上三角矩阵(又称上Hessenberg矩阵),如果当i>j+1时,aij=0.任何实矩阵A=[aij]n n都可经过相似变换变为拟上三角矩阵,这个过程称为拟上三角化。对实矩阵A施行拟上三角化的具体算法如下:

记A

(1)

=A- rur-urpr

(5)继续

此算法执行完毕后,得拟上三角矩阵A(n-似与A。

a11c1

A(n-1)

2)

a(1,n-n-2

1)a(1,n-n-1

n-a(1n

1)

1)

相

!!!!

n-2)a(n-2,n-2

1)a(n,n-n-1

n-a(nn

1)

=(1)

=A,并记A

(r)

=A=[aij]n

(r)

n

cn-

2

对于r=1,2,!n-2,执行

收稿日期:2002-04-11

基金项目:本课题受到北京建筑工程学院博士启动基金资助。

:),,,

第4期 刘长河 寿玉亭等:大型实对称矩阵特征值的数值解法

59

如果A为实对称阵,A(n-1)是对称的三对角阵:

a11c1

A(n-1)=

!

c1

2)a(22

Rn中的标准正交基q1,q2,!,qn,算法如下: ProcedureMGS(b1,b2,!,bn)

fori=1,2,!ndo bi:=bi;

(1)

c2 !an-(n-1)

1,n-11)a(n,n-n-1

!cn-2

an,n-n-a(nn

(n-

1)

1

1)

endforforj=1,2,!,n qj:=b(jj)/#b(jj)#; fori=j+1,!,n

(j)

b(ij+1):=b(ij)-qj(q*jbi);

(2)

由于相似矩阵具有相同的特征值,因此,要求矩阵A的特征值,只需求矩阵A

(n-1)

的特征值即可。

在此过程中,乘法运算的次数大致为O(n3)。

endforendfor

Davidson算法还可用来求矩阵A的多个特征值。设l n,求矩阵A的l个最大(或最小)特征及其相应的特征向量的算如下。其中,m是一个给定的数,用以控制矩阵Vk的最大维数。

选择初始标准正交矩阵V1:=[v1,v2,!,vl,] Rnxl;fork=1,!,do

1 计算矩阵Wk:=AVk;

2 计算Rayleigh矩阵Hk:=VkWk;

3 计算Hk的l个最大(或最小)特征值 k,i,及其相应的特征向量yk,i,(1%i%l);4 计算l个Ritz向量xk,i:=Vkyk,i,(1%i%l);

5 计算l个残量rk,i:= k,i-Wkyk,i,(1%i%l);

如果#rk,i#<!,结束;否则,

6 计算l个新的方向向量tk,i:=ck,irk,i,(1%i%l);

7 ifdim(Vk)%m-lthen Vk+1:=MGS(Vk,tk,l!,tk,l) else

Vk+1:=MGS(Vk,l!xk,l,tk,l,tk,l); endifendfor

计算结束时的 k,i和xk,i分别是所要求的l个最大(或最小)特征值及相应的特征向量。算法中的Ck,i被称为proconditioner.可取Ck,i=diag( k,i,l,

!,

-1

i,),,i,,M-aj,j|

t

2 Davidson算法

Davidson算法最初用来求解稀疏矩阵A的最大(或最小)特征值及相应的特征向量。算法如下:选择一个单位向量v1:V1=[v1];

for k=1,!, do

t

计算迭代矩阵HK:=VkAVk;

计算矩阵Hk的最大(或最小)特征值与特征向量( K,yk);

计算相应的Ritz向量:xk:=VKyk 计算残量:rk:=( kI-A)xk; 如果#rk#<!,结束;否则, 计算向量tK+1:=( kI-D)-1rk; 将[Vk;tk+1]标准正交化,记为Vk+1;endfor

当Vk的维数变得过大时,可以用最后一个Ritz向量作为初始向量,重新开始。计算结束时的 k,xK分别是所要求的最大(或最小)特征值及相应的特征向量。

本文中,[!]表示矩阵,# #表示Euclidean范数,向量x={x1,x2,!xn}T的范数#x#=

2

xi.!1

为一事先给定的非常小的数,用以控制程序的结束。计算矩阵Hk的最大(或最小)特征值与特征向量( K,yk)可用带双步位移的QR方法(参阅文献[4])。标准正交化过程采用改进的Gram-Schmidt正交化方法[1]。具体如下:

b,b,

60

-1

的Ck,j=( 可参阅文献[5]。k,j-M)

北京建筑工程学院学报 第18卷

参考文献:

[1] B.N.PARLETT.Thesymmetriceigenvalueproblem.Prentice-Hall[M] EnglewoodCliffs,NJ,1980.

[2] 张凯院,徐仲 数值代数[M] 西安:西北工业大学出版社,

2000.

3 大型稀疏矩阵特征值的计算

设A是一个大型实对称阵,现求其l个最大(或最小)特征值。这里,l n.由第一部分知,经过拟三角化,可把A化为三对角矩阵。由于n很大,此三对角矩阵就是一个稀疏矩阵且与A相似,再利用Davidson算法即可求出此稀疏矩阵的l个最大(或最小)特征值即A的l个最大(或最小)特征值。

[3] M.CROUZEIX,B.PHILIPPE,ANDM.SADKANE.TheDavid

sonMethod[J] PUT.1994,(15)PP62-76.[4] 颜庆津 数值分析[M] 北京:北京航空航天大学出版社,

2000.

[5] R.B.MORGANANDD.S.SCOTT.GeneralizationsofDavid

son&smethodforcomputingeigenvaluesofsparsesymmetricma trices[J] put.(1986),7,PP817-825.

MethodsofComputingtheEigenvaluesof

LargeRealandSymmetricMatrix

LiuChanghe ShouYuting MaLongyou DaiXiwu LiuShixiang

(Dept.ofBasicSciences,Beijing100044)

Abstract:Inthispaper,DavidsonMethod,amethodofcomputingtheeigenvaluesoflarge、binedwiththemethodofquasi-tridiagonalizingthismethodcanbeusedtocomputetheeigenvaluesoflarge、realandsymmetricmatrix.

Keywords:eigenvalue;tridiagonalmatrix;DavidsonMethod