并行计算-矩阵特征值计算 -

9 矩阵特征值计算在实际的工程计算中，经常会遇到求n 阶方阵 A 的特征值(Eigenvalue)与特征向量(Eigenvector)的问题。

对于一个方阵A，如果数值λ使方程组Ax=λx即(A-λI n )x=0 有非零解向量(Solution Vector)x，则称λ为方阵A的特征值，而非零向量x为特征值λ所对应的特征向量，其中I n 为n阶单位矩阵。

由于根据定义直接求矩阵特征值的过程比较复杂，因此在实际计算中，往往采取一些数值方法。

本章主要介绍求一般方阵绝对值最大的特征值的乘幂(Power)法、求对称方阵特征值的雅可比法和单侧旋转(One-side Rotation)法以及求一般矩阵全部特征值的QR 方法及一些相关的并行算法。

1.1 求解矩阵最大特征值的乘幂法1.1.1 乘幂法及其串行算法在许多实际问题中，只需要计算绝对值最大的特征值，而并不需要求矩阵的全部特征值。

乘幂法是一种求矩阵绝对值最大的特征值的方法。

记实方阵A的n个特征值为λi i=(1,2, …,n)，且满足：│λ1 │≥│λ2 │≥│λ3 │≥…≥│λn │特征值λi 对应的特征向量为x i 。

乘幂法的做法是：①取n维非零向量v0 作为初始向量；②对于k=1,2, …,做如下迭代：直至uk+1 ∞-uku k =Av k-1 v k = u k /║u k ║∞<ε为止，这时v k+1 就是A的绝对值最大的特征值λ1 所对应的特征向∞量x1 。

若v k-1 与v k 的各个分量同号且成比例，则λ1 =║u k ║∞；若v k-1 与v k 的各个分量异号且成比例，则λ1 = -║u k ║∞。

若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时间，则因为一轮计算要做一次矩阵向量相乘、一次求最大元操作和一次规格化操作，所以下述乘幂法串行算法21.1 的一轮计算时间为n2+2n=O(n2 )。

算法21.1 单处理器上乘幂法求解矩阵最大特征值的算法输入：系数矩阵A n×n ，初始向量v n×1 ，ε输出：最大的特征值m axBeginwhile (│diff│>ε) do(1)for i=1 to n do(1.1)sum=0(1.2)for j= 1 to n dosum=sum+a[i,j]*x[j]end for(1.3)b[i]= sumend for(2)max=│b[1]│(3)for i=2 to n doif (│b[i]│>max) then max=│b[i]│ end ifend for(4)for i=1 to n dox[i] =b[i]/maxend for(5)diff=max-oldmax , oldmax=maxend whileEnd1.1.2 乘幂法并行算法乘幂法求矩阵特征值由反复进行矩阵向量相乘来实现，因而可以采用矩阵向量相乘的数据划分方法。

设处理器个数为p，对矩阵A 按行划分为p 块，每块含有连续的m 行向量，这里m=⎡n/ p⎤，编号为i的处理器含有A的第i m 至第(i+1)m-1 行数据，(i=0,1, …,p-1)，初始向量v被广播给所有处理器。

各处理器并行地对存于局部存储器中a 的行块和向量v 做乘积操作，并求其m 维结果向量中的最大值localmax，然后通过归约操作的求最大值运算得到整个n 维结果向量中的最大值max 并广播给所有处理器，各处理器利用max 进行规格化操作。

最后通过扩展收集操作将各处理器中的m维结果向量按处理器编号连接起来并广播给所有处理器，以进行下一次迭代计算，直至收敛。

具体算法框架描述如下：算法21.2 乘幂法求解矩阵最大特征值的并行算法输入：系数矩阵A n×n ，初始向量v n×1 ，ε输出：最大的特征值m axBegin对所有处理器m y_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 dosum=sum+a[i,j]*x[j]end for(1.3)b[i]= sumend for(2)localmax=│b[0]│/*对所计算的特征向量的相应分量求新旧值之差的最大值l ocalmax */(3)for i=1 to m-1 doif (│b[i]│>localmax) then localmax=│b[i]│ end ifend for(4)用A llreduce 操作求出所有处理器中l ocalmax 值的最大值m ax并广播到所有处理器中(5)for i=0to m-1 do /*对所计算的特征向量归一化*/b [i ] =b [i ]/maxend for (6)用 A llgather 操作将各个处理器中计算出的特征向量的分量的新值组合并广播到 所有处理器中(7)diff=max -oldmax, oldmax =max end while End若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时间，一轮迭代的计算时间为 m n +2m ；一轮迭代中，各处理器做一次归约操作，通信量为 1， 一次扩展收集操作，通信量为 m ，则通信时间为 4t s (p - 1) + (m + 1)t w ( p - 1) 。

因此乘幂法的一轮并行计算时间为T p = mn + 2m + 4t s (MPI 源程序请参见所附光盘。

p - 1) + (m + 1)t w ( p - 1) 。

1.2 求对称矩阵特征值的雅可比法1.2.1 雅可比法求对称矩阵特征值的串行算法若矩阵A =[a ij ]是n 阶实对称矩阵，则存在一个正交矩阵U ，使得U T AU =D 其中 D 是一个对角矩阵，即⎡λ1⎢ D = ⎢0⎢ ⎢ ⎢⎣0λ20 ⎤⎥ 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 ,θ) T R (p ,q ,θ)=I n不妨记B= R (p ,q ,θ)T AR (p ,q ,θ)，如果将右端展开，则可知矩阵B 的元素与矩阵A 的元素之 间有如下关系：b pp = a pp cos 2θ+a qq sin 2θ+a pq sin2θ b qq = a pp sin 2θ+a qq cos 2θ-a pq sin2θb pq = ((a qq -a pp )sin2θ)/2+a pq cos2θ b qp = b pq b pj = a pj cos θ+a qj sin θ b qj = -a pj sin θ +a qj cos θ b ip = a ip cos θ+a iq sin θ b iq = -a ip sin θ +a iq cos θ b ij = a ij其中 1 ≤ i , j ≤ n 且i ,j ≠ p ,q 。

因为A 为对称矩阵，R 为正交矩阵，所以B 亦为对称矩阵。

若 要求矩阵B 的元素b pq =0，则只需令 ((a qq -a pp )sin2θ)/2+a pq cos2θ=0，即：pqqqtg 2θ =-a pq (a q q - a pp ) 2在实际应用时，考虑到并不需要解出 θ，而只需要求出 s in2θ，sin θ 和 c os θ 就可以了。

如果限制 θ 值在-π/2 < 2θ ≤ π/2，则可令m = - a pq , n = 1( a- app ) , w = sgn( n ) m容易推出： sin 2θ = w , 2sin θ = 2(1 + w ,1 - w2 )cos θ = m 21 - sin2 θ+ n 2利用sin2θ，sin θ和cos θ的值，即得矩阵B 的各元素。

矩阵A 经过旋转变换，选定的非主对角元素a pq 及a qp （一般是绝对值最大的）就被消去，且其主对角元素的平方和增加了 2a2，而非主对角元素的平方和减少了 2a 2 pq ，矩阵元素总的平方和不变。

通过反复选取主元素a pq ， 并作旋转变换，就逐步将矩阵A 变为对角矩阵。

在对称矩阵中共有(n 2-n )/2 个非主对角元素要 被消去, 而每消去一个非主对角元素需要对 2n 个元素进行旋转变换, 对一个元素进行旋转 变换需要 2 次乘法和 1 次加法，若各取一次乘法运算时间或一次加法运算时间为一个单位时 间，则消去一个非主对角元素需要 3 个单位运算时间，所以下述算法 21.3 的一轮计算时间 复杂度为 (n 2-n )/2*2n *3=3n 2(n-1)=O(n 3)。

算法 21.3 单处理器上雅可比法求对称矩阵特征值的算法输入：矩阵A n ×n ，ε 输出：矩阵特征值 E igenvalueBegin(1)while (max >ε) do(1.1) max =a[1,2] (1.2)for i =1 to n dofor j = i +1 to n doif (│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 )， si 2=w ，si 1=w /sqrt(2(1+ sqrt(1-w *w )))，co1=sqrt(1-si 1*si 1)， b [p ,p ]= a [p ,p ]*co1*co1+ a [q ,q ]*si 1*si 1+ a [p ,q ]*si 2， b [q ,q ]= a [p ,p ]*si 1*si 1+ a [q ,q ]*co1*co1- a [p ,q ]*si 2， b [q , p ]=0，b [p ,q ]=0(1.4)for j =1 to n doif ((j ≠ p ) and ( j ≠ q )) then (i)b [p ,j ]= a [p ,j ]*co1+ a [q ,j ]*si 1 (ii)b [q ,j ]= -a [p ,j ]*si 1 + a [q ,j ]*co1 end if end for(1.5)for i =1 to n doEndif ((i ≠ p ) and (i ≠ q )) then(i)b [i , p ]= a [i ,p ]*co1+ a [i ,q ]*si 1 (ii)b [i , q ]= - a [i ,p ]*si 1+ a [i ,q ]*co1 end if end for(1.6)for i=1 to n dofor j=1 to n do a [i,j ]=b [i,j ] end for end for end while (2)for i =1 to n doEigenvalue [i ]= a [i , i ] end for1.2.2 雅可比法求对称矩阵特征值的并行算法串行雅可比算法逐次寻找非主对角元绝对值最大的元素的方法并不适合于并行计算。