基于ADMM迭代法的PageRank问题的求解_吴文军.pdf

收稿日期:2022 09 26基金项目:国家自然科学基金青年基金(11901024);福建省自然科学基金面上资助项目(2020J01166，2021J01661)作者简介:吴文军(1998)，男，江西省抚州市人，福建师范大学数学与统计学院硕士研究生，主要从事数值代数研究通信作者:唐嘉(1983)，女，湖南省湘潭市人，理学博士，福建师范大学数学与统计学院副教授，硕士生导师，主要从事数值代数及其应用研究基于 ADMM 迭代法的 Pageank 问题的求解吴文军，唐嘉(福建师范大学 数学与统计学院，福建 福州 350007)摘要:针对 Pageank 问题导出的线性方程组，首先对方程组的系数矩阵进行了 LU 分解，提出了一种基于交替方向乘子法(ADMM)形式的迭代算法，用于求解该线性方程组的最小二乘解;然后，证明了所提出的算法的收敛性;最后，数值结果表明了该算法的可行性关键词:Pageank;ADMM;LU 分解中图分类号:O242 2文献标识码:A文章编号:1673 1670(2023)02 0008 060引言随着互联网及其技术的蓬勃发展，网络搜索引擎已经成为检索信息最受欢迎的工具 谷歌的 Pag-eank 算法成为该领域重要的网页排名算法，它通过网页的链接结构来量化每个页面的重要性，在最近几年引起了相当大的关注1 479由谷歌设计的 Pageank 算法2 是一个网页内容中立的排序函数，被应用于许多领域，如 Web 的超链接结构和用马尔可夫链建模的图，神经和社会网络分析等3 Pageank 问题就是求解 Google 矩阵 A 的首特征值 1 所对应的特征向量，即 Pageank 向量 x满足Ax=x(1)其中 A=P+(1 )veT，(0，1)是阻尼因子，v是元素和为 1 的向量，e=(1，1，1)T，P 是列随机矩阵x 是一个非负向量，满足x1=1 又因为x1=1，所以式(1)的 Pageank 问题可以改写为线性方程(I P)x=(1 )v(2)其中 In n是一个单位矩阵，且 eTx=1 对于求解线性系统(2)，学者们提出了很多迭代方法 如内外迭代(inner-outer，IO)法4 353、幂内外(PIO)迭代法5 23、广义内外(GIO)迭代方法1 483和多步幂内外(MPIO)迭代方法6 91 近年来并没有出现矩阵分解来求解 Pageank 问题，对此笔者考虑利用LU 分解，之后将问题转化为最小二乘问题来求解，对最 小 二 乘 问 题 利 用 的 是 交 替 方 向 乘 子 法(ADMM)来求解本文的结构如下:在第 1 节中，阐述几个重要的定义;在第 2 节中，简要介绍交替方向乘子法(ADMM)和所提出的算法;第 3 节详细证明所提出算法的收敛性;最后，在第 4 节中给出几个数值例子来说明该方法的可行性1预备知识定义 17 51对任意的 x，yD，0，1，有x+(1 )yD，则称 D 为凸集定义 27 54对任意的 x，yn，有 f(y)f(x)+f(x)，y x，则称连续可微函数 fn 为在 n上的凸函数定义 37 56对任意的 x，yn，0，有f(y)f(x)+f(x)，y x)+2y x22，则称连续可微函数 fn 为在 n上的强凸函数注 18 设 fDn 是二次连续可微函数，f 是强凸函数当且仅当2f(x)I 是对称半正定的 其中 0，2f(x)是 f(x)的 Hessian 阵注 27 57设 fDn 是强凸函数，对任意的 x，yD，0，有以下两个不等式成立:第 38 卷第 2 期2023 年 4 月平顶山学院学报Journal of Pingdingshan UniversityVol 38 No 2Apr 2023f(x)f(y)f(x)，xy+2x y22，f(x)f(y)，x y x y222改进的交替方向乘子法2 1ADMM交替方向乘子法(ADMM)由 Gabay、Mercier、Glowinski 和 Marrocco 于 1970 年首次提出 该方法具有可分解性，与对偶分解法非常相似，并具有类似于乘子法的收敛性 ADMM 在分布式计算中的简单实现非常适合于大规模凸优化问题9 4 其最显著的优点是能够对变量进行分离处理，并且充分利用目标函数的可分离结构，求解收敛速度快，易于工程实现1983 年，Gabay、Glowinski 和 Tallec10 展示了ADMM 的全局收敛性 Yang Junfeng 和 Yuan Xi-aoming 在文献 11 中证明了 ADMM 收敛于两个迭代步的不精确解，因而适用于迭代方法 此外，最近还发现，当其中有一个函数是非凸时，文献 12 讨论了 ADMM 算法的收敛性 贾泽慧等13 证明了一类非凸优化问题 ADMM 算法的全局收敛性2 2算法将 Pageank 问题转化为方程组(2)，对于其系数矩阵 I P 是非奇异的 M-矩阵的情形，可对I P 进行 LU 分解，之后将其转化为如下的最小二乘解问题:minx12LUx (1 )v22(3)令 Ux=y，则式(3)的等价形式如下:minx，y12Ly (1 )v22，s t Ux y=0(4)问题(4)的增广拉格朗日函数为L(x，y，)=12Ly (1 )v22T(Ux y)+2Ux y22(5)其中 n，0 是惩罚参数计算增广拉格朗日函数(5)的平稳点，可得xL(x，y，)=UT(Ux y)UT=0UTUx=UT(+y)，(6)yL(x，y，)=LTLy (1 )LTv+y Ux=0(LTL+I)y=(1 )LTv +Ux (7)因为 0，所以矩阵 UTU 和 LTL+I 是对称正定的，故 x=(UTU)1UT(+y)和 y=(LTL+I)1(1 )LTv +Ux)是 L(x，y，)的稳定点 因此，利用 ADMM 算法求解问题(5)，其 ADMM的迭代格式9 14为xk+1=arg min Lx(x，yk，k)，(8)yk+1=arg min Ly(xk+1，y，k)，(9)k+1=k(Uxk+1 yk+1)(10)根据式(6)和式(7)，可以将 ADMM 的迭代更新为式(11):xk+1=(UTU)1UT(k+yk)，yk+1=(LTL+I)1(1 )LTv k+Uxk+1)，k+1=k(Uxk+1 yk+1)(11)算法 1(用于求解 Pageank 问题(2)的LUADMM 迭代法)步骤一给定初始的向量 v，x0，y0，矩阵 P，迭代收敛精度，参数，令 k=0步骤二对 I P 进行 LU 分解，得到 L 和 U矩阵步骤三计算xk+1=(UTU)1UT(k+yk)，yk+1=(LTL+I)1(1 )LTv k+Uxk+1)，k+1=k(Uxk+1 yk+1)步骤四若 Pxk+1+(1 )v xk+12，x*=xk+1xk+11，结束 否则，k=k+1，回到步骤三3算法的全局收敛性分析利用以下的引理和定理，讨论算法的收敛性引理 1增广拉格朗日函数(5)是关于 x 和 y的强凸函数，其参数为 min(UTU)和 min(LTL)，同时满足L(x+x，y，)L(x，y，)xTxL(x，y，)+min(UTU)2x2，(12)L(x，y+y，)L(x，y，)yTyL(x，y，)+min(LTL)2y2(13)证明根据注 1 和 2，需证明增广拉格朗日函数的 Hessian 阵是关于 x 和 y 对称半正定的，xL(x，y，)=UT(Ux y)UT2xL(x，y，L)=UTU，(14)yL(x，y，)=LTLy (1 )LTv+y Ux2yL(x，y，)=LTL+I(15)由式(14)和式(15)可得2xL(x，y，)min(UTU)I=9第 2 期吴文军，唐嘉:基于 ADMM 迭代法的 Pageank 问题的求解(UTU min(UTU)I)，(16)2yL(x，y，)min(LTL)I=LTL min(LTL)I+I(17)因为 I P 是非奇异的 M-矩阵，L 和 U 矩阵是经过 LU 分解得到的，所以 L 和 U 矩阵是可逆的，所以 UTU 是对称正定的，UTU min(UTU)I 是对称半正定矩阵 因此，对于 0，增广拉格朗日函数是关于 x 的强凸函数，其参数为 min(UTU)同样地，LTL 是对称正定的，LTL min(LTL)I 是对称半正定矩阵 显然，对于 0，LTL min(LTL)I+I 是对称半正定的，因此 L 是关于 y 的强凸函数，其参数为 min(LTL)定理 1若存在 x*，y*，*n使得 Ux*y*=0，并且对于任意的 x，yn，(x*，y*)是L(x，y，)上的最小点，则(x*，y*)是式(4)的最优解证明由约束条件 Ux y=0，可得 Ux y22=0，T(Ux y)=0 因此问题(4)可以看成min(x，y)，Ux=y12Ly (1 )v22=min(x，y)，Ux=y12Ly (1 )v22T(Ux y)+2Ux y2=min(x，y)，Ux=yL(x，y，)(18)接下来，忽略约束条件 Ux=y 来扩大 L(x，y，*)最小值的集合，因此不等式(19)是成立的:min(x，y)，Ux=yL(x，y，*)min(x，y)L(x，y，*)(19)此外，由于对于任意的 x，y，n，(x*，y*)是 L(x，y，)上的最小点，所以下面等式成立:min(x，y)，Ux=yL(x，y，*)=12Ly*(1 )v22T(Ux y)+2Ux y22(20)最后，由 Ux*y*=0 和式(20)，可知对任意的 x，yn，不等式(21)成立:min(x，y)，Ux=y12Ly (1 )v2212Ly*(1 )v22(21)由于 Ux*y*=0，(x*，y*)是问题(4)的可行解，因此不等式(21)表明(x*，y*)是问题(4)的最优解此外，还需要考虑问题(4)和它的拉格朗日对偶问题(22):maxminx，yL(x，y，)(22)假设(x*，y*)是一个最优解，*是一个最优对偶解 现需要找到问题(5)的充分必要的最优性KKT 条件 显然，这两个问题的互补松弛和对偶可行性条件是需要同时满足的 此外，初始可行性和稳定条件如下初始可行性:Ux*y*=0;稳定条件:x12Ly (1 )v()22(x*，y*)=0=UT*，y12Ly (1 )()v(x*，y*)=LT(Ly*(1 )v)=*(23)因此 KKT 条件如下:UT*=0(24)引理 2若 Un n，x，y，n，0，则有2Ux y 1221222=2Ux y22 T(Ux y)证明2Ux y 1221222=2y Ux+1221222=2y Ux22+T(y Ux)[phone]=2Ux y22 T(Ux y)定理 2若 k 是有界序列，而且k=0k+1 k22，(25)则由算法生成的序列(xk，yk，k)的极限点满足KKT 最优性条件(24)证明由引理 1 可得，L 是关于 x 和 y 的强凸函数 令 1=min(UTU)，2=min(LTL)，可得L(xk，yk，k)L(xk+1，yk，k)(xk xk+1)T(xL(xk+1，yk，k)+12xk xk+1I22(26)由于 xk+1使得xL(x，y，)xk+1=0，可得L(xk，yk，k)L(xk+1，yk，k)12xk xk+122，(27)同理可得L(xk+1，yk，k)L(xk+1，yk+1，k)01平顶山学院学报2023 年22yk yk+122，(28)由式(11)可得 Uxk+1 yk+1=1(k+1k)，所以L(xk+1，yk+1，k)L(xk+1，yk+1，k+1)=(k+1)T(Uxk+1 yk+1)(k)T(Uxk+1 yk+1)=1(k+1 k)22(29)通过式(27)、式(28)、式(29)，可以得到L(xk，yk，k)L(xk+1，yk+1，k+1)=L(xk，yk，k)L(xk+1，yk，k)+L(xk+1，yk，k)L(xk+1，yk+1，k)+L(xk+1，yk+1，k)L(xk+1，yk+1，k+1)12xk xk+122+22yk yk+1221(k+1 k)22(30)此外，对任意的 x，y，kn，L(x，y，k)是有界的 又因为 k 是有界的，故由引理 2 可得L(x，y，k)=12Ly (1 )v22(k)T(Ux y)+2Ux y22=12Ly (1 )v22+2Ux y 1k2212k22，(31)对不等式(30)进行累加求和，得到k=012xk xk+122+22yk yk+12()2k=01(k+1 k)22(32)因为k=0k+1 k22，所以k+1 k0并且k=012xk xk+122+22yk yk+12()2 (33)由式(33)可知，xk+1 xk0，yk+1 yk0 再由式(11)，可以得到UTU(xk+1 xk)=UT(yk+k)UTUxk=UTk UT(Uxk yk)，(34)(LTL+I)(yk+1 yk)=Uxk+1 k+(1 )LTv LTLyk yk=(Uxk+1 yk+1)+(yk+1 yk)LT(Lyk(1 )v)k，(35)k+1 k=(Uxk+1 yk+1)(36)由于 k+1 k0，根据式(36)可知，Uxk+1yk+10 将其带入式(34)、式(35)，可得UTk0，(37)LT(Lyk(1 )v)+k0(38)因为 k 是有界的，所以(xk，yk)也有界(x，y)是(xk，yk)的极限点，所以存在子序列(xki，yki)收敛于(x，y)k 有界意味着存在子序列 ki 收敛于 因此，(x，y，)是由算法生成的序列(xk，yk，k)的一个极限点，对此满足Ux y=0，UT=0，LT(Ly (1 )v)=(39)因此该极限点满足 KKT 最优性条件4数值实验数值实验在双核处理器(1 60 GHz，8 GBAM)上的 MATLAB 2018a 中进行 IT 表示迭代次数，CPU 表示运行时间，终止准则为 ES 1012，定义为(1 )v (I P)xk2(1 )v2 将本文提出的算法与 IO4 353、PIO5 23、MPIO6 91进行对比，在这三个实验中，内迭代次数都为 2，而在 MPIO 算法中取m=3 进行实验表1 列举了3 个测试矩阵，可从 https:/sparsetamu edu/获得 为了公平起见，初始向量 x0=v，取=0 5，阻尼因子取 =0 90，0 95，0 99，在 LU-ADMM 中取 =105表 1测试矩阵的性质矩阵大小非零元稠密度Harvard500500 5002 6360 054Minnesota2 642 2 6426 6060 095Califonia9 664 9 66416 1500 001 7从表 2、表 3、表 4 可以看出，通过比较这四种算法的迭代次数 IT 和 ES，得出本文的 LUADMM算法是可行的 并且可以明显地观察到，的取值不会影响到 LUADMM 算法的迭代次数，所以采用本文提出的算法时，并不需要考虑 的取值，无论 是否无限地接近 1 通过对比三个测试矩阵 Har-vard500、Minnesota 和 Califonia，在调用 LUADMM算法时，它们的迭代步数在 6 到 8 之间 由此可见，测试矩阵的大小并不会影响该算法的迭代次数 但是该算法明显的缺点是 CPU 时间过长11第 2 期吴文军，唐嘉:基于 ADMM 迭代法的 Pageank 问题的求解表 2Harvard500 矩阵的测试结果参数IOPIOMPIOLUADMM0 90ITCPUES1140 052 98 56 1013690 053 98 88 1013500 052 47 13 101370 267 81 45 10130 95ITCPUES2130 031 49 05 10131280 034 29 69 1013920 033 79 10 101380 243 554 37 10130 99ITCPUES5530 041 59 94 10133320 037 39 88 10132380 034 29 10 1013100 256 55 83 1013表 3Minnesota 矩阵的测试结果参数IOPIOMPIOLUADMM0 90ITCPUES1250 017 19 98 1013760 016 09 47 1013550 015 77 88 101360 263 42 56 10150 95ITCPUES2580 021 59 54 10131560 020 09 20 10131120 019 78 48 101360 250 85 48 10150 99ITCPUES1 2660 060 39 96 10137610 055 69 78 10135440 050 89 71 101370 257 02 26 1014表 4Califonia 矩阵的测试结果参数IOPIOMPIOLUADMM0 90ITCPUES640 032 09 63 1013390 033 57 73 1013280 031 87 29 101360 618 11 08 10140 95ITCPUES860 045 48 44 1013520 039 77 04 1013370 037 67 81 101370 611 65 85 10140 99ITCPUES1180 046 78 42 1013710 043 87 03 1013500 044 09 87 101370 607 17 31 1014通过图 1 和图 2 可以得到，无论是 Minnesota矩阵还是 Califonia 矩阵，它们的迭代次数随着 的值变大而增加 由此可得，LUADMM 算法的迭代次数与 的取值有关，取值越大迭代次数越多结合表 2、表 3、表 4 和图 1、图 2 可知，本文提出的 LUADMM 算法的迭代次数不会因为矩阵的变化而变化，而是随着 的取值变化而变化图 1Minnesota 矩阵在 取不同值时的迭代次数21平顶山学院学报2023 年图 2Califonia 矩阵在 取不同值时的迭代次数5结论本文 主 要 是 使 用 一 种 交 替 方 向 乘 子 法(ADMM)来 求 解 Pageank 问 题，为 找 到 满 足ADMM 的情况，先对 Pageank 问题导出的线性方程组的矩阵进行 LU 分解，将其转化为最小二乘问题，之后利用 ADMM 的迭代格式来求解 Pageank问题 在适当的条件下证明了所提方法的收敛性实验证明了本文所提出的方法是有效的，并且是可行的参考文献:1 TIAN Z L，LIU Y，ZHANG Y，et al The general inner-out-er iteration method based on regular splittings for the Pag-eank problemJ Applied Mathematics and Computa-tion，2019，356 2 PAGE L，BIN S，MOTWAMI，et al The Pageank cita-tion ranking:bringing order to the web Palo Alto:Stanford University，1998 3 SHEPELYANSKY D L，ZHIOV O V Towards Google ma-trix of brainJ Physics Letters A，2010，374(31/32):3206 3209 4 GLEICH D F，GAY A P，GEIF C，et al An inner-outeriteration for computing PageankJ SIAM Journal onScientific Computing，2010，32(1)5 GU C Q，XIE F，ZHANG K A two-step matrix splitting it-eration for computing PageankJ Journal of Computa-tional and Applied Mathematics，2015，278(C)6 WEN C，HUANG T Z，SHEN Z L A note on the two-stepmatrix splitting iteration for computing Pageank J Jour-nal of Computational and Applied Mathematics，2017，315(C)7NESTEOV Y Introductory lectures on convex optimiza-tion M Dordrecht:Kluwer Academic Publishers，2004 8 BOYD S，VANDENBEGHE L Convex optimization M Cambridge:Cambridge University Press，2004 9 BOYD S，PAIKH N，CHU E，et al Distributed optimiza-tion and statistical learning via the alternating directionmethod of multipliersJ Foundations and Trends in Ma-chine Learning，2011，3(1)10 GLOWINSKI，TALLEC P L Augmented Lagrangianmethods for the solution of V ariational problems Madison:University of Wisconsin，1987 11 YANG J F，YUAN X M An inexact alternating directionmethod for trace norm regularized least squares problem Nanjing:Nanjing University，2010 12THEMELIS A，PATINOS P Douglas-achford splittingand ADMM for nonconvex optimization:tight convergenceresults J SIAM Journal on Optimization，2020，30(1):149 181 13 JIA Z H，HUANG J，WU Z M An incremental aggrega-ted proximal ADMM for linearly constrained nonconvex op-timization with application to sparse logistic regressionproblems J Journal of Computational and Applied Math-ematics，2021，390:113384(责任编辑:王彦江)31第 2 期吴文军，唐嘉:基于 ADMM 迭代法的 Pageank 问题的求解