本节书摘来华章计算机出版社《R的极客理想——高级开发篇 A》一书中的第2章，第2.2节，作者：张丹 更多章节内容可以访问云栖社区“异步社区”公众号查看。
问题如何用R语言实现PageRank算法？
引言Google搜索，早已成为我每天必用的工具，我无数次惊叹它搜索结果的准确性。同时，我也在做Google的SEO，推广自己的博客。经过几个月尝试，我的博客PR到2了，外链也有几万个。总结下来，还是感叹PageRank的神奇。笔者认为PageRank是改变互联网的算法！2.2.1 PageRank算法介绍PageRank是Google专有的算法，用于衡量特定网页相对于搜索引擎索引中的其他网页而言的重要程度。它由Larry Page和Sergey Brin在20世纪90年代后期发明。PageRank实现了将链接价值概念作为排名因素。PageRank是让链接来“投票”。一个页面的“得票数”由所有链向它的页面的重要性来决定，指向一个页面的超链接相当于对该页投一票。一个页面的PageRank是由所有链向它的页面（链入页面）的重要性经过递归算法得到的。一个有较多链入的页面会有较高的等级，相反如果一个页面没有任何链入页面，那么它没有等级。简单一句话概括，就是从许多优质的网页链接过来的网页，必定还是优质网页。PageRank的计算基于以下两个基本假设，一是数量假设，即如果一个页面节点接收到的其他网页指向的入链数量越多，那么这个页面越重要；二是质量假设，即指向页面A的入链质量不同，质量高的页面会通过链接向其他页面传递更多的权重，所以越是质量高的页面指向页面A，页面A越重要。要提高PageRank有3个要点，也就是我们做SEO时必须要考虑的问题：（1）反向链接数：反向链接数越多，指向页面的权重越高。（2）反向链接是否来自PageRank较高的页面：反向链接页面的权重越高，指向页面分配的权重就越高。（3）反向链接源页面的链接数：反向链接页面的链接数越多，则反向链接源页面的权重越高。循环解释前面两条。2.2.2 PageRank算法原理在初始阶段，网页通过链接关系构建起有向图，每个页面设置相同的PageRank值，通过若干轮的递归计算，会得到每个页面所获得的最终PageRank值。随着每一轮的计算进行，网页当前的PageRank值会不断得到更新。在一轮更新页面PageRank得分的计算中，每个页面将其当前的PageRank值平均分配到本页面所包含的出链上，这样每个链接即获得了相应的权值。而每个页面将所有指向本页面的入链所传入的权值求和，即可得到新的PageRank得分。当每个页面都获得了更新后的PageRank值，就完成了一轮PageRank计算。1.?算法原理PageRank算法建立在随机冲浪者模型上，其基本思想是：网页的重要性排序是由网页间的链接关系所决定的，算法是依靠网页间的链接结构来评价每个页面的等级和重要性，一个网页的PR值涉及指向它的链接网页数，还涉及指向它的网页的其他网页本身的重要性。PageRank具有两大特性，一是PR值的传递性，即网页A指向网页B时，网页A的PR值会部分传递给网页B；二是重要性的传递性，即一个重要的网页比一个不重要网页，传递的权重要更多。2.?计算公式PageRank算法公式为（2.1）：
（2.1）其中PR(pi)是pi页面的PageRank值，n是所有页面的数量，pi是不同的网页p1、p2、p3，M(i)是pi链入网页的集合，L(j)是pj链出网页的数量，d是阻尼系数，任意时刻，用户到达某页面后并继续向后浏览的概率。通过经验值，Google设置d=0.85，则1-d=0.15表示用户停止点击随机跳到新URL的概率，d的取值范围是03.?构造实例：以4个页面的数据为例我们以4个页面的数据为例，构造一个简单的PageRank实例模型，如图2-2所示。在图2-2中，ID=1的页面链向页面2、3、4，所以一个用户从ID=1的页面跳转到2、3、4的概率各为1/3，ID=2的页面链向页面3和4，所以一个用户从ID=2的页面跳转到3,4的概率各为1/2，ID=3的页面链向页面4，所以一个用户从ID=3的页面跳转到4的概率各为1，ID=4的页面链向页面2，所以一个用户从ID=4的页面跳转到2的概率各为1。下面我们开始构建PageRank的数据模型和转移矩阵。第一步，通过网页的链接关系，构造出邻接表：
链接源页面 链接目标页面
1 2,3,4 2 3,4 3 4 4 2第二步，通过邻接表，构建邻接矩阵（方阵），其中列表示源页面，行表示目标页面。
[,1] [,2] [,3] [,4] [1,] 0 0 0 0 [2,] 1 0 0 1 [3,] 1 1 0 0 [4,] 1 1 1 0第三步，把邻接矩阵，转换为概率矩阵（转移矩阵）。
[,1] [,2] [,3] [,4] [1,] 0 0 0 0 [2,] 1/3 0 0 1 [3,] 1/3 1/2 0 0 [4,] 1/3 1/2 1 0这样，通过网页的链接关系，我们就构造出了概率转移矩阵。2.2.3 R语言单机算法实现下面我们用R语言实现从邻接表到转移矩阵的构建过程，创建数据文件page.csv，用于表示网页的邻接表。
1,2 1,3 1,4 2,3 2,4 3,4 4,2分别用下面3种方式实现PageRank模型，即未考虑阻尼系统的情况，考虑阻尼系统的情况，直接用R的矩阵特征值函数计算。1.?未考虑阻尼系统的情况R语言实现如下：
> adjacencyMatrix<-function(pages){ # 构建邻接矩阵 + n<-max(apply(pages,2,max)) + A <- matrix(0,n,n) + for(i in 1:nrow(pages)) A[pages[i,]$dist,pages[i,]$src]<-1 + A +} > probabilityMatrix<-function(G){ # 变换概率矩阵 + cs <- colSums(G) + cs[cs==0] <- 1 + n <- nrow(G) + A <- matrix(0,nrow(G),ncol(G)) + for (i in 1:n) A[i,] <- A[i,] + G[i,]/cs + A +} > eigenMatrix<-function(G,iter=100){ # 递归计算矩阵特征值 + iter<-10 + n<-nrow(G) + x <- rep(1,n) + for (i in 1:iter) x <- G %*% x + x/sum(x) +}结果解读如下：ID=1的页面，PR值是0，因为没有指向ID=1的页面。ID=2的页面，PR值是0.4，权重最高，因为页面1和4都指向页面2，页面4权重较高，并且页面4只有一个链接指向页面2，权重传递没有损失。ID=3的页面，PR值是0.19，虽有页面1和2都指向了页面3，但是页面1和2还指向其他页面，权重被分散了，所以ID=3的页面PR并不高。ID=4的页面，PR值是0.39，权重很高，因为被页面1，2，3都指向了。从上面的结果，我们发现ID=1的页面，PR值是0，那么ID=1的页面，就不能向其他页面输出权重了，计算就会不合理！所以，增加阻尼系数d，修正没有链接指向的页面，保证页面的最小PR值>0。2.?考虑阻尼系统的情况这里增加函数dProbabilityMatrix。
dProbabilityMatrix<-function(G,d=0.85){ # 变换概率矩阵，考虑阻尼系数d的情况
+ cs <- colSums(G) + cs[cs==0] <- 1 + n <- nrow(G) + delta <- (1-d)/n + A <- matrix(delta,nrow(G),ncol(G)) + for (i in 1:n) A[i,] <- A[i,] + d*G[i,]/cs + A +} # 运行程序 > pages<-read.table(f?ile="page.csv",header=FALSE,sep=",") # 读数据文件到内存 > names(pages)<-c("src","dist") # 设置数据的header > A<-adjacencyMatrix(pages);A # 构建邻接矩阵 [,1] [,2] [,3] [,4] [1,] 0 0 0 0 [2,] 1 0 0 1 [3,] 1 1 0 0 [4,] 1 1 1 0 > G<-dProbabilityMatrix(A);G # 概率转移矩阵，考虑阻尼系数d的情况 [,1] [,2] [,3] [,4] [1,] 0.0375000 0.0375 0.0375 0.0375 [2,] 0.3208333 0.0375 0.0375 0.8875 [3,] 0.3208333 0.4625 0.0375 0.0375 [4,] 0.3208333 0.4625 0.8875 0.0375 > q<-eigenMatrix(G,100);q # PageRank值 [,1] [1,] 0.0375000 [2,] 0.3738930 [3,] 0.2063759 [4,] 0.3822311 增加阻尼系数后，ID=1的页面，就有值了PR(1)=(1-d)/n=(1-0.85)/4=0.0375，即无外链页面的最小值。 3.?直接用R的特征值计算函数 上面的算法，我们做的都是递归计算，通过多次循环得到矩阵的特征值。R语言本身提供了直接计算特征值的函数，我们可以直接调用，从而简化代码的复杂度。下面增加计算特征值函数calcEigenMatrix。 > calcEigenMatrix<-function(G){ # 直接计算矩阵特征值 + x <- Re(eigen(G)$vectors[,1]) + x/sum(x) + } # 运行程序 > pages<-read.table(f?ile="page.csv",header=FALSE,sep=",") # 读数据文件到内存 > names(pages)<-c("src","dist") # 设置数据的header > A<-adjacencyMatrix(pages);A # 构建邻接矩阵 [,1] [,2] [,3] [,4] [1,] 0 0 0 0 [2,] 1 0 0 1 [3,] 1 1 0 0 [4,] 1 1 1 0 > G<-dProbabilityMatrix(A);G # 通过特征值函数计算，概率转移矩阵 [,1] [,2] [,3] [,4] [1,] 0.0375000 0.0375 0.0375 0.0375 [2,] 0.3208333 0.0375 0.0375 0.8875 [3,] 0.3208333 0.4625 0.0375 0.0375 [4,] 0.3208333 0.4625 0.8875 0.0375 > q<-calcEigenMatrix(G);q # PageRank值 [1] 0.0375000 0.3732476 0.2067552 0.3824972直接计算矩阵特征值，可以有效地减少循环的操作，提高程序运行效率。2.2.4 R语言分步式算法实现使用PageRank的场景大都是基于海量数据的，通常会用Hadoop的解决方案，通过分步式并行算法计算出每个页面的得分。利用R语言，我可以先构建出分步式算法的原型，保证分步式算法的结果与单机算法的结果是一致的，从而总结出分步式算法的计算公式，再用Hadoop的MapReduce进行重写，从而达到生产环境的上线要求。1.?PageRank算法分步式原理PageRank的分步式算法原理，简单来讲，就是通过矩阵计算实现并行化。我们直接从邻接矩阵开始说起，对于Hadoop来说需要把邻接矩阵的列按数据行存储，而R语言可以直接用矩阵对象表示。
[,1] [,2] [,3] [,4] [1,] 0.0375000 0.0375 0.0375 0.0375 [2,] 0.3208333 0.0375 0.0375 0.8875 [3,] 0.3208333 0.4625 0.0375 0.0375 [4,] 0.3208333 0.4625 0.8875 0.0375接下来，迭代求矩阵特征值，按照MapReduce的设计思路，把R语言的算法也分为map过程和reduce过程。map过程主要用于数据分组，reduce过程主要用于数据计算，让map和reduce循环迭代，从而计算出矩阵特征值。Hadoop的MapReduce计算过程如图2-3所示。对于Hadoop的map过程，输入为邻接矩阵和pr值；输出中key为pr的行号，value为邻接矩阵和pr值的乘法求和公式。对于Hadoop的reduce过程，输入中key为pr的行号，value为邻接矩阵和pr值的乘法求和公式；输出中key为pr的行号，value为计算的结果，即pr值。2.?R语言程序模拟我们通过R语言来实现这个程序，把矩阵计算分解为map和reduce两个过程。在map过程中，把数据按一定规则切分成多块，海量数据分布在Hadoop的多个存储节点上用于模拟。本节使用的数据集仅包括了4个网页，我把数据切分成2组，模拟存储在2个Hadoop节点上。
> map <- function(S0, node = "a") { # map过程函数 + S <- apply(S0, 2, function(x) x/sum(x)) + if (node == "a") + S[, 1:2] else S[, 3:4] + }``` 在reduce过程中，数据从多个map节点汇总到reduce节点，进行数据合并计算。通过R语言，我直接在reduce()函数内部完成了迭代计算的过程。而对完全的分步式的MapReduce来说，每次迭代都要包括完成的map过程和reduce过程。reduce <- function(A, B, a = 0.85, niter = 100) { # reduce过程函数
n <- nrow(A)q <- rep(1, n)Ga <- a A + (1 - a)/n (A[A != 1] = 1)Gb <- a B + (1 - a)/n (B[B != 1] = 1)for (i in 1:niter) { # 迭代计算qa <- as.matrix(q[1:ncol(A)])qb <- as.matrix(q[(ncol(A) + 1):n])q <- Ga %% qa + Gb %% qb}as.vector(q/sum(q)) # 标准化PR值} 运行程序，进行100次迭代，计算矩阵特征值（PR值），比较计算结果，发现这与单机的计算结果是一致的。S0 <- t(matrix(c(0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0),nrow = 4))
# 初始数据矩阵A <- map(S0, "a") # 矩阵数据分组B <- map(S0, "b")reduce(A, B) # 矩阵数据计算[1] 0.0375000 0.3732476 0.2067552 0.3824972
3.?矩阵的计算过程 对于每次迭代计算的PR值，我们可以打印出来，从而清楚地看到每次迭代的结果。把reduce()函数的代码稍加修改，增加原始PR值输出和标准化的PR值输出。reduce <- function(A, B, a = 0.85, niter = 100) {
n <- nrow(A)q <- rep(1, n)Ga <- a A + (1 - a)/n (A[A != 1] = 1)Gb <- a B + (1 - a)/n (B[B != 1] = 1)for (i in 1:niter) {qa <- as.matrix(q[1:ncol(A)])qb <- as.matrix(q[(ncol(A) + 1):n])q <- Ga %% qa + Gb %% qb}print(q) # 原始PR值print(as.vector(q/sum(q)))
}第1次迭代的计算结果。
reduce(A, B, niter=1) # 原始PR值
[,1][1,] 0.1500000
[2,] 1.2833333[3,] 0.8583333[4,] 1.7083333[1] 0.0375000 0.3208333 0.2145833 0.4270833 # 标准化的PR值
第1次迭代时，原始PR值的矩阵计算公式。0.0375000 0.0375 0.0375 0.0375 1 0.1500000.3208333 0.0375 0.0375 0.8875 * 1 = 1.2833330.3208333 0.4625 0.0375 0.0375 1 0.8583330.3208333 0.4625 0.8875 0.0375 1 1.708333
第2次迭代的计算结果。reduce(A, B, niter=2)
[,1][1,] 0.1500000
[2,] 1.6445833[3,] 0.7379167[4,] 1.4675000[1] 0.0375000 0.4111458 0.1844792 0.3668750第2次迭代时，原始PR值的矩阵计算公式。0.0375000 0.0375 0.0375 0.0375 0.150000 0.1500000.3208333 0.0375 0.0375 0.8875 * 1.283333 = 1.64458330.3208333 0.4625 0.0375 0.0375 0.858333 0.73791670.3208333 0.4625 0.8875 0.0375 1.708333 1.4675000
第10次迭代计算结果。reduce(A, B, niter=10)
[,1][1,] 0.1500000
[2,] 1.4955721[3,] 0.8255034[4,] 1.5289245[1] 0.0375000 0.3738930 0.2063759 0.3822311
第10次迭代时，标准化PR值的矩阵计算公式。0.150000 0.03750001.4955721 / (0.15+1.4955721+0.8255034+1.5289245) = 0.37389300.8255034 0.20637591.5289245 0.3822311
在了解PageRank的原理后，使用R语言构建PageRank模型，是非常容易的。实际应用中，我们也愿意用比较简单的方式建模、验证后，再用其他更底层的语言实现企业应用！ 在我的博客中还有一些进阶的内容，介绍了如何用MapReduce分步式算法来实现PageRank模型，请参考文章《PageRank算法并行实现》（http://blog.fens.me/algorithm-pagerank-mapreduce/）。