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

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, …,做如下迭代：

直至u

k+1 ∞

-

u

k

u 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 ax

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 行向量，这里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 ax

Begin

对所有处理器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 do

sum=sum+a[i,j]*x[j]

end for

(1.3)b[i]= sum

end for

(2)localmax=│b[0]│/*对所计算的特征向量的相应分量

求新旧值之差的最大值l ocalmax */

(3)for i=1 to m-1 do

if (│b[i]│>localmax) then localmax=│b[i]│ end if

end for

(4)用A llreduce 操作求出所有处理器中l ocalmax 值的最大值m ax

并广播到所有处理器中

(5)for i=0to m-1 do /*对所计算的特征向量归一化*/

b [i ] =b [i ]/max

end 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

λ2

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 ,θ) 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，即：