求解PageRank 問題的Arnoldi 松弛兩步分裂算法

顧傳青,付友花, 王金波

(1. 上海大學(xué) 理學(xué)院, 上海200444;2. 保密通信重點(diǎn)實(shí)驗(yàn)室, 成都610041)

PageRank 算法是著名的超鏈分析算法[1], 是當(dāng)今搜索引擎使用最重要的排序算法之一.PageRank 問題就是求解Google 矩陣首特征值1 所對(duì)應(yīng)的特征向量. Google 矩陣的定義如下:

式中, α ∈(0,1)是阻尼因子, E = veT, e = (1,1,··· ,1)T∈Rn, v = e/n, P是一個(gè)列隨機(jī)矩陣. PageRank 向量x 是Google 矩陣的主特征向量, 滿足

為解決PageRank 問題, 最經(jīng)典的算法是易于計(jì)算的冪法[1]. 但當(dāng)阻尼因子α 接近1 時(shí),冪法收斂的速度非常慢. 為了改進(jìn)算法的收斂速度, Gleich 等[2]利用Richardson 迭代法提出了內(nèi)外迭代法. Gu 等[3]將冪法和內(nèi)外迭代法相結(jié)合,提出了兩步分裂迭代(power-inner-outer,PIO)算法. 隨后，Xie 等[4]在PIO 算法基礎(chǔ)上提出了松弛兩步分裂(relaxed PIO, RPIO)算法.用于PageRank 問題的Arnoldi 算法[5]是一種重啟的Arnoldi 算法[6]. Wu 等[7]將冪法和重啟的Arnoldi 算法相結(jié)合, 提出了深度重啟的Arnoldi 算法. 本工作在文獻(xiàn)[3]的基礎(chǔ)上, 在原有的PIO 算法中加入一個(gè)新的松弛參數(shù), 并且運(yùn)用深度重啟的Arnoldi 算法來加速算法的收斂性, 得到了一個(gè)Arnoldi 松弛兩步分裂(Arnoldi-RPIO)算法. 與原有的內(nèi)外迭代法[2]和PIO 算法[3]相比, 本工作給出的數(shù)值算例顯示了新算法的有效性.

1 內(nèi)外迭代法和深度重啟的Arnoldi 算法

1.1 內(nèi)外迭代法

根據(jù)式(1)和(2), 特征值問題可以轉(zhuǎn)化為求線性系統(tǒng)

的解, 其中eTx = 1. 考慮到問題在阻尼因子α 越小時(shí)越容易求解PageRank, 線性系統(tǒng)(3)可以寫成如下的等價(jià)方程:

并通過不動(dòng)點(diǎn)迭代求解式(4), 即

式(5)可稱為外迭代. 但是, 求解系數(shù)矩陣I -βP 的線性系統(tǒng)仍然比較困難. 文獻(xiàn)[2]提出了使用Richardson 迭代法來計(jì)算xk+1, 將式(5)中的右端項(xiàng)記為f, 即

定義內(nèi)線性系統(tǒng)

然后使用Richardson 內(nèi)迭代:

外迭代的停止條件為

內(nèi)迭代的停止條件為

1.2 PIO 算法

在內(nèi)外迭代法的基礎(chǔ)上, Gu 等[3]提出了PIO 算法來求解PageRank 問題. PIO 算法的迭代格式如下: 給定一個(gè)初始向量x(0), 計(jì)算

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

1.3 深度重啟的Arnoldi 算法

考慮如下的特征值問題:

式中, A ∈Cn×n, (λ,x)是矩陣A 的特征對(duì). 給定一個(gè)‖v1‖2= 1 的初始向量, 由Arnoldi 過程生成一組Krylov 子空間

的一組正交基v1,v2,··· ,vm. 在實(shí)際求解中,可以用修正的Gram-Schmidt 算法來得到Krylov子空間[8]. 由Arnoldi 過程可以得到

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

深度重啟的Arnoldi 算法的具體步驟如下.

步驟1 給定Krylov 子空間維數(shù)m, 所要求的特征對(duì)個(gè)數(shù)p ≤m, 控制精度tol 和單位初始向量v1.

步驟2 運(yùn)行m 步關(guān)于A 和v1的Arnoldi 過程, 得到計(jì)算Hm的所有特征對(duì)i = 1,2,··· ,m. 將Hm的特征值按模數(shù)遞減排序, 從中選出p 個(gè)最大的特征對(duì), 轉(zhuǎn)向步驟6.

步驟3 將vp+1作為初始值運(yùn)行m-p 步Arnoldi 過程, 得到. 計(jì)算Hm的所有特征對(duì)(), i = 1,2,··· ,m. 將Hm的特征值按模數(shù)遞減排序, 從中選出前p 個(gè)特征對(duì).

步驟4 檢驗(yàn)收斂性. 對(duì)于得到的最大的Ritz 值λ1和Ritz 向量y1計(jì)算其殘量. 若滿足計(jì)算精度, 則選取x1=Vmy1作為PageRank 向量的近似解, 算法停止, 否則繼續(xù).

步驟5 將選出的p 個(gè)特征向量yi, i=1,2,···, p, 正交化為2 模單位向量, 得到m×p 矩陣Wp= [w1,w2,··· ,wp]. 如果特征向量yi是復(fù)向量, 需要實(shí)部和虛部分開存取, 構(gòu)造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 算法的基礎(chǔ)上引入了一個(gè)新的參數(shù)γ, 并使用Arnoldi 算法進(jìn)行深度重啟來加速算法收斂.

Arnoldi-RPIO 算法的具體步驟如下.

步驟1 給定初始化向量v, 選取Krylov 子空間的維數(shù)m, 欲保留的特征向量數(shù)p, 外迭代的收斂精度τ, 內(nèi)迭代的收斂精度η, 參數(shù)α, β, 控制內(nèi)外迭代與Arnoldi 算法階段的轉(zhuǎn)換參數(shù)α1和α2, maxit, restart=0, 外迭代的誤差r =1, 內(nèi)迭代的誤差d=1,d0=d.

步驟2 運(yùn)行1.3 節(jié)算法2～3 次. 若是第1 次運(yùn)行, 則運(yùn)行步驟2～7, 否則運(yùn)行步驟3～7.如果殘差范數(shù)滿足收斂精度, 則算法停止, 否則繼續(xù).

步驟3 將由深度重啟的Arnoldi 算法得到的PageRank 向量的近似解x1作為內(nèi)外迭代法的初始向量.

end if r ＜τ, 算法停止, 否則轉(zhuǎn)向步驟2.

關(guān)于Arnoldi-RPIO 算法有如下幾點(diǎn)說明: ①參數(shù)α1, α2, restart, maxit 是用來控制內(nèi)外迭代法和深度重啟的Arnoldi 算法的轉(zhuǎn)換; ②定義rcurr為一次內(nèi)迭代后的殘差范數(shù), rpre為一次內(nèi)迭代前的殘差范數(shù), ratio = rcurr/rpre為內(nèi)迭代前后相鄰殘差范數(shù)比; 同時(shí)定義ratio1為內(nèi)迭代中當(dāng)前步的殘差范數(shù)和上一步的殘差范數(shù)之比(這里給出ratio 和ratio1的定義是為了保證算法的運(yùn)算效率); ③通過比較每一次進(jìn)入內(nèi)外迭代前后的殘差范數(shù)比ratio = rcurr/rpre與α1的大小, 來決定是否執(zhí)行restart=restart+1; 通過maxit 控制restart 的次數(shù), 如果restart＞maxit, 則中止內(nèi)外迭代法階段, 進(jìn)入深度重啟的Arnoldi 算法階段, 否則繼續(xù)運(yùn)行內(nèi)外迭代法.

在實(shí)際求解中, α1和α2的選擇非常關(guān)鍵. 由于ratio1=d/d0, ratio=r/r0, 所以這兩個(gè)參數(shù)很自然地小于α (冪法的收斂速率是α). 本工作中選擇α1=α-0.1, α2=α-0.1. 事實(shí)上,如果令γ =1, 則Arnoldi-RPIO 算法就等價(jià)于文獻(xiàn)[7]提出的Arnoldi 算法.

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

RPIO 算法和深度重啟的Arnoldi 算法的收斂性分別在文獻(xiàn)[4, 9-10]中被證明. 不妨假設(shè)A 是一個(gè)對(duì)角化矩陣, 并且x1和子空間Km的距離定義為

式中, Lm-1表示次數(shù)小于m-1 的多項(xiàng)式集, 并且p(λ1)=1, σ(A)表示A 的特征值集. 不失一般性, 假設(shè)A 的特征值降序排列:

引理1[9]假設(shè)A 是一個(gè)對(duì)角化矩陣, 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,

ω 是正規(guī)化因子. 在Arnoldi-RPIO 算法的下一步循環(huán)中, v*1將被作為Arnoldi 過程的初始向量, 這樣就可以得到一個(gè)新的Krylov 子空間:

下面的定理證明了本工作給出的新算法的收斂性.

因?yàn)?/p>

可以推出

假設(shè)πi是P 的一個(gè)特征值,可知是(I -βP)-1的一個(gè)特征值. 由

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

對(duì)于i=2,3,··· ,n, 由于|λi|≤α, γ ≥1, 可以推出

從而得到如下兩個(gè)關(guān)系式:

現(xiàn)在令p(λ)=q(λ)/q(1), 則p(1)=1, 從而得到

因此, 證明了

4 數(shù)值算例

本工作選擇Stanford-Berkeley 矩陣作為測試矩陣, 給出數(shù)值算例來展示Arnoldi-RPIO算法的有效性, 并將其與內(nèi)外迭代法和PIO 算法進(jìn)行對(duì)比. 數(shù)值算例是在內(nèi)存為4 GB, 處理器為2.53 GHz Intel(R)Core(TM)i3 M CPU 的電腦上用MATLAB(R2014a)程序進(jìn)行的. 這里, Mat-Vec 表示矩陣向量乘積的次數(shù), CPU 表示運(yùn)算時(shí)間, 單位為s.

為保證實(shí)驗(yàn)的公平性, 阻尼因子α 的取值分別為0.990, 0.993, 0,995, 0.998, β =0.5, 內(nèi)外終止條件均為η =10-2, τ =10-8. 在Arnoldi-RPIO 算法中, 參數(shù)α1=α-0.1, α2=α-0.1.為了描述Arnoldi-RPIO 算法的有效性, 定義

表1 列出了測試矩陣的特征, 包括矩陣的維數(shù)(n)、非零元個(gè)數(shù)(nn)和稠密度(den), 并定義表2 和圖1 列出了3 種算法的數(shù)值結(jié)果.

表1 測試矩陣的特征Table 1 Characteristics of test matrix

表2 關(guān)于Stanford-Berkeley 矩陣3 種算法的數(shù)值結(jié)果Table 2 Numerical results of three algorithms for Stanford-Berkeley matrix

圖1 關(guān)于Stanford-Berkeley 矩陣3 種算法的收斂性分析, τ =10-8Fig.1 Convergence analysis of three algorithms for Standford-Berkeley matrices, τ =10-8

圖1 描述了阻尼因子α 取不同值時(shí)3 種算法的收斂軌跡. 算例中Arnoldi-RPIO 算法選取的參數(shù)為γ = 6,m = 8,p = 4, maxit=40. 觀察表2 的數(shù)值結(jié)果發(fā)現(xiàn), Arnoldi-RPIO 算法在Mat-Vec 和CPU 兩個(gè)方面的表現(xiàn)都是最好的. 對(duì)于Stanford-Berkeley 矩陣, 數(shù)值結(jié)果顯示: 當(dāng)α = 0.998 時(shí), Arnoldi-RPIO 算法耗時(shí)250.101 4 s 達(dá)到的收斂精度, PIO 算法需耗時(shí)575.583 5 s 才能達(dá)到; 并且阻尼因子α 取值越接近于1 時(shí), 本工作給出的Arnoldi-RPIO 算法的計(jì)算效果就越好.