# DNSMOS: A Non-Intrusive Perceptual Objective Speech Quality metric to
evaluate Noise Suppressors

###### Abstract

Human subjective evaluation is the “gold standard” to evaluate speech quality
optimized for human perception. Perceptual objective metrics serve as a proxy
for subjective scores. The conventional and widely used metrics require a
reference clean speech signal, which is unavailable in real recordings.
Previous no-reference approaches correlate poorly with human ratings and are
not widely adopted in the research community. One of the biggest use cases of
these perceptual objective metrics is to evaluate noise suppression
algorithms. This paper introduces a multi-stage self-teaching based perceptual
objective metric that is designed to evaluate noise suppressors. The proposed
method generalizes well in challenging test conditions with a high correlation
to human ratings.

Index Terms—  Speech, Perceptual Speech Quality, Objective Metric, Deep Noise
Suppressor, Metric.

##  1 Introduction

Subjective evaluation of speech quality is the most reliable way to evaluate
Speech Enhancement (SE) methods [1]. However, subjective tests are not
scalable as they require a considerable number of listeners, the process is
laborious, time-consuming, and expensive. Conventional objective speech
quality metrics such as Perceptual Evaluation of Speech Quality (PESQ) [2],
Perceptual Objective Listening Quality Analysis (POLQA) [3] and Signal to
Distortion Ratio (SDR) are widely used to evaluate Speech Enhancement (SE)
algorithms optimized for human perception. Some of these metrics are designed
to predict the subjective Mean Opinion Score (MOS) obtained using the
Telecommunication Standardization Sector of the International
Telecommunication Union (ITU-T) Recommendation P.800 [4]. There are also the
intrusive metrics requiring a reference clean speech to compute the quality
score. However, they are shown to correlate poorly with human rating when used
for SE tasks that involve perceptually invariant transformations [1]. Also,
intrusive metrics cannot be used to evaluate real recordings when clean
reference is unavailable in realistic scenarios.

ITU-T Recommendation P.563 is a non-intrusive technique and can directly
operate on the degraded signal [5]. However, it was developed for narrow-band
applications and works well on limited impairment types. Recently, Deep Neural
Networks (DNN) based approaches have been proposed to estimate the speech
quality scores [6, 7, 8]. Some of these learning-based approaches use other
objective metrics as the ground truth to train their speech quality predictor.
A handful of published methods use MOS obtained using P.800 as the ground
truth to train their models. In [9], the authors trained the model to identify
the Just Noticeable Difference (JND). MOS predictors trained on actual human
ratings are more reliable than the ones trained to predict other objective
metrics like PESQ or POLQA. The accuracy and robustness of the learned models
depend on the quality of the human labels and also the quantity and diversity
of the audio clips. There is no large scale human-labeled data set publicly
available to train a model that can robustly generalize to a variety of audio
impairments, especially for SE task. Though the learning-based approaches have
a higher correlation to human ratings when tested on a limited and small test
set, none of these speech quality metrics are as widely used as some of the
intrusive objective metrics. Most of the publications in the area of noise
suppression and SE test their methods on synthetic noisy speech test sets
using intrusive speech quality metrics such as PESQ, POLQA and SDR.

In this paper, we have developed a robust objective perceptual speech quality
metric called DNSMOS. It can be used to stack rank different Deep Noise
Suppression (DNS) methods based on MOS estimates with great accuracy and hence
the name DNSMOS. A simple Convolutional Neural Network (CNN) based model is
trained using the ground truth human ratings obtained using ITU-T P.808 [10,
11]. P.808 is an online subjective testing framework that is highly
reproducible and shown to stack rank noise suppression models with high
accuracy when each model is tested as an average over a statistically
significant number of clips. However, P.808 ratings per clip can be noisy due
to various factors. The absolute MOS per clip might vary for the same clip in
different P.808 runs due to the inherent noise in human rating [12]. The
number of ratings per clip varies due to spam removal. Also, all the spammers
might not be filtered due to the limitations of spam filtering logic in P.808
implementation. We train a multi-stage self-teaching model inspired by
continual life long learning [13] to learn in the presence of noisy labels. We
show that the multi-stage training significantly improves the correlation and
generalizes remarkably well to other impairment types that are very different
from what was used in the training set. We found the DNSMOS very useful in our
audio/speech research as it generalizes well to a variety of speech
impairments. Hence, we are providing the DNSMOS as an Azure service for other
researchers to use. The details of the API are at www.microsoft.com/en-
us/research/dns-challenge/dnsmos.

##  2 Data and subjective ratings

We used 600 noisy speech test clips comprised of a combination of synthetic
and real recordings. The synthetic test set is a combination of both non-
reverberant and reverberant clips. The real recordings were captured in a
variety of noise types and Signal to Noise Ratio (SNR) and target levels. The
test set is comprised of over 100 noise types and speakers. More details about
the creation of these test sets can be found in [11]. These 600 clips were
processed by over 200 different noise suppression methods. Some of these
methods improved the overall speech quality, but they also introduced a
variety of artifacts as a consequence of over or under suppression of noise,
resulting in speech and noise distortions. The speech quality ratings of the
processed clips varied from very poor (MOS=1) to excellent (MOS=5). The
distribution of the MOS scores in the training data is shown in Figure 1a. The
scores are highly skewed with most ratings populated in the range 3<MOS<4 and
fewer ratings in both the tails. The distributions in Figure 1a also shows
that the current state of the art noise suppression techniques are not good
enough to shift the overall distribution above MOS of 4. Figure 1b shows the
distribution of the number of ratings per clip in the training set. The
skewness in the distribution is due to the fact that some of the ratings were
removed in the process of removing spam ratings. Also, we used 10 ratings per
clip in some of the experiments and 5 ratings per clip in others. Figure 1c
shows the distribution of standard deviation per clip. There are a significant
number of clips with a standard deviation of over 1 MOS, which shows that the
audio clips are highly subjective in nature. Increasing the number of ratings
per clip will reduce the standard deviation, but also increases the cost.

The subjective human ratings are obtained in several P.808 runs conducted over
several months. Multiple noise suppression methods are compared in each P.808
run. Each P.808 run included the best performing noise suppressor, original
noisy speech, and a couple of methods with intermediate perceptual quality
from previous runs as anchors. Hence, some of the clips were rated multiple
times. In total, we have about 120,000 audio clips with associated MOS scores
as ground truth. The average length of each audio clip was about 9 seconds. In
[10], we show that P.808 is highly reproducible with high SRCC between runs
when we average the ratings across clips processed by a given model. However,
the absolute values of the scores vary between P.808 runs as the raters differ
and humans tend to be inconsistent in perceptual tasks due to various biases.

(a) Distribution of MOS ratings

(b) Distribution of number of votes

(c) Distribution of standard deviation per clip

Fig. 1: Statistics of the training set

##  3 DNSMOS

###  3.1 Learning in the presence of label biases

The variation in the absolute values between P.808 runs might pose a challenge
in training models that learn to predict MOS scores when the training data
includes the same audio clip with different ratings. We know from the
literature that DNN models can generalize well in spite of the noise for
classification tasks [14, 15]. The variation in the absolute scores from
different runs can be treated as noise in the labels. Also, the clips with
fewer ratings can be treated as noise due to higher standard deviation. There
are other unknown human biases that play a role in the variation of ratings.
We initially took a direct approach where we trained a deep model by including
all the available data without any modification. The model tends to learn the
weighted average of the ratings across runs and also tries to average out
other biases. This plain approach gave a good correlation with human ratings
when tested on the data that was very similar to the training set. However,
the correlation dropped when tested on a more challenging test set that was
very different from the training data.

We improved the generalization and accuracy of the predictor by training a few
stages of self-teaching models. In [16], a sequence of self-teachers were used
to improve the generalizability of a classification task with binary cross-
entropy loss in noisy and weakly labeled conditions. We use a similar approach
for the regression problem at hand with Mean Squared Error (MSE) as the loss
function. Unlike knowledge distillation where a larger teacher model is used
to produce soft labels to train a smaller student model, we use the same model
architecture for both teacher and the student, and hence the name self-
teacher. The student model uses the weighted average of the predictions from
one or multiple teacher models and the original human MOS ratings. The primary
teacher model $\textbf{{M}}_{0}$ is trained until the loss saturates using the
original ground truth human MOS ratings r. The student model
$\textbf{{M}}_{s}$ at stage s, is trained using the new target given by,

$\hat{\textbf{{r}}}_{s}=\alpha_{0}\textbf{{r}}+\sum_{i=0}^{s-1}\alpha_{i+1}\hat{\textbf{{r}}}_{i}$
(1)

where $\hat{\textbf{{r}}}_{i}$ is the prediction of model $\textbf{{M}}_{i}$
and $\sum_{i=0}^{s}\alpha_{i}=1$. In this work, we restrict s to 2. With some
assumptions to the noise distribution, it can be theoretically proven that the
$\textbf{{M}}_{s}$ is better than $\textbf{{M}}_{0}$ for $\textit{s}>0$. Since
the distribution of noise in our case is complex, we show through empirical
experiments that the predicted MOS of $\textbf{{M}}_{i}$ is of higher accuracy
than $\textbf{{M}}_{0}$ for appropriately chosen values of $\alpha$.

###  3.2 Features

Recently, researchers have seen success in learning features within the model
for tasks such as SE [17], speech, and music synthesis [18] and to learn
acoustic models [19]. They show that using the time-domain waveform requires a
larger model trained on a larger and diverse data set to ensure
generalization. The ground truth MOS scores are obtained for audio clips with
an average length of 9 secs sampled at 16 kHz. This leads to a very large
input dimension if we are treating it as a vector and the model requires many
layers to compress and extract input features. We used Log power Mel
spectrogram as input feature as it correlates well with human perception and
is proven to work very well for analyzing speech quality [7]. For spectral
features, we used a frame size of 20 ms with a hop length of 10 ms and 120 Mel
frequency bands. The input features are then converted to dB scale.

###  3.3 Prediction model

For predicting the MOS scores, we explored different configurations of CNN
based models. The architecture for the best performing model is shown in Table
1, which is the architecture for all M. The input to the model is log power
Mel Spectrogram with 120 Mel bands computed over a clip of length 9 secs
sampled at 16 kHz with a frame size of 20 ms and hop length of 10 ms. This
results in an input dimension of 900 x 120. The model was trained with a batch
size of 32 using the Adam optimizer and MSE loss function until the loss
saturated. We experimented by adding batch normalization layers after every
Conv layer in Table 1. However, adding batch normalization reduces the
prediction accuracy of low volume clips. Humans tend to give lower ratings to
the clips with low amplitudes. We want the model to capture the variations in
the target levels of the data. Hence, we avoid any kind of feature
normalization. We also explored different network architectures including CNN
followed by LSTM. The model in Table 1 generalized the best and was of least
complexity.

Table 1: DNSMOS Prediction Model Layer | Output dimension
---|---
Input | 900 x 120 x 1
Conv: 32, (3 x 3), ‘ReLU’ | 900 x 120 x 32
MaxPool: (2 x 2), Dropout(0.3) | 450 x 60 x 32
Conv: 32, (3 x 3), ‘ReLU’ | 450 x 60 x 32
MaxPool: (2 x 2), Dropout(0.3) | 225 x 30 x 32
Conv: 32, (3 x 3), ‘ReLU’ | 225 x 30 x 32
MaxPool: (2 x 2), Dropout(0.3) | 112 x 15 x 32
Conv: 64, (3 x 3), ‘ReLU’ | 112 x 15 x 64
GlobalMaxPool | 1 x 64
Dense: 64, ‘ReLU’ | 1 x 64
Dense: 64, ‘ReLU’ | 1 x 64
Dense: 1 | 1 x 1
Table 2: Correlation of DNSMOS with other widely used objective metrics | PESQ | SDR | POLQA | DNSMOS ($\textbf{{M}}_{0}$)
---|---|---|---|---
PCC | 0.78 | 0.23 | 0.79 | 0.93
SRCC | 0.82 | 0.25 | 0.84 | 0.94

##  4 Experimental Results

###  4.1 Evaluation metric

Pearson Correlation Coefficient (PCC) or MSE between the predictions of the
developed objective metric and the ground truth human ratings is commonly used
to measure the accuracy of the model [7, 8]. From the earlier discussion, we
know that P.808 correlation is highly repeatable between runs when averaged
across a set of clips per condition, which can be formed by grouping clips
enhanced by a particular SE model or based on other criteria like SNR or
reverb RT60 times. The PCC computed on the average of ratings per group across
different runs is >0.9. We also found that PCC computed on the same clips but
from two different P.808 runs is only about 0.5 due to the high rating noise
per clip.

Hence, for stack ranking different noise suppressors we evaluate by computing
the average of ratings across the entire test set for each model. Therefore,
we compute SRCC and PCC between averaged human ratings and averaged DNSMOS per
model. SRCC gives us the stack ranking accuracy of various SE models.

###  4.2 Results

####  4.2.1 DNSMOS vs other objective speech quality metrics

Table 2 shows the PCC and SRCC between human ratings and widely used PESQ,
POLQA, SDR and DNSMOS ($\textbf{{M}}_{0}$) computed on the Interspeech DNS
challenge blind test results [11]. The results show that SDR correlates poorly
with MOS. PESQ and POLQA reasonable PCC and SRCC. But the primary DNSMOS model
($\textbf{{M}}_{0}$) correlates significantly better than other objective
metrics and is highly reliable in stack ranking the SE models. The scatter
plots in Figure. 2 shows that DNSMOS aligns with MOS better than SDR, PESQ,
and POLQA. Both PESQ and POLQA have a cluster on the left top of the graph
indicating that they tend to penalize certain artifacts more than we humans
do.

(a) MOS vs PESQ

(b) MOS vs SDR

(c) MOS vs POLQA

(d) MOS vs DNSMOS

Fig. 2: Scatter plots of different objective metrics vs MOS Fig. 3: SRCC
between MOS and $\textbf{{M}}_{1}$ for different $\alpha_{0}$ values  Table 3:
Per category SRCC between ($\textbf{{M}}_{0}$) and MOS on DNS Challenge 2
results

$\textbf{{M}}_{i}$($\alpha_{0}$,$\alpha_{1}$,$\alpha_{2}$) | English | Non-English | Tonal | Emotional | Singing | Overall
---|---|---|---|---|---|---
$\textbf{{M}}_{0}$(1 ,- ,-) | 0.78 | 0.67 | 0.51 | 0.69 | 0.36 | 0.65
$\textbf{{M}}_{1}$(0.8, 0.2, -) | 0.81 | 0.86 | 0.81 | 0.58 | 0.61 | 0.90
$\textbf{{M}}_{1}$(0.6, 0.4, -) | 0.90 | 0.85 | 0.82 | 0.72 | 0.58 | 0.88
$\textbf{{M}}_{1}$(0.4, 0.6, -) | 0.93 | 0.92 | 0.87 | 0.67 | 0.57 | 0.88
$\textbf{{M}}_{1}$(0.2, 0.8, -) | 0.90 | 0.85 | 0.82 | 0.72 | 0.58 | 0.88
$\textbf{{M}}_{2}$(-, 0.5, 0.5) | 0.90 | 0.86 | 0.81 | 0.75 | 0.54 | 0.90
$\textbf{{M}}_{2}$(0.6, 0.2, 0.2) | 0.91 | 0.85 | 0.82 | 0.78 | 0.72 | 0.89
$\textbf{{M}}_{2}$(0.8, 0.1, 0.1) | 0.88 | 0.85 | 0.82 | 0.68 | 0.69 | 0.88

####  4.2.2 Generalizability of DNSMOS

We initially analyzed the effect of $\alpha_{0}$ on $\textbf{{M}}_{1}$ to
understand the weightage given to the teacher model $\textbf{{M}}_{0}$ vs
human ratings. Figure. 3 shows the SRCC between MOS and $\textbf{{M}}_{1}$ for
different values of $\alpha_{0}$ computed on one of the P.808 runs with 12
closely stack ranked models. We see in Figure. 3 that the SRCC is higher for
$\alpha_{0}=0.8$.

The noisy and the processed clips comprised of a similar distribution as that
of the training set. In order to test the generalization capability of
training multiple stages of self-teaching models, we evaluated on the ICASSP
2021 DNS Challenge 2 final submissions [20], which comprised of the blind test
set processed by 18 different noise suppressors. The blind set included
emotional speech, singing, utterances in English and non-English languages
including a few tonal languages. The DNSMOS models were trained on the data
which was predominantly in English with very little or no emotional content
such as laughter, crying, yelling, anger, etc. The training set did not
include singing and other non-English languages as well. Table 3 shows the
SRCC between the DNSMOS models at various stages and P.808 MOS for each
category in the blind set. The SRCC of the primary model $\textbf{{M}}_{0}$ is
quite low as it does not generalize well for categories other than English.
$\textbf{{M}}_{1}$ significantly improved SRCC from 0.65 to 0.9 for
$\alpha_{0}=0.8$. The optimal values of $\alpha_{0}$ were different for each
category in table 3. We also experimented with $\textbf{{M}}_{2}$ for a couple
of settings of ($\alpha_{0}$,$\alpha_{1}$,$\alpha_{2}$). For
$\textbf{{M}}_{2}$(-, 0.5, 0.5) with $\textbf{{M}}_{1}$(0.2, 0.8, -) for stage
2 gave the best results. The other settings did not show much improvement over
$\textbf{{M}}_{1}$.

##  5 Conclusion and Future work

DNSMOS is a robust speech quality metric designed to stack rank noise
suppressors with great accuracy. The multi-stage self-teaching approach
significantly improved accuracy. In the future, we would like to deepen our
theoretical understanding of how the choice of $\alpha_{i}$ influences the
accuracy of $\textbf{{M}}_{i}$ for different datasets. This will help in
training DNSMOS on other impairment types such as network distortions, codec
artifacts, and reverberation without the need for large scale data with human
labels.

## References

  * [1] Chandan KA Reddy, Ebrahim Beyrami, Jamie Pool, Ross Cutler, Sriram Srinivasan, and Johannes Gehrke,  “A scalable noisy speech dataset and online subjective test framework,”  Proc. Interspeech 2019, pp. 1816–1820, 2019.
  * [2] A. W. Rix, J. G. Beerends, M. P. Hollier, and A. P. Hekstra,  “Perceptual evaluation of speech quality (pesq)-a new method for speech quality assessment of telephone networks and codecs,”  in 2001 IEEE International Conference on Acoustics, Speech, and Signal Processing. Proceedings (Cat. No.01CH37221), 2001, vol. 2, pp. 749–752 vol.2.
  * [3] John Beerends, Christian Schmidmer, Jens Berger, Matthias Obermann, Raphael Ullmann, Joachim Pomy, and Michael Keyhl,  “Perceptual objective listening quality assessment (polqa), the third generation itu-t standard for end-to-end speech quality measurement part ii-perceptual model,”  AES: Journal of the Audio Engineering Society, vol. 61, pp. 385–402, 06 2013.
  * [4] “Itu-t recommendation p.800: Methods for subjective determination of transmission quality, itu-t recommendation p.800,”  Feb 1998.
  * [5] “Itu-t, recommendation p.563: Single-ended method for objective speech quality assessment in narrowband telephony applications,”  2004\.
  * [6] Xuan Dong and Donald S Williamson,  “An attention enhanced multi-task model for objective speech assessment in real-world environments,”  in ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2020, pp. 911–915.
  * [7] Hannes Gamper, Chandan KA Reddy, Ross Cutler, Ivan J Tashev, and Johannes Gehrke,  “Intrusive and non-intrusive perceptual speech quality assessment using a convolutional neural network,”  in 2019 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA). IEEE, 2019, pp. 85–89.
  * [8] A. R. Avila, H. Gamper, C. Reddy, R. Cutler, I. Tashev, and J. Gehrke,  “Non-intrusive speech quality assessment using neural networks,”  in ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2019, pp. 631–635.
  * [9] Pranay Manocha, Adam Finkelstein, Zeyu Jin, Nicholas J Bryan, Richard Zhang, and Gautham J Mysore,  “A differentiable perceptual audio metric learned from just noticeable differences,”  arXiv preprint arXiv:2001.04460, 2020.
  * [10] Babak Naderi and Ross Cutler,  “An open source implementation of itu-t recommendation p.808 with validation,” 2020.
  * [11] Chandan KA Reddy, Vishak Gopal, Ross Cutler, Ebrahim Beyrami, Roger Cheng, Harishchandra Dubey, Sergiy Matusevych, Robert Aichner, Ashkan Aazami, Sebastian Braun, et al.,  “The interspeech 2020 deep noise suppression challenge: Datasets, subjective testing framework, and challenge results,”  arXiv preprint arXiv:2005.13981, 2020.
  * [12] ITUT Rec,  “P. 1401, methods, metrics and procedures for statistical evaluation, qualification and comparison of objective quality prediction models,”  International Telecommunication Union, Geneva, Switzerland, 2012\.
  * [13] German I Parisi, Ronald Kemker, Jose L Part, Christopher Kanan, and Stefan Wermter,  “Continual lifelong learning with neural networks: A review,”  Neural Networks, vol. 113, pp. 54–71, 2019.
  * [14] Chandan KA Reddy, Ross Cutler, and Johannes Gehrke,  “Supervised classifiers for audio impairments with noisy labels,”  Proc. Interspeech 2019, pp. 2573–2577, 2019.
  * [15] David Rolnick, Andreas Veit, Serge J. Belongie, and Nir Shavit,  “Deep learning is robust to massive label noise,”  ArXiv, vol. abs/1705.10694, 2017.
  * [16] Anurag Kumar and Vamsi Krishna Ithapu,  “A sequential self teaching approach for improving generalization in sound event recognition,”  .
  * [17] Alexandre Defossez, Gabriel Synnaeve, and Yossi Adi,  “Real time speech enhancement in the waveform domain,” 2020.
  * [18] Kou Tanaka, Hirokazu Kameoka, Takuhiro Kaneko, and Nobukatsu Hojo,  “Wavecyclegan2: Time-domain neural post-filter for speech waveform generation,”  arXiv preprint arXiv:1904.02892, 2019.
  * [19] Y. Hoshen, R. J. Weiss, and K. W. Wilson,  “Speech acoustic modeling from raw multichannel waveforms,”  in 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2015, pp. 4624–4628.
  * [20] Chandan K A Reddy, Harishchandra Dubey, Vishak Gopal, Ross Cutler, Sebastian Braun, Hannes Gamper, Robert Aichner, and Sriram Srinivasan,  “Icassp 2021 deep noise suppression challenge,”  arXiv preprint arXiv:2009.06122, 2020.


# A Parallel PageRank Algorithm For Undirected Graph

Qi Zhang  Rongxia Tang  Zhengan Yao  Zan-Bo Zhang  Department of Mathematics,
Sun Yat-sen University, China  School of Statistics and Mathematics, Guangdong
University of Finance and Economics, China

###### Abstract

As a measure of vertex importance according to the graph structure, PageRank
has been widely applied in various fields. While many PageRank algorithms have
been proposed in the past decades, few of them take into account whether the
graph under investigation is directed or not. Thus, some important properties
of undirected graph—symmetry on edges, for example—is ignored. In this paper,
we propose a parallel PageRank algorithm specifically designed for undirected
graphs that can fully leverage their symmetry. Formally, our algorithm extends
the Chebyshev Polynomial approximation from the field of real function to the
field of matrix function. Essentially, it reflects the symmetry on edges of
undirected graph and the density of diagonalizable matrix. Theoretical
analysis indicates that our algorithm has a higher convergence rate and
requires less computation than the Power method, with the convergence rate
being up to 50% higher with a damping factor of $c=0.85$. Experiments on six
datasets illustrate that our algorithm with 38 parallelism can be up to 39
times faster than the Power method.

###### keywords:

PageRank, undirected graph, Chebyshev Polynomial, parallel

††journal: Applied Mathematics and Computation

##  1 Introduction

PageRank[1] is used to measure the importance of vertices according to graph
structure. Generally, a vertex has higher importance if it is linked by more
vertices or the vertices linking to it have higher importance themselves.
Despite originating from search engine, PageRank has been applied in many
fields such as social network analysis, chemistry, molecular biology and so
on[2, 3, 4].

In the past decades, plenty of PageRank algorithms had been proposed, among
which Power method[5] and Monte Carlo(MC) method[6] were two main approaches.
Existing algorithms[7, 8, 9, 10] have advantages respectively, however, most
of them do not distinguish the graph under investigation is directed or not.
For example, Power method treats any graph as directed, while undirected graph
widely exists in applications such as power grids, traffic nets and so on.
Compared with directed graph, the most significant property of undirected
graph is the symmetry on edges. Hence, simply treating undirected graph as
directed omitted this property. A few algorithms based on the MC method
benefited from symmetry. For example, Sarma[11] accelerated MC method on
undirected graph by borrowing ideas from distributed random walk[12], Luo[13,
14] proposed Radar Push on undirected graph and decreased the variance of
result. Thus, taking advantage of symmetry is a feasible solution of designing
PageRank algorithm for undirected graph.

PageRank vector can be obtained by computing $(\bm{I}-c\bm{P})^{-1}\bm{p}$,
where $\bm{I}$ is identity matrix, $\bm{P}$ is the possibility transition
matrix of the graph and $\bm{p}$ is the personalized vector. The state-of-the-
art PageRank algorithm, i.e., Forward Push(FP)[15], approximates
$(\bm{I}-c\bm{P})^{-1}\bm{p}$ by $\sum\limits_{i=0}^{k}(c\bm{P})^{i}\bm{p}$.
From the perspective of polynomial approximation,
$\sum\limits_{i=0}^{k}(cx)^{i}$ might not be the optimal approximating
polynomial of $(1-cx)^{-1}$ as polynomial with higher degree is required to
obtain an approximation with less error. Thus, constructing the optimal
polynomial is a natural idea.

Chebyshev Polynomial[16] is a set of orthogonal polynomials with good
properties, by which the optimal uniform approximating polynomial of any real,
continuous and square integrable function can be constructed. Directly
extending Chebyshev Polynomial approximation from the field of real function
to the field of matrix function is always infeasible. However, if $\bm{p}$ can
be linearly represented by $\bm{P}$’s eigenvectors, then
$(\bm{I}-c\bm{P})^{-1}\bm{p}$ can be converted to the linear combination of
$\bm{P}$’s eigenvectors, where the coefficients are functions of the
corresponding eigenvalues. By applying the Chebyshev Polynomial approximation
in the coefficients, a new approach of computing PageRank can be obtained. To
this end, two conditions required that (1)$\bm{P}$ has enough linearly
independent eigenvectors and (2)$\bm{P}$’s eigenvalues are all real. The first
condition is sufficed by the density of diagonalisable matrix. For undirected
graph, the second condition can be ensured by the symmetry. Motivated by this,
a parallel PageRank algorithm specially for undirected graph is proposed in
this paper. The contributions are as below.

  1. (1)

We present a solution of computing PageRank via the Chebyshev Polynomial
approximation. Based on the recursion relation of Chebyshev Polynomial,
PageRank vector can be calculated iteratively.

  2. (2)

We proposes a parallel PageRank algorithm, i.e., the Chebyshev Polynomial
Approximation Algorithm(CPAA). Theoretically, CPAA has higher convergence rate
and less computation amount compared with the Power method.

  3. (3)

Experimental results on six data sets suggest that CPAA markedly outperforms
the Power method, when damping factor $c=0.85$, CPAA with 38 parallelism
requires only 60% iteration rounds to get converged and can be at most 39
times faster.

The remaining of this paper are as follows. In section 2, PageRank and
Chebyshev Polynomial are introduced briefly. In section 3, the method of
computing PageRank via Chebyshev Polynomial approximation is elaborated.
Section 4 proposes CPAA and the theoretical analysis as well. Experiments on
six data sets are conducted in section 5. A summarization of this paper is
presented in Section 6.

##  2 Preliminary

In this section, PageRank and Chebyshev Polynomial are introduced briefly.

###  2.1 PageRank

Given graph $G=(V,E)$, where $V=\\{v_{1},v_{2},\cdots,v_{n}\\}$,
$E=\\{(v_{i},v_{j}):i,j=1,2,\cdots,n\\}$ and $|E|=m$. Let
$\bm{A}=(a_{ij})_{n\times n}$ denote the adjacency matrix where $a_{ij}$ is 1
if $(v_{j},v_{i})\in E$ and 0 else. Let $\bm{P}=(p_{ij})_{n\times n}$ denote
the probability transition matrix where

$p_{ij}=\left\\{\begin{array}[]{ll}a_{ij}/\sum\limits_{i=1}^{n}{a_{ij}},&if\quad\sum\limits_{i=1}^{n}{a_{ij}\neq
0},\\\ 0,&else.\end{array}\right.$

Let $\bm{d}=(d_{1},d_{2},...,d_{n})^{T}$ where $d_{i}$ is 1 if $v_{i}$ is
dangling vertex and 0 else. Let $\bm{p}=(p_{1},p_{2},\cdots,p_{n})^{T}$ denote
the $n$-dimensional probability distribution. Let
$\bm{P}^{{}^{\prime}}=\bm{P}+\bm{p}\bm{d}^{T}$ and
$\bm{P}^{{}^{\prime\prime}}=c\bm{P}^{{}^{\prime}}+(1-c)\bm{p}\bm{e}^{T}$,
where $c\in(0,1)$ is damping factor and $\bm{e}=(1,1,...,1)^{T}$. Denote by
$\bm{\pi}$ the PageRank vector. According to[17], let
$\bm{p}=\frac{\bm{e}}{n}$, one of PageRank’s definitions is

$\small\bm{\pi}=\bm{P}^{{}^{\prime\prime}}\bm{\pi},\pi_{i}>0,\sum\limits_{i=1}^{n}\pi_{i}=1.$
(1)

Based on random walk model, a random walk starts at random vertex, with
probability $c$ walks according to graph and with probability $1-c$
terminates, then $\bm{\pi}$ is the probability distribution of random walk
terminating at each vertex.

If $G(V,E)$ is undirected graph, then $\bm{d}=\bm{0}$ and
$\bm{P}=\bm{P}^{{}^{\prime}}$. Let $\bm{D}=diag\\{d_{1},d_{2},\cdots,d_{n}\\}$
where $d_{i}=\sum\limits_{j=1}^{n}{a_{ij}}$, then $d_{i}>0$ and
$\bm{P}=\bm{A}\bm{D}^{-1}$. From Formula (1), we have

$(\bm{I}-c\bm{P})\bm{\pi}=(1-c)\bm{p}$.

Since $||\bm{P}||_{2}=1$, $\bm{I}-c\bm{P}$ is invertible, thus

$\bm{\pi}=(1-c)(\bm{I}-c\bm{P})^{-1}\bm{p}$.

In fact, $(1-c)(\bm{I}-c\bm{P})^{-1}\bm{p}$ is the algebraic form of Forward
Push. With initial mass distribution $\bm{p}$, each vertex reserves $1-c$
proportion of mass receiving from its source vertices and evenly pushes the
remaining $c$ proportion to its target vertices, then $\bm{\pi}$ is the final
mass distribution with

$\small\pi_{i}=\frac{\bm{e}^{T}_{i}(\bm{I}-c\bm{P})^{-1}\bm{p}}{\bm{e}^{T}(\bm{I}-c\bm{P})^{-1}\bm{p}},$
(2)

where $\bm{e}_{i}$ is $n$-dimensional vector with the $i_{th}$ element is 1
and 0 others.

###  2.2 Chebyshev Polynomial

Chebyshev Polynomial is a set of polynomials defined as

$T_{n}(x)=\cos{(n\arccos{x})}$, $n=0,1,2,\cdots$.

$T_{n}(x)$ satisfies:

  1. $(1)$

$|T_{n}(x)|\leq 1$.

  2. $(2)$

$\frac{2}{\pi}\int_{-1}^{1}\frac{T_{m}(x)T_{n}(x)}{\sqrt{1-x^{2}}}dx=\left\\{\begin{array}[]{ll}2,&if\quad
m=n=0,\\\ \delta_{mn},&else.\end{array}\right.$

  3. $(3)$

$T_{n+1}(x)=2xT_{n}(x)-T_{n-1}(x)$.

Chebyshev Polynomial is usually utilized in constructing the optimal uniform
approximating polynomial. For any given function $f(x)\in C(a,b)$ and
$\int_{a}^{b}f^{2}(x)dx<\infty$, it follows that

$f(x)=\frac{c_{0}}{2}+\sum\limits_{k=1}^{\infty}c_{k}T_{k}(x)$,

where
$c_{k}=\frac{2}{b-a}\int_{a}^{b}\frac{1}{\sqrt{1-(\frac{2}{b-a}x-\frac{a+b}{b-a})^{2}}}f(x)T_{k}(\frac{2}{b-a}x-\frac{a+b}{b-a})dx$,
$\frac{c_{0}}{2}+\sum\limits_{k=1}^{n}c_{k}T_{k}(x)$ is the optimal uniform
approximating polynomial of $f(x)$. For

$f(x)=(1-cx)^{-1},x\in(-1,1)$,

where $c\in(0,1)$ is constant, we have
$c_{k}=\frac{2}{\pi}\int_{0}^{\pi}\frac{\cos{(kt)}}{1-c\cos{t}}dt$,
$\\{c_{k}:k=0,1,2,\cdots\\}$ decreases monotonously to 0, and

$\lim\limits_{M\to\infty}\sum\limits_{k=M}^{\infty}c_{k}=\lim\limits_{M\to\infty}\frac{2}{\pi}\int_{0}^{\pi}\frac{\sum\limits_{k=M}^{\infty}\cos{(kt)}}{1-c\cos{t}}dt=0$.

Thus, $\forall\epsilon>0,\exists M_{1}\in\mathbb{N}^{+},s.t.$, $\forall
M>M_{1}$,

$\small\max\limits_{x\in(-1,1)}|f(x)-(\frac{c_{0}}{2}+\sum\limits_{k=1}^{M}c_{k}T_{k}(x))|<\epsilon.$
(3)

##  3 Computing PageRank by Chebyshev Polynomial approximation

In this section, Formula (3) is extended from the field of real function to
the filed of matrix function. We mainly prove that, $\forall\epsilon>0$,
$\exists M_{1}\in\mathbb{N}^{+}$, $s.t.$, $\forall M>M_{1}$,

$\small||(\bm{I}-c\bm{P})^{-1}\bm{p}-(\frac{c_{0}}{2}\bm{p}+\sum\limits_{k=1}^{M}c_{k}T_{k}(\bm{P})\bm{p})||_{2}<\epsilon.$
(4)

To this end, we firstly give two lemmas as below.

######  Lemma 1.

[18] Diagonalizable matrices are Zariski dense.

Lemma 1 implies that, for any given matrix $\bm{P}\in\mathbb{R}^{n\times n}$,
there exists diagonalizable matrix sequence
$\left\\{\bm{P}_{t}:t=0,1,2,\cdots\right\\}$, where
$\bm{P}_{t}\in\mathbb{R}^{n\times n}$ and $||\bm{P}_{t}||_{2}=1$, $s.t.$,
$\bm{P}_{t}\rightrightarrows\bm{P}$. Thus, we have

$\lim\limits_{t\to\infty}||\bm{P}_{t}-\bm{P}||_{2}=0$.

Since both $\sum\limits_{k=0}^{\infty}(c\bm{P}_{t})^{k}$ and
$\sum\limits_{k=1}^{M}c_{k}T_{k}(\bm{P}_{t})$ are uniformly convergent with
respect to $t$, $\forall\epsilon>0$, $\exists N_{1}\in\mathbb{N}^{+}$, $s.t.$,
$\forall t>N_{1}$,

$||\sum\limits_{k=0}^{\infty}(c\bm{P})^{k}\bm{p}-\sum\limits_{k=0}^{\infty}(c\bm{P}_{t})^{k}\bm{p}||_{2}<\epsilon$,

and

$||\sum\limits_{k=1}^{M}c_{k}T_{k}(\bm{P})\bm{p}-\sum\limits_{k=1}^{M}c_{k}T_{k}(\bm{P}_{t})\bm{p}||_{2}<\epsilon$.

Since
$(\bm{I}-c\bm{P})^{-1}\bm{p}=\sum\limits_{k=0}^{\infty}(c\bm{P})^{k}\bm{p}$,
Formula (4) can be hold by proving

$\small\lim\limits_{t\to\infty}||\sum\limits_{k=0}^{\infty}(c\bm{P}_{t})^{k}\bm{p}-(\frac{c_{0}}{2}\bm{p}+\sum\limits_{k=1}^{M}c_{k}T_{k}(\bm{P}_{t})\bm{p})||_{2}<\epsilon.$
(5)

######  Lemma 2.

The whole eigenvalues of $\bm{P}$ are real.

###### Proof.

We first prove that $\bm{D}-c\bm{A}$ is positive definite matrix.
$\forall\bm{x}\in\mathbb{R}^{n}$ and $\bm{x}\neq\bm{0}$, it follows that

$\begin{aligned}
\bm{x}^{T}(\bm{D}-c\bm{A})\bm{x}=&\sum\limits_{i=1}^{n}{d_{i}x_{i}^{2}}-c\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}{a_{ij}x_{i}x_{j}}\\\
=&\frac{1}{2}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}{a_{ij}(x_{i}^{2}+x_{j}^{2}-2cx_{i}x_{j})}.\end{aligned}$

Since $x_{i}^{2}+x_{j}^{2}-2cx_{i}x_{j}>0,i,j=1,2,\cdots,n$, $\bm{D}-c\bm{A}$
is positive definite.

Let $\bm{B}=\bm{D}-c\bm{A}$, then

$(\bm{B}\bm{P}^{T})^{T}=(\bm{A}-c\bm{A}\bm{D}^{-1}\bm{A})^{T}=\bm{A}-c\bm{A}\bm{D}^{-1}\bm{A}=\bm{B}\bm{P}^{T}$.

$\forall\bm{y}\in\mathbb{C}^{n}$ and $\bm{y}\neq\bm{0}$, assuming
$\bm{y}=\bm{y}_{1}+\mathrm{i}\bm{y}_{2}$, where
$\bm{y}_{1},\bm{y}_{2}\in\mathbb{R}^{n}$, we have

$\begin{aligned}
\bm{y}^{*}\bm{B}\bm{y}&=(\bm{y}_{1}-\mathrm{i}\bm{y}_{2})^{T}\bm{B}(\bm{y}_{1}+\mathrm{i}\bm{y}_{2})\\\
&=\bm{y}_{1}^{T}\bm{B}\bm{y}_{1}+\bm{y}_{2}^{T}\bm{B}\bm{y}_{2}+\mathrm{i}(\bm{y}_{1}^{T}\bm{B}\bm{y}_{2}-\bm{y}_{2}^{T}\bm{B}\bm{y}_{1})\\\
&=\bm{y}_{1}^{T}\bm{B}\bm{y}_{1}+\bm{y}_{2}^{T}\bm{B}\bm{y}_{2}>0.\end{aligned}$

Denote by $\lambda$ the eigenvalue of $\bm{P}^{T}$, by $\bm{\chi}$ the
corresponding eigenvector, then

$\lambda\bm{\chi}^{*}\bm{B}\bm{\chi}=\bm{\chi}^{*}\bm{B}\bm{P}^{T}\bm{\chi}=\bm{\chi}^{*}(\bm{B}\bm{P}^{T})^{T}\bm{\chi}=(\bm{P}^{T}\overline{\bm{\chi}})^{T}\bm{B}\bm{\chi}=\overline{\lambda}\bm{\chi}^{*}\bm{B}\bm{\chi}.$

Since $\bm{\chi}^{*}\bm{B}\bm{\chi}>0$, we have $\overline{\lambda}=\lambda$,
i.e., $\lambda$ is real. ∎

Since $\bm{P}_{t}$ is diagonalizable, there exist $n$ linearly independent
eigenvectors. Denote by $\bm{\chi}_{l}(t)$ the eigenvector of $\bm{P}_{t}$, by
$\lambda_{l}(t)=\lambda_{l}^{(1)}(t)+\mathrm{i}\lambda_{l}^{(2)}(t)$ the
corresponding eigenvalue, where $||\bm{\chi}_{l}(t)||_{2}=1$,
$\lambda_{l}^{(1)}(t),\lambda_{l}^{(2)}(t)\in\mathbb{R}$ and $l=1,2,\cdots,n$,
we have $|\lambda_{l}(t)|\leq 1$ and

$||\bm{P}\bm{\chi}_{l}(t)-\lambda_{l}(t)\bm{\chi}_{l}(t)||_{2}=||\bm{P}\bm{\chi}_{l}(t)-\bm{P}_{t}\bm{\chi}_{l}(t)||_{2}\leq||\bm{P}-\bm{P}_{t}||_{2}$.

Given $l$, since $\mathbb{C}$ is complete and $\\{\lambda_{l}(t)\\}$ is
bounded, there exists convergent subsequence
$\\{\lambda_{l}(t_{k})\\}\subseteq\\{\lambda_{l}(t)\\}$, let
$\lambda_{l}=\lim\limits_{k\to\infty}\lambda_{l}(t_{k})$. Similarly, since
$\mathbb{R}^{n}$ is complete and $\\{\bm{\chi}_{l}(t_{k})\\}$ is bounded,
there exists convergent subsequence
$\\{\bm{\chi}_{l}(t_{k_{s}})\\}\subseteq\\{\bm{\chi}_{l}(t_{k})\\}$, let
$\bm{\chi}_{l}=\lim\limits_{s\to\infty}\bm{\chi}_{l}(t_{k_{s}})$. For
convenience, we still mark $t_{k_{s}}$ as $t$ in the following. Since

$\lim\limits_{t\to\infty}||\bm{P}\bm{\chi}_{l}(t)-\lambda_{l}(t)\bm{\chi}_{l}(t)||_{2}\leq\lim\limits_{t\to\infty}||\bm{P}-\bm{P}_{t}||_{2}=0$,

it follows that

$\bm{P}\bm{\chi}_{l}=\lim\limits_{t\to\infty}\bm{P}\bm{\chi}_{l}(t)=\lim\limits_{t\to\infty}\lambda_{l}(t)\bm{\chi}_{l}(t)=\lambda_{l}\bm{\chi}_{l}$,

i.e, $\lambda_{l}$ is the eigenvalue of $\bm{P}$ and $\bm{\chi}_{l}$ is the
corresponding eigenvector. According to Lemma 2, $\lambda_{l}$ is real, thus
for arbitrary $l=0,1,\cdots,n$,

$\lim\limits_{t\to\infty}\lambda_{l}^{(2)}(t)=0.$ (6)

For any given $t$, $\\{\bm{\chi}_{l}(t):l=1,2,\cdots,n\\}$ is basis of
$\mathbb{R}^{n}$, by which $\bm{p}$ can be linearly represented. Assuming

$\bm{p}=\sum\limits_{l=1}^{n}\alpha_{l}(t)\bm{\chi}_{l}(t)$,

where $\alpha_{l}(t)\in\mathbb{R}$ and $|\alpha_{l}(t)|<1$, then

$\sum\limits_{k=0}^{\infty}(c\bm{P}_{t})^{k}\bm{p}=\sum\limits_{k=0}^{\infty}(c\bm{P}_{t})^{k}\sum\limits_{l=1}^{n}\alpha_{l}(t)\bm{\chi}_{l}(t)=\sum\limits_{l=1}^{n}\alpha_{l}(t)\sum\limits_{k=0}^{\infty}(c\lambda_{l}(t))^{k}\bm{\chi}_{l}(t)$.

Since both $\sum\limits_{k=0}^{\infty}(c\lambda_{l}(t))^{k}$ and
$\sum\limits_{k=0}^{\infty}(c\lambda_{l}^{(1)}(t))^{k}$ are uniformly
convergent with respect to $t$, we have

$\begin{aligned}
&||\sum\limits_{l=1}^{n}\alpha_{l}(t)\sum\limits_{k=0}^{\infty}(c\lambda_{l}(t))^{k}\bm{\chi}_{l}(t)-\sum\limits_{l=1}^{n}\alpha_{l}(t)\sum\limits_{k=0}^{\infty}(c\lambda^{(1)}_{l}(t))^{k}\bm{\chi}_{l}(t)||_{2}\\\
=&|\sum\limits_{l=1}^{n}\alpha_{l}(t)(\sum\limits_{k=0}^{\infty}(c\lambda_{l}(t))^{k}-\sum\limits_{k=0}^{\infty}(c\lambda^{(1)}_{l}(t))^{k})|\cdot||\bm{\chi}_{l}(t)||_{2}\\\
=&|\sum\limits_{l=1}^{n}\alpha_{l}(t)(\sum\limits_{k=0}^{\infty}(c(\lambda^{(1)}_{l}(t)+\mathrm{i}\lambda^{(2)}_{l}(t)))^{k}-\sum\limits_{k=0}^{\infty}(c\lambda^{(1)}_{l}(t))^{k})|,\end{aligned}$

from Formula (6), $\forall\epsilon>0$, $\exists N_{2}\in\mathbb{N}^{+}$,
$s.t$, $\forall t>N_{2}$,

$||\sum\limits_{l=1}^{n}\alpha_{l}(t)\sum\limits_{k=0}^{\infty}(c\lambda_{l}(t))^{k}\bm{\chi}_{l}(t)-\sum\limits_{l=1}^{n}\alpha_{l}(t)\sum\limits_{k=0}^{\infty}(c\lambda^{(1)}_{l}(t))^{k}\bm{\chi}_{l}(t)||_{2}<\epsilon$,

i.e.,

$||\sum\limits_{k=0}^{\infty}(c\bm{P}_{t})^{k}\bm{p}-\sum\limits_{l=1}^{n}\alpha_{l}(t)\sum\limits_{k=0}^{\infty}(c\lambda^{(1)}_{l}(t))^{k}\bm{\chi}_{l}(t)||_{2}<\epsilon$.

Since
$(1-c\lambda^{(1)}_{l}(t))^{-1}=\sum\limits_{k=0}^{\infty}(c\lambda^{(1)}_{l}(t))^{k}$,
it follows that

$\small||\sum\limits_{k=0}^{\infty}(c\bm{P}_{t})^{k}\bm{p}-\sum\limits_{l=1}^{n}\alpha_{l}(t)(1-c\lambda^{(1)}_{l}(t))^{-1}\bm{\chi}_{l}(t)||_{2}<\epsilon.$
(7)

Approximating $f(\lambda^{(1)}_{l}(t))=(1-c\lambda^{(1)}_{l}(t))^{-1}$ by
Chebyshev Polynomial, from Formula (3), $\forall\epsilon>0$, $\exists
M_{1}\in\mathbb{N}^{+}$, $s.t.$, $\forall M>M_{1}$,

$||\sum\limits_{l=1}^{n}\alpha_{l}(t)(1-c\lambda_{l}^{(1)}(t))^{-1}\bm{\chi}_{l}(t)-\sum\limits_{l=1}^{n}\alpha_{l}(t)(\frac{c_{0}}{2}+\sum\limits_{k=1}^{M}c_{k}T_{k}(\lambda_{l}^{(1)}(t)))\bm{\chi}_{l}(t)||_{2}<\epsilon$,

i.e.,

$||\sum\limits_{l=1}^{n}\alpha_{l}(t)(1-c\lambda_{l}^{(1)}(t))^{-1}\bm{\chi}_{l}(t)-(\frac{c_{0}}{2}\bm{p}+\sum\limits_{k=1}^{M}c_{k}\sum\limits_{l=1}^{n}T_{k}(\lambda_{l}^{(1)}(t))\alpha_{l}(t)\bm{\chi}_{l}(t))||_{2}<\epsilon$.

$T_{k}(x)$ is a polynomial of degree $k$, assuming

$T_{k}(x)=\sum\limits_{i=0}^{k}a_{i}x^{i}$,

where the coefficient $a_{i}\in\mathbb{R}$, we have

$\begin{aligned}
T_{k}(\bm{P}_{t})\bm{p}&=\sum\limits_{i=0}^{k}a_{i}\bm{P}^{i}_{t}\bm{p}\\\
&=\sum\limits_{i=0}^{k}a_{i}\bm{P}^{i}_{t}\sum\limits_{l=1}^{n}\alpha_{l}(t)\bm{\chi}_{l}(t)\\\
&=\sum\limits_{l=1}^{n}\alpha_{l}(t)\sum\limits_{i=0}^{k}a_{i}\bm{P}^{i}_{t}\bm{\chi}_{l}(t)\\\
&=\sum\limits_{l=1}^{n}\alpha_{l}(t)\sum\limits_{i=0}^{k}a_{i}(\lambda^{(1)}_{l}(t)+\mathrm{i}\lambda^{(2)}_{l}(t))^{i}\bm{\chi}_{l}(t).\end{aligned}$

From Formula (6), $\forall k$ and $\forall\epsilon>0$, $\exists
N_{3}\in\mathbb{N}^{+}$, $s.t.$, $\forall t>N_{3}$,

$||\sum\limits_{l=1}^{n}\alpha_{l}(t)\sum\limits_{i=0}^{k}a_{i}(\lambda^{(1)}_{l}(t)+\mathrm{i}\lambda^{(2)}_{l}(t))^{i}\bm{\chi}_{l}(t)-\sum\limits_{l=1}^{n}\alpha_{l}(t)\sum\limits_{i=0}^{k}a_{i}(\lambda_{l}^{(1)}(t))^{i}\bm{\chi}_{l}(t)||_{2}<\epsilon$,

i.e.,

$||T_{k}(\bm{P}_{t})\bm{p}-\sum\limits_{l=1}^{n}T_{k}(\lambda_{l}^{(1)}(t))\alpha_{l}(t)\bm{\chi}_{l}(t)||_{2}<\epsilon$.

Thus, $\forall M$ and $\forall\epsilon>0$, $\exists N_{4}\in\mathbb{N}^{+}$,
$s.t.$, $\forall t>N_{4}$,

$||\sum\limits_{l=1}^{n}\alpha_{l}(t)(1-c\lambda_{l}^{(1)}(t))^{-1}\bm{\chi}_{l}(t)-(\frac{c_{0}}{2}\bm{p}+\sum\limits_{k=1}^{M}c_{k}T_{k}(\bm{P}_{t})\bm{p})||_{2}<\epsilon$.

From Formula (7), $\forall\epsilon>0$, $\exists M_{1}\in\mathbb{N}^{+}$,
$s.t.$, $\forall M>M_{1}$, $\exists N\in\mathbb{N}^{+}$, $s.t.$, $\forall
t>N$,

$||\sum\limits_{k=0}^{\infty}(c\bm{P}_{t})^{k}\bm{p}-(\frac{c_{0}}{2}\bm{p}+\sum\limits_{k=1}^{M}c_{k}T_{k}(\bm{P}_{t})\bm{p})||_{2}<\epsilon$,

i.e.,

$\lim\limits_{t\to\infty}||\sum\limits_{k=0}^{\infty}(c\bm{P}_{t})^{k}\bm{p}-(\frac{c_{0}}{2}\bm{p}+\sum\limits_{k=1}^{M}c_{k}T_{k}(\bm{P}_{t})\bm{p})||_{2}=0$.

Thus Formula (5) is hold, so as Formula (4).

It should be noted that, Formula (4) requires the whole eigenvalues of
$\bm{P}$ are real, thus it might not be hold on some directed graphs. From the
perspective of polynomial approximation, FP approximates $f(x)=(1-cx)^{-1}$ by
taking $\\{1,x,x^{2},\cdots\\}$ as the basis, while Formula (4) approximates
$f(x)$ by taking Chebyshev Polynomial $\\{T_{k}(x):k=1,2,\cdots\\}$ as basis.
Since Chebyshev Polynomial are orthogonal, Formula (4) is more efficient.

##  4 Algorithm and analysis

In this section, a parallel PageRank algorithm for undirected graph is
proposed and then the theoretical analysis is given.

###  4.1 Algorithm designing

From Formula (4), PageRank vector can be obtained by calculating
$\\{c_{k}T_{k}(\bm{P})\bm{p}:k=0,1,2.\cdots\\}$:

  1. (1)

$\\{c_{k}:k=0,1,2,\cdots\\}$ relies on damping factor $c$ only, thus we can
calculate and store them beforehand.

  2. (2)

$\\{T_{k}(\bm{P})\bm{p}:k=0,1,2,\cdots\\}$ can be calculated iteratively by

$T_{k+1}(\bm{P})\bm{p}=2\bm{P}T_{k}(\bm{P})\bm{p}-T_{k-1}(\bm{P})\bm{p}$,

it is sufficed for each vertex reserving the $(k-1)_{th}$ and $k_{th}$
iteration results to generate the $(k+1)_{th}$ iteration result.

  3. (3)

During the same iteration round, each vertex calculates independently, thus
the computing can be done in parallel.

Restricted by the computing resource, assigning exclusive thread for each
vertex is in feasible, we can invoke $K$ independent threads and assigning
vertices to them. The Chebyshev Polynomial Approximation Algorithm(CPAA) is
given as Algorithm 1.

Algorithm 1 Chebyshev Polynomial Approximation Algorithm

1:

2:$K$: The parallelism.

3:$M$: The upper bound of iteration rounds.

4:$\\{c_{0},c_{1},\cdots,c_{M}\\}$: The coefficient of Chebyshev Polynomial
approximation.

5:$\bm{\pi}$: PageRank vector.

6:Each vertex $v_{i}$ maintains
$T_{i},T_{i}^{{}^{\prime}},T_{i}^{{}^{\prime\prime}}$ and
$\overline{\pi}_{i}$.

7:Assign vertices to $K$ calculating threads, denote by $S_{j}$ the set of
vertices belonging to thread $j$.

8:Initially set $T_{i}=1$ and $\overline{\pi}_{i}=\frac{c_{0}}{2}T_{i}$.

9:

10:for $k=1;k\leq M;k++$ do

11: for all $j$ do $\triangleright$ [Do in parallel.]

12: if $k=1$ then

13: for $u\in S_{j}$ do

14: $T_{u}^{{}^{\prime}}=0$;

15: for $v_{i}\in N(u)$ do$\triangleright$ [$N(u)$ is the set of $u$’s
adjacency vertices.]

16: $T_{u}^{{}^{\prime}}=T_{u}^{{}^{\prime}}+\frac{T_{i}}{deg(v_{i})}$;

17: end for

18: $\overline{\pi}_{u}=\overline{\pi}_{u}+c_{1}*T_{u}^{{}^{\prime}}$;

19: end for

20: else

21: for $u\in S_{j}$ do

22: $T^{{}^{\prime\prime}}_{u}=0$;

23: for $v_{i}\in N(u)$ do

24:
$T_{u}^{{}^{\prime\prime}}=T_{u}^{{}^{\prime\prime}}+\frac{T_{i}^{{}^{\prime}}}{deg(v_{i})}$;

25: end for

26: $T_{u}^{{}^{\prime\prime}}=2T_{u}^{{}^{\prime\prime}}-T_{u}$;

27: $\overline{\pi}_{i}=\overline{\pi}_{i}+c_{k}*T_{u}^{{}^{\prime\prime}}$;

28: end for

29: end if

30: end for

31: for all $j$ do $\triangleright$ [Do in parallel.]

32: for $u\in S_{j}$ do

33: $T_{u}=T_{u}^{{}^{\prime}}$;

34: $T_{u}^{{}^{\prime}}=T_{u}^{{}^{\prime\prime}}$;

35: end for

36: end for

37:end for

38:Calculate PageRank vector $\bm{\pi}$ by
$\pi_{i}=\frac{\overline{\pi}_{i}}{\sum\limits_{i=1}^{n}\overline{\pi}_{i}}$.

From the perspective of vertex, CPAA is a loop of two stages, i.e., mass
generating and mass accumulating. Initially, each vertex holds 1 mass and 0
accumulating mass. During the $k_{th}$ iteration, at generating stage, each
vertex generates the $(k+1)_{th}$ iteration mass by reading the $k_{th}$
iteration mass from its adjacency vertices and subtracting its own
$(k-1)_{th}$ iteration mass; at accumulating stage, each vertex adds the
$(k+1)_{th}$ iteration mass weighted by $c_{k+1}$ to the accumulating mass. In
the end, the proportion of accumulating mass of one vertex is its PageRank
value. Although the mass each vertex holding may change at each iteration, the
total mass of the graph is constant at $n$.

From the whole, the total accumulating mass $S$ of the graph is constant and
falls into two parts, i.e., the accumulated and unaccumulated. CPAA is a
process that converting the unaccumulated mass to the accumulated mass.
Initially, the accumulated mass is 0 and the unaccumulated mass is $S$. While
CPAA running, the accumulated mass increases and the unaccumulated mass
decreases, specifically, at the $k_{th}$ iteration, the accumulated mass
increases by $c_{k}n$ and the unaccumulated mass decreases by the same amount.
When the accumulated mass is $S$ and the unaccumulated mass is 0, CPAA gets
converged. The distribution of accumulated mass is the PageRank vector.

It should be noted that, the dependency between the $k_{th}$ iteration and
$(k+1)_{th}$ iteration can not be eliminated, paralleling only exists in the
same iteration round.

###  4.2 Algorithm analysis

We analysis CPAA from three aspects, convergence rate, error and computation
amount.

####  4.2.1 Convergence rate

Convergence rate reflects the rapidity of unaccumulated mass decreasing to 0.
Let

$S_{k}=\left\\{\begin{array}[]{ll}n\frac{c_{0}}{2},&if\quad k=0,\\\
n\frac{c_{0}}{2}+n\sum\limits_{i=1}^{k}c_{i},&else.\end{array}\right.$

$S_{k}$ is accumulated mass at the end of the $k_{th}$ iteration, and
$S=\lim\limits_{k\to\infty}S_{k}$. Let

$\sigma_{k}=\frac{S-S_{k}}{S-S_{k-1}},k=1,2,3,\cdots$.

Then $1-\sigma_{k}$ is the ratio of unaccumulated mass decreasing at the
$k_{th}$ iteration.

######  Proposition 1.

$\sigma_{k}=\frac{c^{2}-(2-c)(1-\sqrt{1-c^{2}})}{c^{2}-c(1-\sqrt{1-c^{2}})}$,
$c\in(0,1)$, $k=1,2,3,\cdots$.

###### Proof.

Firstly, let $B_{k}=\sum\limits_{i=k}^{\infty}c_{k}$, then $B_{k}$ is
absolutely convergent, and

$1-\sigma_{k}=\frac{c_{k}}{B_{k}}$.

since $c_{k}=\frac{2}{\pi}\int_{0}^{\pi}\frac{\cos{kt}}{1-c\cos{t}}dt$, it
follows the recurrence relation

$c_{k-1}+c_{k+1}=\frac{2}{\pi}\int_{0}^{\pi}\frac{2\cos{kt}\cos{t}}{1-c\cos{t}}dt=\frac{2}{c}c_{k}$.

According to

$\begin{aligned} &c_{k-1}+c_{k+1}=\frac{2}{c}c_{k},\\\
&c_{k}+c_{k+2}=\frac{2}{c}c_{k+1},\\\ &\cdots,\end{aligned}$

we have

$B_{k}+c_{k-1}+B_{k}-c_{k}=\frac{2}{c}B_{k}$,

thus $B_{k}=\frac{c_{k-1}-c_{k}}{\frac{2}{c}-2}$ and
$\frac{c_{k}}{B_{k}}=\frac{\frac{2}{c}-2}{\frac{c_{k-1}}{c_{k}}-1}$.

Since $c_{0}=\frac{2}{\sqrt{1-c^{2}}}$,
$c_{1}=\frac{2}{c}(\frac{1-\sqrt{1-c^{2}}}{\sqrt{1-c^{2}}})$,
$c_{2}=\frac{2}{c^{2}}\frac{2(1-\sqrt{1-c^{2}})-c^{2}}{\sqrt{1-c^{2}}}$, it
follows that

$\frac{c_{0}}{c_{1}}=\frac{c_{1}}{c_{2}}=\frac{c}{1-\sqrt{1-c^{2}}}$.

According to $\frac{c_{0}}{c_{1}}+\frac{c_{2}}{c_{1}}=\frac{2}{c}$, we have

$\frac{c_{0}}{c_{1}}=\frac{c_{1}}{c_{2}}=\cdots=\frac{c_{k-1}}{c_{k}}=\frac{c}{1-\sqrt{1-c^{2}}}$.

Thus,
$\sigma_{k}=\frac{c^{2}-(2-c)(1-\sqrt{1-c^{2}})}{c^{2}-c(1-\sqrt{1-c^{2}})}$.
∎

Proposition(1) shows that $\sigma_{k}$ relates with $c$ only. For any given
$c$, $\sigma_{k}$ is a constant. For convenience, we mark $\sigma_{k}$ as
$\sigma_{c}$ in the following. While CPAA running, the unaccumulated mass
decreases by $1-\sigma_{c}$ proportion per iteration. Referring to the Power
method, $\sigma_{c}$ indicates the convergence rate of CPAA. Since

$\frac{\sigma_{c}}{c}=\frac{c-(\frac{2}{c}-1)(1-\sqrt{1-c^{2}})}{c^{2}-c(1-\sqrt{1-c^{2}})}<1$,
$c\in(0,1)$,

CPAA has higher convergence rate compared with the Power method. Figure 1
shows the relation between $\sigma_{c}$ and $c$. When $c=0.85$,
$\sigma_{c}=0.5567$, which implies that CPAA only takes 65% iteration rounds
of Power method to get converged.

Figure 1: Convergence rate

####  4.2.2 Error

It is hard to estimate the error of each vertex, we discuss the relation
between the relative error and the upper bound of iteration rounds from the
whole. Denote by

$ERR_{M}=1-\frac{S_{M}}{S}$

the relative error, then we have

$S_{M}=\frac{c_{0}}{2}+\sum\limits_{k=1}^{M}c_{k}=\frac{c_{0}}{2}+\frac{1-\beta^{M}}{1-\beta}\beta
c_{0}$,

thus,

$\small ERR_{M}=1-\frac{S_{M}}{S}=\frac{2\beta^{M+1}}{1+\beta},$ (8)

where $\beta=\frac{1-\sqrt{1-c^{2}}}{c}$. Figure 2 shows that $ERR_{M}$
decreases exponentially with respect to $M$, when $c=0.85$, $ERR_{M}$ can be
less than $10^{-4}$ within 20 iteration rounds. It should be noted that the
error estimation above mentioned is very rough as taking no graph structure
into consideration.

Figure 2: Relative error

####  4.2.3 Computation amount

While CPAA running, for vertex $v$, there are $deg(v)$ multiplications and
$deg(v)$ additions at generating stage, and $1$ addition at accumulating
stage. Therefore, the computation consists of $m$ multiplications and $m+2n$
additions per iteration. The total computation are $Mm$ multiplications and
$M(m+2n)$ additions.

##  5 Experiments

In this section, we experimentally demonstrate CPAA. Firstly, the hardware and
measurements are introduced briefly. Then, the convergence of CPAA is
illustrated. At last, the comparison of CPAA with other algorithms is
elaborated.

###  5.1 Experimental setting

All algorithms in this experiment are implemented with C++ on serves with
Intel(R) Xeon(R) Silver 4210R CPU 2.40GHz 40 processors and 190GB memory. The
operation system is Ubuntu 18.04.5 LTS. Six data sets are illustrated in table
1, where $n$, $m$ and $deg=\frac{m}{n}$ represent the number of vertices, the
number of edges and the average degree respectively. The CPU time consumption
$T$, iteration rounds $k$ and max relative error $ERR=\max\limits_{v_{i}\in
V}\frac{|\overline{\pi}_{i}-\pi_{i}|}{\pi_{i}}$ are taken to estimate the
algorithms, where the true PageRank value $\pi_{i}$ is obtained by Power
method at the $210_{th}$ iteration. The damping factor $c=0.85$.

Data sets | $n$ | $m$ | $deg$
---|---|---|---
NACA0015 | 1039183 | 6229636 | 5.99
delaunay-n21 | 2097152 | 12582816 | 6.0
M6 | 3501776 | 21003872 | 6.0
NLR | 4163763 | 24975952 | 6.0
CHANNEL | 4802000 | 85362744 | 17.78
kmer-V2 | 55042369 | 117217600 | 2.13
Table 1: Data sets

###  5.2 Convergence

The relation of $ERR$, $T$ and $k$ on six data sets is illustrated in Figure
3.

  1. (1)

The blue lines show that $ERR$ has a negative exponential relation with $k$,
which is in consistent with Formula (8). With the increasing of iteration
rounds, the result gets closer to the true value. The max relative error can
be less than $10^{-4}$ within 20 iterations.

  2. (2)

The red lines show that $T$ has a positive linear relation with $k$, thus, It
also implies that the computation of each iteration is constant. On kmer-V2,
the max relative error is lower than $10^{-2}$ within 15 seconds.

  3. (3)

The blue lines almost stay horizontal after 50 iteration rounds, which seems
CPAA has limitation in precision. We believe it is not a algorithm flaw, but
caused by the insufficiency of C++’s DOUBLE data type, whose significant digit
number is 15. The computer considers $1$ is equal to $1+10^{-16}$, so that,
$c_{k}$ with $k>50$ can not change the result.

Figure 3: $k$ versus $ERR$ and $T$

###  5.3 Comparison with other algorithms

We compare CPAA with Power method (SPI) [1, 19], parallel Power method (MPI)
[20] and an improvement of FP, i.e., IFP1 [21]. Both MPI, IFP1 and CPAA
execute with 38 parallelism. The relation of $ERR$ and $T$ are illustrated in
Figures 4. Table 2 shows iteration rounds and CPU time consumption when
$ERR<10^{-3}$.

  1. (1)

The green, cyan and magenta lines are lower than the blue lines, which
indicates that both MPI, IFP1 and CPAA are much faster than SPI. For example,
on kmer-V2, it takes near 900 seconds for SPI to get $ERR<10^{-3}$, while the
remaining algorithms cost less than 90 seconds. It implies that
parallelization is an effective solution to accelerate PageRank computing.

  2. (2)

The magenta lines are lower than the green and cyan lines, which implies that
CPAA outperforms MPI and IFP1. For example, on kmer-V2, it takes about 67
seconds for MPI to get $ERR<10^{-3}$, while CPAA costs less than 24 seconds.
Since the whole algorithms are executed with the same parallelism, we believe
the advantages are owe to the higher convergence rate of CPAA.

  3. (3)

On almost all six data sets, Power method takes 20 iterations to get
$ERR<10^{-3}$, while CPAA only takes 12 iterations, and can be at most 39
times faster than SPI. It implies that CPAA takes 60% iterations of Power
method to get converged, which is consistent with the theoretical analysis.

Figure 4: $T$ versus $ERR$ Data sets | SPI | MPI | ITA | CPAA
---|---|---|---|---
$k$ | $T$ | $k$ | $T$ | $k$ | $T$ | $k$ | $T$
NACA0015 | 20 | 13.541 | 20 | 0.887 | - | 1.390 | 12 | 0.604
delaunay-n21 | 20 | 20.175 | 20 | 1.854 | - | 2.921 | 12 | 1.111
M6 | 20 | 46.771 | 20 | 3.130 | - | 5.097 | 12 | 1.999
NLR | 20 | 57.876 | 20 | 3.467 | - | 5.983 | 12 | 2.379
CHANNEL | 20 | 102.828 | 20 | 6.819 | - | 9.689 | 12 | 5.845
kmer-V2 | 40 | 899.406 | 40 | 67.068 | - | 70.346 | 15 | 23.597
Table 2: Iterations rounds and CPU time consumption when $ERR<10^{-3}$

##  6 Conclusion

How to compute PageRank more efficiently is an attractive problem. In this
paper, a parallel PageRank algorithm specially for undirected graph called
CPAA is proposed via Chebyshev Polynomial approximation. Both the theoretical
analysis and experimental results indicate that CPAA has higher convergence
rate and less computation amount compared with the Power method. Moreover,
CPAA implies a more profound connotation, i.e., the symmetry of undirected
graph and the density of diagonalizable matrix. Based on the latter, the
problem of computing PageRank is converted to the problem of approximating
$f(x)=(1-cx)^{-1}$. This is a innovative approach for PageRank computing. In
the future, some other orthogonal polynomials—Laguerre polynomial, for
example—can be taken into consideration.

## References

  * [1] S. Brin, L. Page, Reprint of: The anatomy of a large-scale hypertextual web search engine, Computer networks 56 (18) (2012) 3825–3833.
  * [2] D. F. Gleich, Pagerank beyond the web, Siam Review 57 (3) (2015) 321–363.
  * [3] S. Brown, A pagerank model for player performance assessment in basketball, soccer and hockey, arXiv preprint arXiv:1704.00583 (2017).
  * [4] Ying, Hou, Zhou, Dou, Wang, Shao, The changes of central cultural cities based on the analysis of the agglomeration of literati’s footprints in tang and song dynasties, Journal of Geo-information Science 22 (5) (2020) 945–953.
  * [5] T. Haveliwala, S. Kamvar, D. Klein, C. Manning, G. Golub, Computing pagerank using power extrapolation, Tech. rep., Stanford (2003).
  * [6] K. Avrachenkov, N. Litvak, D. Nemirovsky, N. Osipova, Monte carlo methods in pagerank computation: When one iteration is sufficient, SIAM Journal on Numerical Analysis 45 (2) (2007) 890–904.
  * [7] S. Kamvar, T. Haveliwala, G. Golub, Adaptive methods for the computation of pagerank, Linear Algebra and its Applications 386 (2004) 51–65.
  * [8] S. D. Kamvar, T. H. Haveliwala, C. D. Manning, G. H. Golub, Extrapolation methods for accelerating pagerank computations, in: Proceedings of the 12th international conference on World Wide Web, 2003, pp. 261–270.
  * [9] G. Wu, Y. Wei, A power-arnoldi algorithm for computing pagerank, Numerical Linear Algebra with Applications 14 (7) (2007) 521–546.
  * [10] Y. Zhu, S. Ye, X. Li, Distributed pagerank computation based on iterative aggregation-disaggregation methods, in: Proceedings of the 14th ACM international conference on Information and knowledge management, 2005, pp. 578–585.
  * [11] A. D. Sarma, A. R. Molla, G. Pandurangan, E. Upfal, Fast distributed pagerank computation, in: International Conference on Distributed Computing and Networking, Springer, 2013, pp. 11–26.
  * [12] A. Das Sarma, D. Nanongkai, G. Pandurangan, P. Tetali, Distributed random walks, Journal of the ACM (JACM) 60 (1) (2013) 1–31.
  * [13] S. Luo, Distributed pagerank computation: An improved theoretical study, Proceedings of the AAAI Conference on Artificial Intelligence 33 (2019) 4496–4503.
  * [14] S. Luo, Improved communication cost in distributed pagerank computation–a theoretical study, in: International Conference on Machine Learning, PMLR, 2020, pp. 6459–6467.
  * [15] R. Andersen, F. R. K. Chung, K. J. Lang, Local graph partitioning using pagerank vectors, in: 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS 2006), 21-24 October 2006, Berkeley, California, USA, Proceedings, 2006.
  * [16] D. I. Shuman, P. Vandergheynst, P. Frossard, Chebyshev polynomial approximation for distributed signal processing, in: Distributed Computing in Sensor Systems, 7th IEEE International Conference and Workshops, DCOSS 2011, Barcelona, Spain, 27-29 June, 2011, Proceedings, 2011.
  * [17] Berkhin, Pavel, A survey on pagerank computing, Internet Mathematics 2 (1) (2005) 73–120.
  * [18] R. J. de la Cruz, P. Saltenberger, Density of diagonalizable matrices in sets of structured matrices defined from indefinite scalar products, arXiv preprint arXiv:2007.00269 (2020).
  * [19] L. Page, S. Brin, R. Motwani, T. Winograd, The pagerank citation ranking: Bringing order to the web., Tech. rep., Stanford InfoLab (1999).
  * [20] N. T. Duong, Q. A. P. Nguyen, A. T. Nguyen, H.-D. Nguyen, Parallel pagerank computation using gpus, in: Proceedings of the Third Symposium on Information and Communication Technology, 2012, pp. 223–230.
  * [21] Q. Zhang, R. Tang, Z. Yao, J. Liang, Two parallel pagerank algorithms via improving forward push, arXiv preprint arXiv:2302.03245 (2023).