预处理变形共轭梯度法并行求解矩阵的Moore-Penrose逆曹方颖;吕全义【摘要】提出了一种求解Moore-Penrose逆的并行预处理变形共轭梯度法,将求解Moore-Penrose 逆转化求解矩阵方程极小范数解或极小范数最小二乘解的问题.给出了两种预处理方法.一种方法是给出预处理矩阵是可逆对角矩阵,然后并行求解预处理矩阵方程；另一种方法是给出预处理矩阵是严格对角占优矩阵,该方法提出了迭代法的预处理模式,构造并行迭代求解预处理矩阵方程的迭代格式,进而使用变形共轭梯度法并行求解.通过数值试验,这两种预处理方法与直接使用变形共轭梯度法相比较,第二种方法有效提高了收敛速度,而且具有很好的并行性.【期刊名称】《纺织高校基础科学学报》【年(卷),期】2013(026)001【总页数】6页(P137-142)【关键词】并行算法;预处理变形共轭梯度法;预处理矩阵方程;Moore-Penrose逆【作者】曹方颖;吕全义【作者单位】西北工业大学应用数学系,陕西西安710129;西北工业大学应用数学系,陕西西安710129【正文语种】中文【中图分类】O246广义逆矩阵是矩阵理论中的一项重要的发展，特别自近五十年代以来，广义逆矩阵的理论和计算方法的研究取得了长足的发展，并在数值分析，数学规划，数理统计，控制论，博弈论和网络理论等领域发挥着重要的应用［1］.求解Moore-Penrose 逆的算法最终转化为求解线性矩阵方程.这种线性矩阵方程是无解的，因此求解线性矩阵方程的极小范数最小二乘解.对于求解线性矩阵方程，虽然串行算法已经比较完善，但随着计算规模的增大，求解线性矩阵方程的存储需求和计算量快速增加，单台处理机往往难以承受，并行计算势在必行.近年来，一般求解线性矩阵方程的直接法［2］主要有LU分解法，Cholesky分解法，QR分解法等；迭代法［2］有Jacobi迭代法，Gauss-Seidel迭代法，JGS迭代法，SOR迭代法，SSOR迭代法，SAOR迭代法，参数迭代法等.关于矩阵方程的最小二乘解或矩阵Moore-Penrose逆的研究已经取得了许多成果［3-5］，但对于预处理与并行计算的研究还很少见.对于求解大型线性矩阵方程以共轭梯度法和变形共轭梯度法为主，其具有存储量少，计算量少和适合并行计算等优点.文献［2］研究了关于求解矩阵方程极小范数解或极小范数最小二乘解的变形共轭梯度法.变形共轭梯度法作为最基本的Krylov子空间方法，易于并行化.变形共轭梯度法的收敛速度与系数矩阵的条件数紧密相关，条件数愈小，收敛性愈好，该算法可以在很少的几步就会获得高精度的近似解.但当系数矩阵的条件数很大时，收敛速度就很慢.于是出现了预处理变形共轭梯度法［6-7］（简称PMCG法），它是通过适当的预处理方法引入预处理矩阵M，使矩阵的特征值分布更为集中，降低矩阵条件数，以达到提高收敛速度的目的.本文针对矩阵方程极小范数解或极小范数最小二乘解问题提出了两种并行解决方案.第一种方案是为了降低求解预条件矩阵方程AM-1＝的计算复杂度，选取预处理矩阵M为可逆对角矩阵，预处理矩阵方程在求解过程中只在相邻处理机间有通信.第二种方案是为了降低求解预条件矩阵方程AM-1 ＝的计算复杂度，对预处理矩阵方程采取若干次迭代方法来求解，而且在求解过程中只在相邻处理机间有通信.两种方法都有效的改变了原矩阵方程的条件数又具有很好的并行性.最后在HP2600集群上进行了数值试验，并与传统的变形共轭梯度法的计算结果进行了比较.文中（ATA）⊗I表示矩阵（ATA）与I的Kronecker积（X）表示将矩阵X按行拉直构成的列向量.定义同型矩阵A与B的内积为［A，B］＝tr（ATB），由此导出矩阵的Frobenius范数1 预处理方法1.1 方案一求解矩阵A的Moore-Penrose逆就是求解线性矩阵方程AX＝I极小范数解或极小范数最小二乘解，AX＝I正规线性矩阵方程为其中 A为m×n矩阵，X为n×m矩阵.将矩阵方程（1）通过按行拉直转化为线性方程组设λ1，λ2，…，λn是 ATA的特征值，由文献［1］的结论可知矩阵（ATA）⊗I的特征值为λ1，…，λ1，λ2，…，λ2，…，λn，…，λn，可见矩阵（ATA）⊗I与矩阵ATA的条件数相同.若矩阵A的行数m≤n，则取M为n阶可逆的对角矩阵为这样，矩阵方程（1）可转化为只要使M-1ATAM-1的条件数小于ATA，则采用变形共轭梯度法求解矩阵方程（3）的速度将得到提高.对于方程（3）首先计算AM-1＝，设MX＝Y.则方程（3）可转化为求解若m＞n，可将矩阵A的Moore-Penrose逆转化成AT的Moore-Penrose逆采用同样预处理矩阵即可.1.2 方案二选取预条件矩阵M使得AM-1尽量接近矩阵A的等价标准形，如果不考虑AM-1的精确解，采用迭代算法计算AM-1的近似矩阵，只要迭代法适合并行计算，那么对于M的选取就简单的多.由于Jacobi迭代算法适合并行计算，所以采用该算法，又Jacobi迭代算法收敛的充分条件为系数矩阵为严格对角占优，故若m＜n，选取预处理矩阵M为n阶可逆且含矩阵A的带状部分，即每行的半带宽ri的选取要求矩阵M是严格对角占优矩阵，若矩阵A主对角元有零元，在矩阵M中取1.若m＞n，同样将矩阵A的Moore-Penrose逆转化成AT的Moore-Penrose逆即可.在求解AM-1的近似值时，设＝AM-1采用的迭代格式为其中 M＝D-L-U，D，-L，-U分别是矩阵M的对角矩阵，严格下三角矩阵与严格上三角矩阵.求解线性矩阵方程AX＝I转化为求解设MX＝Y，则式（5）可转化为求解将文献［2］矩阵方程的变形共轭梯度法应用于方程（4）和（6），推出矩阵方程（4）和（6）变形共轭梯度算法为步骤1 任给初始矩阵Y（0），计算步骤2 对k＝0，1，…，计算步骤3 计算 MX＝ Y（k＋1）.显然，方案一是方案二的一种特殊情况.即方案二中，初始矩阵＝O且l＝1的情况.由于这样选取M可以保证在求解预处理矩阵方程过程中只需要相邻处理机进行通讯.在方案二中，不需要准确计算AM-1，而是利用上面的迭代算法迭代l次，这样既降低了矩阵的条件数，又减少了计算量.2 算法的并行实现为叙述方便，设处理机台数为p，p整除m且m＝pl，p整除n且n＝ps，pi （i ＝1，2，…，p）表示第i台处理机，myid表示当前处理机，mod表示整除.记×l列块子矩阵.由于矩阵的存储比较复杂，所以在本文中采用了3种矩阵相乘的并行算法，分别是［8-9］方法1：行列划分矩阵相乘；方法2：列行划分矩阵相乘；方法3：列列划分矩阵相乘.（1）存储方式.将 Ai，Bi均存储于第pi （i＝1，2，…，p）台处理机中.（2）预处理过程.对于方案2，任选初始矩阵（0），其中用方法1计算用方法2计算矩阵乘积＝ Q（0）.（3）循环过程.（ⅰ）各进程计算Y＝Y（k）转到（4）；否则，继续计算.（ⅱ）各进程计算αk ＝［R（k），R（k）］／［Z（k），Z（k）］.（ⅲ）各进程计算（ⅳ）用方法2计算DY（k＋1），各进程计算.各进程计算到［R（k＋1），R （k＋1）］.（ⅴ）用方法3计算（ⅵ）置k：＝k＋1，返回到（ⅰ）.（4）求解X过程.选初始矩阵X（0）＝O，迭代计算其中用方法3计算T＝（L＋U）X（i），各进程计算 Wi＝Ti＋Yi，再用方法3计算X（i＋1）＝D-1 W.3 数值算例与结果分析3.1 算例在并行机上分别采用MCG算法与本文给出的PMCG算法计算2个算例，进行了比较.在迭代过程中，初始矩阵Y（0），X（0）为零矩阵，终止条件ε＝10-10. 例1 设矩阵A＝（aij）m×n，其中MoorePenrose A的-逆矩阵.这里取m＝1 000，n＝1 200.求解结果见表1.表1 例1m＝1 000，n＝1 200的计算结果算法 p 10 12 14 16变形共轭梯度法T ／s 3 211.932 6 2 731.052 3 2 413.631 3 2 132.829 2 K 1 207 1 206 1 206 1 206 S 11.760 8 13.307 5 15.059 5 E 0.98 0.95 0.94 E 0.17 0.17 0.16 0.16本文算法方案一T／s 3 152.603 2 2 704.986 4 2 393.522 3 2 067.931 6 K 1 190 1 191 1 190 1 183 S 11.555 0 13.171 4 15.245 E 0.96 0.94 0.95 E 0.17 0.16 0.16 0.16本文算法方案二l＝1 T／s 3 388.900 4 2 890.228 6 2 541.558 1 2 201.466 2 K 1 190 1 191 1 190 1 183 S 11.73 13.33 15.39 E 0.98 0.95 0.96 E 0.16 0.15 0.15 0.15续表1 例1m＝1 000，n＝1 200的计算结果算法 p 10 12 14 16本文算法方案二l＝2 T／s 530.000 8 442.980 9 382.557 3 343.035 6 K 181 181 181 181 S11.96 13.85 15.45 E 0.99 0.98 0.96 E 1 0.99 0.98 0.96本文算法方案二l＝3 T／s 2 635.418 7 2 226.532 3 1 953.942 3 1 707.247 0 K 914 915 914 914 S 11.83 13.49 15.44 E 0.99 0.96 0.96 E 0.20 0.20 0.19 0.19例2 设矩阵 A＝（aij）m×n，其中求A的Moore-Penrose逆矩阵.这里取m＝100，n＝300，求解结果见表2.表2 例2m＝100，n＝300的计算结果算法 p 0 2 4 6变形共轭梯度法T／s 924.698 3 520.135 1 282.739 4 219.731 4 K 7 781 7 874 7 718 7 718 S 1.78 3.27 4.40 E 0.89 0.82 0.73 E 0.11 0.10 0.09 0.08本文算法方案一T／s 407.9622 250.5793 143.814 1 84.209 6 K 3 427 3 953 3 954 3 092 S 1.62 2.84 4.84E 0.81 0.71 0.81 E 0.25 0.20 0.18 0.20本文算法方案二l＝1 T／s 474.3510 283.651 6 154.900 9 87.411 8 K 3 427 3 953 3 954 3 092 S 1.67 3.06 5.43 E 0.84 0.77 0.90 E 0.22 0.18 0.17 0.20本文算法方案二l＝2 T／s 357.948 3 182.376 0 103.425 8 74.069 5 K 2 586 2 515 2 577 2 520 S 1.96 3.46 4.83 E 0.98 0.87 0.81 E 0.29 0.28 0.25 0.23本文算法方案二l＝3 T／s 102.549 6 51.875 0 28.800 7 21.062 7 K 766 730 737 734 S 1.98 3.56 4.87 E 0.99 0.89 0.81 E 1 0.99 0.89 0.81注1 表1，2中p表示处理机台数，T（s）表示时间，K表示迭代次数，S表示加速比，E表示相对并行效率表示绝对并行效率.注2 由于一台处理机的计算时间较长，因此采用多台计算机的计算时间，并且使用多台处理机中最小的计算时间作为基准来计算加速比和并行效率.3.2 结果分析（1）本文算法明显提高了迭代速度，提高了变形共轭梯度法的效率，具有良好的并行性.（2）例1中，与变形共轭梯度法相比较，本文算法方案一与方案二中l＝1的迭代次数没有明显减少，说明原线性矩阵方程的条件数没有明显变化，预处理算法基本失效，这是因为方案一的预处理矩阵近似数量矩阵.在本文算法方案二中，当内迭代次数l＝2时效果最好，迭代速度已经达到很快并且并行效率很高，优于l＝1的情况，即优于传统的Jacobi预处理方法.（3）例2的结果可以看出，与变形共轭梯度法相比较，方案一与方案二的迭代次数明显减少，大大地提高了收敛速度.致谢：感谢西北工业大学高性能计算研究与发展中心的大力支持！【相关文献】［1］王松佳，杨振海.广义逆矩阵及其应用［M］.北京：北京工业大学出版社，2006.［2］张凯院，徐仲.数值代数［M］.2版修订本.北京：科学出版社，2010.［3］袁仕芳，廖安平，雷渊.矩阵方程AXB＋CYD＝E的对称极小范数最小二乘解［J］.计算数学，2007，29（2）：203-216.［4］安晓虹，徐仲，陆全，等.行满秩Toeplitz 型矩阵Moore-Penrose逆的快速算法［J］.计算机工程与应用，2010，46（30）：1-4.［5］李书连，张凯院.多变量LME一种异类约束最小二乘解的迭代算法［J］.纺织高校基础科学学报，2011，24（3）：326-332.［6］张宝琳，谷同祥，莫则尧.数值并行原理与方法［M］.北京：国防工业出版社，1999.［7］胡家赣.解线性代数方程组的迭代解法［M］.北京：科学出版社，1999：173-201.［8］陈国良，安虹，陈俊，等.并行算法实践［M］.北京：高等教育出版社，2004.［9］李晓梅，吴建平.数值并行算法与软件［M］.北京：科学出版社，2007.。