并行计算-矩阵特征值计算--
9矩阵特征值计算 在实际的工程计算中,常常会遇到求 n 阶方阵 A 的特征值(Eigenvalue)与特征向量(Eigenvector)的问题。对于一个方阵 A,假如数值 λ 使方程组 Ax=λx 即 (A-λIn )x=0 有非零解向量(Solution Vector)x,则称λ为方阵A的特征值,而非零向量x为特征 值λ所对应的特征向量,其中In 为n阶单位矩阵。 由于依据定义干脆求矩阵特征值的过程比较困难,因此在实际计算中,往往实行一些 数值方法。本章主要介绍求一般方阵肯定值最大的特征值的乘幂(Power)法、求对称方阵特 征值的雅可比法和单侧旋转(One-side Rotation)法以及求一般矩阵全部特征值的 QR 方法及 一些相关的并行算法。 1.1 求解矩阵最大特征值的乘幂法 1.1.1 乘幂法及其串行算法 在很多实际问题中,只须要计算肯定值最大的特征值,而并不须要求矩阵的全部特征值。 乘幂法是一种求矩阵肯定值最大的特征值的方法。记实方阵A的n个特征值为λi i=(1,2, …,n), 且满意: │λ1 │≥│λ2 │≥│λ3 │≥…≥│λn │ 特征值λi 对应的特征向量为xi 。乘幂法的做法是:①取n维非零向量v0 作为初始向量;②对于k=1,2, …,做如下迭代: 直至 u k +1 ¥ - uk uk =Avk-1vk = uk /║uk ║∞ ε) 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 的第 im 至第(i+1)m-1 行数据,(i=0,1, …,p-1),初 始向量 v 被广播给全部处理器。 各处理器并行地对存于局部存储器中 a 的行块和向量 v 做乘积操作,并求其 m 维结果 向量中的最大值 localmax,然后通过归约操作的求最大值运算得到整个 n 维结果向量中的最 大值 max 并广播给全部处理器,各处理器利用 max 进行规格化操作。最终通过扩展收集操 作将各处理器中的 m 维结果向量按处理器编号连接起来并广播给全部处理器,以进行下一 次迭代计算,直至收敛。详细算法框架描述如下: 算法 21.2乘幂法求解矩阵最大特征值的并行算法 输入:系数矩阵An×n ,初始向量v n×1 ,ε 输出:最大的特征值 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, 一次扩展收集操作,通信量为 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=[aij ]是n阶实对称矩阵,则存在一个正交矩阵U,使得 UTAU=D 其中 D 是一个对角矩阵,即 él1 ê D = ê0 êM ê êë0 0 L l2 L MO 0L 0 ù ú 0 ú M ú ú ln úû 这里λi (i=1,2,…,n)是A的特征值,U的第i列是与λi 对应的特征向量。雅可比算法主要是通过 正交相像变换将一个实对称矩阵对角化,从而求出该矩阵的全部特征值和对应的特征向量。 因此可以用一系列的初等正交变换逐步消去A的非对角线元素,从而使矩阵A对角化。设初 等正交矩阵为R(p,q,θ),其(p,p)( q,q)位置的数据是cosθ,(p, q)( q, p)位置的数