并行计算矩阵特征值计算.docx
《并行计算矩阵特征值计算.docx》由会员分享，可在线阅读，更多相关《并行计算矩阵特征值计算.docx（64页珍藏版）》请在冰豆网上搜索。
并行计算矩阵特征值计算
9矩阵特征值计算
在实际的工程计算中，经常会遇到求n阶方阵A的特征值（Eigenvalue）与特征向量（Eigenvector）的问题。
对于一个方阵A，如果数值λ使方程组
Ax=λx
即（A-λIn）x=0有非零解向量（SolutionVector）x，则称λ为方阵A的特征值，而非零向量x为特征值λ所对应的特征向量，其中In为n阶单位矩阵。
由于根据定义直接求矩阵特征值的过程比较复杂，因此在实际计算中，往往采取一些数值方法。
本章主要介绍求一般方阵绝对值最大的特征值的乘幂（Power）法、求对称方阵特征值的雅可比法和单侧旋转（One-sideRotation）法以及求一般矩阵全部特征值的QR方法及一些相关的并行算法。
1.1求解矩阵最大特征值的乘幂法
1.1.1乘幂法及其串行算法
在许多实际问题中，只需要计算绝对值最大的特征值，而并不需要求矩阵的全部特征值。
乘幂法是一种求矩阵绝对值最大的特征值的方法。
记实方阵A的n个特征值为λii=（1,2,…,n），且满足：
│λ1│≥│λ2│≥│λ3│≥…≥│λn│
特征值λi对应的特征向量为xi。
乘幂法的做法是：
①取n维非零向量v0作为初始向量；②对于k=1,2,…,做如下迭代：
直至u
k+1∞-uk
uk=Avk-1vk=uk/║uk║∞
<ε为止，这时vk+1就是A的绝对值最大的特征值λ1所对应的特征向
∞
量x1。
若vk-1与vk的各个分量同号且成比例，则λ1=║uk║∞；若vk-1与vk的各个分量异号且成比例，则λ1=-║uk║∞。
若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时间，则因为一轮计算要做一次矩阵向量相乘、一次求最大元操作和一次规格化操作，所以下述乘幂法串行算法21.1的一轮计算时间为n2+2n=O（n2）。
算法21.1单处理器上乘幂法求解矩阵最大特征值的算法
输入：
系数矩阵An×n，初始向量vn×1，ε
输出：
最大的特征值max
Begin
while（│diff│>ε）do
（1）fori=1tondo
（1.1）sum=0
（1.2）forj=1tondosum=sum+a[i,j]*x[j]endfor
（1.3）b[i]=sum
endfor
（2）max=│b[1]│（3）fori=2tondo
if（│b[i]│>max）thenmax=│b[i]│endifendfor
（4）fori=1tondo
x[i]=b[i]/max
endfor
（5）diff=max-oldmax,oldmax=max
endwhile
End
1.1.2乘幂法并行算法
乘幂法求矩阵特征值由反复进行矩阵向量相乘来实现，因而可以采用矩阵向量相乘的数据划分方法。
设处理器个数为p，对矩阵A按行划分为p块，每块含有连续的m行向量，这里m=⎡n/p⎤，编号为i的处理器含有A的第im至第（i+1）m-1行数据，（i=0,1,…,p-1），初始向量v被广播给所有处理器。
各处理器并行地对存于局部存储器中a的行块和向量v做乘积操作，并求其m维结果向量中的最大值localmax，然后通过归约操作的求最大值运算得到整个n维结果向量中的最大值max并广播给所有处理器，各处理器利用max进行规格化操作。
最后通过扩展收集操作将各处理器中的m维结果向量按处理器编号连接起来并广播给所有处理器，以进行下一次迭代计算，直至收敛。
具体算法框架描述如下：
算法21.2乘幂法求解矩阵最大特征值的并行算法
输入：
系数矩阵An×n，初始向量vn×1，ε
输出：
最大的特征值max
Begin
对所有处理器my_rank（my_rank=0,…,p-1）同时执行如下的算法:
while（│diff│>ε）do/*diff为特征向量的各个分量新旧值之差的最大值*/
（1）fori=0tom-1do/*对所存的行计算特征向量的相应分量*/
（1.1）sum=0
（1.2）forj=0ton-1do
sum=sum+a[i,j]*x[j]
endfor
（1.3）b[i]=sum
endfor
（2）localmax=│b[0]│/*对所计算的特征向量的相应分量求新旧值之差的最大值localmax*/
（3）fori=1tom-1do
if（│b[i]│>localmax）thenlocalmax=│b[i]│endifendfor
（4）用Allreduce操作求出所有处理器中localmax值的最大值max
并广播到所有处理器中
（5）fori=0tom-1do/*对所计算的特征向量归一化*/
b[i]=b[i]/max
endfor
（6）用Allgather操作将各个处理器中计算出的特征向量的分量的新值组合并广播到所有处理器中
（7）diff=max-oldmax,oldmax=max
endwhileEnd
若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时
间，一轮迭代的计算时间为mn+2m；一轮迭代中，各处理器做一次归约操作，通信量为1，
一次扩展收集操作，通信量为m，则通信时间为4ts（
p-1）+（m+1）tw（p-1）。
因此乘幂法的
一轮并行计算时间为Tp=mn+2m+4ts（
MPI源程序请参见所附光盘。
p-1）+（m+1）tw（p-1）。
1.2求对称矩阵特征值的雅可比法
1.2.1雅可比法求对称矩阵特征值的串行算法
若矩阵A=[aij]是n阶实对称矩阵，则存在一个正交矩阵U，使得
UTAU=D
其中D是一个对角矩阵，即
⎡λ1
⎢
D=⎢0
⎢
⎢
⎢⎣0
0
λ2
0
0⎤
⎥
0⎥
⎥
⎥
λn⎥⎦
这里λi（i=1,2,…,n）是A的特征值，U的第i列是与λi对应的特征向量。
雅可比算法主要是通过
正交相似变换将一个实对称矩阵对角化，从而求出该矩阵的全部特征值和对应的特征向量。
因此可以用一系列的初等正交变换逐步消去A的非对角线元素，从而使矩阵A对角化。
设初等正交矩阵为R（p,q,θ），其（p,p）（q,q）位置的数据是cosθ，（p,q）（q,p）位置的数据分别是-sinθ和sinθ（p≠q），其它位置的数据和一个同阶数的单位矩阵相同。
显然可以得到：
R（p,q,θ）TR（p,q,θ）=In
不妨记B=R（p,q,θ）TAR（p,q,θ），如果将右端展开，则可知矩阵B的元素与矩阵A的元素之间有如下关系：
bpp=appcos2θ+aqqsin2θ+apqsin2θbqq=appsin2θ+aqqcos2θ-apqsin2θ
bpq=（（aqq-app）sin2θ）/2+apqcos2θbqp=bpq
bpj=apjcosθ+aqjsinθbqj=-apjsinθ+aqjcosθbip=aipcosθ+aiqsinθbiq=-aipsinθ+aiqcosθbij=aij
其中1≤i,j≤n且i,j≠p,q。
因为A为对称矩阵，R为正交矩阵，所以B亦为对称矩阵。
若
要求矩阵B的元素bpq=0，则只需令（（aqq-app）sin2θ）/2+apqcos2θ=0，即：
tg2θ=
-apq
（aqq-app）2
在实际应用时，考虑到并不需要解出θ，而只需要求出sin2θ，sinθ和cosθ就可以了。
如果限制θ值在-π/2<2θ≤π/2，则可令
qq
m=-apq,
n=1（a
-app）,
w=sgn（n）m
容易推出：
sin2θ=w,
2
sinθ=
2（1+
w,
1-w2）
cosθ=
m2
1-sin2θ
+n2
利用sin2θ，sinθ和cosθ的值，即得矩阵B的各元素。
矩阵A经过旋转变换，选定的非主
pq
对角元素apq及aqp（一般是绝对值最大的）就被消去，且其主对角元素的平方和增加了2a2，
而非主对角元素的平方和减少了2a2pq，矩阵元素总的平方和不变。
通过反复选取主元素apq，并作旋转变换，就逐步将矩阵A变为对角矩阵。
在对称矩阵中共有（n2-n）/2个非主对角元素要被消去,而每消去一个非主对角元素需要对2n个元素进行旋转变换,对一个元素进行旋转变换需要2次乘法和1次加法，若各取一次乘法运算时间或一次加法运算时间为一个单位时间，则消去一个非主对角元素需要3个单位运算时间，所以下述算法21.3的一轮计算时间复杂度为（n2-n）/2*2n*3=3n2（n-1）=O（n3）。
算法21.3单处理器上雅可比法求对称矩阵特征值的算法
输入：
矩阵An×n，ε
输出：
矩阵特征值Eigenvalue
Begin
（1）while（max>ε）do（1.1）max=a[1,2]（1.2）fori=1tondo
forj=i+1tondo
if（│a[i,j]）│>max）thenmax=│a[i,j]│，p=i，q=jendifendfor
endfor
（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）forj=1tondo
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
endifendfor
（1.5）fori=1tondo
End
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
endifendfor
（1.6）fori=1tondoforj=1tondo
a[i,j]=b[i,j]
endforendfor
endwhile
（2）fori=1tondo
Eigenvalue[i]=a[i,i]
endfor
1.2.2雅可比法求对称矩阵特征值的并行算法
串行雅可比算法逐次寻找非主对角元绝对值最大的元素的方法并不适合于并行计算。
因此，在并行处理中，我们每次并不寻找绝对值最大的非主对角元消去，而是按一定的顺序将A中的所有上三角元素全部消去一遍，这样的过程称为一轮。
由于对称性，在一轮中，A的下三角元素也同时被消去一遍。
经过若干轮，可使A的非主对角元的绝对值减少，收敛于一个对角方阵。
具体算法如下：
设处理器个数为p，对矩阵A按行划分为2p块，每块含有连续的m/2行向量，记m=⎡n/p⎤，这些行块依次记为A0,A1,…,A2p-1，并将A2i与A2i+1存放与标号为i的处理器中。
每轮计算开始，各处理器首先对其局部存储器中所存的两个行块的所有行两两配对进行旋转变换，消去相应的非对角线元素。
然后按图21.1所示规律将数据块在不同处理器之间
传送，以消去其它非主对角元素。
图1.1p=4时的雅可比算法求对称矩阵特征值的数据交换图
这里长方形框中两个方格内的整数被看成是所移动的行块的编号。
在要构成新的行块配对时，只要将每个处理器所对应的行块按箭头方向移至相邻的处理器即可，这样的传送可以在行块之间实现完全配对。
当编号为i和j的两个行块被送至同一处理器时，令编号为i的行块中的每行元素和编号为j的行块中的每行元素配对，以消去相应的非主对角元素，这样在所有的行块都两两配对之后，可以将所有的非主对角元素消去一遍，从而完成一轮计算。
由图
1.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号处理器发送来的后数据块作为自己的后数据块
endif
（2）if（（my-rank=p-1）and（my-rankmod2≠0））then/*处理器p-1且其为奇数*/（2.1）fori=m/2tom-1do/*将后数据块保存在缓冲区buffer中*/
forj=0ton-1do
buffer[i-m/2,j]=a[i,j]
endforendfor
（2.2）fori=m/2tom-1do
buf[i-m/2]=b[i]
endfor
（2.3）fori=0tom/2-1do/*将前数据块移至后数据块的位置上*/
forj=0ton-1do
a[i+m/2,j]=a[i,j]
endforendfor
（2.4）fori=0tom/2-1do
b[i+m/2]=b[i]
endfor
（2.5）接收p-2号处理器发送的数据块作为自己的前数据块
（2.6）将buffer中的后数据块发送给编号为p-2的处理器
endif
（3）if（（my-rank=p-1）and（my-rankmod2=0））then/*处理器p-1且其为偶数*/（3.1）将后数据块发送给编号为p-2的处理器
（3.2）fori=0tom/2-1do/*将前数据块移至后数据块的位置上*/
forj=0ton-1do
a[i+m/2,j]=a[i,j]
endforendfor
（3.3）fori=0tom/2-1do
b[i+m/2]=b[i]
endfor
（3.4）接收p-2号处理器发送的数据块作为自己的前数据块
endif
（4）if（（my-rank≠p-1）and（my-rank≠0））then/*其它的处理器*/（4.1）if（my-rankmod2=0）then/*偶数号处理器*/
（i）将前数据块发送给编号为my_rank+1的处理器
（ii）将后数据块发送给编号为my_rank-1的处理器
（ii）接收编号为my_rank-1的处理器发送的数据块作为自己的前数据块
（iv）接收编号为my_rank+1的处理器发送的数据块作为自己的后数据块
else/*奇数号处理器*/
（v）fori=0tom-1do/*将前后数据块保存在缓冲区buffer中*/
forj=0ton-1do
buffer[i,j]=a[i,j]
endforendfor
（vi）fori=0tom-1do
buf[i]=b[i]
endfor
（vii）接收编号为my_rank-1的处理器发送的数据块作为自己的前数据块
End
endif
（viii）接收编号为my_rank+1的处理器发送的数据块作为自己的后数据块
（ix）将存于buffer中的前数据块发送给编号为my_rank+1的处理器
（x）将存于buffer中的后数据块发送给编号为my_rank-1的处理器
endif
各处理器并行地对其局部存储器中的非主对角元素aij进行消去，首先计算旋转参数并对第i行和第j行两行元素进行旋转行变换。
然后通过扩展收集操作将相应的旋转参数及第i列和第j列的列号按处理器编号连接起来并广播给所有处理器。
各处理器在收到这些旋转参数和列号之后，按0,1,…,p-1的顺序依次读取旋转参数及列号并对其m行中的第i列和第j列元素进行旋转列变换。
经过一轮计算的2p-1次数据交换之后，原矩阵A的所有非主对角元素都被消去一次。
此时，各处理器求其局部存储器中的非主对角元素的最大元localmax，然后通过归约操作的求最大值运算求得将整个n阶矩阵非主对角元素的最大元max，并广播给所有处理器以决定是否进行下一轮迭代。
具体算法框架描述如下：
算法21.5雅可比法求对称矩阵特征值的并行算法
输入：
矩阵An×n，ε
输出：
矩阵A的主对角元素即为A的特征值
Begin
对所有处理器my_rank（my_rank=0,…,p-1）同时执行如下的算法:
（a）fori=0tom-1do
b[i]=myrank*m+i/*b记录处理器中的行块数据在原矩阵A中的实际行号*/
endfor
（b）while（│max│>ε）do/*max为A中所有非对角元最大的绝对值*/
（1）fori=my_rank*mto（my_rank+1）*m-2do
/*对本处理器内部所有行两两配对进行旋转变换*/
forj=i+1to（my_rank+1）*m-1do
（1.1）r=imodm,t=jmodm/*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）forv=0ton-1do/*对两个主行其余元素的旋转变换*/
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
endif
endfor
（iii）forv=0ton-1do
if（（v≠i）and（v≠j））then
a[r,v]=br[v]
endifendfor
（iv）forv=0tom-1do
/*对两个主列在本处理器内的其余元素的旋转变换*/
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
endifendfor
（v）forv=0tom-1do
if（（v≠r）and（v≠t））then
a[v,i]=bi[v]
endifendfor
（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
endif
（1.3）将所有处理器temp1中的旋转参数及主行的行号按处理器编号连接起来并广播给所有处理器,存于temp2中
（1.4）current=0
（1.5）forv=1topdo
/*根据temp2中的其它处理器的旋转参数及主行的行号对相关的列在本处理器的部分进行旋转变换*/
（i）Compute:
s1=temp2[（v-1）*4+0],c1=temp2[（v-1）*4+1]，
i1
升级会员 