《矩阵特征值计算.docx》由会员分享，可在线阅读，更多相关《矩阵特征值计算.docx（31页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、9矩阵特征值计算在实际的工程计算中，经常会遇到求n阶方阵A的特征值(Eigenvalue)与特征向量(Eigenvector)的问题。对于一个方阵4如果数值X使方程组Ax=kx即(4“)x=0有非零解向量(Solution Vector则称/为方阵4的特征值，而非零向量x为特征值7所对应的特征向量，其中1“为阶单位矩阵。由于根据定义直接求矩阵特征值的过程比较复杂，因此在实际计算中，往往采取一些数值方法。本章主要介绍求一般方阵绝对值最大的特征值的乘幕(Power)法、求对称方阵特征值的雅可比法和单侧旋转(One-side Rotation)法以及求一般矩阵全部特征值的QR方法及一些相关的并行算法
2、。1.1 求解矩阵最大特征值的乘幕法1.1.1 乘幕法及其串行算法在许多实际问题中，只需要计算绝对值最大的特征值,而并不需要求矩阵的全部特征值。乘暴法是一种求矩阵绝对值最大的特征值的方法。记实方阵/的个特征值为九/=(1,2,且满足：|Ai|a2|a3|-|a|特征值不对应的特征向量为X*乘基法的做法是：取维非零向量V。作为初始向量；对于Q1,2,，做如下迭代：以=4味一|V*=M*/|w*|x直至)do(l)for i=l to n do(1.1)sum=Q(1.(2) r j=1 to h dosum=sum+aij*x/end for(1.(3) i= sumend for(2)max=
3、|/1(3)for i=2 to n doif (| bi|max) then max= bi end ifend for(4)for i= to n doxz=bi/maxend for(5)dff=/nax-oldmax, oldmax=maxend whileEnd1.1.2乘塞法并行算法乘幕法求矩阵特征值由反复进行矩阵向量相乘来实现，因而可以采用矩阵向量相乘的数据划分方法。设处理器个数为p,对矩阵A按行划分为p块，每块含有连续的m行向量，这里m编号为i的处理器含有力的第加至第(汁1)办1行数据，0=0,1,p-1),初始向量v被广播给所有处理器。各处理器并行地对存于局部存储器中的行块和
4、向量v做乘积操作，并求其用维结果向量中的最大值localmax,然后通过归约操作的求最大值运算得到整个维结果向量中的最大值max并广播给所有处理器，各处理器利用max进行规格化操作。最后通过扩展收集操作将各处理器中的“维结果向量按处理器编号连接起来并广播给所有处理器，以进行下一次迭代计算，直至收敛。具体算法框架描述如下：算法21.2乘幕法求解矩阵最大特征值的并行算法输入：系数矩阵4x“，初始向量输出：最大的特征值maxBegin对所有处理器my_rank(my_rank=O,pl)同时执行如卜.的算法：while (|diff|)do /* diff为特征向量的各个分量新旧值之差的最大值*/(
5、l)for z=0 to w-1 do /*对所存的行计算特征向量的相应分量*/(1.(1) W=0(1.(2) r /=0 to w-1 dosum=sum+aij*x/end for(3)hi=sumend for(2)localmax= h0|/*对所计算的特征向量的相应分量求新旧值之差的最大值localmax */(3)for z=l to m- doif (| bi|!ocalmax) then localmax= bi end ifend for(4)用Allreduce操作求出所有处理器中/oca加at值的最大值max并广播到所有处理器中(5)for z=Oto m- do /*
6、对所计算的特征向量归一化*/hi=bi/max end for(6)用Allgather操作将各个处理器中计算出的特征向量的分量的新值组合并广播到所有处理器中(7)di ff=max-oldmax, oldmax=max end whileEnd若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时间，一轮迭代的计算时间为m+2m；一轮迭代中，各处理器做一次归约操作，通信量为1,一次扩展收集操作，通信量为机，则通信时间为4G(而-1)+(加+1)*5-1)。因此乘嘉法的一轮并行计算时间为Tp =加+2加+44(后-1)+(加+ l)，w(p-l)。MPI源程序请参见所附光
7、盘。1.2求对称矩阵特征值的雅可比法1.2.1 雅可比法求对称矩阵特征值的串行算法若矩阵4=%是阶实对称矩阵，则存在一个正交矩阵U,使得utau=d其中。是一个对角矩阵，即0 A 00 A, A 0D= lM M 0 M00 A这里Z(i=l,2,是4的特征值，U的第i列是与尢对应的特征向量。雅可比算法主要是通过正交相似变换将一个实对称矩阵对角化，从而求出该矩阵的全部特征值和对应的特征向量。因此可以用系列的初等正交变换逐步消去4的非对角线元素，从而使矩阵4对角化。设初等正交矩阵为R(p,q,9),其(p,p)( g,q)位置的数据是cos。，(p, g)( g, p)位置的数据分别是-sin。
8、和sin6(p彳必，其它位置的数据和一个同阶数的单位矩阵相同。显然可以得到：R(p,q,0) TR(p,q,e)=I不妨记B= R(p,q,0)TAR(p,q,0),如果将右端展开，则可知矩阵B的元素与矩阵A的元素之间有如下关系：bpp = OppCOS2火ggSin之% apgSin26bqq =cos2-aP9sin20bpg =(4g -%p)sin20)/2+a* cos206卯=bMbpj=40cos 火 arsinebqj =-apjSind +a/COS。biIrCos火asin9bi(j=-a扪sin。+acosebij= ay其中IgiJW且ijrp,夕。因为力为对称矩阵，H
9、为正交矩阵，所以3亦为对称矩阵。若要求矩阵B的元素bpg =0,则只需令(aw -ap/,)sin20)/2+aw cos2=0,即：(aqq -app)2在实际应用时，考虑到并不需要解出仇而只需要求出sin2优sin。和cos。就可以了。如果限制6值在力/2e) do(1.(1) war=al,2(1.(2) r z=l to n dofor /= i+1 to n doif (| tz/j/)|max) then max = aij, p=i, q=j end if end forend for(1.(3) mpute:m=- ap,q, n=(aq,q- ap,p)/2, vv=sgn(
10、)*加乐qrt(机*机+*)， s，2=w, si 1=w/sqrt(2(1+ sqrt(1- w*w), col=sqrt(l-5zl*5zl), bp,p= ap,p*col*col+ap9q*si29bq,q= ap,p*5/l*5zl+ a*col*col- ap,q*si2, W%p=0, bp,q=0(1.(4) r7=1 to n doif(/Mp) and(D) then(i-=4pJ*col+ aqjsil(ii)贴力=-apJYsi+ a”*col end ifend for(1.(5) r Z=1 to n doif(z * p) and (/+ q) then(i)ft
11、z, p= az,p*col+ ai9q*silq=- ai,p*s+ ai,q*co end if end for(1.(6) r i=l to n doforj=/ to n doend for end for end while(2)for i=l to n doEigenvaluei= ai,/ end forEnd(1.(7) 比法求对称矩阵特征值的并行算法串行雅可比算法逐次寻找非主对角元绝对值最大的元素的方法并不适合于并行计算。因此，在并行处理中，我们每次并不寻找绝对值最大的非主对角元消去，而是按一定的顺序将”中的所有上三角元素全部消去一遍，这样的过程称为一轮。由于对称性，在一轮中
12、，/的下三角元素也同时被消去-一遍。经过若干轮，可使4的非主对角元的绝对值减少，收敛于一个对角方阵。具体算法如下：设处理器个数为p,对矩阵A按行划分为2p块，每块含有连续的m/2行向量，记 m这些行块依次记为山川，并将血，与4方+/存放与标号为，的处理器中。每轮计算开始，各处理器首先对其局部存储器中所存的两个行块的所有行两两配对进行旋转变换，消去相应的非对角线元素。然后按图21.1所示规律将数据块在不同处理器之间传送，以消去其它非主对角元素。图1.1 p=4时的雅可比算法求对称矩阵特征值的数据交换图这里长方形框中两个方格内的整数被看成是所移动的行块的编号。在要构成新的行块配对时，只要将每个处理
13、器所对应的行块按箭头方向移至相邻的处理器即可，这样的传送可以在行块之间实现完全配对。当编号为/和j的两个行块被送至同一处理器时，令编号为i的行块中的每行元素和编号为J的行块中的每行元素配对，以消去相应的非主对角元素，这样在所有的行块都两两配对之后，可以将所有的非主对角元素消去一遍，从而完成一轮计算。由图1可以看出，在轮计算中，处理器之间要2p-l次交换行块。为保证不同行块配对时可以将原矩阵X的非主对角元素消去，引入变量6记录每个处理器中的行块数据在原矩阵力中的实际行号。在行块交换时;变量方也跟着相应的行块在各处理器之间传送。处理器之间的数据块交换存在如下规律：0号处理器前一个行块(简称前数据块
14、，后一个行块简称后数据块)始终保持不变，将后数据块发送给1号处理器，作为1号处理器的前数据块。同时接收1号处理器发送的后数据块作为自己的后数据块。p-l号处理器首先将后数据块发送给编号为p-2的处理器，作为p-2号处理器的后数据块。然后将前数据块移至后数据块的位置上，最后接收p-2号处理器发送的前数据块作为自己的前数据块。编号为my_rank的其余处理器将前数据块发送给编号为my_rank+l的处理器，作为 my_rank+l号处理器的前数据块。将后数据块发送给编号为my_rank-l的处理器，作为 my_rank-l号处理器的后数据块。为了避免在通信过程中发生死锁，奇数号处理器和偶数号处理器
15、的收发顺序被错开。使偶数号处理器先发送后接收；而奇数号处理器先将数据保存在缓冲区中，然后接收偶数号处理器发送的数据，最后再将缓冲区中的数据发送给偶数号处理器。数据块传送的具体过程描述如下：算法21.4雅可比法求对称矩阵特征值的并行过程中处理器之间的数据块交换算法输入：矩阵4的行数据块和向量6的数据块分布于各个处理器中输出：在处理器阵列中传送后的矩阵4的行数据块和向量b的数据块Begin对所有处理器my_rank(my_rank=O,p-1)同时执行如下的算法：/*矩阵/和向量b为要传送的数据块*/(l)if(my-rank=O) then /*0号处理器*/(2.(1) 数据块发送给1号处理器
16、(2.(2) 1号处理器发送来的后数据块作为自己的后数据块 end if(2)if (my-rank=p-1) and ( my-rank mod 2丰0) then /*处理器 p-1且其为奇数*/(2.1)for i=m/2 to m- do /*将后数据块保存在缓冲区历娘厂中*/for j=0 to m-1 dobuJfferi-m/2j=aiJend forend for(2.(3) r i-ml2 to m- dobufi-m/2=hiend for(2.(4) r /=0 to m/2- do /*将前数据块移至后数据块的位置上*/for j=0 to w-l doai+m/2j=
17、aij end forend for(2.(5) r z=0 to m/2-1 dobi+m/2=biend for(2.(6) p-2号处理器发送的数据块作为自己的前数据块(2.(7) uffer中的后数据块发送给编号为p-2的处理器end if(3)if (my-rank=p-l) and ( my-rank mod 2=0) then /*处理器 p-1且其为偶数*/(3.(1) 数据块发送给编号为p-2的处理器(3.(2) r i=0 to m/2-1 do /*将前数据块移至后数据块的位置上*/for j=0 to n-1 doai+m/2j=aij end forend for(3
18、.(3) r z=0 to /w/2-1 dobi+m/2=biend for(3.(4) p2号处理器发送的数据块作为自己的前数据块end if(4)if (my-rankand ( my-rank，0) then /*其它的处理器*/(4.(1) (my-rank mod 2=0) then /*偶数号处理器*/将前数据块发送给编号为myank+l的处理器(ii)将后数据块发送给编号为my_rank-l的处理器(ii)接收编号为my_rank-l的处理器发送的数据块作为自己的前数据块(iv)接收编号为my_rank+l的处理器发送的数据块作为自己的后数据块else /*奇数号处理器*/(v
19、)for z=0 to m- do /*将前后数据块保存在缓冲区buffer中*/for7=0 to n- do bufferij=aijend forend for(vi)for z=0 to m- dobuAi=biend for(vii)接收编号为my_rank-l的处理器发送的数据块作为自己的前数据块(viii)接收编号为my_rank+l的处理器发送的数据块作为自己的后数据块(ix)将存于buffer中的前数据块发送给编号为my_rank+l的处理器(x)将存于6叨泊中的后数据块发送给编号为my_rank-l的处理器 end ifend ifEnd各处理器并行地对其局部存储器中的非主
20、对角元素与进行消去，首先计算旋转参数并对第/行和第j行两行元素进行旋转行变换。然后通过扩展收集操作将相应的旋转参数及第 i列和第j列的列号按处理器编号连接起来并广播给所有处理器。各处理器在收到这些旋转参数和列号之后，按0,1,p-1的顺序依次读取旋转参数及列号并对其m行中的第/列和第；列元素进行旋转列变换。经过一轮计算的2/1次数据交换之后，原矩阵/的所有非主对角元素都被消去一次。此时，各处理器求其局部存储器中的非主对角元素的最大元/。8加4X,然后通过归约操作的求最大值运算求得将整个阶矩阵非主对角元素的最大元并广播给所有处理器以决定是否进行下一-轮迭代。具体算法框架描述如下：算法21.5雅可
21、比法求对称矩阵特征值的并行算法输入：矩阵4*”，输出：矩阵/的主对角元素即为/的特征值Begin对所有处理器my_rank(my_rank=O,pl)同时执行如F的算法：(a)for z=0 to m- dobi=myrank*/w+f /* b记录处理器中的行块数据在原矩阵/中的实际行号*/ end for(b)while (| max e) do /* max为A中所有非对角元最大的绝对值*/(1)for z=my_rank*/n to (my_rank+do/*对本处理器内部所有行两两配对进行旋转变换*/forj=z+l to (my_rank+ do(l.(1) z mod m , t
22、=j mod m /*i, j为进行旋转变换行(称为主行)的实际行号,r, t为它们在块内的相对行号*/(l.(2) (ary+0) then /*对四个主元素的旋转变换*/(i)Compute:户。&J，g=(上月-。卜/)/2, h=sgn(g)*J7sqrV户g*g),sin2=h , sin 1=h/sqrt(Q*(l +s夕%(1*), cos 1=sqrt(-sin 11),hpp= ar9i*cos 1*cos 1+atj*sinl *szn i+arj*sin2， bqq= ar,z*cos 1*cos 1*sin2,bpq=0, bqp=0(ii)for v=0 to n-1
23、 do /*对两个主行其余元素的旋转变换*/if (v */) and ( v #j) thenbrv= ar,v*cos+ at,vsinat9v=-ar,v* sin+ a/,v* cos Iend ifend for(iii)for v=0 to h-1 doif (v*/) and (/) thenar,v=brvend if end for(iv)for v=0 to w-1 do/*对两个主列在本处理器内的其余元素的旋转变换*/if ( v r)and (。)thenbiv= av9i*cosl + ayjsinlavj=- av,i* sin+ HM* cosl end if e
24、nd for(v)for v=0 to Idoif (v * r) and ( v /) then av,z= bivend ifend for(vi)Compute:ar,ibpp , atj=bqq , arj=bpq , at,ibqp,/*用temp保存本处理器主行的行号和旋转参数*/ tempi0=5/w 1, temp 11=cos 1,temp 12=(float)z ,/ezwpl3=(float)/ else(vii)Compute:temp 10=0, tempi=0temp 12=0, temp3=0end if(L3)将所有处理器化肥1中的旋转参数及主行的行号按处理器编
25、号连接起来并广播给所有处理器，存于tempi中(1 A)current=0(1.5)for v=l to p do/*根据tempi中的其它处理器的旋转参数及主行的行号对相关的列在本处理器的部分进行旋转变换*/(i)Compute:s 1=temp2(v-1)*4+0, c 1=temp2(v-1)*4+1/z l=(int)/emp2(v-1)*4+2,j 1=(iiit)/efwp2(v-1)*4+3(ii)if(sl、cl、/I 中有一不为0)thenif (my-rank * current) thenfor z=0 to m- doziz=azj*cl + azjl*sa zjl=
26、az9il*sl +azl/l*cl end for for z=0 to m- dodfz,/l= ziz end for end if end if(i i i )current=current+1end forend forend for(2)for counter= to 2p-2 do/*进行2p-2次处理器间的数据交换，并对交换后处理器中所有行两两配对进行旋转变换*/(2.(1) ta_exchange()/*处理器间的数据交换*/(2.(2) r i=0 to m/2- dofor j=m/2 to m- do(1) if (HA河川彳0) then /*对四个主元素的旋转变换*
27、/Compute:户叩力/，g=(如力刃卜叩刖)/2, A=sg(g)3w7(/*#g*g)，sin2=h, sinl=h/sqrt(2*(+sqrt(1,cos =sqrt(-sin 1*sin 1),bpp= ai,bicos*cos+ aj,bj*sinlsinl-ai9bf*sin2, bqq=叩，如* sin 1*sin 1+a/,6/* cos 1cos 1-az,ft/*szw2, bpq=。, bqp=0for v=0 to /?-l do /*对两个主行其余元素的旋转变换*/ if (v丰 bi) and (h/) thenbrv= ai,v*cos+ a/,v*sin a/
28、,v=-ai,v* sinl + tz/,v* cos end ifend for for v=0 to n- doif (v声 hi) and (力/) then ai,v=brvend if end for for v=0 to m- do/*对本处理器内两个主列的其余元素旋转变换*/ if (v i) and ( vy) thenbiv= ay9 bi*cosl + av9 bj*sinav, b/=- av,6/* sin+ av,/* cos end ifend for for v=0 to /w-1 do if (v 声/) and ( v/y) then。匕 ft/=Wvend
29、if end forCompute:ai, bi=bpp , aj,6/=6gg , ai, bU=bpq , aj, bi=hc/p/*用temp保存本处理器主行的行号和旋转参数*/temp0=sin, temp 11=cos 1,temp 12=(float)/, temp3=(float)/?/elseCompute:temp 10=0,temp 11=0,temp2=0,/ezwpl3=0end if(ii)将所有处理器tempi中的旋转参数及主行的行号按处理器编号连接起来并广播给所有处理器，存于tempi中()current=Q(iv)for v=l top do/*根据2切2中的其
30、它处理器的旋转参数及主行的行号对相关的列在本处理器的部分进行旋转变换*/Compute:s 1=templ (v-1)*4+0,c 1=temp2(v-1)*4+1,/I=(int)/ewp2(v-l )*4+2,yl =(int)/ezp2(v-1)*4+3 if (si、cl、中有一不为0)then if (my-rank * current) thenfor z=0 to m- dozzz=az,zl*cl + azjl*s a zj=- az,isl + azj*cend forfor z=0 to m- do azyi= zizend forend ifend ifcurrent=c
31、urrent-1end forend forend forend for(3)Data-exchange()/*进行一轮中的最后一次处理器间的数据交换，使数据回到原来的位置*/(4)localmax=maxlocalmax为本处理器中非对角元最大的绝对值*/(5)for i=0 to tn- dofor j=0 to n-1 doif (/n*my-rank+z)* j) thenif (| aij| localmax) then localmax= end ifend ifend for end for(6)通过Allreduce操作求出所有处理器中localmax的最大值为新的max值 e
32、nd whileEnd在上述算法中，每个处理器在一轮迭代中要处理2 p个行块对，由于每一个行块对含有 mH行，因而对每一个行块对的处理要有(加/2)2=毋/4个行的配对，即消去加2/4个非主对角元素.每消去一个非主对角元素要对同行m个元素和同列n个元素进行旋转变换.由于按行划分数据块，对同行n个元素进行旋转变换的过程是在本处理器中进行的.而对同列n 个元素进行旋转变换的过程则分布在所有处理器中进行.但由于所有处理器同时在为其它处理器的元素消去过程进行列的旋转变换，对每个处理器而言，对列元素进行旋转变换的处理总量仍然是个元素。因此，一个处理器每消去一个非主对角元素共要对2个元素进行旋转变换。而对
33、一个元素进行旋转变换需要2次乘法和1次加法，若取一次乘法运算时间或一次加法运算时间为一个单位时间，则其需要3个单位运算时间，即一轮迭代的计算时间为7i =3X2p X2 x/4=33仍在每轮迭代中，各个处理器之间以点对点的通信方式相互错开交换数据2p1次，通信量为加+闭扩展收集操作(-1)/(20)次，通信量为4,另外有归约操作1次，通信量为1,从而不难得出上述算法求解过程中的总通信时间为：T2=4/v + m(n +1)/M.(4p-2)+n(n-)/p +2上诉-l)+2n(n-l)/p +1/K,(p -1)因此雅可比算法求对矩阵特征值的一轮并行计算时间为7；=1+72。我们的大量实验结
34、果说明，对于阶方阵，用上述算法进行并行计算，一般需要不超过O(/Og)轮就可以收敛。MPI源程序请参见章末附录。1.3求对称矩阵特征值的单侧旋转法1.3.1 单侧旋转法的算法描述求解对称方阵特征值及特征向量的雅可比算法在每次消去计算前都要对非主对角元素选择最大元，导致非实际计算开销很大。在消去计算时，必须对方阵同时进行行、列旋转变换，这称之为双侧旋转(Two-side Rotation)变换。在双侧旋转变换中，方阵的行、列方向都有数据相关关系，使得整个计算中的相关关系复杂，特别是在对高阶矩阵进行特征值求解时，增加了处理器间的通信开销。针对传统雅可比算法的缺点，可以将双侧旋转改为单侧旋转，得出一
35、种求对称矩阵特征值的快速算法。这一算法对矩阵仅实施列变换，使得数据相关关系仅在同列之间，因此方便数据划分，适合并行计算，称为单侧旋转法(One-side Rotation)o 若/为一对称方阵，入是4的特征值，x是Z的特征向量，则有：Ax=)jCo又A%,所以 zJEXXXr,这说明了好是49的特征值，x是的特征向量，即4%的特征值是N 的特征值的平方，并且它们的特征向量相同。我们使用18.7.1节中所介绍的Givens旋转变换对A进行一系列的列变换，得到方阵0使其各列两两正交，即这里/为表示正交化过程的变换方阵。因=dbg(%,%,6“)为阶对角方阵，可见这里叫,力，力是矩阵4的特征值，/的
36、各列是49的特征向量。由此可得：8.=2/(i=l,2,)。设。的第i列为，则 qq：=?，|%|2=扇?=何=。因此在将A进行列变换得到各列两两正交的方阵。后，其各列的谱范数的启即为A的特征值的绝对值。记/的第i列为%, AV=Q,所以Avi = qj。又因为修为力的特征向量，所以4%=入科，即推得入陷=%。因此无的符号可由向量%及均的对应分量是否同号判别，实际上在具体算法中我们只要判断它们的第一个分量是否同号即可。若相同，则无取正值，否则取负值。求对称矩阵特征值的单侧旋转法的串行算法如下：算法21.6求对称矩阵特征值的单侧旋转法输入：对称矩阵4,精度e输出：矩阵4的特征值存于向量8中Beg
37、in(l)while (pe )p=0for i=l to n dofor7=z+l to n do(1.1)5w/w0=0, suml=Oy sum2=0(1.2)for A=1 to n dosum0= sumO+ak,iak,j sum1= sum+ ak,i* ak,i sum2= sum2+ ak9j* ak,j end for(1.3)if (15z/w0| e ) then(i) aa=2*sum0 bb=sum-sum2 rr=sqrt(aa*aa-hb*bh)(ii)if (仍K) then c=sqrt(S6+zr)/(2*zr)5=aa/(2*rr*c) elses=sq
38、rt(zrbb)/(2*/r) c=aa/(2*rr*s) end if(iii) for k=l to n dotempk=c*ak,i+sakj akj=-s*akj+c*akj end for(iv) for k= to n dotempk= c*ekj+s*ekj ekj=-s*ek,i-c*ekj end for(v) for k= to n doak,i= tempkend for(vi) for k= to n doek9i= temp 1kend for(vii) if (15w/n0|p) then p= stun0| end if end ifend for end for end while(2)for z=l to n dosu=0fory=l to n dosu=su+aj,i* aj,i end for8/=sqrtMif (e0j*a0je) then /*各个处理器并行进行对i和j列正交化*/ aa=2*s?w0 bb=s?l -sum2 rr=sqrt(aa*aa+bb*bb)if (bb=0) thenc=sqrt(/4-rr)/(2*