## A Hessenberg-type Algorithm for Computing PageRank Problems

Xian-Ming Gu , 1 Siu-Long Lei,* , 2 Ke Zhang , 3 Zhao-Li Shen , 4 Chun Wen , 5 Bruno Carpentieri 6

Abstract: PageRank is a widespread model for analysing the relative relevance of nodes within large graphs arising in several applications. In the current paper, we present a cost-effective Hessenberg-type method built upon the Hessenberg process for the solution of difficult PageRank problems. The new method is very competitive with other popular algorithms in this field, such as Arnoldi-type methods, especially when the damping factor is close to 1 and the dimension of the search subspace is large. The convergence and the complexity of the proposed algorithm are investigated. Numerical experiments are reported to show the efficiency of the new solver for practical PageRank computations.

Keywords: PageRank vector, Hessenberg process, Arnoldi process, Krylov subspace method, Ritz values.

AMS classifications: 65F15, 65F10, 65Y20.

## 1 Introduction

The PageRank model was originally introduced by S. Brin and L. Page in 1999 [1] to develop fast web search engines, and then studied and enhanced in a vast number of research papers (see, e.g., [2-7]). The model provides a powerful network centrality measure to identify the most important nodes within large graphs arising in several applications fields, such as in chemistry, bioinformatics, neuroscience, bibliometrics, and others [8]. In the original Web problem, the PageRank algorithm determines the ranking of each Web page by computing the stationary probability vector of a random walking process on the Web link graph, which is a directed graph representing the linking structure of the Web [7,9]. The Web link graph is a binary matrix G ∈ N n × n , where n denotes the number of pages, such that G i, j ( ) = 1 when page j has a link pointing to page i , and G i, j ( ) = 0 otherwise. From a linear algebra viewpoint, the algorithm finds the vector x that satisfies

$$A x = x, \ \| x \| _ { 1 } = 1, \ x > 0,$$

1 School of Economic Mathematics/Institute of Mathematics, Southwestern University of Finance and Economics, Chengdu, Sichuan 611130, P.R. China and Bernoulli Institute for Mathematics, Computer Science and Artificial Intelligence, University of Groningen, Nijenborgh 9, P.O. Box 407, 9700 AK Groningen, The Netherlands. E-mail: guxianming@live.cn, guxm@swufe.edu.cn .

2 Corresponding author: Department of Mathematics, Faculty of Science and Technology, University of Macau, Avenida da Universidade, Macao, P.R. China. E-mail: sllei@um.edu.mo .

3 College of Arts and Sciences, Shanghai Maritime University, Shanghai 201306, P.R. China. E-mail: xznuzk123@126.com . 4 Department of Applied Mathematics, College of Science, Sichuan Agricultural University, Yaan, Sichuan 625014, P.R. China. E-mail: szlxiaoyao@163.com .

5 School of Mathematical Sciences, University of Electronic Science and Technology of China, Chengdu, Sichuan 611731, P.R. China. E-mail: wchun17@163.com .

6 Facolt` di Scienze e Tecnologie informatiche, Libera Universit` di Bolzano, Dominikanerplatz 3 - piazza Domenicani, a a 3 Italy - 39100, Bozen-Bolzano. E-mail: Bruno.Carpentieri@unibz.it .

that is it finds the principal unit positive eigenvector x [9] of the Google matrix

$$A = \alpha ( P + v d ^ { \top } ) + ( 1 - \alpha ) v e ^ { \top }.$$

In Eq. (1.2), matrix P ∈ R n × n is called the transition matrix with respect to the random walking process, and is defined as

 The damping factor 0 &lt; α &lt; 1 defines the probability that random Web surfers choose a random link from the page they are visiting [25]. The teleporting vector v = [ v , 1 · · · , v n ] /latticetop ∈ R n × 1 ( v ≥ 0 and ‖ v ‖ 1 = 1) defines the probability v i that the Web surfer jumps to an external page i . Finally, d ∈ N n × 1 is such that d ( ) = 1 if page i i has no hyperlink and 0 otherwise, and e = [1 1 , , . . . , 1] /latticetop ∈ N n × 1 .

$$P ( i, j ) = \begin{cases} \frac { 1 } { \sum _ { k = 1 } ^ { n } G ( k, j ) }, & \text{if $G(i,j)=1$}, \\ 0, & \text{otherwise}. \end{cases} \\ ) < \alpha < 1 \text{ defines the probability that random Web surfers choose a random link from} \\ \text{int} \in \mathbb { R } \ \text{The telernot} \ n = \mathbb { R } ^ { n \times 1 } ( \mathbb { r } > 0 \, \text{and} \ \| \mathbb { R } \ \|. = 1 \}$$

The value of the damping factor α plays an important role in the PageRank model. Theoretically, it represents an upper bound 0 &lt; λ | 2 | ≤ α &lt; 1 for the second largest eigenvalue, λ 2 , of A . Further properties of the Google matrix can be found in [2,4-7,10]. For low values of α (e.g., α = 0 85), . λ 2 is well separated from the largest eigenvalue of A , which is λ 1 = 1, ensuring rapid convergence of the power algorithm applied to problem (1.1). On the other hand, convergence tends to slow down noticeably when α is very close to 1, requiring more robust algorithms than the simple power method. Some computational approaches proposed in the literature include Monte Carlo methods [11], adaptive algorithms [3, 12], extrapolation techniques [2, 7, 13], singular value decompositions [19-21] reordering [14, 15] and innerouter solution strategies [16].

A significant amount of work has been devoted in the last years to the use of Krylov subspace methods based on the Arnoldi decomposition [17,18] for large PageRank computations, mainly due to their memory efficiency and attractive inherent parallelism. Golub and Greif extended the refined Arnoldi procedure to PageRank by forcing a relevant shift to be 1, being able to circumvent the drawbacks due to complex arithmetic and showing overall very good algorithmic efficiency [22]. Many techniques have attempted to combine the conventional Arnoldi method and the power algorithm to produce faster solvers, e.g., the Power-Arnoldi [23-25], the Arnoldi-Extrapolation [26] and the Arnoldi-Inout [27] methods. In the technique proposed in [28], the weighted least squares problem is changed adaptively according to the component of the residual. Then, the generalized Arnoldi method is used to compute the approximate PageRank vector. However, when the dimension of the Krylov subspace is large, Arnoldi-based solvers tend to become very expensive in terms of memory and computational costs; on the other hand, if the dimension of the Krylov subspace is low, they sometimes fail to accelerate the basic power method, especially when the damping factor is high [22-25]. Similarly to the restarted GMRES algorithm [31], they may stagnate in many circumstances [32].

Motivated by costs considerations, other work developed PageRank solvers based on the Bi-Lanczos orthogonalization procedure [18, pp. 139-145] (see e.g. [29, 30]) instead of Arnoldi. In this paper, we look in particular at the Hessenberg reduction process [33-36] that was introduced by K. Hessenberg in 1940 [33], and revived recently to establish a number of cost-effective Krylov subspace solvers for sparse matrix systems, due to its lower arithmetic and storage requirements. The method has been extended to compute the characteristic polynomial of matrices [33, 34, 37], to solve general nonsymmetric linear systems [17, 35, 36, 38, 39], including systems with multiple right-hand sides [40-44] and multi-shifted coefficient matrices [45-47], other types of matrix equations [42, 48, 49], the action of a matrix function

f ( A ) v [47] and other related problems [51]. Theoretical and numerical studies have investigated the mathematical properties of the Hessenberg process, especially in relation with the more conventional Arnoldi procedure. The Arnoldi method was first introduced in 1951 as a means of reducing a dense matrix A into a Hessenberg form by unitary transformations, whereas the Hessenberg process applies similarity transformations [50] and is more suitable for parallel computing. In his paper, Arnoldi hinted that the eigenvalues of the Hessenberg matrix obtained after k /lessmuch n steps, where n is the size of A , could provide accurate approximations of some eigenvalue of A . It was later discovered that this strategy can lead to efficient techniques for approximating eigenvalues of large sparse matrices. In the current work, following a similar development, we modify the Hessenberg process to establish a new family of eigenvalue solvers. We combine the new solvers with the refined and explicitly restarted techniques introduced in [19, 22] to compute realistic PageRank problems. Finally, we analyse their convergence behavior and computational complexity.

The rest of this paper is organized as follows. In Section 2, the Hessenberg process is introduced and a novel family of eigenvalue solvers based on this procedure is described. Moreover, theoretical aspects of such eigenvalue solvers are highlighted in comparison with the classical Arnoldi-like methods. In Section 3, we derive the Hessenberg-type method with explicit restarting and refined techniques for computing PageRank. Both the convergence behavior and the computational cost of the proposed method are discussed. Numerical results in Section 4 show the effectiveness of the proposed algorithm, also against other popular PageRank algorithms. In Section 5, we present some conclusions arising from our study.

## 2 The Hessenberg process with applications to eigenvalue computations

In this section, we briefly review the Hessenberg procedure that is at the basis of our development. We recall some fundamental properties of the algorithm and then we describe a Hessenberg-based projection technique for computing eigenvalues of large nonsymmetric matrices. Our theoretical analyses demonstrates the feasibility of the method, showing some computational advantages over the more conventional Arnoldi procedure.

## 2.1 The Hessenberg process

The Hessenberg process is an an oblique projection technique that reduces a given nonsymmetric matrix A ∈ R n × n to a Hessenberg form [34, 45, pp. 377-381]. Originally, the method was described as a way to compute the characteristic polynomial of a matrix [33]. The basic procedure is presented in Algorithm 1, where a pivoting strategy is included to ensure numerical stability.

Let L k = [ l 1 , . . . , l m ] denote a matrix, ¯ H m = [ h i,j ] be an upper Hessenberg matrix and H m the submatrix obtained from ¯ H m by deleting its last row. Finally, we denote by P /latticetop k = [ e p 1 , e p 2 , . . . , e p n ] where the scalars p i 's (for i = 1 , . . . , n ) are defined in Algorithm 1. After k steps of Algorithm 1, the following matrix equation can be easily established,

$$\ A L _ { k } & = L _ { k + 1 } \bar { H } _ { k } \\ & = L _ { k } H _ { k } + h _ { k + 1, k } l _ { k + 1 } e _ { k } ^ { \top },$$

and P k L k is lower trapezoidal [35, 39]. Unlike Arnoldi, the Hessenberg procedure with pivoting is not

## Algorithm 1 The Hessenberg procedure with pivoting strategy

Input: a general (nonsymmetric) square matrix A ∈ R n × n , an initial vector q 0 and an integer m Output: L k = [ l 1 , . . . , l m ] ∈ R n × m and H m ∈ R ( m +1) × m

/negationslash

Happy breakdown

- 1: Initialize the permutation vector p = [1 2 , , . . . , n ] /latticetop and determine i 0 such that | ( v ) i 0 | = ‖ v ‖ ∞ 2: Compute β = ( v ) i 0 , then l 1 = v /β and p (1) ↔ p ( i 0 ), where ' ↔ ' is used to swap contents [35]. 3: for j = 1 2 , , . . . , m , do 4: Compute u = A l j 5: for i = 1 2 , , . . . , j , do 6: h i,j = ( u ) p ( ) i 7: u = u -h i,j l i 8: end for 9: if ( j &lt; n and u = 0 ) then 10: Determine i 0 ∈ { j +1 , ..., n } such that | ( u ) p ( i 0 ) | = ‖ ( u ) p ( j +1): p ( n ) ‖ ∞ ; 11: h j +1 ,j = ( u ) p ( i 0 ) ; l j +1 = u /h j +1 ,j ; p ( j +1) ↔ p ( i 0 ) 12: else 13: h j +1 ,j = 0; Stop ✄ 14: end if 15: end for

guaranteed to be backward stable in finite precision arithmetic [50]. However, the backward error is reported to be small for most practical problems [36,39,45]. This is also confirmed by our computational experiences. We did not observe noticeable numerical instabilities due to the non-orthogonality of the Krylov basis in our numerical experiments; see also [45] for a discussion of this topic.

## 2.2 Approximation of eigenpairs based on the Hessenberg process

Methods to approximate eigenpairs of a large nonsymmetric matrix A usually compute them from the Hessenberg decomposition of A given by Eq. (2.4). The upper Hessenberg matrix H m can be seen as the projection of A onto the Krylov subspace

$$\mathcal { K } _ { m } ( A, v ) = \text{span} \{ v, A v, \dots, A ^ { m - 1 } v \},$$

and the columns of the matrix L m are a basis of K m ( A, v ). Under certain conditions, the eigenvalues of H m converge to the eigenvalues of A [18,34,50]. Various eigenvalue solvers are built upon this idea, differing each other mainly on the type of projection technique that is used to derive the decomposition (2.4), for example the Arnoldi process, the Bi-Lanczos procedure and the Induced Dimension Reduction (IDR) strategy [52, 53]. The approximate eigenpairs of A are retrieved in the form ( θ , i x ( ) i = L m y ( ) i ), where ( θ , i y ( ) i ) are eigenpairs of the small dimensional matrix H m , such that

$$H _ { m } y ^ { ( i ) } = \theta _ { i } y ^ { ( i ) } \quad \text{with } \| y ^ { ( i ) } \| = 1, \ \ i = 1, 2, \cdots, m.$$

A bound on the residual error can be established directly from Eq. (2.4), by writing

$$A x ^ { ( i ) } - \theta _ { i } x ^ { ( i ) } & = A L _ { m } y ^ { ( i ) } - \theta _ { i } L _ { m } y ^ { ( i ) } \\ & = h _ { m + 1, m } l _ { m + 1 } e _ { m } ^ { \top } y ^ { ( i ) }.$$

If we denote as [ y ( ) i ] m the m th component of the vector y ( ) i , then we obtain or, if we normalize the vector l m ,

$$\Phi _ { 1 } \Phi _ { 2 } \Phi _ { 3 } \Phi _ { 4 } \Phi _ { 5 } \Phi _ { 6 } \Phi _ { 7 } \Phi _ { 8 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phiconner \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi_ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 }\Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9} \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi }$$

This analysis is in line with the results described in [19].

For the ' west0479 ' problem, a real-valued 479-by-479 sparse matrix that has both real and complex eigenvalues, in Fig. 1 we plot the eigenvalues of A computed by the MATLAB command eig , and those of the Hessenberg matrix produced by the IDR( s = 4) projection technique [52], the Sonneveld pencil [52,53], the Arnoldi and the Hessenberg procedures. We clearly see that the Hessenberg procedure can estimate the exterior Ritz values very well, in some cases even slightly more accurately than the Arnoldi procedure. The condition number of the Krylov basis matrix L m is an effective metric to determine the accuracy of the method used. Fig. 2 illustrates that for the Hessenberg process, this condition number does not grow significantly when the dimension of the Krylov subspace increases. The numerical error of the Hessenberg decomposition often remains small 1 . This observation is also exemplified by the stochastic analysis presented in [52, Section 3.3], that can produce (numerical) evidence that the Hessenberg process can be as efficient as the Arnoldi process for practical eigenvalue computations.

$$\| A x ^ { ( i ) } - \theta _ { i } x ^ { ( i ) } \| \leq | h _ { m + 1, m } | \left | [ y ^ { ( i ) } ] _ { m } \right |. & & ( 2. 8 ) \\ \text{e with the results described in } [ 19 ]. & \\'9'" problem, a real-valued 479-by-479 sparse matrix that has both real and complex$$

## 2.3 Relation between the Arnoldi and Hessenberg decompositions

In this subsection, we provide more theoretical background supporting the use of the Hessenberg process to approximate effectively eigenpairs of a given nonsymmetric matrix A . The starting point of our analysis is a comparison between the Hessenberg decompositions computed by the Arnoldi and by the Hessenberg procedures. After m steps of the Arnoldi method applied to A , starting with an initial vector v 0 and assuming no breakdown, the following Hessenberg decomposition is derived:

$$A V _ { m } = V _ { m } H _ { m } + h _ { m + 1, m } \upsilon _ { m + 1 } e _ { m } ^ { \top } = V _ { m + 1 } \bar { H } _ { m }.$$

On the other hand, after m steps of the Hessenberg procedure applied to A with same initial vector l 0 , the resulting matrix factorization writes as

$$A L _ { m } = L _ { m } H _ { m } ^ { ( h ) } + h _ { m + 1, m } ^ { ( h ) } l _ { m + 1 } e _ { m } ^ { \top } = L _ { m + 1 } \hat { H } _ { m } ^ { ( h ) }. & & ( 2. 1 0 )$$

Differently from Arnoldi, however, the columns of L m in Eq. 2.10 are not mutually orthogonal. By computing the reduced QR factorization of L m +1 ,

$$L _ { m + 1 } & = Q _ { m + 1 } R _ { m + 1 }, & ( 2. 1 1 )$$

we can establish the following relation between Eq. (2.10) and Eq. (2.11):

$$A Q _ { m } = Q _ { m + 1 } R _ { m + 1 } \hat { H } _ { m } ^ { ( h ) } R _ { m } ^ { - 1 }.$$

Due to the uniqueness of the Arnoldi decomposition (see Section 3.3 in [52]), by comparing Eq. (2.9) and Eq. (2.12) we conclude that Q m +1 = V m +1 and

$$\bar { H } _ { m } = R _ { m + 1 } \hat { H } _ { m } ^ { ( h ) } R _ { m } ^ { - 1 }.$$

The above result is summarized in the following proposition:

1 Here we give its executable codes in the website: https://github.com/Hsien-Ming-Ku/PageRank-Hessenberg .

Fig. 1: Plots of the Ritz values generated by the IDR( s = 4) factorization, the Sonneveld pencil, the Arnoldi and the Hessenberg procedures.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-0f321c53-a92f-4029-8391-4b8bb96ae138.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-a9291fa8-66c5-4319-ac0c-9a62fd6e133e.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-3ad8917f-a0eb-4021-9f2a-fd6dd68cf289.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-c28e4562-8e56-41fe-ad03-94f299669bb4.jpg)

## Proposition 2.1. It follows that

$$H _ { m } \ = \ R _ { m } H _ { m } ^ { ( h ) } R _ { m } ^ { - 1 } + \frac { h _ { m + 1, m } ^ { ( h ) } } { r _ { m, m } } \tilde { r } e _ { m } ^ { \top },$$

$$\frac { h _ { m + 1, m } ^ { ( h ) } } { r _$$

where ˜ = [ r r i,m +1 ] m i =1 is the vector containing the first m components of the ( m +1) th column of R m +1 .

In fact, Eqs. (2.14)-(2.15) can be also found in [38, 52]. According to Eq. (2.15), in exact arithmetic both procedures produce upper Hessenberg matrices with the same eigenvalues. If the Arnoldi process terminates successfully (i.e., a happy breakdown situation h m +1 ,m = 0 occurs), so does the Hessenberg procedure ( h ( h ) m +1 ,m = 0). On the other hand, a direct consequence of Proposition 2.1 is that the Ritz values produced by the Arnoldi and by the Hessenberg processes are not the same. The condition number of the Krylov basis matrix L m , which is the same as the condition number of matrix R m , gives a clear

Fig. 2: The quality of different basis matrices; Left: the error of different Hessenberg decompositions; Right: the condition number of basis matrices generated by different Hessenberg decompositions.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-ef413832-548c-481d-a3f5-2b93501794c2.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-5ad49a50-a3ae-4bb7-aa6b-e0ea24c3cf7c.jpg)

indication of the accuracy of the eigenvalues of H ( h ) m compared to those resulting from the Arnoldi process. In conclusion, it can be expected that the Hessenberg process can produce feasible approximations of eigenpairs of large nonsymmetric matrices matrix.

## 3 A Hessenberg-type algorithm for computing PageRank

In this section, we will propose a Hessenberg-based algorithm to compute the PageRank vector, that is the positive unit eigenvector corresponding to the largest eigenvalue of the Google matrix. Golub and Greif suggested that the explicitly restarted Arnoldi process for computing eigenvalues and eigenvectors should be implemented in complex arithmetic, thus it is not suitable (it needs to be refined) to compute the PageRank vector [22]. Moreover, since the Hessenberg process is similar in nature to Arnoldi, except that it produces a non-orthogonal basis of the Krylov subspace, we follow a similar development to the refined Arnoldi method for PageRank problems. In other words, the solver described in this section may be called the refined Hessenberg method for PageRank. However, we do not approximate the eigenvectors of A from those of H m (the so-called Ritz-like vectors). Instead, we compute the refined Ritz-like vectors, i.e., the singular vectors associated with the smallest singular values of A -θ I i [19, 22], where { θ i } m i =1 are named the Ritz-like values; cf. Eq. (2.6). The Hessenberg-based method enjoys similar numerical properties to the Arnoldi-based variant: an effective separation of the eigenvectors is ensured, complex arithmetic is avoided by using a shift equal to 1 (due to the fact that the largest eigenvalue of the Google matrix is known), the smallest singular value converges more smoothly to zero than the largest Ritz value to 1 [22, Section 3]. The Hessenberg-type method for computing the PageRank vector is presented in Algorithm 2. The following convergence result can be established after each cycle of m iterations of Algorithm 2.

Theorem 3.1. Let Q m = [ q 1 , q 2 , · · · , q m ] be the matrix obtained from running m -steps of either the Arnoldi or the Hessenberg procedures applied to A starting from an initial vector q 0 , then the Hessenberg

## Algorithm 2 The Hessenberg-type algorithm for computing PageRank

Input: the Google matrix A ∈ R n × n , initial guess q 0 and m

Output: the approximate PageRank vector q

- 1: Choose q 0 and determine i 0 such that | ( q 0 ) i 0 | = ‖ q 0 ‖ ∞ , m is the dimension of search subspace
- 2: Compute β = ( q 0 ) i 0 , then q 0 = q 0 /β
- 3: for /lscript = 1 2 , , . . . , until convergence do
- 4: Compute L m and H m +1 ,m via running Algorithm 1 with the inputs A, m, q 0
- 5: Compute the singular value decomposition: H m +1 ,m -[ I m ; 0 ] = U Γ V /latticetop
- 6: Compute q = L m m v
- 7: Compute r = σ m L m +1 u m
- 8: if ‖ r ‖ 1 / ‖ q ‖ 1 &lt; tol then
- 9: Stop and exit
- 10: else
- 11: Set q 0 = q
- 12: end if
- 13: end for

matrix decomposition can be uniformly written as

$$A Q _ { m } & = Q _ { m } H _ { m } + h _ {$$

Denote as σ m the smallest singular value of H m +1 ,m -[ I m ; 0] , then v m at Line 6 of Algorithm 2 is the corresponding right singular vector, and Q m m v is the approximate PageRank vector. The residual vector at each restarting cycle can be computed as r = σ m Q m +1 u m .

Proof . According to Eq. (3.16) and Algorithm 2, it follows that

$$\begin{array} {$$

where q m is an approximate PageRank vector. Thus, the assertion is verified.

✷

Below, we give the norms of the residual vectors computed by the Arnoldi and by the Hessenberg procedures, respectively:

/negationslash and

$$\| r \| _ { 2 } = \begin{cases} \sigma _ { m }$$

$$\| r \| _ { 1 } = \sigma _ { m } \| Q _ {$$

/negationslash

Although the 2-norm of the Arnoldi residual vectors is much cheaper than the 2-norm of the Hessenberg residual vectors, it should be noted that for PageRank computations it is generally recommended to use the 1-norm; refer, e.g., to [2, 25, 26]. Therefore, the computational complexity of the stopping criterion (at line 8 of Algorithm 2) is almost the same for both methods [24,28,30].

Before we end this section, we provide estimates on the storage requirement and the computational complexity of the new algorithm, also compared against other popular methods. Table 1 shows the

Table 1: Memory requirement for running k iterations of different PageRank algorithms.

| Algorithm   | dim( n )                    | dim( k )   | total                    |
|-------------|-----------------------------|------------|--------------------------|
| Power       | x k , x k - 1               | -          | 2 n                      |
| QE-Power    | k , x k - 1 , x k - 2 , x k | -          | 4 n                      |
| Inner-Outer | x , y , f                   | -          | 3 n                      |
| Arnoldi     | Q k , w                     | H k 1      | ( k +1) n + k 2 / 2 +2 k |
| Hessenberg  | L k , u                     | H k        | ( k +1) n + k 2 / 2 +2 k |
| A-Arnoldi   | Q k , w , r                 | H k        | ( k +2) n + k 2 / 2 +2 k |

1 To minimize the memory requirements, both L k and ¯ H k can be written into the same array as A . Hence, the Hessenberg process requires slightly less storage than the Arnoldi's procedure; refer to [35,39] for discussions on this issue.

memory required in addition to A for running k iterations of the power method (referred to as Power in the table), the power method with quadratic extrapolation (called as QE-power ), the Arnoldi-type method (abbreviated as Arnoldi ), the adaptively accelerated Arnoldi method (called as A-Arnoldi ) and the Hessenberg-type method (abbreviated as Hessenberg ). Here, w x u , , and r are intermediate working vectors used at the k -th step, and Q k denotes the k orthonormal vectors in the modified GramSchmidt process. Analogously, L k denotes the n × k non-orthonormal matrix in the variant of the LU -like factorization process.

Table 2 shows the computational workloads required to execute one cycle of each different iterative algorithm. Here, N z represents the number of nonzero entries of matrix A . In fact, both Arnoldi and A-Arnoldi for computing PageRank enjoy the similar pseudo-code of Algorithm 2, the only difference is to choose the Hessenberg, Arnoldi, or generalized Arnoldi process at Line 4 of Algorithm 2. It implies that we need to compare the cost of the Hessenberg process with both the Arnoldi and generalized Arnoldi procedures. We can see that one cycle of the Hessenberg-type method is cheaper than for Arnoldi and for the generalized Arnoldi methods. Thus its use can be computationally attractive for large PageRank computation. The convergence performance of the Hessenberg method is also superior to Arnoldi algorithms (i.e., Arnoldi and A-Arnoldi ), as proved numerically in the next section. Besides matrix-vector multiplications, also the computation of vector norms and SAXPY 2 (which stands for 'SinglePrecision A X Plus Y' and is a combination of the scalar multiplication and vector addition) operations · determines the total computational cost of these three algorithms. Overall, when m increases, the cost of each cycle increases too but the total number of iterations decreases. The optimal value of the restart parameter that minimizes the total solution time remains problem dependent, and this issue will be examined in our numerical experiments section.

2 See the details from https://developer.nvidia.com/blog/six-ways-saxpy/ .

Table 2: Computational cost of one cycle of different algorithms for computing PageRank.

| operation                                             | Arnoldi     | Hessenberg   | A-Arnoldi        |
|-------------------------------------------------------|-------------|--------------|------------------|
| matrix-vector product (see Line 4 of Algorithm 1)     | 2 mN z      | 2 mN z       | 2 mN z           |
| inner product (e.g., see Line 6 of Arnoldi )          | m ( m +1) n | 0            | 3 m ( m +1) n/ 2 |
| ( u ) p ( i ) (see Line 6 of Algorithm 1)             | 0           | mn           | 0                |
| SAXPY : x + α y (see Line 7 of Algorithm 1)           | m ( m +1) n | m ( m +1) n  | m ( m +1) n      |
| ‖ u ‖ 2 or ‖ u ‖ ∞ (e.g., see Line 10 of Algorithm 1) | 2 mn        | mn           | 3 mn             |
| vector scaling (see Line 11 of Algorithm 1)           | mn          | mn           | mn               |

## 4 Numerical experiments

In this section, numerical experiments are reported to illustrate the efficiency of the Hessenberg-based PageRank algorithm presented in this paper also against other popular PageRank algorithms, that are the conventional power method including its variants with quadratic extrapolation [2] and with linear extrapolation [13], the Arnoldi-based PageRank method introduced in [22], and the adaptively accelerated Arnoldi method [28] . 3 The performance of these methods were assessed in terms of number of matrixvector products (or, equivalently, number of iteration steps for the first three algorithms) and elapsed CPU time (in seconds) required to achieve convergence to a prescribed accuracy. Unless otherwise stated, the stopping criterion used in our runs was

$$\frac { \| A q - q \| _ { 1 } } { \| q \| _ { 1 } } < t o l \ = 1 0 ^ { - 8 },$$

and all the algorithms were started from the initial vector q 0 = e / ‖ e ‖ 1 , where e = [1 , . . . , 1] T . According to Theorem 3.1, the cost of implementing the above stopping criterion can be alleviated for both Arnoldiand Hessenberg-type methods, since A q -q = σ m Q m +1 u m and the computation of σ m Q m +1 u m is actually cheaper than that of A q -q , when m is not large. In our experiments with the method denoted as QE-power , the quadratic extrapolation technique was applied every five iterations, following the observations made in [2]. The experiments were run in MATLAB R2017b (64 bit) on a computer equipped with Intel Core i5-8250U processor (CPU 1.60 ∼ 1.80 GHz), 8 GB of RAM using double precision floating point arithmetic with machine epsilon set equal to 10 -16 .

The matrix problems used in our runs are obtained from the SuiteSparse Matrix Collection , which is available online at https://sparse.tamu.edu/ . In Table 3, we describe the characteristics of our test matrices, including number of rows ( n ), number of nonzeros ( N z ), number of zero columns ( zcol ), average nonzeros of every row ( aN z ) and density ( den ) which is defined as

$$d e n = \frac { N _ { z } } { n \times n } \times 1 0 0.$$

Here, the number of zero columns corresponds to the number of dangling nodes. The largest problem in our set has size 5 , 363 260 and 79 023 142 nonzeros. , , ,

3 Numerical results with the IDR( s )-based PageRank method are omitted due to its unsatisfactory performance for large values of s and m . However, the MATLAB code of the IDR(s)-based PageRank method is still included in our GitHub repository: https://github.com/Hsien-Ming-Ku/PageRank-Hessenberg for testing purposes.

Table 3: The characteristic of test matrices.

| Matrix ID   | Matrix name      | Size      | Nonzeros   | zcol      |   aN z | den              |
|-------------|------------------|-----------|------------|-----------|--------|------------------|
| I           | soc-Slashdot0902 | 82,168    | 948,464    | 3,727     | 11.543 | 1 . 405 × 10 - 2 |
| II          | amazon0312       | 400,727   | 3,200,440  | 12,353    |  7.987 | 1 . 993 × 10 - 3 |
| III         | amazon-2008      | 735,323   | 5,158,388  | 88,557    |  7.015 | 9 . 540 × 10 - 4 |
| IV          | wiki-Talk        | 2,394,385 | 5,021,410  | 2,246,783 |  2.097 | 8 . 759 × 10 - 5 |
| V           | ljournal-2008    | 5,363,260 | 79,023,142 | 545,626   | 14.734 | 2 . 747 × 10 - 4 |

## 4.1 Choice of the restart value m

First, we investigate the effect of the restart parameter m on the convergence of the Arnoldi (A-P) , A-Arnoldi (GA-P) and Hessenberg (H-P) methods in terms of number of iterations and elapsed CPU time, since this parameter may noticeably affect the performance of Krylov subspace-based methods. The results are presented numerically in Tables 4-6. In Figs. 3-4, for the test matrix ' soc-Slashdot0902 ' we plot the curves showing the total CPU time versus m for different damping factors and tol 's values.

Table 4: Number of iterations required by the Arnoldi-, GArnoldi- and Hessenberg-based algorithms with different restart numbers. (problem soc-Slashdot0902 and tol = 10 -7 ).

| α    | m = 6   | m = 6   | m = 6   | m = 7   | m = 7   | m = 7   | m = 8   | m = 8   | m = 8   | m = 9   | m = 9   | m = 9   | m = 10   | m = 10   | m = 10   |
|------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|----------|----------|----------|
|      | A-P     | GA-P    | H-P     | A-P     | GA-P    | H-P     | A-P     | GA-P    | H-P     | A-P     | GA-P    | H-P     | A-P      | GA-P     | H-P      |
| 0.85 | 5       | 5       | 5       | 4       | 4       | 4       | 3       | 3       | 4       | 3       | 3       | 3       | 3        | 3        | 3        |
| 0.90 | 6       | 6       | 7       | 5       | 5       | 5       | 4       | 4       | 4       | 4       | 4       | 4       | 3        | 3        | 4        |
| 0.95 | 9       | 8       | 9       | 7       | 7       | 8       | 5       | 5       | 6       | 5       | 5       | 6       | 4        | 4        | 4        |
| 0.99 | 16      | 15      | 21      | 15      | 13      | 15      | 11      | 10      | 11      | 9       | 8       | 10      | 7        | 7        | 8        |

Table 5: Number of iterations required by the Arnoldi-, GArnoldi- and Hessenberg-type algorithms with different restart numbers (problem soc-Slashdot0902 and tol = 10 -8 ).

| α    | m = 6   | m = 6   | m = 6   | m = 7   | m = 7   | m = 7   | m = 8   | m = 8   | m = 8   | m = 9   | m = 9   | m = 9   | m = 10   | m = 10   | m = 10   |
|------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|----------|----------|----------|
|      | A-P     | GA-P    | H-P     | A-P     | GA-P    | H-P     | A-P     | GA-P    | H-P     | A-P     | GA-P    | H-P     | A-P      | GA-P     | H-P      |
| 0.85 | 6       | 5       | 6       | 5       | 4       | 5       | 4       | 4       | 4       | 4       | 4       | 4       | 3        | 3        | 3        |
| 0.90 | 7       | 7       | 7       | 6       | 5       | 6       | 5       | 5       | 5       | 4       | 4       | 5       | 4        | 4        | 4        |
| 0.95 | 11      | 9       | 10      | 8       | 7       | 9       | 6       | 6       | 7       | 6       | 5       | 6       | 5        | 5        | 5        |
| 0.99 | 18      | 19      | 24      | 17      | 15      | 16      | 14      | 11      | 12      | 11      | 10      | 12      | 9        | 8        | 9        |

According to the results reported in Tables 4-6, the number of iterations required to converge by these three algorithms tends to decrease for higher restart numbers m , especially for larger damping factors. This behaviour is expected because larger search spaces may provide better approximations. On the other hand, the total solution time of the three methods is not significantly reduced. As mentioned in Section 3, and explained in Tables 1-2, the storage requirement and the computational cost of one Arnoldi and Hessenberg cycles increase with m . However, it should be noted that Hessenberg is more cost effective than both Arnoldi and A-Arnoldi for larger m . In our numerical experiments, we choose

the restart numbers equal to m = 8 10 due to memory constraints , since the number of iterations and the , 4 total elapsed CPU time are still acceptable. It may be worth investigating techniques that can effectively reduce the dimension of the Krylov subspace for the Hessenberg method, e.g., by optimizing the choice of the starting vector [52] and by utilizing vector extrapolations [13], but this analyses is beyond the scope of this study.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-2c06a56a-b809-46b9-be15-1a5aff0d8550.jpg)

m

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-99308990-ee26-4167-a8b0-381e1435edf1.jpg)

m

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-54239943-8df9-4069-94f2-f927da58678c.jpg)

m

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-51959a0d-c905-4df0-9dd9-3e15eca05dc4.jpg)

m

Fig. 3: Plot of the elapsed CPU time in seconds versus the restart number (i.e., m ) for the test problem ' soc-Slashdot0902 ' using tol = 10 -7 .

## 4.2 Effect of damping factors on the CPU time and the number of iterations

For the five matrix problems listed in Table 3, we report on the number of matrix-vector products ( Mvp in short) and the elapsed CPU time of the power method, the power methods with quadratic extrapolation

4 In fact, the restart number is often chosen in the PageRank literature as m ≤ 10 [27,30,31].

Fig. 4: Plot of the elapsed CPU time in seconds versus the restart number m for the test problem ' soc-Slashdot0902 ' using tol = 10 -8 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-a4e9edd7-3e64-450d-80bd-d01c8ab52c99.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-49194f76-087b-4796-bc87-4433536abea4.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-46d4ee43-f54f-45d8-86b8-b8b7723fca3c.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00235v2-figure-99dae7fb-c210-47c1-92ce-c5fcd05afaae.jpg)

and with linear extrapolation, the Arnoldi-type method, the adaptively accelerated Arnoldi method and the Hessenberg-type method for various values of the damping factor α ranging from 0 85 to 0 99. . .

We can see from the results of Table 6 that the power method accelerated by quadratic extrapolation outperforms the conventional power method and its linearly extrapolated variant, while in most cases our Hessenberg-based solver is the fastest method in terms of elapsed CPU time, with the only exception for matrix 'IV' using α = 0 99. It can be observed that . Arnoldi is more cost-effective than A-Arnoldi at equal number of Mvp , especially for large problems. This behaviour is in agreement with the cost analysis presented in Table 2. Except these few cases, the A-Arnoldi method is still attractive to consider. On the other hand, one observes from Table 6 that the numerical behaviors of the Arnoldi , A-Arnoldi and Hessenberg algorithms relies on the choice of m and α . For example, when m is small, say m = 8, these three algorithms are only slightly better than the Power , Power-Tan and QE-Power methods. However, as m and α increase, their improvements become gradually more significant. Unlike the Arnoldi A-Arnoldi , ,

14

Table 6: Matrix-vector products and CPU time in seconds versus damping factors

| ID   | α    |      | Power   | Power-Tan   |       | QE-Power   | QE-Power   | Arnoldi   | Arnoldi   | Arnoldi   | Arnoldi   | A-Arnoldi   | A-Arnoldi   | A-Arnoldi   | A-Arnoldi   | Hessenberg   | Hessenberg   | Hessenberg   | Hessenberg   |
|------|------|------|---------|-------------|-------|------------|------------|-----------|-----------|-----------|-----------|-------------|-------------|-------------|-------------|--------------|--------------|--------------|--------------|
|      | α    |      |         |             |       |            |            | m =       | m =       | m =       | m =       | m = 8       | m = 8       | m = 10      | m = 10      | m =          | m =          | m =          | m =          |
|      | α    | Mvp  | CPU     | Mvp         | CPU   | Mvp        | CPU        | Mvp       | CPU       | Mvp       | CPU       | Mvp         | CPU         | Mvp         | CPU         | Mvp          | CPU          | Mvp          | CPU          |
| I    | 0.85 | 82   | 0.216   | 83          | 0.228 | 52         | 0.168      | 32        | 0.152     | 30        | 0.136     | 32          | 0.122       | 30          | 0.132       | 32           | 0.122        | 30           | 0.118        |
| I    | 0.90 | 125  | 0.324   | 126         | 0.333 | 52         | 0.159      | 40        | 0.186     | 40        | 0.187     | 40          | 0.159       | 40          | 0.163       | 40           | 0.144        | 40           | 0.146        |
| I    | 0.95 | 241  | 0.628   | 243         | 0.634 | 98         | 0.298      | 48        | 0.192     | 50        | 0.206     | 48          | 0.195       | 50          | 0.202       | 56           | 0.187        | 50           | 0.176        |
| I    | 0.99 | 824  | 2.089   | 832         | 2.105 | 169        | 0.504      | 112       | 0.419     | 90        | 0.350     | 88          | 0.327       | 80          | 0.308       | 96           | 0.319        | 90           | 0.303        |
| II   | 0.85 | 80   | 1.032   | 80          | 1.035 | 66         | 1.014      | 40        | 0.951     | 40        | 0.943     | 40          | 0.996       | 40          | 1.089       | 48           | 0.854        | 40           | 0.726        |
| II   | 0.90 | 120  | 1.553   | 121         | 1.562 | 98         | 1.522      | 56        | 1.253     | 40        | 1.187     | 48          | 1.166       | 40          | 1.292       | 56           | 0.993        | 40           | 0.918        |
| II   | 0.95 | 234  | 3.033   | 235         | 3.041 | 152        | 2.422      | 80        | 1.779     | 56        | 1.681     | 72          | 1.732       | 70          | 1.802       | 88           | 1.554        | 80           | 1.393        |
| II   | 0.99 | 1071 | 13.80   | 1076        | 14.16 | 391        | 6.196      | 224       | 4.917     | 190       | 4.404     | 176         | 4.212       | 160         | 4.095       | 216          | 3.798        | 250          | 4.387        |
| III  | 0.85 | 79   | 1.782   | 80          | 1.821 | 64         | 1.779      | 40        | 1.633     | 40        | 1.704     | 40          | 1.761       | 40          | 1.896       | 40           | 1.315        | 40           | 1.321        |
| III  | 0.90 | 119  | 2.782   | 120         | 2.839 | 95         | 2.659      | 48        | 1.981     | 50        | 2.304     | 48          | 2.130       | 50          | 2.376       | 56           | 1.828        | 50           | 1.627        |
| III  | 0.95 | 235  | 5.248   | 236         | 5.289 | 152        | 4.215      | 72        | 2.889     | 70        | 3.034     | 72          | 3.152       | 70          | 3.283       | 72           | 2.345        | 80           | 2.632        |
| III  | 0.99 | 1082 | 24.85   | 1084        | 25.02 | 411        | 11.79      | 208       | 8.428     | 180       | 7.779     | 152         | 6.817       | 150         | 7.107       | 176          | 5.843        | 190          | 6.358        |
| IV   | 0.85 | 67   | 4.798   | 67          | 4.819 | 52         | 4.727      | 32        | 4.321     | 30        | 4.403     | 32          | 4.683       | 30          | 4.808       | 32           | 3.327        | 40           | 4.264        |
| IV   | 0.90 | 104  | 7.344   | 104         | 7.464 | 52         | 4.770      | 40        | 5.352     | 40        | 5.838     | 40          | 5.948       | 40          | 6.266       | 40           | 4.217        | 40           | 4.299        |
| IV   | 0.95 | 214  | 15.11   | 214         | 15.27 | 84         | 7.327      | 56        | 7.413     | 60        | 8.736     | 48          | 7.108       | 60          | 9.532       | 64           | 6.768        | 50           | 5.335        |
| IV   | 0.99 | 1066 | 20.21   | 1068        | 75.96 | 175        | 15.57      | 156       | 20.21     | 150       | 21.47     | 112         | 16.24       | 110         | 17.08       | 160          | 16.69        | 150          | 15.93        |
| V    | 0.85 | 80   | 36.59   | 79          | 36.47 | 61         | 30.61      | 48        | 28.53     | 50        | 30.88     | 48          | 29.79       | 50          | 32.38       | 56           | 29.44        | 50           | 26.87        |
| V    | 0.90 | 122  | 54.77   | 122         | 55.12 | 92         | 45.32      | 64        | 38.26     | 60        | 36.66     | 64          | 39.21       | 60          | 38.65       | 72           | 37.62        | 60           | 31.82        |
| V    | 0.95 | 247  | 110.9   | 244         | 109.5 | 148        | 71.61      | 104       | 61.85     | 90        | 55.92     | 96          | 59.05       | 90          | 57.84       | 104          | 54.59        | 90           | 47.97        |
| V    | 0.99 | 1153 | 794.7   | 1061        | 529.1 | 630        | 322.6      | 352       | 214.7     | 240       | 145.8     | 216         | 132.3       | 190         | 121.4       | 248          | 127.7        | 220          | 115.9        |

and Hessenberg , the Power , Power-Tan and QE-Power methods are simple and their main computational cost is the evaluation of matrix-vector products. These characteristics often make them still feasible for computing PageRank when the damping factor α is not large.

In addition, it is interesting to mention that the Hessenberg-type method often needs more Mvp 's for convergence than both Arnoldi and A-Arnoldi , whereas the total CPU time of Hessenberg is still less. This is because the Hessenberg uses the cheap similarity transformations to reduce the large matrix into the Hessenberg form, whereas the latter two methods use expensive (weighted) unitary transformations.

## 5 Conclusions

In this paper, we proposed a novel approach for computing the PageRank problem. The proposed method has lower computational cost than both Arnoldi and A-Arnoldi to find the approximate PageRank vector, thus it can afford to use higher dimensional Krylov subspaces. Extensive numerical experiments are reported to illustrate the efficiency of the proposed method also compared to other stateof-the-art matrix solvers for this problem class, especially when the damping factor is large. Hence, we conclude that the Hessenberg method as well as the Arnoldi and A-Arnoldi methods can be useful computational tools for practical large-scale PageRank computations.

Future research will focus on the theory of the Hessenberg process and the convergence of the Hessenberg-type algorithm is still required to be further analyzed. In addition, it is interesting to study how to reduce the restart number m and improve the convergence speed of our methods. Moreover, the proposed method can be extended to compute the more general Markov chains [8,29], e.g., in ProteinRank and CiteRank.

## Acknowledgements

The authors would like to thank Prof. Zhongxiao Jia for his comments about the strategy used in the refined Arnoldi algorithm. Meanwhile, the authors are grateful to Dr. Reinaldo Astudillo (ASML Holding N.V.) for his kind suggestions about executing the IDR-based Hessenberg decompositions used in Section 2.2. This research is supported by NSFC (11601323 and 11801463), the Applied Basic Research Program of Sichuan Province (2020YJ0007) and the research grants MYRG2018-00025-FST, MYRG2020-00208FST from University of Macau. The last author is member of the Gruppo Nazionale per il Calcolo Scientifico (GNCS) of the Istituto Nazionale di Alta Matematica (INdAM) and his work was partially supported by INdAM-GNCS under Progetti di Ricerca 2020.

## References

- [1] L. Page, S. Brin, R. Motwani, T. Winograd, The PageRank citation ranking: Bringing order to the web, Technical Report No. 1999-66 , Stanford InfoLab., Jan. 29, 1999, 17 pages. Available online at: http://ilpubs.stanford.edu:8090/422/ .
- [2] S.D. Kamvar, T.H. Haveliwala, C.D. Manning, G.H. Golub, Extrapolation methods for accelerating PageRank computations, in: WWW '03 Proceedings of the 12th international conference on World Wide Web , Budapest, Hungary, May 20-24, 2003, ACM New York, NY (2003): 261-270. DOI: 10.1145/775152.775190.

- [3] S. Kamvar, T. Haveliwala, G. Golub, Adaptive methods for the computation of PageRank, Linear Algebra Appl. , 386 (2004): 51-65.
- [4] A.N. Langville, C.D. Meyer, Deeper inside PageRank, Internet Math. , 1(3) (2005): 335-380.
- [5] A.N. Langville, C.D. Meyer, A survey of eigenvector methods of web information retrieval, SIAM Rev. , 47(1) (2005): 135-161.
- [6] P. Berkhin, A survey on PageRank computing, Internet Math. , 2(1) (2005): 73-120.
- [7] A.N. Langville, C.D. Meyer, Google's PageRank and Beyond: The Science of Search Engine Rankings , Princeton University Press, Princeton, NJ (2006).
- [8] D.F. Gleich, PageRank beyond the web, SIAM Rev. , 57(3) (2015): 321-363.
- [9] K. Bryan, T. Leise, The 25 000 000 000 eigenvector: The linear algebra behind Google, , , , SIAM Rev. , 48(3) (2006): 569-581.
- [10] A. Cicone, S. Serra-Capizzano, Google PageRanking problem: The model and the analysis, J. Comput. Appl. Math. , 234(11) (2010): 3140-3169.
- [11] K. Avrachenkov, N. Litvak, D. Nemirovsky, N. Osipova, Monte Carlo methods in PageRank computation: when one iteration is sufficient, SIAM J. Numer. Anal. , 45(2) (2007): 890-904.
- [12] W. Liu, G. Li, J. Cheng, Fast PageRank approximation by adaptive sampling, Knowl. Inf. Syst. , 42(1) (2015): 127-146.
- [13] X. Tan, A new extrapolation method for PageRank computations, J. Comput. Appl. Math. , 313 (2017): 383-392.
- [14] A.N. Langville, C.D. Meyer, A reordering for the PageRank problem, SIAM J. Sci. Comput. , 27(6) (2006): 2112-2120.
- [15] Y. Lin, X. Shi, Y. Wei, On computing PageRank via lumping the Google matrix, J. Comput. Appl. Math. , 224(2) (2009): 702-708.
- [16] D.F. Gleich, A.P. Gray, C. Greif, T. Lau, An inner-outer iteration for computing PageRank, SIAM J. Sci. Comput. , 32(1) (2010): 349-371.
- [17] M. Heyouni, H. Sadok, On a variable smoothing procedure for Krylov subspace methods, Linear Algebra Appl. , 268 (1998): 131-149.
- [18] Y. Saad, Numerical Methods for Large Eigenvalue Problems (Revised Ed.), SIAM, Philadelphia, PA (2011).
- [19] Z. Jia, Refined iterative algorithms based on Arnoldi's process for large unsymmetric eigenproblems, Linear Algebra Appl. , 259 (1997): 1-23.
- [20] Z. Jia, Polynomial characterizations of the approximate eigenvectors by the refined Arnoldi method and an implicitly restarted refined Arnoldi algorithm, Linear Algebra Appl. , 287(1-3) (1999): 191-214.

- [21] Z. Jia, A refined subspace iteration algorithm for large sparse eigenproblems, Appl. Numer. Math. , 32(1) (2000): 35-52.
- [22] G.H. Golub, C. Greif, An Arnoldi-type algorithm for computing page rank, BIT , 46(4) (2006): 759-771.
- [23] G. Wu, Y. Wei, A Power-Arnoldi algorithm for computing PageRank, Numer. Linear Algebra Appl. , 14(7) (2007): 521-546.
- [24] G.-J. Yin, J.-F. Yin, On Arnoldi method accelerating PageRank cmputations, in: Web Information Systems and Mining. WISM 2010 (F.-L. Wang, Z. Gong, X. Luo, J. Lei, eds.) , Lecture Notes in Computer Science, vol 6318, Springer, Berlin, Heidelberg (2010): 378-385. DOI: 10.1007/978-3-642-16515-3 47.
- [25] G. Wu, Y. Zhang, Y. Wei, Accelerating the Arnoldi-type algorithm for the PageRank problem and the ProteinRank problem, J. Sci. Comput. , 57(1) (2013): 74-104.
- [26] G. Wu, Y. Wei, An Arnoldi-extrapolation algorithm for computing PageRank, J. Comput. Appl. Math. , 234(11) (2010): 3196-3212.
- [27] C. Gu, W. Wang, An Arnoldi-Inout algorithm for computing PageRank problems, J. Comput. Appl. Math. , 309 (2017): 219-229.
- [28] J.-F. Yin, G.-J. Yin, M. Ng, On adaptively accelerated Arnoldi method for computing PageRank, Numer. Linear Algebra Appl. , 19(1) (2012): 73-85.
- [29] R.W. Freund, M. Hochbruck, On the use of two QMR algorithms for solving singular systems and applications in Markov chain modeling, Numer. Linear Algebra Appl. , 1(4) (1994): 403-420.
- [30] K. Teramoto, T. Nodera, A note on Lanczos algorithm for computing PageRank, in: Forging Connections between Computational Mathematics and Computational Geometry (K. Chen, A. Ravindran, eds.) , Springer Proceedings in Mathematics &amp; Statistics, Vol. 124, Springer, Cham, Switzerland (2016): 25-33. DOI: 10.5176/2251-1911 CMCGS14.15 3.
- [31] G. Wu, Y.-C. Wang, X.-Q. Jin, A preconditioned and shifted GMRES algorithm for the PageRank problem with multiple damping factors, SIAM J. Sci. Comput. , 34(5) (2012): A2558-A2575.
- [32] G. Wu, Y. Wei, Arnoldi versus GMRES for computing pageRank: A theoretical contribution to google's pageRank problem, ACM Trans Inf. Syst. , 28(3) (2010): 11. DOI: 10.1145/1777432.1777434.
- [33] K. Hessenberg, BehandLung Linearer Eigenwertaufgaben Mit Hilfe Der Hamilton-Cayleyschen Gleichung , Numerische Verfahren, Bericht 1, Institut f¨r Praktische Mathematik (IPM), Technische u Hochschule Darmstadt (1940). The scanned report and a biographical sketch of Karl Hessenberg's life are available at http://www.hessenberg.de/karl1.html .
- [34] J.H. Wilkinson, The Algebraic Eigenvalue Problem , Clarendon Press, Oxford, UK (1965).
- [35] H. Sadok, CMRH: A new method for solving nonsymmetric linear systems based on the Hessenberg reduction algorithm, Numer. Algorithms , 20(4) (1999), pp. 303-321.

- [36] D. Stephens, ELMRES: An Oblique Projection Method To Solve Sparse Non-Symmetric Linear Systems (Ph.D dissertation) , Florida Institute of Technology, Melbourne, USA (1999). Available online at http://ncsu.edu/hpc/Documents/Publications/gary\_howell/stephens.pdf .
- [37] A.S. Householder, F.L. Bauer, On certain methods for expanding the characteristic polynomial, Numer. Math. , 1(1) (1959): 29-37.
- [38] H. Sadok, D.B. Szyld, A new look at CMRH and its relation to GMRES, BIT , 52(2) (2012): 485-501.
- [39] M. Heyouni, H. Sadok, A new implementation of the CMRH method for solving dense linear systems, J. Comput. Appl. Math. , 213(2) (2008): 387-399.
- [40] K. Zhang, C. Gu, Flexible global generalized Hessenberg methods for linear systems with multiple right-hand sides, J. Comput. Appl. Math. , 263 (2014): 312-325.
- [41] M. Heyouni, The global Hessenberg and CMRH methods for linear systems with multiple right-hand sides, Numer. Algorithms , 26(4) (2001): 317-332.
- [42] M. Heyouni, A. Essai, Matrix Krylov subspace methods for linear systems with multiple right-hand sides, Numer. Algorithms , 40(2) (2005): 137-156.
- [43] S. Amini, F. Toutounian, M. Gachpazan, The block CMRH method for solving nonsymmetric linear systems with multiple right-hand sides, J. Comput. Appl. Math. , 337 (2018): 166-174.
- [44] S. Amini, F. Toutounian, Weighted and flexible versions of block CMRH method for solving nonsymmetric linear systems with multiple right-hand sides, Comput. Math. Appl. , 76(8) (2018): 2011-2021.
- [45] X.-M. Gu, T.-Z. Huang, G. Yin, B. Carpentieri, C. Wen, L. Du, Restarted Hessenberg method for solving shifted nonsymmetric linear systems, J. Comput. Appl. Math. , 331 (2018): 166-177.
- [46] X.-M. Gu, T.-Z. Huang, B. Carpentieri, A. Imakura, K. Zhang, L. Du, Efficient variants of the CMRH method for solving a sequence of multi-shifted non-Hermitian linear systems simultaneously, J. Comput. Appl. Math. , 375 (2020): 112788. DOI: 10.1016/j.cam.2020.112788.
- [47] Z. Ramezani, F. Toutounian, Extended and rational Hessenberg methods for the evaluation of matrix functions, BIT , 59(2) (2019): 523-545.
- [48] M. Addam, M. Heyouni, H. Sadok, The block Hessenberg process for matrix equations, Electron. Trans. Numer. Anal. , 46 (2017): 460-473.
- [49] M. Heyouni, F. Saberi-Movahed, A. Tajaddini, On global Hessenberg based methods for solving Sylvester matrix equations, Comput. Math. Appl. , 77(1) (2019): 77-92.
- [50] P.A. Businger, Reducing a matrix to Hessenberg form, Math. Comp. , 23(108) (1969): 819-821.
- [51] M. Heyouni, Newton Generalized Hessenberg method for solving nonlinear systems of equations, Numer. Algorithms , 21(1-4) (1999): 225-246.
- [52] R. Astudillo, M.B. van Gijzen, A restarted Induced Dimension Reduction method to approximate eigenpairs of large unsymmetric matrices, J. Comput. Appl. Math. , 296 (2016): 24-35.

- [53] M.H. Gutknecht, J.-P.M. Zemke, Eigenvalue computations based on IDR, SIAM J. Matrix Anal. Appl. , 34(2) (2013): 283-311.
Xian-Ming Gu Siu-Long Lei Ke Zhang Zhao-Li Shen Chun Wen Bruno Carpentieri
A Hessenberg-type Algorithm for Computing PageRank Problems
PageRank is a widespread model for analysing the relative relevance of nodes
within large graphs arising in several applications. In the current paper, we
present a cost-effective Hessenberg-type method built upon the Hessenberg
process for the solution of difficult PageRank problems. The new method is very
competitive with other popular algorithms in this field, such as Arnoldi-type
methods, especially when the damping factor is close to $1$ and the dimension
of the search subspace is large. The convergence and the complexity of the
proposed algorithm are investigated. Numerical experiments are reported to show
the efficiency of the new solver for practical PageRank computations.