Acceleration of an adaptive generalized Arnoldi method for computing PageRank

1.   Introduction

Using the hyperlink structure of web pages, Google's PageRank becomes one of the most successful methods for measuring the importance of each page ^[1]. From the viewpoint of numerical computations, the core of PageRank problems can be regarded as the problem of solving a dominant eigenvector of the Google matrix $A$:

where $x\in\mathbb{R}^{n}$ is the PageRank vector, $\alpha \in (0, 1)$ is a damping factor, $e = [1, 1, \cdots, 1]^{\rm T} \in \mathbb{R}^{n}$, $v = e/n$, $P\in\mathbb{R}^{n\times n}$ is a column-stochastic matrix, see ^[2] for details.

As an iterative method based on matrix-vector products, the power method is widely used for computing PageRank ^[1,3]. However, when the damping factor $\alpha$ is close to 1, the power method suffers from slow convergence such that some accelerated techniques are developed. For example, based on the inner-outer iteration method proposed by Gleich et al. ^[4], Tian et al. ^[5] developed a general inner-outer iteration method for solving PageRank problems. Using the trace of the Google matrix $A$, Tan ^[6] introduced an extrapolation strategy and presented the power method with extrapolation process based on trace (PET) for improving the computation of PageRank problems.

On the other hand, Krylov subspace methods based on the Arnoldi process have been applied to compute PageRank problems. Golub and Greif ^[7] proposed an Arnoldi-type method by using the singular value decomposition (SVD), where the known largest eigenvalue 1 is considered as a shift such that the computation of the largest Ritz value is avoided. Wu and Wei ^[8] developed a Power-Arnoldi algorithm by periodically combining the power method with the thick restarted Arnoldi algorithm ^[9]. Hu et al. ^[10] proposed a variant of the Power-Arnoldi algorithm by employing the PET method.

Recently, the idea of introducing weighted inner products into an Arnoldi process has successfully been applied to many academic fields ^[11,12]. Yin et al. ^[13] proposed an adaptive generalized Arnoldi (GArnoldi) method for computing PageRank by applying a weighted inner product into an Arnoldi-type method. Wen et al. ^[14] developed an adaptive Power-GArnoldi algorithm by making use of the power method and the adaptive GArnoldi method together. Motivated by these works, with the aim of accelerating the adaptive GArnoldi method, a new method is proposed by periodically knitting the PET method with the adaptive GArnoldi method for PageRank problems. The new method is denoted as GArnoldi-PET method. Convergence performance of our proposed method is studied in detail, and numerical results are used to show its feasibility and effectiveness.

The remainder of this paper is organized as follows. In Section 2, we briefly introduce the PET method and the adaptive GArnoldi method for PageRank problems. In Section 3, we propose the GArnoldi-PET method and discuss its convergence. In Section 4, numerical results and comparisons are reported. Finally, conclusions are given in Section 5.

2.   Previous work

In this section, we give simple introductions of the PET method and the adaptive GArnoldi method for computing PageRank.

2.1.   The PET method for computing PageRank

Here, we first give the algorithmic version of the PET method for PageRank problems as follows, see ^[6] for more details.

Now, some illustrations of Algorithm 1 are given as follows.

● In step 1, the parameter $\mu$ is the trace of the Google matrix $A$.

● In step 2, the power method is run $m_1$ steps, which means the extrapolation technique is not employed to the power method in each iteration, but is used every $m_1$ power iterations.

● In step 3, we can see that the extrapolation strategy based on $x^{(m_1 - 1)}, x^{(m_1)}$ and $\mu$ is easy to implement as given in step 3.1.

2.2.   The adaptive GArnoldi method for computing PageRank

As described in ^[13], let $G = (g_{ij})$ be an $n\times n$ symmetric positive definite (SPD) matrix, then the GArnoldi process based on a weighted inner product is presented as Algorithm 2.

In Algorithm 2, $(\cdot, \cdot)_G$ is a $G$-inner product defined as $(x, y)_G = x^{\rm T}Gy, \forall x\in \mathbb{R}^n, y\in \mathbb{R}^n$, and $\|\cdot\|_G$ is a $G$-norm defined by

where $Q\in \mathbb{R}^{n\times n}$ is an orthogonal matrix, $D = \text{diag}\{d_1, d_2, \cdots, d_n\}$ is a $n\times n$ diagonal matrix with $d_i > 0$, $i = 1, 2, \cdots, n$, and $G = Q^{\rm T}DQ$ is a diagonalized decomposition of $G$. Let $e_m \in \mathbb{R} ^{m}$ be the $m$-th co-ordinate vector, then the GArnoldi process has the following relations ^[13]

where the matrix $V_k = [\widetilde{v}_1, \widetilde{v}_2, \cdots, \widetilde{v}_k]\ (k = m, m+1)$ is an $n\times k$ $G$-orthogonal matrix, $H_m = (h_{ij})$ is an $m\times m$ Hessenberg matrix.

From Algorithm 2, it is obvious that different SPD matrices $G$ will lead to different GArnoldi methods. Since every SPD matrix can be diagonalized, for simplicity, we let $G = \text{diag}\{d_1, d_2, \cdots, d_n\}$, $d_i > 0$, $i = 1, 2, \cdots, n$. It is seen that in each outer iteration of the GArnoldi method, we hope to find a vector $v$ satisfying $\text{min}\|Av - v\|_G$, where $v$ is taken from a Krylov subspace $\mathcal {K}_m(A, v_1) = \text{span}(v_1, Av_1, \cdots, A^{m-1}v_1)$. Denote $r = Av - v \triangleq [r_1, r_2, \cdots, r_n]^{\rm T}$, it has $\text{min}\|Av - v\|_G = \text{min}\sqrt{\sum^n_{i = 1}d_i r_i^2}$, which leads to a weighted least squares problem where $d_i$ is actually the weight for the $i$-th component of residual $r_i$, $i = 1, 2, \cdots, n$. In order to speed up the computation of PageRank problems, Yin et al. ^[13] changed the weights adaptively according to the changing of the current residual corresponding to the approximate PageRank vector. Therefore, one choice of the matrix $G$ is that

where $r$ is the residual vector computed by the last calculation, and $\sum_{i = 1}^n d_i = 1$. And the algorithmic version of the adaptive GArnoldi algorithm for computing PageRank is presented as Algorithm 3.

Note that, in step 5 of Algorithm 3, $s_m$ and $u_m$ denote the right and left singular vector of $H_{m+1, m} - [I; 0]^{\rm T}$ associated with the minimal singular value $\sigma_m$, respectively.

3.   The GArnoldi-PET method for computing PageRank

In order to accelerate the computation of PageRank problems, we develop a new method by combining the PET method with the adaptive GArnoldi method. The new method is called GArnoldi-PET method. Here we first describe the construction of the GArnoldi-PET method, and then discuss its convergence.

3.1.   The GArnoldi-PET method

As described in the subsection 2.1, based on the trace of the Google matrix $A$, an extrapolation strategy has been presented to speed up the convergence of the power method. Numerical experiments in ^[6] have illustrated that the PET method has a faster convergence than the power method when the damping factor $\alpha$ is close to 1. On the other hand, since the Arnoldi method is more computationally intense than applying the same number of iterations of the power method ^[8], thus it is natural to consider using the extrapolation strategy based on trace and the adaptive GArnoldi method together.

Similar to the construction of the Power-Arnoldi algorithm ^[8], the mechanism of our proposed method can be presented as follows: Given a unit positive vector $x^{(0)}$, and an approximate PageRank vector is obtained by iterating the adaptive GArnoldi method (Algorithm 3) for a few times (e.g., 2–3 times). If this approximate PageRank vector does not satisfy our prescribed tolerance, then we run the PET method to obtain another approximate vector with the resulting vector as the initial guess. If this approximate PageRank vector still can not satisfy our accuracy, then we return to Algorithm 3 with the new approximation as the starting vector. Repeating the above procedure until the described accuracy is achieved.

There is a problem about how to control the conversion between the PET method and the adaptive GArnoldi method. Many strategies have been developed to deal with this problem. Here, as given in ^[8], three parameters $\beta$, $restart$, $maxit$ are used to control the procedure. Let $\tau^{(curr)}$ be the residual norm of the current iteration, and $\tau^{(prev)}$ be the residual norm of the previous iteration. Computing $ratio = \tau^{(curr)}/\tau^{(prev)}$, if $ratio > \beta$, then $restart = restart + 1$. If $restart \geq maxit$, then we terminate the PET method and trigger the adaptive GArnoldi method. The specific implementation of the GArnoldi-PET method is given as follows.

Now, some remarks about the GArnoldi-PET method are given as follows.

● As shown in the step 3.16, the matrix $G$ is adaptively changed according to the current residual.

● According to the construction of the GArnoldi-PET method, it is natural to treat the PET method as an accelerated technique for the adaptive GArnoldi method.

● In each iteration of the GArnoldi-PET method, the storage requirements are approximately $m+1$ length$-n$ vectors in the adaptive GArnoldi method and two vectors in the PET method. Its main computational cost consists of $m$ matrix-vector products, $\frac{m(m+1)}{2}$ inner products in the adaptive GArnoldi method and one matrix-vector product in the PET method.

3.2.   Convergence analysis of the GArnoldi-PET method

Here we discuss the convergence analysis of the GArnoldi-PET algorithm. Particularly, we focus on the procedure when turning from the PET method to the adaptive GArnodli method.

Assume $\sigma(A)$ denote the set of eigenvalues of the Google matrix $A$, and its eigenvalues are arranged as $1 = |\lambda _1| > \left| {{\lambda _2}} \right| \ge \cdots \ge \left| {{\lambda _n}} \right|$. Let $\mathcal {L}_{m-1}$ be the set of polynomials of degree not exceeding $m-1$, and $\mathcal {K}_m(A, v_1)$ be a Krylov subspace. If $(\lambda_i, \varphi_i), i = 1, 2, \cdots, n$ are the eigenpairs of $A$, and $(\widetilde{\lambda}_j, \widetilde{y}_j), j = 1, 2, \cdots, m$ are the eigenpairs of $H_m$, then $\widetilde{\lambda}_j$ is often used to approximate $\lambda_j$, and $\widetilde{\varphi}_j = V_m \widetilde{y}_j$ is applied to approximate $\varphi_j$ in the standard Arnoldi method. However, instead of using Ritz vectors $\widetilde{\varphi}_j$ as approximate eigenvectors, Jia ^[15] proposed a new strategy such that for each $\widetilde{\lambda}_j$, a unit norm vector $\widetilde{u}_j \in \mathcal {K}_m(A, v_1)$ satisfying the condition

is used to approximate $\varphi_j$. Here $\widetilde{u}_j$ is called a refined approximate eigenvector corresponding to $\lambda_j$. The convergence of the refined Arnoldi method is given as follows.

Theorem 1 ^[15]. Assume that $v_1 = \sum_{i = 1}^n \gamma _ix_i$ with respect to the eigenbasis $\{x_i \}_{i = 1, 2, \cdots, n}$ in which $\|x_i\|_2 = 1, i = 1, 2, \cdots, n$ and $\gamma_i \ne 0$, let $S = [x_1, x_2, \cdots, x_n]$, and

Then

where $\widetilde{u}_j$ is a refined approximate eigenvector as above, $\sigma_{max}(S)$ and $\sigma_{min}(S)$ are the largest and smallest singular value of the matrix $S$, respectively.

Before analyzing the convergence of the GArnoldi-PET algorithm, some useful conclusions are introduced as follows.

Lemma 1 ^[14]. Let $G = \text{diag}\{d_1, d_2, \cdots, d_n\}$, $d_i > 0, 1\leq i\leq n$, be a diagonal matrix. Then for any vector $x\in\mathbb{R}^{n}$, according to the definitions of the $G$-norm and the 2-norm, it has

Lemma 2 ^[10]. Let $v_1$ be the initial vector for the PET method, which is from the previous adaptive GArnoldi method. Then the PET iteration in Algorithm 4 produces the vector

where $k \geq maxit$, $\eta$ is the normalizing factor, $\mu$ is the trace of the matrix $A$, $m_1$ is a given number, $T$ is called as the iterative matrix and $I$ is an $n\times n$ identity matrix.

Theorem 2 ^[16]. Assume that the spectrum of the column-stochastic matrix $P$ is $\{1, \lambda_2, \cdots, \lambda_n\}$, then the spectrum of the matrix $A = \alpha P + (1-\alpha)ve^{\rm T}$ is $\{1, \alpha\lambda_2, \cdots, \alpha\lambda_n\}$, where $0 < \alpha < 1$, $v$ is a vector with nonnegative elements such that ${e^{\rm T}}v = 1$.

Theorem 3 ^[17]. Let $P$ be an $n \times n$ column-stochastic matrix. Let $\alpha$ be a real number such that $0 < \alpha < 1$. $E$ is the $n \times n$ rank-one column-stochastic matrix $E = v{e^{\rm T}}$, where $e$ is the $n$-vector of all ones and $v$ is an $n$-vector whose elements are all non-negative and sum to 1, $A = \alpha P + (1-\alpha)ve^{\rm T}$ is the $n \times n$ column-stochastic matrix, then its dominant eigenvalue ${\lambda _1} = 1$, $\left| {{\lambda _2}} \right| \le \alpha$.

In the next cycle of the GArnoldi-PET method, $v_1^{new}$ will be used as the initial vector for an $m$-step GArnoldi process, so that the new Krylov subspace

is constructed. The following theorem shows the convergence of the GArnoldi-PET method.

Theorem 4. Assume that $v_1 = \sum_{i = 1}^n \gamma _ix_i$ with respect to the eigenbasis $\{x_i\}_{i = 1, 2, \cdots, n}$ in which $\|x_i\|_2 = 1, i = 1, 2, \cdots, n$ and $\gamma_1 \ne 0$. Let $G = \text{diag}\{d_1, d_2, \cdots, d_n\}$, $d_i > 0, i = 1, 2, 3, \cdots, n$, $S = [x_1, x_2, x_3, \cdots, x_n]$, and

Then

where $u$ is taken from the Krylov subspace $\mathcal {K}_m(A, v_1)$, $\mu = 1+\alpha \left(\frac{l}{n}-1\right)$ with the number of dangling nodes $l$, $k \geq maxit$, $\sigma_{min}(S)$ is the smallest singular value of the matrix $S$.

Proof. Let $\{1, \pi_2, \cdots, \pi_n\}$ be the eigenvalue set of the matrix $P$, then $\{1, \lambda_2 = \alpha\pi_2, \cdots, \alpha_n = \alpha\pi_n\}$ is the eigenvalue set of the matrix $A$ by Theorem 2. According to (3.3), we have

Let $\varphi_i\; (i = 1, 2, \cdots, n)$ be eigenvalues of the iterative matrix $T$. Then, we have

Since $\mu = 1+\alpha \left(\frac{l}{n}-1\right)$ with the number of dangling nodes $l$, it has $1-\mu = \alpha \left(1-\frac{l}{n}\right)\geq 0$. On the other hand, according to Theorem 3 and $1 = |\lambda _1| > \left| {{\lambda _2}} \right| \ge \cdots \ge \left| {{\lambda _n}} \right|$, it has $|\lambda_i| \leq \alpha$ and

Hence, we obtain $|\lambda_i-\mu +1|\leq \alpha-\mu+1, i = 2, \cdots, n$, and

For any $u \in \mathcal {K}_m(A, v_1^{new})$, there exists $q(x) \in \mathcal {L}_{m-1}$ such that

where we used the facts that $\lambda_1 = 1$, $Ax_i = \lambda_ix_i, Tx_i = \varphi_ix_i, i = 1, 2, \cdots, n$. According to (3.2) and (3.5), for the numerator of (3.6), it has

For the denominator of (3.6), it has

Combining (3.5), (3.7) and (3.8) into (3.6), we have

where $\frac{\alpha^{{m_{\;1\;}} - 1} \;(\alpha-\mu + 1)}{2-\mu} < 1$. Let $p(\lambda) = q(\lambda)/q(1)$, where $p(1) = 1$, then we get

4.   Numerical experiments

In this section, we test the effectiveness of the GArnoldi-PET method and compare it with the PET method (Algorithm 1) ^[6], the Power-Arnoldi algorithm (called as PA) ^[8] and the adaptive GArnoldi method (called as A-Arnoldi) ^[13] in terms of the number of matrix-vector products (Mv) and the computing time in seconds (CPU). All the numerical results are obtained by using MATLAB 2018b on the Windows 10 (64 bit) operating system with 1.7 GHz Intel(R) Core(TM) i5 CPU and RAM 4.00 GB.

In Table 1, we list the characteristics of test matrices including the matrix size ($n$), the number of nonzero elements ($nnz$), the number of dangling nodes ($numd$) and the density ($den$) which is defined by $den = \frac{nnz}{n\times n} \times 100$.

For a fair comparison, all algorithms use the same initial guess $x^{(0)} = v = e/n$ with $e = [1, 1, \cdots, 1]^{\rm T}$. The tolerance is chosen as $tol = 10^{-8}$. The values of the damping factor $\alpha$ are 0.990, 0.993, 0.995 and 0.997, respectively. The parameter $\beta = \alpha - 0.1$. We run the thick restarted Arnoldi procedure, with the number of approximate eigenpairs $p$, two times per cycle in the Power-Arnoldi method. Similarly, we run the adaptive GArnoldi procedure two times per cycle in the GArnoldi-PET method. In addition, for describing the efficiency of the GArnoldi-PET method, we define

Example 1. The first test matrix is the wb-cs-stanford matrix, which contains 9914 pages, 36854 links and 2861 dangling nodes. It is available from https://sparse.tamu.edu/Gleich/wb-cs-stanford. In this example, we set the parameters $m = 5, p = 3, maxit = 6$ and $m_1 = 40$. Numerical results of the PET method, the adaptive GArnoldi method, the Power-Arnoldi algorithm and the GArnoldi-PET algorithm are reported in Table 2. Figure 1 plots the convergence history of the four methods with different values of $\alpha$.

From Table 2, we can see that the Power-Arnoldi algorithm works better than the PET method and the adaptive GArnoldi method in terms of the number of matrix-vector products and the computing time. However, the GArnoldi-PET algorithm performs the best. For example, when $\alpha = 0.997$, the Power-Arnoldi algorithm needs 0.1473 seconds to reach the desired accuracy, while the GArnoldi-PET algorithm only uses 0.1038 seconds, and the speedup is $29.53\%$.

From Figure 1, it is easy to find that the GArnoldi-PET algorithm has a faster convergence speed than the PET method and Power-Arnoldi algorithm, even though its iteration counts are slightly inferior to the adaptive GArnoldi method. Obviously, only the number of iterations can not describe the whole story.

Example 2. The second test matrix is the web-Stanford matrix, which contains 281903 nodes, 2312497 links and 172 dangling nodes. It is available from https://sparse.tamu.edu/SNAP/web-Stanford. In this example, we choose the parameters $m = 5, p = 3, maxit = 12$ and $m_1 = 35$. Numerical results of the PET method, the adaptive GArnoldi method, the Power-Arnoldi algorithm and the GArnoldi-PET algorithm are given in Table 3. Figure 2 depicts the convergence of the four methods with different values of $\alpha$

From Table 3, it observes that the GArnoldi-PET algorithm outperforms the other three methods in terms of the number of matrix-vector products and the computing time. Although the speedup is only $6.57\%$ relative to the Power-Arnoldi algorithm when $\alpha = 0.993$. However, when $\alpha$ increases, e.g., $\alpha = 0.997$, the Power-Arnoldi algorithm needs 13.7626 seconds to reach the desired accuracy, the GArnoldi-PET algorithm only takes 10.7224 seconds, and the speedup becomes $22.09\%$.

From Figure 2, it shows that the GArnoldi-PET algorithm converges faster than the PET method and the Power-Arnoldi algorithm. When $\alpha$ is close to one, e.g., $\alpha = 0.995$ and $\alpha = 0.997$, the iteration counts of the GArnoldi-PET algorithm are less than those of the adaptive GArnoldi method. This suggests that our new algorithm has some potential.

Example 3. The last test matrix is the wikipedia-20051105 matrix, which contains 1634989 nodes, 19753078 links and 72556 dangling nodes. It is available from https://sparse.tamu.edu/Gleich/wikipedia-20051105. In this example, we make the parameters $m = 5, p = 3, maxit = 8$ and $m_1 = 50$. Numerical results of the PET method, the adaptive GArnoldi method, the Power-Arnoldi algorithm and the GArnoldi-PET algorithm are listed in Table 4. Figure 3 shows the convergence curves of the four methods with different values of $\alpha$.

From Table 4, we also see that the GArnoldi-PET algorithm makes great improvements on the PET method, the adaptive GArnoldi method and the Power-Arnoldi algorithm in terms of the number of matrix-vector products and the computing time. For a large damping factor such as $\alpha = 0.997$, the Power-Arnoldi algorithm takes 56.5771 seconds to reach the desired accuracy, while the GArnoldi-PET algorithm takes 41.9286 seconds to achieve the same accuracy, and the speedup is $25.89\%$.

From Figure 3, we again find that the GArnoldi-PET algorithm converges faster than the PET method, the adaptive GArnoldi method and the Power-Arnoldi algorithm. For different values of $\alpha$, the iteration counts of the GArnoldi-PET algorithm are the least.

5.   Conclusions

In this paper, by combining the PET method with the adaptive GArnoldi method, we propose a new method called as GArnoldi-PET method for accelerating the computation of PageRank problems. Its construction and theoretical analysis can be found in Section 3. Numerical results in Section 4 show that our proposed method is quite efficient and better than the existing methods, especially when the damping factor is close to 1. However, much research still needs further study. For example, determining the optimal choice of the parameters, or considering to use some preconditioning strategies as given in ^[18,19,20] for the GArnoldi method.

Acknowledgments

This research is supported by Science Challenge Project (TZ2016002-TZZT2019-B1.4), the Key Fund Project of Sichuan Provincial Department of Education (17za0003) and the Applied Basic Research Project of Sichuan Province (2020YJ0007). We would like to thank the anonymous referees for their valuable comments and suggestions.

Conflict of interest

All authors declare that there is no conflict of interest in this paper.