求解PageRank 问题的Arnoldi 松弛两步分裂算法

顾传青,付友花, 王金波

(1. 上海大学 理学院, 上海200444;2. 保密通信重点实验室, 成都610041)

PageRank 算法是著名的超链分析算法[1], 是当今搜索引擎使用最重要的排序算法之一.PageRank 问题就是求解Google 矩阵首特征值1 所对应的特征向量. Google 矩阵的定义如下:

式中, α ∈(0,1)是阻尼因子, E = veT, e = (1,1,··· ,1)T∈Rn, v = e/n, P是一个列随机矩阵. PageRank 向量x 是Google 矩阵的主特征向量, 满足

为解决PageRank 问题, 最经典的算法是易于计算的幂法[1]. 但当阻尼因子α 接近1 时,幂法收敛的速度非常慢. 为了改进算法的收敛速度, Gleich 等[2]利用Richardson 迭代法提出了内外迭代法. Gu 等[3]将幂法和内外迭代法相结合,提出了两步分裂迭代(power-inner-outer,PIO)算法. 随后，Xie 等[4]在PIO 算法基础上提出了松弛两步分裂(relaxed PIO, RPIO)算法.用于PageRank 问题的Arnoldi 算法[5]是一种重启的Arnoldi 算法[6]. Wu 等[7]将幂法和重启的Arnoldi 算法相结合, 提出了深度重启的Arnoldi 算法. 本工作在文献[3]的基础上, 在原有的PIO 算法中加入一个新的松弛参数, 并且运用深度重启的Arnoldi 算法来加速算法的收敛性, 得到了一个Arnoldi 松弛两步分裂(Arnoldi-RPIO)算法. 与原有的内外迭代法[2]和PIO 算法[3]相比, 本工作给出的数值算例显示了新算法的有效性.

1 内外迭代法和深度重启的Arnoldi 算法

1.1 内外迭代法

根据式(1)和(2), 特征值问题可以转化为求线性系统

的解, 其中eTx = 1. 考虑到问题在阻尼因子α 越小时越容易求解PageRank, 线性系统(3)可以写成如下的等价方程:

并通过不动点迭代求解式(4), 即

式(5)可称为外迭代. 但是, 求解系数矩阵I -βP 的线性系统仍然比较困难. 文献[2]提出了使用Richardson 迭代法来计算xk+1, 将式(5)中的右端项记为f, 即

定义内线性系统

然后使用Richardson 内迭代:

外迭代的停止条件为

内迭代的停止条件为

1.2 PIO 算法

在内外迭代法的基础上, Gu 等[3]提出了PIO 算法来求解PageRank 问题. PIO 算法的迭代格式如下: 给定一个初始向量x(0), 计算

直到向量序列{xk}收敛, 其中0 ＜β ＜α ＜1.

1.3 深度重启的Arnoldi 算法

考虑如下的特征值问题:

式中, A ∈Cn×n, (λ,x)是矩阵A 的特征对. 给定一个‖v1‖2= 1 的初始向量, 由Arnoldi 过程生成一组Krylov 子空间

的一组正交基v1,v2,··· ,vm. 在实际求解中,可以用修正的Gram-Schmidt 算法来得到Krylov子空间[8]. 由Arnoldi 过程可以得到

式中, Vm= (v1,v2,··· ,vm)是Krylov 子空间的一组正交基, Hm= {hi,j}m×m∈Cm×m是一个m×m 的上Hessenberg 矩阵, em= (0,0,··· ,0,1)T,Hm∈C(m+1)×m是一个如下的上Hessenberg 矩阵:

深度重启的Arnoldi 算法的具体步骤如下.

步骤1 给定Krylov 子空间维数m, 所要求的特征对个数p ≤m, 控制精度tol 和单位初始向量v1.

步骤2 运行m 步关于A 和v1的Arnoldi 过程, 得到计算Hm的所有特征对i = 1,2,··· ,m. 将Hm的特征值按模数递减排序, 从中选出p 个最大的特征对, 转向步骤6.

步骤3 将vp+1作为初始值运行m-p 步Arnoldi 过程, 得到. 计算Hm的所有特征对(), i = 1,2,··· ,m. 将Hm的特征值按模数递减排序, 从中选出前p 个特征对.

步骤4 检验收敛性. 对于得到的最大的Ritz 值λ1和Ritz 向量y1计算其残量. 若满足计算精度, 则选取x1=Vmy1作为PageRank 向量的近似解, 算法停止, 否则继续.

步骤5 将选出的p 个特征向量yi, i=1,2,···, p, 正交化为2 模单位向量, 得到m×p 矩阵Wp= [w1,w2,··· ,wp]. 如果特征向量yi是复向量, 需要实部和虚部分开存取, 构造m×p阶矩阵.

步骤6 通过在Wp的底部增加一行0 形成(m+1)×p 的矩阵Wp= [Wp;0]. 然后令是(m +1) × (m +1) 阶单位矩阵Im+1的最后一列. 显然(m+1)×(p+1)的矩阵Wp+1是正交矩阵.

步骤7 使用旧的Vm+1和Hm形成新的Vm+1和Vm+1Wp+1. 然后令, 返回步骤3.

2 深度重启的Arnoldi-RPIO 算法

本工作提出利用深度重启的Arnoldi-RPIO 算法来求解PageRank 问题. 迭代格式如下:

Arnoldi-RPIO 算法在PIO 算法的基础上引入了一个新的参数γ, 并使用Arnoldi 算法进行深度重启来加速算法收敛.

Arnoldi-RPIO 算法的具体步骤如下.

步骤1 给定初始化向量v, 选取Krylov 子空间的维数m, 欲保留的特征向量数p, 外迭代的收敛精度τ, 内迭代的收敛精度η, 参数α, β, 控制内外迭代与Arnoldi 算法阶段的转换参数α1和α2, maxit, restart=0, 外迭代的误差r =1, 内迭代的误差d=1,d0=d.

步骤2 运行1.3 节算法2～3 次. 若是第1 次运行, 则运行步骤2～7, 否则运行步骤3～7.如果残差范数满足收敛精度, 则算法停止, 否则继续.

步骤3 将由深度重启的Arnoldi 算法得到的PageRank 向量的近似解x1作为内外迭代法的初始向量.

end if r ＜τ, 算法停止, 否则转向步骤2.

关于Arnoldi-RPIO 算法有如下几点说明: ①参数α1, α2, restart, maxit 是用来控制内外迭代法和深度重启的Arnoldi 算法的转换; ②定义rcurr为一次内迭代后的残差范数, rpre为一次内迭代前的残差范数, ratio = rcurr/rpre为内迭代前后相邻残差范数比; 同时定义ratio1为内迭代中当前步的残差范数和上一步的残差范数之比(这里给出ratio 和ratio1的定义是为了保证算法的运算效率); ③通过比较每一次进入内外迭代前后的残差范数比ratio = rcurr/rpre与α1的大小, 来决定是否执行restart=restart+1; 通过maxit 控制restart 的次数, 如果restart＞maxit, 则中止内外迭代法阶段, 进入深度重启的Arnoldi 算法阶段, 否则继续运行内外迭代法.

在实际求解中, α1和α2的选择非常关键. 由于ratio1=d/d0, ratio=r/r0, 所以这两个参数很自然地小于α (幂法的收敛速率是α). 本工作中选择α1=α-0.1, α2=α-0.1. 事实上,如果令γ =1, 则Arnoldi-RPIO 算法就等价于文献[7]提出的Arnoldi 算法.

3 Arnoldi-RPIO 算法的收敛性分析

RPIO 算法和深度重启的Arnoldi 算法的收敛性分别在文献[4, 9-10]中被证明. 不妨假设A 是一个对角化矩阵, 并且x1和子空间Km的距离定义为

式中, Lm-1表示次数小于m-1 的多项式集, 并且p(λ1)=1, σ(A)表示A 的特征值集. 不失一般性, 假设A 的特征值降序排列:

引理1[9]假设A 是一个对角化矩阵, Arnoldi 算法的初始向量v1可以展开为v1=是特征基, 且‖xi‖2= 1, i = 1,2,··· ,n,, γi/= 0, 则如下的不等式成立:

式中, Pm是子空间Km(A,v1)上的正交投影,

将由深度重启的Arnoldi 算法得到的v1作为RPIO 算法的初始向量. RPIO 算法得到向量v*1=ωGkv1, 其中k ≥maxit,

ω 是正规化因子. 在Arnoldi-RPIO 算法的下一步循环中, v*1将被作为Arnoldi 过程的初始向量, 这样就可以得到一个新的Krylov 子空间:

下面的定理证明了本工作给出的新算法的收敛性.

因为

可以推出

假设πi是P 的一个特征值,可知是(I -βP)-1的一个特征值. 由

和λ1=1, 得到Gx1=x1,Gkx1=x1. 再令

对于i=2,3,··· ,n, 由于|λi|≤α, γ ≥1, 可以推出

从而得到如下两个关系式:

现在令p(λ)=q(λ)/q(1), 则p(1)=1, 从而得到

因此, 证明了

4 数值算例

本工作选择Stanford-Berkeley 矩阵作为测试矩阵, 给出数值算例来展示Arnoldi-RPIO算法的有效性, 并将其与内外迭代法和PIO 算法进行对比. 数值算例是在内存为4 GB, 处理器为2.53 GHz Intel(R)Core(TM)i3 M CPU 的电脑上用MATLAB(R2014a)程序进行的. 这里, Mat-Vec 表示矩阵向量乘积的次数, CPU 表示运算时间, 单位为s.

为保证实验的公平性, 阻尼因子α 的取值分别为0.990, 0.993, 0,995, 0.998, β =0.5, 内外终止条件均为η =10-2, τ =10-8. 在Arnoldi-RPIO 算法中, 参数α1=α-0.1, α2=α-0.1.为了描述Arnoldi-RPIO 算法的有效性, 定义

表1 列出了测试矩阵的特征, 包括矩阵的维数(n)、非零元个数(nn)和稠密度(den), 并定义表2 和图1 列出了3 种算法的数值结果.

表1 测试矩阵的特征Table 1 Characteristics of test matrix

表2 关于Stanford-Berkeley 矩阵3 种算法的数值结果Table 2 Numerical results of three algorithms for Stanford-Berkeley matrix

图1 关于Stanford-Berkeley 矩阵3 种算法的收敛性分析, τ =10-8Fig.1 Convergence analysis of three algorithms for Standford-Berkeley matrices, τ =10-8

图1 描述了阻尼因子α 取不同值时3 种算法的收敛轨迹. 算例中Arnoldi-RPIO 算法选取的参数为γ = 6,m = 8,p = 4, maxit=40. 观察表2 的数值结果发现, Arnoldi-RPIO 算法在Mat-Vec 和CPU 两个方面的表现都是最好的. 对于Stanford-Berkeley 矩阵, 数值结果显示: 当α = 0.998 时, Arnoldi-RPIO 算法耗时250.101 4 s 达到的收敛精度, PIO 算法需耗时575.583 5 s 才能达到; 并且阻尼因子α 取值越接近于1 时, 本工作给出的Arnoldi-RPIO 算法的计算效果就越好.