从线性代数角度理解 PageRank 算法细节 (包含基于 NumPy 的代码实现)_综合

前言

近期在研究自然语言处理 (NLP) 领域中的关键词抽取 (Keyword Extraction) 任务，遇到一个经典算法 TextRank ¹，其中心思想便是借鉴于 PageRank 算法。

PageRank 算法² 最早由 Google 提出，用于对搜索返回的网页进行重要性排序。该算法的思想很直观：对于一个网页，当链接到它的其他网页的重要性分数总和越高，它就越重要；这同时包含了两种情况：(1) 链接到它的网页的数量较多；(2) 链接到它的网页的重要性较高 (A page has high rank if the sum of the ranks of its backlinks is high. This covers both the case when a page has many backlinks and when a page has a few highly ranked backlinks)。目前 PageRank 虽不是 Google 用来排序网页的唯一算法，但是其最早使用的，也是最著名的网页排序算法。

虽然 PageRank 算法在理解上和实现上都比较简单，但有些细节还是值得好好推敲一下。

一些概念

这里先介绍一些描述算法时需要用到的概念术语：

悬挂节点 (dangling nodes)：在网络分析中，悬挂节点指没有传出链接 (outgoing links) 的节点；一个不包含任何外部链接的网页就是网络中的一个悬挂节点。本文会用“悬挂网页”指网络中属于悬挂节点的网页。
随机矩阵 (stochastic matrix) ：用于描述一个马尔科夫链的转变的非负实方阵，又称概率矩阵、转移矩阵 (transition matrix)。左随机矩阵 (left stochastic matrix) 每一列元素的和均为 1，而右随机矩阵 (right stochastic matrix) 中的每一行元素的和为 1。注意：随机矩阵不是指随机生成的矩阵。
强连通图 (strongly connected graph)：是指在有向图 G 中，如果对于每一对 $v_i, v_j$ ( $vi≠vjv_i \neq v_j$ )，从 $v_i$ 到 $v_j$ 和从 $v_j$ 到 $v_i$ 都存在路径 (这条路径中间可以存在多个其它节点)。
前向链接 (forward links) 与 反向链接 (backlinks)：如果网页 $i$ 中有指向网页 $j$ 的链接，那么网页 $j$ 就属于网页 $i$ 的前向链接，而网页 $i$ 属于网页 $j$ 的反向链接。
正矩阵 (positive matrix)：当一个矩阵内的所有元素均大于 0 时，该矩阵即为正矩阵。注意不要与正定矩阵相混淆。

PageRank 算法简易版本

算法思想

根据前言那一章中提到的算法思想，可以先简单地把网页 $i$ 的分数 (重要性分数的简写，下同) 定义为它的反向链接的分数的总和。但当一个网页包含前向链接越多时，我们希望它传递给这些前向网页的分数越少。

所以当网页 $j$ 包含 $n_j$ 个前向链接时，我们令这 $n_j$ 个网页从网页 $j$ 获取的传递分数为 $x_j/n_j$ ，而不是 $x_j$ . 令 $L_i$ 表示网页 $i$ 的反向链接的集合，即该集合中的每个网页都包含指向网页 $i$ 的链接，那么网页 $i$ 的分数为：
$xi=∑j∈Lixjnj(1)x_i = \sum_{j\in L_i} \frac{x_j}{n_j} \tag{1}$

下图为文献³ 中的一个例子，其中矩形表示网页，有向箭头表示链接。

a simple web

图 1. 一个简单的网络拓扑结构

先用 $n×nn\times n$ 的邻接矩阵 $A\mathbf{A}$ 来表示该网络结构， $n$ 为网页总数，该网络结构是一个强连通图，当网页 $j$ 是 $i$ 的反向链接时， $A[i][j]=1\mathbf{A}[i][j] = 1$ ，否则 $A[i][j]=0\mathbf{A}[i][j] = 0$ ；令 $A\mathbf{A}$ 中所在列向量和不为 0 的元素除以该列的和，得到一个左随机矩阵 $S\mathbf{S}$ ，显然该矩阵所有元素是非负的：
$A=[0011100011011100],S=[00112130001312012131200]\mathbf{A} = \begin{bmatrix} 0 & 0 & 1 & 1\\ 1 & 0 & 0 & 0\\ 1 & 1 & 0 & 1\\ 1 & 1 & 0 & 0\\ \end{bmatrix}, \quad \mathbf{S} = \begin{bmatrix} 0 & 0 & 1 & \frac{1}{2} \\ \frac{1}{3} & 0 & 0 & 0\\ \frac{1}{3} & \frac{1}{2} & 0 & \frac{1}{2}\\ \frac{1}{3} & \frac{1}{2} & 0 & 0\\ \end{bmatrix}$

这时我就可将公式 (1) 写成矩阵形式：
$Sx=x(2)\mathbf{S}\mathbf{x}=\mathbf{x} \tag{2}$

其中 $x=[x1,x2,x3,x4]T\mathbf{x}=[x_1,x_2,x_3,x_4]^T$ . 所以，求 $x\mathbf{x}$ 就相当于求矩阵 $S\mathbf{S}$ 在特征值等于 1 时的特征向量。而上述的矩阵 $S\mathbf{S}$ 确实有特征值等于 1，且该特征值对应的特征向量均为 $12, 4, 9, 16]^T$ 的常数倍(注意，一特征值对应的特征向量与一常数相乘仍为该特征值的特征向量)，也说明 $S\mathbf{S}$ 在特征值为 1 时的特征空间的维数为 1。我们希望所有网页的重要性分数相加等于 1，那么可将上述向量缩放使其元素总和为 1，得到 $x1=1231,x2=431,x3=931,x4=1631x_1=\frac{12}{31}, x_2=\frac{4}{31}, x_3=\frac{9}{31}, x_4=\frac{16}{31}$ ，这些就是该例子中四个网页的 PageRank 值。

其实对于不存在悬挂节点的网络，其对应的矩阵 $S\mathbf{S}$ 必有等于 1 的特征值。为什么？首先，一个网络如果不存在悬挂节点，那么 $A\mathbf{A}$ 的每一列的和均大于 0，那么 $S\mathbf{S}$ 必然是一个左随机矩阵，接下来证明左随机矩阵均有等于 1 的特征值：

令 $S\mathbf{S}$ 为一个 $n×nn\times n$ 的左随机矩阵，而 $e\mathbf{e}$ 是一个元素均为 1 的 $n$ 维列向量。对于一个方阵，它与它自身的转置有着相同的特征值，所以 $S\mathbf{S}$ 与 $ST\mathbf{S}^T$ 的特征值相同。由于 $S\mathbf{S}$ 是左随机矩阵，所以有 $STe=e\mathbf{S}^Te=e$ ，也就是 $ST\mathbf{S}^T$ 与 $S\mathbf{S}$ 均有等于 1 的特征值。

虽然，这里介绍的 PageRank 简易版本能够很好地计算图 1 所示的例子，但它还是有一定的缺陷。

简易版本的缺陷

对于非全连通的网络存在不唯一的 PageRank 值

我们使用 $V1(S)V_1{(\mathbf{S})}$ 来表示左随机矩阵 $S\mathbf{S}$ 在特征值为 1 时的特征空间。当 $dim(V1(S))=1\text{dim}(V_1{(\mathbf{S}))} = 1$ 时，说明 $S\mathbf{S}$ 在特征值为 1 时的特征向量均线性相关，令其中一个特征向量为 $a\mathbf{a}$ ，那其他特征向量均可由 $kak\mathbf{a}$ 表示 ( $k$ 为任意常数)，此时只有一个特征向量 $x\mathbf{x}$ 满足 $∑ixi=1\sum_i x_i = 1$ ，也就是当 $k=1∑iaik=\frac{1}{\sum_i a_i}$ 时.

two disconnected webs

图 2. 两个不相连通的子网络

对于具有强连通图结构的网络，其对应的左随机矩阵 $S\mathbf{S}$ 满足 $dim(V1(S))=1\text{dim}(V_1{(\mathbf{S}))} = 1$ ，相关证明可见本文最后一章“补充证明”。此时我们能获得唯一的 PageRank 结果，即 $S\mathbf{S}$ 对应于特征值 1 的唯一满足 $∑ixi=1\sum_i x_i = 1$ 条件的特征向量 $x\mathbf{x}$ 。但网络结构并非都是强连通的。以图 2 为例，该网络包含两个不相连通的子网络，所以该网络不是强连通的，该网络对应的左随机矩阵 $S\mathbf{S}$ 为：
$S=[010001000000011200101200000]\mathbf{S} = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & \frac{1}{2} \\ 0 & 0 & 1 & 0 & \frac{1}{2} \\ 0 & 0 & 0 & 0 & 0 \\ \end{bmatrix}$
此时 $dim(V1(S))=2\text{dim}(V_1{(\mathbf{S}))} = 2$ ，该特征空间的一组基向量为： $x=[1/2,1/2,0,0,0]T,y=[0,0,1/2,1/2,0]T\mathbf{x} = [1/2, 1/2, 0, 0, 0]^T, \mathbf{y} = [0,0,1/2,1/2,0]^T$ ,这时特征值为 1 时满足 $∑ixi=1\sum_i x_i = 1$ 条件的特征向量 $x\mathbf{x}$ 并不唯一，我们并不清楚应该应用哪个特征向量作为 PageRank 的结果。其实，如果一个网络由 $r$ 个非相互连通的子网络 $W_1,...,W_r$ 组成，那么其对应的 $S\mathbf{S}$ 满足 $dim(V1(S))≥r\text{dim}(V_1(\mathbf{S})) \geq r$ . 这里我不再进行证明，感兴趣的读者可以参考文献³。

无法应对悬挂节点

当一个网络存在悬挂节点，那其对应的邻接矩阵 $A\mathbf{A}$ 必存在和为 0 的列向量，所以其对应的 $S\mathbf{S}$ 不是左随机矩阵，也无法按照公式 (2) 获取所有网页的 PageRank 值。

PageRank 算法通用版本

应对非强连通网络

在不存在悬挂节点的情况下，要如何应对非强连通网络的情况呢？令 $N\mathbf{N}$ 表示一个元素均为 $1 / n$ 的 $n×nn\times n$ 矩阵，显然它是一个左随机矩阵，其对应的网络结构是强连通的，所以 $dim(V1(N))=1\text{dim}(V_1(\mathbf{N})) = 1$ 。然后，将 $S\mathbf{S}$ 替换成 $M\mathbf{M}$ :

$M=(1?d)S+dN\mathbf{M} = (1-d) \mathbf{S} + d \mathbf{N}$
其中 $\leq m \leq 1$ ， $M\mathbf{M}$ 为 $S\mathbf{S}$ 和 $N\mathbf{N}$ 的加权平均。 $S\mathbf{S}$ 和 $N\mathbf{N}$ 均是左随机矩阵，所以对于 $m∈[0,1]m\in [0,1]$ ， $M\mathbf{M}$ 均为左随机矩阵，证明很简单，可参考“补充证明”一章；对于 $\in (0,1]$ ， $M\mathbf{M}$ 是正矩阵 (所有元素均大于 0)，说明在其对应的网络结构中，所有节点都直接相连，所以其对应的网络结构也是强连通的 (注意：强连通仅要求间接相连)，此时 $dim(V1(M))=1\text{dim}(V_1(\mathbf{M})) = 1$ 。所以，当 $\in (0,1]$ 时， $M\mathbf{M}$ 可替代公式 (2) 中的 $S\mathbf{S}$ 以应对非强连通网络。所以公式 (2) 可替换成：
$x=(1?d)Sx+dn(3)\mathbf{x} = (1-d) \mathbf{S} \mathbf{x} + d \mathbf{n} \tag{3}$
其中 $n\mathbf{n}$ 是一个元素均为 $1 / n$ 的列向量；当 $∑ixi=1\sum_i x_i = 1$ 时，有 $Nx=n\mathbf{N} \mathbf{x} = \mathbf{n}$ .

引用维基百科里的描述，我们可以这样理解上述的做法（以下提到的“没有外部链接的网页”即为本文中悬挂网页）：

这个方程式引入了随机浏览者（random surfer）的概念，即假设某人在浏览器中随机打开某些页面并点击了某些链接。为了便于理解，这里假设上网者不断点击网页上的链接直到进入一个没有外部链接的网页，此时他会随机浏览其他的网页（可以与之前的网页无关）。为了处理那些“没有外部链接的页面”（这些页面就像“黑洞”一样吞噬掉用户继续向下浏览的概率）所带来的问题，我们假设：这类页面链接到集合中所有的网页（不管它们是否相关），使得这类网页的 PageRank 值将被所有网页均分。对于这种残差概率（residual probability），我们引入阻尼系数 $d$ （damping factor），并声明 $d = 0.85$ ，其意义是：任意时刻，用户访问到某页面后继续访问下一个页面的概率，相对应的 $1 ? d = 0.15$ 则是用户停止点击，随机浏览新网页的概率。 $d$ 的大小由一般上网者使用浏览器书签功能的频率的平均值估算得到。

对于一个存在悬挂节点的网络，使用 $M\mathbf{M}$ 替代 $S\mathbf{S}$ ，能对悬挂网页赋予重要性分数 $d / n$ , 但此时的 $M\mathbf{M}$ 仍不属于左随机矩阵 (当 $m > 0$ 时)，所以使用公式 (3) 仍无法解决悬挂节点的问题。

应对悬挂节点

文献²中有提到，Google 的做法就是在计算 PageRank 时直接将悬挂网页删除，也就是计算 $A\mathbf{A}$ 时并不考虑那些悬挂网页。

还有一种做法，就是令悬挂网页为所有其他网页的反向链接，即在计算 $A\mathbf{A}$ 时将所有悬挂节点对应的列向量里的元素均赋值为 $1 / n$ ，之后再计算 $S\mathbf{S}$ 并通过公式 (3) 计算 $x\mathbf{x}$ .

计算各网页的 PageRank 值: 特征向量 $x\mathbf{x}$

有了公式 (3)，那计算机要如何计算 $x\mathbf{x}$ 呢？常用的方法是幂法 (power method)。幂法主要用于近似计算矩阵的主特征值 (dominant eigenvalue) 和主特征向量 (dominant eigenvector)⁴。随机矩阵均有等于 1 的特征值，且其他特征值均小于 1，证明可见文献⁵。所以随机矩阵的主特征值为 1，能够通过幂法获取主特征向量 $x\mathbf{x}$ .

使用幂法：先从一个初始化向量 $x0\mathbf{x}^0$ 开始，使用下面公式迭代计算：
$xt=Mxt?1(4)\mathbf{x}^t = \mathbf{M} \mathbf{x}^{t-1} \tag{4}$
直到 $xt\mathbf{x}^t$ 收敛 ( $?\epsilon$ 为一阈值)：
$∥xt?xt?1∥<?\Vert \mathbf{x}^{t} - \mathbf{x}^{t-1} \Vert < \epsilon$

可证明，当初始化的 $x0\mathbf{x}^0$ 的各元素之和等于 1 时，即满足 $∑ixi0=1\sum_i \mathbf{x}_i^0 = 1$ 时，对于 $\in (0, +\infin)$ ， $xt\mathbf{x}^t$ 均满足 $∑ixit=1\sum_i \mathbf{x}_i^t = 1$ :

由于 $M\mathbf{M}$ 是一个左随机矩阵，即对于 $j = 1, 2, . . ., n$ ，满足 $∑inMij=1\sum_i^n\mathbf{M}_{ij} = 1$ 。所以根据公式 (4) 有：
$∑inxit=∑in∑jnMijxjt?1=∑jnxjt?1∑inMij=∑jnxjt?1\sum_i^n \mathbf{x}_i^t = \sum_i^n \sum_j^n \mathbf{M}_{ij} \mathbf{x}_j^{t-1} = \sum_j^n \mathbf{x}_j^{t-1} \sum_i^n \mathbf{M}_{ij} = \sum_j^n \mathbf{x}_j^{t-1}$
说明与左随机矩阵相乘并不会改变向量的元素之和，所以只要初始 $x0\mathbf{x}^0$ 的元素之和等于 1，那接下来迭代获得的 $xt\mathbf{x}^t$ 的元素之和也会等于 1。

代码实现 (基于 NumPy)

import numpy as npN = 1000 # 网页数量
d = 0.15 # 阻尼系数################ 样例数据生成 #######################
# 随机产生邻接矩阵
A = np.random.randint(0, 2, size=(N, N))
# dangling_nodes = np.random.choice(N, 3) # 随机产生悬挂节点
# A[:, dangling_nodes] = 0
# 获取悬挂节点
danling_nodes = A.sum(axis=0) == 0.0
# normal_nodes = ~danling_nodes
# 处理悬挂节点 (dangling nodes)
A[:, danling_nodes] = 1.0
# data type: int to float
A = A / 1.0 
# 邻接矩阵 -> 左随机矩阵
S = A / A.sum(axis=0).reshape(1, -1) 
# 计算 M
M = (1-d) * S + d * (1 / N)
# 初始化 x 元素均为 1/n，此时 x 的 1-范数等于 1
x = np.full(N, fill_value=1/N, dtype=float)################ 迭代计算 #######################
iters = 0
threshold = 1e-5 # 阈值
while True:x_ = xx = np.dot(M, x)print(f'The current sum of PR: {
      x.sum()}')iters += 1if abs(x - x_).sum() < threshold:print(f'The last iterations: {
      iters}')print(f'The final PRs: {
      x[:10]}')break

运行上面的代码可以得到类似输出：

The current sum of PR: 1.0000000000000007
The current sum of PR: 1.0000000000000007
The current sum of PR: 1.0000000000000009
The current sum of PR: 1.0000000000000009
The last iterations: 4
The final PRs: [0.00099134 0.00101442 0.00096278 0.00098558 0.00097684 0.001024030.00100729 0.00097193 0.00100127 0.00101572]

从此可看出，在迭代过程中，所有网页 PageRank 值的总和一直没变，等于初始值的总和，也就是 1。当我们不对悬挂节点进行处理时，将上述代码进行相应的修改：

################ 样例数据生成 #######################
# 随机产生邻接矩阵
A = np.random.randint(0, 2, size=(N, N))
# 随机产生悬挂节点
dangling_nodes = np.random.choice(N, 3)
A[:, dangling_nodes] = 0
# 获取悬挂节点
dangling_nodes = A.sum(axis=0) == 0.0
normal_nodes = ~dangling_nodes
# 处理悬挂节点 (dangling nodes)
# A[:, danling_nodes] = 1.0
# data type: int to float
A = A / 1.0 
# 邻接矩阵 -> 左随机矩阵
A[:, normal_nodes] /= A[:, normal_nodes].sum(axis=0).reshape(1, -1)
S = A
# 计算 M
M = (1-d) * S + d * (1 / N)
# 初始化 x 元素均为 1/n，此时 x 的 1-范数等于 1
x = np.full(N, fill_value=1/N, dtype=float)

运行，从结果中可以看到随着迭代的进行，所有网页的 PageRank 总和一直在减少，直至为 0：

The current sum of PR: 0.9974500000000005
The current sum of PR: 0.994856218039011
The current sum of PR: 0.9922695031287359
The current sum of PR: 0.9896895490506252
...
The current sum of PR: 0.003835323251689406
The last iterations: 2137
The final PRs: [3.86809739e-06 3.78727177e-06 3.83120473e-06 3.86735285e-063.86316246e-06 3.62640880e-06 3.92647718e-06 3.76577298e-063.54793272e-06 3.77907705e-06]

补充证明

证明：两左随机矩阵的加权平均仍为左随机矩阵

给定两个 $n×nn\times n$ 的左随机矩阵 $S\mathbf{S}$ 和 $N\mathbf{N}$ ，令：
$M=(1?d)S+dN\mathbf{M} = (1-d) \mathbf{S} + d \mathbf{N}$
则对于 $j = 1, 2, . . ., n$ ，有：
$∑iMj=(1?d)∑iSj+d∑iNj=(1?d)+d=1\sum_i \mathbf{M}_j = (1-d) \sum_i \mathbf{S}_j + d \sum_i \mathbf{N}_j = (1-d) + d = 1$
即 $M\mathbf{M}$ 也是左随机矩阵，得证。

证明：对于强连通网络对应的左随机矩阵 $S\mathbf{S}$ ，有 $dim(V1(S))=1\text{dim}(V_1(\mathbf{S})) = 1$

已知矩阵 $S\mathbf{S}$ 对应的是强连通图，根据文献⁶ 中的定理 1.4 或百度百科不可约矩阵词条中的定理 1，可以确定 $S\mathbf{S}$ 是不可约的，这时就可借助 Perron–Frobenius 定理来证明此时 $dim(V1(S))=1\text{dim}(V_1(\mathbf{S})) = 1$ ，具体可参考维基百科：

在这里插入图片描述

有关 Perron–Frobenius 定理的证明可参考文献⁶。

参考源

TextRank: Bringing Order into Texts ??
The PageRank Citation Ranking: Bringing Order to the Web ?? ??
The $25,000,000,000 Eigenvector: The Linear Algebra behind Google ?? ??
《统计学习方法》李航著 ??
Stochastic Matrices ??
Nonnegative and spectral matrix theory(Lecture notes) ?? ??

从线性代数角度理解 PageRank 算法细节 (包含基于 NumPy 的代码实现)

前言