求解PageRank問題的重啟GMRES修正的多分裂迭代法*

肖文可，陳星玎

（北京工商大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，北京 100048）

引 言

Google 搜索引擎以其著名的PageRank 算法成為近年來最成功的搜索引擎之一.Brin 和Page 于1998 年提出了PageRank 算法[1]，該算法通過分析網(wǎng)絡(luò)的鏈接結(jié)構(gòu)獲得網(wǎng)頁的重要程度排名.PageRank 問題就是求解Google 矩陣A的首特征值1 所對應(yīng)的特征向量，即線性系統(tǒng)的解.Google 矩陣A是矩陣P與矩陣veT的凸組合，即

其中α ∈(0,1)是阻尼因子，e=(1,1,···,1)T，v=e/n，矩陣P是網(wǎng)絡(luò)的超鏈接結(jié)構(gòu)[2]所定義的隨機(jī)矩陣，n為矩陣P的維數(shù).如圖1 所示，由4 個網(wǎng)頁構(gòu)成的超鏈接結(jié)構(gòu)，它所定義的隨機(jī)矩陣P4，其元素pij表示從網(wǎng)頁i鏈接到網(wǎng)頁j的概率，元素值都落在[0,1]之間，且滿足每行和為1，

圖1 4 個網(wǎng)頁的超鏈接結(jié)構(gòu)Fig.1 The hyperlink structure of 4 pages

Google 矩陣A是 一個高維不可約、非周期的隨機(jī)矩陣.當(dāng) α越接近于1 時，網(wǎng)絡(luò)矩陣P占比越大，PageRank 問題越接近于真實情況，但越難以求解.

最初求解PageRank 問題(1)的算法是冪法[3].冪法單次迭代計算量小，但收斂速度慢.尤其當(dāng)阻尼因子α接近1 時，冪法幾乎不收斂.隨后產(chǎn)生了冪法的加速方法，如二次外推法[2]、修正的冪外推法[4]等.另一方面，利用eTx=1，由式(1)、(2) 可得

那么，特征值問題(1)可以轉(zhuǎn)化為下面線性方程組的求解，即

其中I表示n階單位矩陣.針對線性方程組(4)，學(xué)者們構(gòu)造了許多迭代方法進(jìn)行求解，如內(nèi)外迭代(inner-outer,IO)法[5]、兩步分裂迭代法[6]、廣義二級分裂迭代法[7]、多步冪法修正的廣義二級分裂迭代法[8]和多分裂迭代(MSI)法[9]等.近年來，基于Krylov 子空間的迭代方法由于其特有的優(yōu)勢被廣泛用于線性方程組的求解.雖然Krylov 子空間方法的收斂速度快，但是這類方法的單次迭代計算量大，存儲空間消耗大，在使用時往往需要與已有的基本迭代法相結(jié)合[10-12].本文將Krylov 子空間方法中的重啟GMRES 方法[13]與MSI 法相結(jié)合，提出了一種重啟GMRES 修正的多分裂迭代(GMSI)方法來求解線性方程組(4).

本文的結(jié)構(gòu)如下：第1 節(jié)簡要介紹了MSI 法和重啟GMRES 方法；第2 節(jié)給出了重啟GMRES 修正的多分裂迭代方法的計算流程；第3 節(jié)給出了新方法的收斂性分析；第4 節(jié)通過數(shù)值實驗驗證了新方法的有效性；最后進(jìn)行總結(jié).

1 MSI 法和重啟GMRES 方法

本節(jié)中，我們將簡要回顧MSI 法和重啟GMRES 方法.

1.1 MSI 法

Gu 等[9]在線性方程組(4)中引入兩個參數(shù) β1和 β2，提出了一種基于多分裂的兩步迭代法，即MSI 法.

MSI 方法 給定初始向量x(0)，對k=0,1,2,···，計算如下迭代：

直至{x(k)}收斂.其中參數(shù) βi滿足0 ≤βi＜α(i=1,2)，，即對每一步迭代解進(jìn)行歸一化.

我們采用IO 法求解線性系統(tǒng)(5).定義如下兩個內(nèi)線性系統(tǒng)：

其中

f1=(α-β1)Px(k)+(1-α)v,

f2=(α-β2)Px(k+1/2)+(1-α)v.

然后，采用Richardson 迭代：

計算線性系統(tǒng)(5)中的x(k+1/2)和x(k+1)，注意y(1

l1)=x(k+1/2)，y(2l2)=x(k+1).

MSI 方法的終止標(biāo)準(zhǔn)如下：假設(shè)η和τ分別代表內(nèi)、外迭代終止指標(biāo)，當(dāng)

時，內(nèi)迭代(8)終止；當(dāng)

時，外迭代(4)終止.

1.2 重啟GMRES 方法

Saad 等[13]提出的GMRES 方法被廣泛應(yīng)用于線性方程組的求解.由于GMRES 方法的時間復(fù)雜度和空間復(fù)雜度都較大，在實際計算中往往采用重啟GMRES 方法[13]，即GMRES(m)，m代表重啟數(shù).重啟GMRES 方法的具體計算步驟將由算法1 給出.

算法1GMRES(m)方法

步1 開始.給定重啟數(shù)m，初始向量x0，迭代終止指標(biāo)τ，計算

r0=(1-α)v-(I-αP)x0, β=‖r0‖2,v1=r0/β.

步2 迭代.對j=1,2,···,m,

步3 計算近似解.xm=x0+Vmym，其中正交基矩陣Vm={vj}1≤j≤m，向量ym是最小二乘問題miny∈Rm‖βe1-y‖2的解，e1=(1,0,···,0)T，={hi j}1≤i≤m+1,1≤j≤m.

步4 重啟.計算rm=(1-α)v-(I-αP)xm，若‖rm‖2＜τ，則算法停止；否則x0←xm，v1=rm/‖rm‖2轉(zhuǎn)至步2.

2 重啟GMRES-MSI 方法

本文把重啟GMRES 方法與MSI 法相結(jié)合，提出了一種求解PageRank 問題(4)的新方法，即重啟GMRES 修正的多分裂迭代法，簡記為GMRES(m)-MSI 方法.該方法的基本思想是：給定一個初始向量x0，首先利用GMRES(m)方法計算一個近似解x1，如果x1未達(dá)到收斂精度，那么將x1作為MSI 方法的初始向量繼續(xù)求解.重復(fù)上述過程，直至達(dá)到收斂精度.算法流程如圖2 所示.

圖2 重啟GMRES-MSI 方法的計算流程圖Fig.2 The calculation flow chart for the restarted GMRES-MSI method

在GMRES(m)-MSI 方法中，用于控制GMRES(m)方法和MSI 方法之間轉(zhuǎn)換的三個參數(shù)分別記為 α1，ν和?.ν表示MSI 方法的當(dāng)前迭代次數(shù)，?表示MSI 方法的最大迭代次數(shù).由于Google 矩陣的第二特征值的模小于等于α[14]，所以 α1的選擇應(yīng)該小于α，取α1=α-0.1.設(shè)rpre與rcur分別表示前后兩步MSI 方法的殘差范數(shù)，定義ζ=rcur/rpre，通過比較ζ和 α1的大小來決定是否執(zhí)行ν=ν+1.若ν ＞?，則終止MSI 方法，進(jìn)入GMRES(m)方法；否則，繼續(xù)運行MSI 方法.具體計算步驟見算法2.

算法2GMRES(m)-MSI 方法

步1 開始.選取初始向量x0，重啟數(shù)m，內(nèi)、外迭代終止指標(biāo)分別為η和τ，參數(shù) β1，β2和 ?，初始誤差取r=1.

步2 運行GMRES(m)方法2 至3 次，若得到的近似解達(dá)到收斂精度，算法停止，否則繼續(xù).

步3 將GMRES(m)方法得到的近似解x1作為MSI 方法的初始向量.

算法停止，否則轉(zhuǎn)向步2.

3 算法的收斂性分析

引理1設(shè)x(k)是MSI 方法的第k步迭代解，那么第k+1步迭代解為x(k+1)=Gx(k)，迭代矩陣

G=(I-β2P)-1[(α-β2)P(I-β1P)-1(A-β1P)+(1-α)veT].

證明因為0 ≤β ＜α ＜1且P是隨機(jī)矩陣，所以I-βiP(i=1,2)是嚴(yán)格對角占優(yōu)矩陣，因此I-βiP(i=1,2)是可逆矩陣.

由于在MSI 方法中采用了歸一化處理，每一步的迭代解滿足eTx(k)=1，那么由式(5)可得

因此，迭代矩陣G為

G=(I-β2P)-1[(α-β2)P(I-β1P)-1(A-β1P)+(1-α)veT].□

將GMRES(m)方法計算的近似解x1作為MSI 方法的初始向量，由引理1 得，經(jīng)過l次MSI 方法，求得近似解x*1為

x*1=Glx1,l≥?.

然后，將x*1作為GMRES(m)方法的初始向量構(gòu)造新的Krylov 子空間Km，即

Km=span{v*1,(I-αP)v*1,···,(I-αP)m-1v*1},

其中

v*1=r*0/‖r*0‖2,r*0=(1-α)v-(I-αP)x*1.

設(shè)Pm是最高次數(shù)不超過m次的全體多項式的集合，不失一般性，假設(shè)矩陣A的特征值按照降序排列為[15]

1=|λ1|＞|λ2|≥··· ≥|λn-1|≥|λn|.

引理2[13]設(shè)I-αP可對角化，且(I-αP)=XΛX-1，這里Λ=diag(s1,s2,···,sn)，si(i=1,2,···,n)是矩陣I-αP的特征值，定義

那么采用GMRES(m)方法后的殘差rm滿足不等式

‖rm‖2≤κ2(X)ε(m)‖r0‖2,

其中

κ2(X)=‖X‖2‖X-1‖2.

引理3[15]設(shè)隨機(jī)矩陣P的特征值為{1,π2,···,πn}，則Google 矩陣A的特征值為{1,απ2,···,απn}.

引理4設(shè)隨機(jī)矩陣P的特征值為{1,π2,···,πn}，則矩陣

A-β1P=(α-β1)P+(1-α)veT

的特征值分別為1-β1和(α-β1)πi(i=2,3,···,n).

證明由于P是隨機(jī)矩陣，所以(1,e)為P的特征對.設(shè)Q=(e,L)是一個非奇異矩陣，Q的第一列為矩陣P的特征向量e.設(shè)那么

則有yTe=1，YTe=0.因此，

因此，YTPL有特征值π2,π3,···,πn(文獻(xiàn)[15]中，定理5.1).對矩陣A-β1P做相似變換，得

因此，矩陣A-β1P的特征值為1-β1和(α-β1)πi(i=2,3,···,n).□

下面，我們給出GMRES(m)-MSI 方法的收斂性分析.

定理1設(shè)I-αP可對角化[16]，則采用GMRES(m)-MSI 方法后的殘差rm滿足不等式

證明由于I-αP可對角化，根據(jù)引理2 有

‖rm‖2≤κ2(X)ε(m)‖r*0‖2.

由引理3 可知，Google 矩陣A的特征值λi=απi(i=2,3,···,n)，矩陣(I-βjP)-1(j=1,2)的特征值為(i=1,2,···,n).由引理4 可知，矩陣A-β1P的特征值為λi-β1πi(i=1,2,···,n).

根據(jù)引理4 的證明，有

(1-α)veTd1=1-αdi=0(i=2,3,···,n)

那么矩陣的特征值為和.

由d1=1-α,π1=1,λ1=1，得?1=1.

對于i=2,3,···,n，有

又因為|πi|≤1(i=2,3,···,n)，有|λi|≤α(i=2,3,···,n)，從而

假設(shè)由GMRES(m)方法計算得到的近似解x1可以表示為其中 γi不全為零，μi(i=1,2,···,n)是矩陣A的一組特征基，且‖μi‖2=1.

由于

及Gμ1=μ1，Glμ1=μ1，有

4 數(shù)值實驗

本節(jié)我們將通過數(shù)值實驗驗證GMRES(m)-MSI 方法的有效性.所有的數(shù)值實驗都是在內(nèi)存8 GB，主頻2.6 GHZ，操作系統(tǒng)Windows 10，處理器為Intel(R) Core(TM)i5-4210M 的計算機(jī)上用MATLAB(2019a)完成的.

我們選取3 個網(wǎng)絡(luò)矩陣進(jìn)行數(shù)值實驗，分別為zkyG 矩陣、Email-Enron 矩陣和web-Stanford 矩陣.zkyG 矩陣記錄了2014 年2 月與中國科學(xué)院主頁http://www.cas.cn 相關(guān)的5 000 個網(wǎng)站的鄰接矩陣，Email-Enron 矩陣記錄了2003 年與美國安然公司有電子郵件來往的36 692 個郵件所構(gòu)成的鄰接矩陣，web-Stanford 矩陣記錄了2002 年與斯坦福大學(xué)主頁https://www.stanford.edu 相關(guān)的281 903 個網(wǎng)站的鄰接矩陣.以上網(wǎng)絡(luò)矩陣都可以從https://sparse.tamu.edu 上獲得，它們的性質(zhì)見表1.表中dimension 代表矩陣的維數(shù)，nonzeros 代表矩陣的非零元個數(shù)，average nonzeros 代表矩陣平均每行非零元的個數(shù).

表1 網(wǎng)絡(luò)矩陣的性質(zhì)Table 1 The nature of the network matrix

當(dāng)阻尼因子α取不同的值時，我們分別比較了內(nèi)外迭代IO 方法、MSI 方法和GMRES(m)-MSI 方法(GMSI)的計算效果.在IO 方法中，選取β=0.5，η=10-2.為了保證數(shù)值實驗的公平，在MSI 方法和GMRES(m)-MSI 方法中同樣取η=10-2.另外，根據(jù)經(jīng)驗取值，在MSI 方法中分別取β1=0.5,β2=0.85(zkyG 矩陣、web-Stanford 矩陣)，以及β1=β2=0.5(Email-Enron 矩陣)，算法會有較好的表現(xiàn).在GMRES(m)-MSI 方法中，選取同樣的 β1和β2，重啟數(shù)m=8，?=10.在這三種方法中，選取的初始向量x0=e/n，迭代終止指標(biāo)τ=10-8.

為了測試GMRES(m)-MSI 方法的加速效果，我們定義相對加速比

其中TMSI代表MSI 方法的運行時間，TGMSI代表GMRES(m)-MSI 方法的運行時間.

表2～4 分別給出了 α取不同值時，這3 種方法在不同測試矩陣上的迭代步數(shù)N和CPU 計算時間T（以秒為單位）.從數(shù)值結(jié)果可見：GMRES(m)-MSI 方法優(yōu)于MSI 方法和IO 方法，并且隨著 α的不斷增大，GMRES(m)-MSI 方法的加速效果越來越明顯.

表2 ZkyG 矩陣Table 2 The zkyG matrix

表3 Email-Enron 矩陣Table 3 The Email-Enron matrix

表4 Web-Stanford 矩陣Table 4 The web-Stanford matrix

圖3～5 分別給出了α取不同值時，3 種方法計算不同測試矩陣的收斂曲線.橫軸N代表迭代步數(shù)，縱軸r代表殘差范數(shù)，即r=‖αz+(1-α)v-x‖2.可以清楚看到，GMRES(m)-MSI 方法的收斂速度是最快的.

圖3 α取不同值時, 3 種方法計算zkyG 矩陣的收斂曲線圖Fig.3 For different α values, the convergence curves of the zkyG matrix with the 3 methods

圖4 α取不同值時, 3 種方法計算Email-Enron 矩陣的收斂曲線圖Fig.4 For different α values, the convergence curves of the Email-Enron matrix with the 3 methods

圖5 α取不同值時, 3 種方法計算web-Stanford 矩陣的收斂曲線圖Fig.5 For different α values, the convergence curves of the web-Stanford matrix with the 3 methods

5 總 結(jié)

本文給出了一種求解PageRank 問題的新方法，即重啟GMRES 方法修正的多分裂迭代法，GMRES(m)-MSI 方法.我們不僅給出了GMRES(m)-MSI 方法的收斂性證明，而且通過數(shù)值實驗驗證了新方法的有效性.數(shù)值結(jié)果表明，當(dāng)阻尼因子α接近于1 時，GMRES(m)-MSI 方法的收斂速度遠(yuǎn)遠(yuǎn)優(yōu)于IO 法和MSI 法.下一步，我們將從理論上研究GMRES(m)-MSI 方法中最優(yōu)參數(shù)的選取原則.