【计算机】矩阵特征值计算

【计算机】矩阵特征值计算 9 矩阵特征值计算 在实际的工程计算中，经常会遇到求n阶方阵A的特征值(Eigenvalue)与特征向量(Eigenvector)的问。对于一个方阵A，如果数值λ使方程组 Ax=λx 即 (A-λI)x=0有非零解向量(Solution Vector)x，则称λ为方阵A的特征值，而非零向量x为n 特征值λ所对应的特征向量，其中I为n阶单位矩阵。 n 由于根据定义直接求矩阵特征值的过程比较复杂，因此在实际计算中，往往采取一些数值。本章主要介绍求一般方阵绝对值最大的特征值的乘幂(Power)法、求对称方阵特征值的雅可比法和单侧旋转(One-side Rotation)法以及求一般矩阵全部特征值的QR方法及一些相关的并行算法。 1.1 求解矩阵最大特征值的乘幂法 1.1.1 乘幂法及其串行算法 在许多实际问题中，只需要计算绝对值最大的特征值，而并不需要求矩阵的全部特征值。乘幂法是一种求矩阵绝对值最大的特征值的方法。记实方阵A的n个特征值为λ ii=(1,2, „,n)，且满足: ?λ???λ???λ??„??λ? 123n 特征值λ对应的特征向量为x。乘幂法的做法是:?取n维非零向量v作为初始向量;?ii0对于k=1,2, „,做如下迭代: u =Av v= u /?u? kk-1 kkk? 直至u,u,ε为止，这时v就是A的绝对值最大的特征值λ所对应的特k+11k，k1,, 征向量x。若v与v的各个分量同号且成比例，则λ=?u?;若v与v的各个分量异号1k-1k1k?k-1k且成比例，则λ= -?u?。若各取一次乘法和加法运算时间、一次除法运算时间、一次比较1k? 运算时间为一个单位时间，则因为一轮计算要做一次矩阵向量相乘、一次求最大元操作和一 22 次规格化操作，所以下述乘幂法串行算法21.1的一轮计算时间为n+2n=O(n)。 算法21.1 单处理器上乘幂法求解矩阵最大特征值的算法 输入:系数矩阵A，初始向量v，ε ××nn n1 输出:最大的特征值max Begin while (?diff?>ε) do (1)for i=1 to n do (1.1)sum=0 (1.2)for j= 1 to n do sum=sum+a[i,j]*x[j] end for (1.3)b[i]= sum end for (2)max=?b[1]? (3)for i=2 to n do if (?b[i]?>max) then max=?b[i]? end if end for (4)for i=1 to n do x[i] =b[i]/max end for (5)diff=max-oldmax , oldmax=max end while End 1.1.2 乘幂法并行算法 乘幂法求矩阵特征值由反复进行矩阵向量相乘来实现，因而可以采用矩阵向量相乘的数据划分方法。设处理器个数为p，对矩阵A按行划分为p块，每块含有连续的m行向量，这里，编号为i的处理器含有A的第im至第(i+1)m-1行数据，(i=0,1, „,p-1)，初m,n/p，， 始向量v被广播给所有处理器。 各处理器并行地对存于局部存储器中a的行块和向量v做乘积操作，并求其m维结果向量中的最大值localmax，然后通过归约操作的求最大值运算得到整个n维结果向量中的最大值max并广播给所有处理器，各处理器利用max进行规格化操作。最后通过扩展收集操 m作将各处理器中的维结果向量按处理器编号连接起来并广播给所有处理器，以进行下一次迭代计算，直至收敛。具体算法框架描述如下: 算法21.2 乘幂法求解矩阵最大特征值的并行算法 输入:系数矩阵A，初始向量v，ε ××nn n1 输出:最大的特征值max Begin 对所有处理器my_rank(my_rank=0,„, p-1)同时执行如下的算法: while (?diff?>ε) do /* diff 为特征向量的各个分量新旧值之差的最大值*/ (1)for i=0 to m-1 do /*对所存的行计算特征向量的相应分量*/ (1.1)sum=0 (1.2)for j= 0 to n-1 do sum=sum+a[i,j]*x[j] end for (1.3)b[i]= sum end for (2)localmax=?b[0]? /*对所计算的特征向量的相应分量 求新旧值之差的最大值localmax */ (3)for i=1 to m-1 do if (?b[i]?>localmax) then localmax=?b[i]? end if end for (4)用Allreduce操作求出所有处理器中localmax值的最大值max 并广播到所有处理器中 (5)for i=0to m-1 do /*对所计算的特征向量归一化 */ b[i] =b[i]/max end for (6)用Allgather操作将各个处理器中计算出的特征向量的分量的新值组合并广播到 所有处理器中 (7)diff=max-oldmax, oldmax=max end while End 若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时间，一轮迭代的计算时间为mn+2m;一轮迭代中，各处理器做一次归约操作，通信量为1， 4t(p,1)，(m，1)t(p,1)。因此乘幂法的一次扩展收集操作，通信量为m，则通信时间为sw 一轮并行计算时间为。 T,mn，2m，4t(p,1)，(m，1)t(p,1)psw MPI源程序请参见所附光盘。 1.2 求对称矩阵特征值的雅可比法 1.2.1 雅可比法求对称矩阵特征值的串行算法 若矩阵A=[a]是n阶实对称矩阵，则存在一个正交矩阵U，使得 ijTUAU=D 其中D是一个对角矩阵，即 ,0？0，，1,,,0？02,,D, ,,？？？？,,00？,,,n，， 这里λ (i=1,2,„,n)是A的特征值，U的第i列是与λ对应的特征向量。雅可比算法主要是通ii 过正交相似变换将一个实对称矩阵对角化，从而求出该矩阵的全部特征值和对应的特征向量。因此可以用一系列的初等正交变换逐步消去A的非对角线元素，从而使矩阵A对角化。设初等正交矩阵为R(p,q,θ)，其(p,p)( q,q)位置的数据是cosθ，(p, q)( q, p)位置的数据分别是-sinθ和sinθ(p ? q)，其它位置的数据和一个同阶数的单位矩阵相同。显然可以得到: TR(p,q,θ) R(p,q,θ)=I n T不妨记B= R(p,q,θ)AR(p,q,θ)，如果将右端展开，则可知矩阵B的元素与矩阵A的元素之间有如下关系: 2222 b = acosθ+asinθ+asin2θ b = asinθ+a cosθ-asin2θ ppppqqpqqqppqqpq b = ((a -a)sin2θ)/2+a cos2θ b= b pqqqpppqqp pq b= acosθ+asinθ b= -asinθ +acosθ pjpjqjqj pjqj b= acosθ+asinθ b= -asinθ +acosθ ipipiqiqipiq b= a ijij 其中1 ? i, j ? n且i,j ? p,q。因为A为对称矩阵，R为正交矩阵，所以B亦为对称矩阵。若要求矩阵B的元素b =0，则只需令 ((a -a)sin2θ)/2+a cos2θ=0，即: pqqqpppq ,apqtg2,, (a,a)2qqpp 在实际应用时，考虑到并不需要解出θ，而只需要求出sin2θ，sinθ和cosθ就可以了。 ? π/2，则可令 如果限制θ值在-π/2 < 2θ 1m ,,,,,ma,n(aa),wsgn(n)pqqqpp222，mn容易推出: w2sin,,, , sin2,,wcos,,1,sin, 22(1，1,w) 利用sin2θ，sinθ和cosθ的值，即得矩阵B的各元素。矩阵A经过旋转变换，选定的 非主对角元素a及a(一般是绝对值最大的)就被消去，且其主对角元素的平方和增加了pqqp 22，而非主对角元素的平方和减少了2a，矩阵元素总的平方和不变。通过反复选取主2apqpq 2元素a，并作旋转变换，就逐步将矩阵A变为对角矩阵。在对称矩阵中共有(n-n)/2个非主pq 对角元素要被消去, 而每消去一个非主对角元素需要对2n个元素进行旋转变换, 对一个元 素进行旋转变换需要2次乘法和1次加法，若各取一次乘法运算时间或一次加法运算时间为 一个单位时间，则消去一个非主对角元素需要3个单位运算时间，所以下述算法21.3的一 223轮计算时间复杂度为 (n-n)/2*2n*3=3n(n-1)=O(n)。 算法21.3 单处理器上雅可比法求对称矩阵特征值的算法 输入:矩阵A，ε ×nn 输出:矩阵特征值Eigenvalue Begin (1)while (max >ε) do (1.1) max=a[1,2] (1.2)for i=1 to n do for j= i+1 to n do if (?a[i,j]) ?>max) then max =?a[i,j] ?，p=i，q=j end if end for end for (1.3)Compute: m=- a[p,q]，n=(a[q,q]- a[p,p])/2，w=sgn(n)*m/sqrt(m*m+n*n)， si2=w，si1=w/sqrt(2(1+ sqrt(1-w*w)))，co1=sqrt(1-si1*si1)， b[p,p]= a[p,p]*co1*co1+ a[q,q]*si1*si1+ a[p,q]*si2， b[q,q]= a[p,p]*si1*si1+ a[q,q]*co1*co1- a[p,q]*si2， b[q, p]=0，b[p,q]=0 (1.4)for j=1 to n do if ((j ? p) and ( j ? q)) then (i)b[p,j]= a[p,j]*co1+ a[q,j]*si1 (ii)b[q,j]= -a[p,j]*si1 + a[q,j]*co1 end if end for (1.5)for i=1 to n do if((i ? p) and (i ? q)) then (i)b[i, p]= a[i,p]*co1+ a[i,q]*si1 (ii)b[i, q]= - a[i,p]*si1+ a[i,q]*co1 end if end for (1.6)for i=1 to n do for j=1 to n do a[i,j]=b[i,j] end for end for end while (2)for i=1 to n do Eigenvalue[i]= a[i, i] end for End 1.2.2 雅可比法求对称矩阵特征值的并行算法 串行雅可比算法逐次寻找非主对角元绝对值最大的元素的方法并不适合于并行计算。因此，在并行处理中，我们每次并不寻找绝对值最大的非主对角元消去，而是按一定的顺序将A中的所有上三角元素全部消去一遍，这样的过程称为一轮。由于对称性，在一轮中，A的下三角元素也同时被消去一遍。经过若干轮，可使A的非主对角元的绝对值减少，收敛于一个对角方阵。具体算法如下: 设处理器个数为p，对矩阵A按行划分为2p块，每块含有连续的m/2行向量，记 ，这些行块依次记为A,A, „,A，并将A与A存放与标号为i的处理器中。 m,n/p012p-12i2i+1，， 每轮计算开始，各处理器首先对其局部存储器中所存的两个行块的所有行两两配对进行旋转变换，消去相应的非对角线元素。然后按图21.1所示规律将数据块在不同处理器之间传送，以消去其它非主对角元素。 01234567 开 始 :(0，1)(2，3)(4，5)(6，7)第一步 :(0，3)(1，5)(2，7)(4，6)第二步 :(0，5)(3，7)(1，6)(2，4)第三步 :(0，7)(5，6)(3，4)(1，2)第四步 :(0，6)(7，4)(5，2)(3，1)第五步 :(0，4)(6，2)(7，1)(5，3)第六步 :(0，2)(4，1)(6，3)(7，5)第七步 :(0，1)(2，3)(4，5)(6，7) 图 错误～文档中没有指定样式的文字。21.1 p=4时的雅可比算法求对称矩阵特征值的数据 交换图 这里长方形框中两个方格内的整数被看成是所移动的行块的编号。在要构成新的行块配对时，只要将每个处理器所对应的行块按箭头方向移至相邻的处理器即可，这样的传送可以 在行块之间实现完全配对。当编号为i和j的两个行块被送至同一处理器时，令编号为i的行块中的每行元素和编号为j的行块中的每行元素配对，以消去相应的非主对角元素，这样在所有的行块都两两配对之后，可以将所有的非主对角元素消去一遍，从而完成一轮计算。由图 错误～文档中没有指定样式的文字。21.1可以看出，在一轮计算中，处理器之间要2p-1次交换行块。为保证不同行块配对时可以将原矩阵A的非主对角元素消去，引入变量b记录每个处理器中的行块数据在原矩阵A中的实际行号。在行块交换时，变量b也跟着相应的行块在各处理器之间传送。 处理器之间的数据块交换存在如下规律: 0号处理器前一个行块(简称前数据块，后一个行块简称后数据块)始终保持不变，将后数据块发送给1号处理器，作为1号处理器的前数据块。同时接收1号处理器发送的后数据块作为自己的后数据块。 p-1号处理器首先将后数据块发送给编号为p-2的处理器，作为p-2号处理器的后数据块。然后将前数据块移至后数据块的位置上，最后接收p-2号处理器发送的前数据块作为自己的前数据块。 编号为my_rank的其余处理器将前数据块发送给编号为my_rank+1的处理器，作为my_rank+1号处理器的前数据块。将后数据块发送给编号为my_rank-1的处理器，作为my_rank-1号处理器的后数据块。 为了避免在通信过程中发生死锁，奇数号处理器和偶数号处理器的收发顺序被错开。使偶数号处理器先发送后接收;而奇数号处理器先将数据保存在缓冲区中，然后接收偶数号处理器发送的数据，最后再将缓冲区中的数据发送给偶数号处理器。数据块传送的具体过程描述如下: 算法21.4 雅可比法求对称矩阵特征值的并行过程中处理器之间的数据块交换算法 输入:矩阵A的行数据块和向量b的数据块分布于各个处理器中 输出:在处理器阵列中传送后的矩阵A的行数据块和向量b的数据块 Begin 对所有处理器my_rank(my_rank=0,„, p-1)同时执行如下的算法: /*矩阵A和向量b为要传送的数据块*/ (1)if (my-rank=0) then /*0号处理器*/ (1.1)将后数据块发送给1号处理器 (1.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-1 do /* 将后数据块保存在缓冲区buffer中*/ for j=0 to n-1 do buffer[i-m/2,j]=a[i,j] end for end for (2.2)for i=m/2 to m-1 do buf [i-m/2] =b[i] end for (2.3)for i=0 to m/2-1 do /*将前数据块移至后数据块的位置上*/ for j=0 to n-1 do a[i+m/2,j]=a[i,j] end for end for (2.4)for i=0 to m/2-1 do b[i+m/2] =b[i] end for (2.5)接收p-2号处理器发送的数据块作为自己的前数据块 (2.6)将buffer中的后数据块发送给编号为p-2的处理器 end if (3)if ((my-rank=p-1) and ( my-rank mod 2=0)) then /*处理器p-1且其为偶数*/ (3.1)将后数据块发送给编号为p-2的处理器 (3.2)for i=0 to m/2-1 do /*将前数据块移至后数据块的位置上*/ for j=0 to n-1 do a[i+m/2,j]=a[i,j] end for end for (3.3)for i=0 to m/2-1 do b[i+m/2] =b[i] end for (3.4)接收p-2号处理器发送的数据块作为自己的前数据块 end if (4)if ((my-rank ? p-1) and ( my-rank ? 0)) then /*其它的处理器*/ (4.1)if (my-rank mod 2=0) then /*偶数号处理器*/ (i)将前数据块发送给编号为my_rank+1的处理器 (ii)将后数据块发送给编号为my_rank-1的处理器 (ii)接收编号为my_rank-1的处理器发送的数据块作为自己的前 数据块 (iv)接收编号为my_rank+1的处理器发送的数据块作为自己的后 数据块 else /*奇数号处理器*/ (v)for i=0 to m-1 do /* 将前后数据块保存在缓冲区buffer中*/ for j=0 to n-1 do buffer[i,j]=a[i,j] end for end for (vi)for i=0 to m-1 do buf[i] =b[i] end for (vii)接收编号为my_rank-1的处理器发送的数据块作为自己的前 数据块 (viii)接收编号为my_rank+1的处理器发送的数据块作为自己的 后数据块 (ix)将存于buffer中的前数据块发送给编号为my_rank+1的处理 器 (x)将存于buffer中的后数据块发送给编号为my_rank-1的处理器 end if end if End 各处理器并行地对其局部存储器中的非主对角元素a进行消去，首先计算旋转参数并ij 对第i行和第j行两行元素进行旋转行变换。然后通过扩展收集操作将相应的旋转参数及第i列和第j列的列号按处理器编号连接起来并广播给所有处理器。各处理器在收到这些旋转参数和列号之后，按0,1,„,p-1的顺序依次读取旋转参数及列号并对其m行中的第i列和第j列元素进行旋转列变换。 经过一轮计算的2p-1次数据交换之后，原矩阵A的所有非主对角元素都被消去一次。此时，各处理器求其局部存储器中的非主对角元素的最大元localmax，然后通过归约操作的求最大值运算求得将整个n阶矩阵非主对角元素的最大元max，并广播给所有处理器以决定是否进行下一轮迭代。具体算法框架描述如下: 算法21.5 雅可比法求对称矩阵特征值的并行算法 输入:矩阵A，ε ×nn 输出:矩阵A的主对角元素即为A的特征值 Begin 对所有处理器my_rank(my_rank=0,„, p-1)同时执行如下的算法: (a)for i=0 to m-1 do b[i] =myrank*m+i /* b记录处理器中的行块数据在原矩阵A中的实际行号*/ end for (b)while (?max?>ε) do /* max 为A中所有非对角元最大的绝对值*/ (1)for i=my_rank*m to (my_rank+1)*m-2 do /*对本处理器内部所有行两两配对进行旋转变换*/ for j=i+1 to (my_rank+1)*m-1 do (1.1)r=i mod m , t=j mod m /*i, j 为进行旋转变换行(称为主行)的 实际行号, r, t为它们在块内的相对行号*/ (1.2)if (a[r,j] ? 0) then /*对四个主元素的旋转变换*/ (i)Compute: f=-a[r,j], g=( a[t,j]- a[r,i])/2, h=sgn(g)*f/sqrt(f*f+g*g) , sin2=h , sin1=h/sqrt(2*(1+sqrt(1-h*h))) , cos1=sqrt(1-sin1*sin1)， bpp= a[r,i]*cos1*cos1+a[t,j]*sin1*sin1+a[r,j]*sin2， bqq= a[r,i]* sin1*sin1+a[t,j]* cos1*cos1-a[r,j]*sin2， bpq=0 , bqp=0 (ii)for v=0 to n-1 do /*对两个主行其余元素的旋转变换*/ if ((v ? i) and ( v ? j)) then br[v] = a[r,v]*cos1 + a[t,v]*sin1 a[t,v] = -a[r,v]* sin1 + a[t,v]* cos1 end if end for (iii)for v=0 to n-1 do if ((v ? i) and ( v ? j)) then a[r,v]=br[v] end if end for (iv)for v=0 to m-1 do /*对两个主列在本处理器内的其余元素的旋转变换*/ if (( v ? r) and ( v ? t)) then bi[v] = a[v,i]*cos1 + a[v,j]*sin1 a[v,j]= - a[v,i]* sin1 + a[v,j]* cos1 end if end for (v)for v=0 to m-1do if ((v ? r) and ( v ? t)) then a[v,i]= bi[v] end if end for (vi)Compute: ， a[r,i]=bpp , a[t,j]=bqq , a[r,j]=bpq , a[t,i]=bqp /*用temp1保存本处理器主行的行号和旋转参数*/ temp1[0]=sin1, temp1[1]=cos1， temp1[2]=(float)i ,temp1[3]= (float)j else (vii)Compute: temp1[0]=0, temp1[1]= 0， temp1[2]= 0 , temp1[3]= 0 end if (1.3)将所有处理器temp1中的旋转参数及主行的行号 按处理器编号连接起来并广播给所有处理器,存于temp2中 (1.4)current=0 (1.5) for v=1 to p do /*根据temp2中的其它处理器的旋转参数及主行的行号对相关的 列在本处理器的部分进行旋转变换*/ (i)Compute: s1=temp2[(v-1)*4+0] , c1=temp2[(v-1)*4+1]， i1=(int)temp2[(v-1)*4+2], j1=(int)temp2[(v-1)*4+3] (ii )if (s1、c1、 i1、 j1中有一不为0) then if (my-rank ? current) then for z=0 to m-1 do zi[z]=a[z,i1]*c1 + a[z,j1]*s1 a[ z,j1]=- a[z,i1]*s1 + a[z,j1]*c1 end for for z=0 to m-1 do a[z,i1]= zi[z] end for end if end if (iii)current=current+1 end for end for end for (2)for counter=1 to 2p-2 do /*进行2p-2次处理器间的数据交换, 并对交换后处理器中所有行两两配对进行旋转变换*/ (2.1)Data_exchange( ) /*处理器间的数据交换*/ (2.2)for i=0 to m/2-1 do for j=m/2 to m-1 do (i) if (a[i,b[j]] ? 0) then /*对四个主元素的旋转变换*/ ?Compute: ， f= -a[i,b[j]]，g=(a[j,b[j]]- a[i,b[i]])/2 h=sgn(g)*f/sqrt(f*f+g*g)， sin2=h, sin1=h/sqrt(2*(1+sqrt(1-h*h)))， cos1=sqrt(1-sin1*sin1)， bpp= a[i,b[i]]*cos1*cos1+ a[j,b[j]]*sin1*sin1+a[i,b[j]]*sin2， bqq= a[i,b[i]]* sin1*sin1+a[j,b[j]]* cos1*cos1-a[i,b[j]]*sin2， bpq=0, bqp=0 ?for v=0 to n-1 do /*对两个主行其余元素的旋转变换*/ if ((v ? b[i]) and ( v ? b[j])) then br[v] = a[i,v]*cos1 + a[j,v]*sin1 a[j,v] = -a[i,v]* sin1 + a[j,v]* cos1 end if end for ?for v=0 to n-1 do if ((v ? b[i]) and ( v ? b[j])) then a[i,v]=br[v] end if end for ?for v=0 to m-1 do /*对本处理器内两个主列的其余元素旋转变换*/ if ((v ? i) and ( v ? j)) then bi[v] = a[v, b[i]]*cos1 + a[v, b[j]]*sin1 a[v, b[j]] = - a[v, b[i]]* sin1 + a[v, b[j]]* cos1 end if end for ?for v=0 to m-1 do if ((v ? i) and ( v ? j)) then a[v, b[i]]=bi[v] end if end for ?Compute: a[i, b[i]]=bpp , a[j, b[j]]=bqq , a[i, b[j]]=bpq , a[j, b[i]]=bqp /*用temp1保存本处理器主行的行号和旋转参数*/ temp1[0]=sin1 , temp1[1]=cos1, temp1[2]=(float)b[i] , temp1[3]= (float)b[j] else ?Compute: temp1[0]=0,temp1[1]= 0, temp1[2]= 0,temp1[3]= 0 end if (ii)将所有处理器temp1中的旋转参数及主行的行号按处理 器编号连接起来并广播给所有处理器,存于temp2中 (iii)current=0 (iv)for v=1 to p do /*根据temp2中的其它处理器的旋转参数及主行的行号 对相关的列在本处理器的部分进行旋转变换*/ ?Compute: s1=temp2[(v-1)*4+0], c1=temp2[(v-1)*4+1], i1=(int)temp2[(v-1)*4+2], j1=(int)temp2[(v-1)*4+3] ?if (s1、c1、 i1、 j1中有一不为0) then if (my-rank ? current) then for z=0 to m-1 do zi[z]=a[z,i1]*c1 + a[z,j1]*s1 a[ z,j1]=- a[z,i1]s1 + a[z,j1]*c1 end for for z=0 to m-1 do a[z,i1]= zi[z] end for end if end if ?current=current+1 end for end for end for end for (3)Data-exchange( ) /*进行一轮中的最后一次处理器间的数据交换,使数据回到原来的位置*/ (4)localmax=max /*localmax为本处理器中非对角元最大的绝对值*/ (5)for i=0 to m-1 do for j=0 to n-1 do if ((m*my-rank+i) ? j) then if (?a[i,j]?> localmax) then localmax=?a[i,j]? end if end if end for end for (6)通过Allreduce操作求出所有处理器中localmax的最大值为新的max值 end while End 在上述算法中, 每个处理器在一轮迭代中要处理2 p个行块对, 由于每一个行块对含有 222=m/4 个行的配对, 即消去m/4 个非主对角m/2行, 因而对每一个行块对的处理要有(m/2) 元素. 每消去一个非主对角元素要对同行n 个元素和同列n 个元素进行旋转变换. 由于按行划分数据块, 对同行n 个元素进行旋转变换的过程是在本处理器中进行的. 而对同列n 个元素进行旋转变换的过程则分布在所有处理器中进行. 但由于所有处理器同时在为其它处理器的元素消去过程进行列的旋转变换，对每个处理器而言，对列元素进行旋转变换的处理总量仍然是n 个元素。 因此，一个处理器每消去一个非主对角元素共要对2n 个元素进行旋转变换。而对一个元素进行旋转变换需要2次乘法和1次加法，若取一次乘法运算时间或一次加法运算时间为一个单位时间，则其需要3个单位运算时间，即一轮迭代的计算时间 23为T,3×2 p ×2n ×m/4,3n/p;在每轮迭代中，各个处理器之间以点对点的通信方式相1 互错开交换数据2p-1次，通信量为mn+m，扩展收集操作n(n-1)/(2p)次，通信量为4，另外有归约操作1次，通信量为1，从而不难得出上述算法求解过程中的总通信时间为: T,[4t，m(n，1)t](4p,2)，[n(n,1)/p，2]t(p,1)，[2n(n,1)/p，1]t(p,1)2swsw T,T，T因此雅可比算法求对矩阵特征值的一轮并行计算时间为。我们的大量实验结果p12 说明，对于n阶方阵，用上述算法进行并行计算，一般需要不超过O(logn)轮就可以收敛。 MPI源程序请参见章末附录。 1.3 求对称矩阵特征值的单侧旋转法 1.3.1 单侧旋转法的算法描述 求解对称方阵特征值及特征向量的雅可比算法在每次消去计算前都要对非主对角元素选择最大元，导致非实际计算开销很大。在消去计算时，必须对方阵同时进行行、列旋转变换，这称之为双侧旋转(Two-side Rotation)变换。在双侧旋转变换中，方阵的行、列方向都有数据相关关系，使得整个计算中的相关关系复杂，特别是在对高阶矩阵进行特征值求解时，增加了处理器间的通信开销。针对传统雅可比算法的缺点，可以将双侧旋转改为单侧旋转，得出一种求对称矩阵特征值的快速算法。这一算法对矩阵仅实施列变换，使得数据相关关系仅在同列之间，因此方便数据划分，适合并行计算，称为单侧旋转法(One-side Rotation)。 T若A为一对称方阵， λ是A的特征值，x是A的特征向量，则有:Ax=λx。又A=A，所以T2TTTAAx=Aλx=λλx，这说明了 λ是AA的特征值，x是AA的特征向量 ，即AA的特征值是A的特征值的平方，并且它们的特征向量相同。 我们使用18.7.1节中所介绍的Givens旋转变换对A进行一系列的列变换，得到方阵Q T T 使其各列两两正交，即AV,Q，这里V为表示正交化过程的变换方阵。 因QQ=(AV)AV=VTTTAAV = diag(δ,δ, „,δ) 为n阶对角方阵，可见这里δ,δ, „,δ是矩阵AA的特征12n12n 2T,,,值，V的各列是AA的特征向量。由此可得: (i=1,2, „,n)。设Q的第i列为q,则 i ii TTq=δ, ||q||= 。因此在将A进行列变换得到各列两两正交的方阵Q后，qqq,,,,iii i2iiii 其各列的谱范数||q||即为A的特征值的绝对值。记V的第i列为v， AV,Q，所以Av = q。i2iii 又因为v 为A的特征向量，所以Av =λv，即推得λv= q。因此λ的符号可由向量 q及viiiiiiiiii 的对应分量是否同号判别，实际上在具体算法中我们只要判断它们的第一个分量是否同号即 可。若相同，则λ取正值，否则取负值。 i 求对称矩阵特征值的单侧旋转法的串行算法如下: 算法21.6 求对称矩阵特征值的单侧旋转法 输入:对称矩阵A，精度e 输出:矩阵A的特征值存于向量B中 Begin (1)while (p>ε) p=0 for i=1 to n do for j=i+1 to n do (1.1 )sum[0]=0, sum[1]=0, sum[2]=0 (1.2 )for k=1 to n do sum[0]= sum[0]+a[k,i]*a[k, j] sum[1]= sum[1]+ a[k,i]* a[k,i] sum[2]= sum[2]+ a[k, j]* a[k, j] end for (1.3 )if (?sum[0]?>ε) then (i) aa=2*sum[0] bb=sum[1]-sum[2] rr=sqrt(aa*aa+bb*bb) (ii) if (bb?0) then c=sqrt((bb+rr)/(2*rr)) s=aa/(2*rr*c) else s=sqrt((rr-bb)/(2*rr)) c=aa/(2*rr*s) end if (iii) for k=1 to n do temp[k]=c*a[k,i]+s*a[k,j] a[k,j]=[-s]*a[k,i]+c*a[k,j] end for (iv) for k=1 to n do temp1[k]= c*e[k,i]+s*e[k,j] e[k,j]=[-s]*e[k,i]+c*e[k,j] end for (v) for k=1 to n do a[k,i]= temp[k] end for (vi) for k=1 to n do e[k,i]= temp1[k] end for (vii) if (?sum[0]?>p) then p=?sum[0] ? end if end if end for end for end while (2)for i=1 to n do su=0 for j=1 to n do su=su+a[j,i]* a[j,i] end for B[j]=sqrt[su] if (e[0,j]*a[0,j]<0) then B[j]= - B[j] endif end for End 上述算法的一轮迭代需进行n(n-1)/2次旋转变换，若取一次乘法或一次加法运算时间为一个单位时间，则一次旋转变换要做3次内积运算而消耗6n个单位时间;与此同时，两列元素进行正交计算还需要12n个单位时间，所以奇异值分解的一轮计算时间复杂度为 3n(n-1)/2*(12 n +6n)= 9n(n-1) n=O(n)。 1.3.2 求对称矩阵特征值的单侧旋转法的并行计算 在求对称矩阵特征值的单侧旋转法的计算中, 主要的计算是矩阵的各列正交化过程.为了进行并行计算, 我们对n阶对称矩阵A按行划分为p块(p为处理器数)，每块含有连续的q行向量，这里q=n/p，第i块包含A的第i×q, „,(i+1)×q-1行向量，其数据元素被分配到第i号处理器上(i=0,1,„,p-1)。E矩阵取阶数为n的单位矩阵I，按同样方式进行数据划分，每块含有连续的q行向量。对矩阵A的每一个列向量而言，它被分成p个长度为q的子向量，分布于p个处理器之中，它们协同地对各列向量做正交计算。在对第i列与第j列进行正交计算时，各个处理器首先求其局部存储器中的q维子向量[i,j]、[i,i]、[j,j]的内积，然后通过归约操作的求和运算求得整个n维向量对[i,j]、[i,i]、[j,j]的内积并广播给所有处理器，最后各处理器利用这些内积对其局部存储器中的第i列及第j列q维子向量的元素做正交计算。算法21.7按这样的方式对所有列对正交化一次以完成一轮运算，重复进行若干轮运算，直到迭代收敛为止。在各列正交后, 编号为0的处理器收集各处理器中的计算结果，由0号处理器计算矩阵的特征值. 具体算法框架描述如下: 算法21.7 求对称矩阵特征值的并行单侧旋转算法 输入:对称矩阵A，精度e 输出:对称矩阵A的特征值存于向量B中 Begin 对所有处理器my_rank(my_rank=0,„, p-1)同时执行如下的算法: (a) while (not convergence) do (1)k=0 (2)for i=0 to n-1 do (2.1)for j=i+1 to n-1 do (i)sum[0]=0, sum[1]=0, sum[2]=0 (ii)计算本处理器所存的列子向量a[*, i]、a[*, j]的内积赋值给 ss[0] (iii)计算本处理器所存的列子向量a[*, i]、a[*, i]的内积赋值给 ss[1] (iv)计算本处理器所存的列子向量a[*, j]、a[*, j]的内积赋值给 ss[2] (v)通过规约操作实现: ?计算所有处理器中ss[0]的和赋给sum [0] ?计算所有处理器中ss[1]的和赋给 sum [1] ?计算所有处理器中ss[2]的和赋给sum[2] ?将sum [0]、sum [1]、sum [2] 的值广播到所有处理器中 e ) then /*各个处理器并行进行对i和j列正交化*/ (vi)if (?sum(0)?> ?aa=2*sum[0] ?bb=sum[1]-sum[2] ?rr=sqrt(aa*aa+bb*bb) ?if (bb>=0) then c=sqrt((bb+rr)/(2*rr)) s=aa/(2*rr*c) else s=sqrt((rr-bb)/(2*rr)) c=aa/(2*rr*s) end if ?for v=0 to q-1 do temp[v]=c*a[v, i]+s*a[v, j] a[v, j]=(-s)*a[v, i]+c*a[v, j] end for ?for v=0 to z-1 do temp1[v]=c*e[v, i]+s*e[v, j] e[v,j]=(-s)*e[v, i]+c*e[v, j] end for ?for v=0 to q-1 do a[v,i]=temp[v] end for ?for v=0 to z-1 do e[v,i]=temp1[v] end for ?k, k+1 end if end for end for end while 0号处理器从其余各处理器中接收子矩阵a和e，得到经过变换的矩阵A和I: /*0号处理器负责计算特征值*/ (b) for i=1 to n do (1)su=0 (2)for j=1 to n do su=su+A[j,i]* A[j,i] end for (3)B[j]=sqrt(su) (4)if (I[0,j]*a[0,j]<0) then B[j]= - B[j] endif end for End 上述算法的一轮迭代要进行n(n-1)/2次旋转变换, 而一次旋转变换需要做3次内积运算，若取一次乘法或一次加法运算时间为一个单位时间，则需要6q个单位时间，另外，还要对四列子向量中元素进行正交计算花费12q个单位时间，所以一轮迭代需要的计算时间=n(n-1)6q;内积计算需要通信，一轮迭代共做归约操作n(n-1)/2次，每次通信量为3，因T1 而通信时间T=。由此得出一轮迭代的并行计算时[2t(p,1)，3t(p,1)]*n(n,1)22sw 间T=T+T。 p12 MPI源程序请参见所附光盘。 1.4 求一般矩阵全部特征值的QR方法 1.4.1 QR方法求一般矩阵全部特征值的串行算法 在18.6节中，我们介绍了对一个n阶方阵A进行QR分解的并行算法，它可将A分解 T成A=QR的形式，其中R是上三角矩阵，Q是正交矩阵，即满足QQ=I。利用QR方法也可求方阵特征值，它的基本思想是: 首先令A= A，并对A进行QR分解，即: A= Q R;然后对A作如下相似变换:A 1111112-1-1= Q A Q= Q Q R Q= RQ，即将RQ相乘得到A;再对A进行QR分解得A = Q R，1111111111122222然后将RQ相乘得A。反复进行这一过程，即得到矩阵序列:A, A, „, A, A ,„，它22312mm+1们满足如下递推关系:A= Q R ;A= R Q (i=1,2, „, m,„)其中Q均为正交阵，R均为上iiii+1iiii ,,A三角方阵。这样得到的矩阵序列或者将收敛于一个以A的特征值为对角线元素的上三i 角矩阵，形如: ,***？*,,1,,,**？*,,2,,, *？*3,,,,？？？,,,n,, 或者将收敛于一个特征值容易计算的块上三角矩阵。 算法21.8 单处理器上QR方法求一般矩阵全部特征值的算法 输入:矩阵A,单位矩阵Q，ε ×nn 输出:矩阵特征值Eigenvalue Begin (1) while ((p>ε) and (count<1000)) do (1.1)p=0 , count = count +1 (1.2)QR_DECOMPOSITION(A,Q,R) /*算法18.11对矩阵A进行QR分解*/ (1.3)MATRIX_TRANSPOSITION(Q) /*由于算法18.11 对A进行QR分解 -1-1-1T-1-1，因而要使用算法18.1对矩阵Q进行转置，得到(Q)=( Q)=Q*/ 只得到Q (1.4)A=MATRIX_PRODUCT(R,Q) /*算法18.5将矩阵RQ相乘，得到新的 A 矩阵 */ (1.5)for i=1 to n do for j=1 to i-1do if (?a[i,j]?> p) then p=?a[i,j]? end if end for end for end while (2) for i=1 to n do Eigenvalue[i]=a[i,i] end for End 显然，QR分解求矩阵特征值算法的一轮时间复杂度为QR分解、矩阵转置和矩阵相乘 算法的时间复杂度之和。 1.4.2 QR方法求一般矩阵全部特征值的并行算法 并行QR分解求矩阵特征值的思想就是反复运用并行QR分解和并行矩阵相乘算法进行 ,,A迭代，直到矩阵序列收敛于一个上三角矩阵或块上三角矩阵为止。具体的并行算法描i 述如下: 算法21.9 QR方法求一般矩阵全部特征值的并行算法 输入:矩阵A,单位矩阵Q，ε ×nn 输出:矩阵特征值Eigenvalue Begin 对所有处理器my_rank(my_rank=0,„, p-1)同时执行如下的算法: (a) while ((p>ε) and (count<1000)) do (1)p=0 , count = count +1 /* 对矩阵A进行QR分解*/ (2)if(my_rank=0) then/*0号处理器*/ (2.1)for j=0 to m-2 do (i)for i=j+1 to m-1 do Turnningtransform( ) /*旋转变换*/ end for (ii)将旋转变换后的A和Q的第j 行传送到第1号处理器 end for (2.2)将旋转变换后的A和Q的第m-1行传送到第1号处理器 end if (3)if ((my_rank>0) and (my_rank<(p-1))) then /*中间处理器*/ (3.1)for j=0 to my_rank*m-1 do (i)接收左邻处理器传送来的A和Q子块中的第j 行 (ii)for i=0 to m-1 do Turnningtransform( ) /*旋转变换*/ end for (iii)将旋转变换后的A和Q子块中的第j 行传送到右邻处理器 end for (3.2)for j=0 to m-2 do (i)z=my_rank*m (ii)for i=j+1 to m-1 do Turnningtransform( ) /*旋转变换*/ end for (iii)将旋转变换后的A和Q子块中的第j 行传送到右邻处理器 end for (3.3)将旋转变换后的A和Q子块中的第m-1行传送到右邻处理器 end if (4)if (my_rank= (p-1)) then /*p-1号处理器*/ (4.1)for j=0 to my_rank*m-1 do (i)接收左邻处理器传送来的A和Q子块中的第j 行 (ii)for i=0 to m-1 do Turnningtransform( ) /*旋转变换*/ end for end for (4.2)for j=0 to m-2 do for i=j+1 to m-1 do f Turnningtransform( ) /*旋转变换*/ end for end for end if p-1号处理器对矩阵Q进行转置 (5) (6)p-1号处理器分别将矩阵R和矩阵Q划分为大小为m*M和M*m的p 块子阵，依次发送给0至p-2号处理器 (7)使用算法18.6对矩阵RQ进行并行相乘得到新的A矩阵 (8)for i=1 to n do /* 求出A中元素的最大绝对值赋予p*/ for j=1 to i-1do if (?a[i,j]?> p) then p=?a[i,j]? end if end for end for end while (b)for i=1 to n do Eigenvalue[i] =a[i,i] end for End 在上述算法的一轮迭代中，实际上是使用算法18.12、18.1、18.6并行地进行矩阵的QR 分解)矩阵的转置)矩阵的相乘各一次，因此, 一轮迭代的计算时间为上述算法的时间复杂 度之和。 MPI源程序请参见所附光盘。 1.5 小结 矩阵的特征值与特征向量在工程技术上应用十分广泛，在诸如系统稳定性问题和数字信号处理的K-L变换中，它都是最基本、常用的算法之一。本章主要讨论了矩阵特征值的几种基本算法，关于该方面的其它讲解与分析读者可参考文献[1]、[2]。 参考文献 [1]. 陈国良 编著(并行算法的设计与分析(修订版)(高等教育出版社，2002.11 [2]. 陈国良,陈崚 编著(VLSI计算理论与并行算法(中国科学技术大学出版社，1991 附录 求对称矩阵特征值的雅可比并行算法MPI源程序 源程序cjacobi.c #include "stdio.h" float sum; #include "stdlib.h" MPI_Status status; #include "mpi.h" #include "math.h" float sgn(float v) #include "string.h" { #define E 0.00001 float f; #define min -1 if (v>=0) f=1; #define intsize sizeof(int) else f=-1; #define charsize sizeof(char) return f; #define floatsize sizeof(float) } /*A是一个N*N的对称矩阵*/ #define A(x,y) A[x*N+y] void read_fileA() /*I是一个N*N的单位矩阵*/ { #define I(x,y) I[x*N+y] int i,j; #define a(x,y) a[x*N+y] FILE *fdA; #define e(x,y) e[x*N+y] #define b(x) b[x] fdA=fopen("dataIn.txt","r"); #define buffer(x,y) buffer[x*N+y] fscanf(fdA,"%d %d", &M, &N); #define buffee(x,y) buffee[x*N+y] if(M != N) { int M,N; puts("The input is error!"); int *b; exit(0); int m,p; } int myid,group_size; m=N/p; if (N%p!=0) m++; float *A,*I; A=(float*)malloc(floatsize*N*m*p); float *a,*e; I=(float*)malloc(floatsize*N*m*p); float max; for(i = 0; i < M; i ++) for(j = 0; j < M; j ++) e(i,j)=I(i,j); fscanf(fdA, "%f", A+i*M+j); } fclose(fdA); } printf("Input of file \"dataIn.txt\"\n"); else printf("%d\t %d\n",M, N); { for(i=0;iE) { for(v=0;vlmax) N,MPI_FLOAT,myid+1, lmax=fabs(a(i,j)); myid, } MPI_COMM_WORLD, &status); /*通过归约操作的求最大值运算求得将 MPI_Recv(&b(m/2),m/2, 整个n阶矩阵非主对角元素的最大 MPI_INT,myid+1,myid,元max，并广播给所有处理器以决 MPI_COMM_WORLD,定是否进行下一轮迭代*/ &status); MPI_Allreduce(&lmax,&max,1, MPI_Send(&(buffer(0,0)), MPI_FLOAT,MPI_MAX, m/2*N,MPI_FLOAT, MPI_COMM_WORLD); myid+1,myid+1, } /* while */ MPI_COMM_WORLD); /*0号处理器收集经过旋转变换的各子矩阵a， MPI_Send(&(buffee(0,0)), 得到非主对角元素全为0的结果矩阵 m/2*N,MPI_FLOAT, A*/ myid+1,myid+1, if (myid==0) MPI_COMM_WORLD); { A=(float*)malloc(floatsize*N*m*p); MPI_Send(&buf[0],m/2, I=(float*)malloc(floatsize*N*m*p); MPI_INT,myid+1, for(i=0;i