简介面向大数据、搜索引擎与复杂网络分析学习者的PageRank计算示例包以谷歌网页链接数据为对象给出三种Python实现基于稀疏矩阵的单机迭代计算、直接调用NetworkX库完成PageRank、在NetworkX中自行编写PageRank逻辑并与库函数对照。三种写法从不同层面覆盖同一算法适合中高级Python开发者理解稀疏矩阵存储、迭代收敛和图计算库的封装差异也可作为课堂实验或面试准备的动手材料。压缩包共4个文件包含3个Python脚本和1个.gz格式的原始网页数据文件整体大小19.38MB脚本可直接运行gz数据解压后即可用于实验目录结构简洁。目前已有824人学习/下载。通过对比三种实现读者能够掌握PageRank公式的稀疏矩阵表达熟悉NetworkX的基本调用和自定义迭代写法并可从耗时与结果一致性角度验证算法的正确性为后续处理更大规模图数据或迁移到分布式计算框架打下基础。1. 用稀疏矩阵啃下 web-Google.txt单机 PageRank 的计算思路第一次拿 web-Google.txt 跑 PageRank 的时候我以为最麻烦的是算法本身结果差点栽在数据读入上。这份数据集记录了 Google 网页之间的超链接关系约 87.5 万个节点、510 万条边如果直接转成稠密矩阵875713 乘以 875713 的浮点数组就能把 16G 内存吃穿。用 Python 的稀疏矩阵做单机计算是性价比最高的路线——不需要分布式环境一台普通笔记本就能把三种 PageRank 实现全部跑完。这篇笔记把预处理、幂迭代、Gauss-Seidel 迭代和稀疏线性求解器四条链路拆开重点放在参数设置和翻车点上适合做大数据课程设计、算法复现或者单纯想搞清楚 PageRank 工程细节的人。2. web-Google.txt 数据预处理逐行解析、出度统计与 CSR 构建2.1 数据格式与读盘策略别用 loadtxt 直接怼web-Google.txt 每行是两个整数表示一条从源网页指向目标网页的超链接例如11342 0表示节点 11342 有一条边指向节点 0。文件头部会有若干以#开头的注释行记录节点数、边数等元信息读入时必须跳过。第一次处理这种规模的数据我习惯先写一个独立脚本做解析与统计不急着算 PageRank。因为后面三种方法都要用到两个基础信息节点总数 N、每个节点的出度。如果直接用numpy.loadtxt加载全部 510 万行它内部会先把所有行读成 Python 对象再转数组内存峰值会冲到好几个 G单机场景下这就是第一个坑。import numpy as np src_list [] dst_list [] max_node 0 with open(web-Google.txt, r) as f: for line in f: line line.strip() if not line or line.startswith(#): continue parts line.split() if len(parts) 2: continue u int(parts[0]) v int(parts[1]) src_list.append(u) dst_list.append(v) max_node max(max_node, u, v) N max_node 1 src np.array(src_list, dtypenp.int32) dst np.array(dst_list, dtypenp.int32) print(fnodes: {N}, edges: {len(src)})这段逻辑很直接src和dst分别存边的起点和终点max_node记录出现的最大节点 IDN 取max_node 1这样能覆盖所有节点。数据类型选int32而不是int64因为 87 万个节点用 32 位整数绰绰有余能省一半数组内存。解析完成后 510 万个int32二元组大约占用 40MB属于可控范围。这里有个容易被忽视的细节web-Google.txt 的节点 ID 是稀疏的也就是说并不是 0 到 N-1 每个 ID 都真实存在。取 max 值加一作为矩阵维度会多出一些空行空列。这些空节点没有出边也没有入边在 PageRank 计算中会被自然当作 dangling 节点处理不影响最终排名结果但会让矩阵维度偏大一点。如果后续想压缩矩阵规模可以建一个 ID 到连续索引的映射不过对 87 万这个量级来说没必要。2.2 从邻接表到稀疏矩阵scipy.sparse 的构建顺序拿到src和dst两个一维数组后下一步是构建稀疏邻接矩阵。这里的核心约束是行表示源节点列表示目标节点A[i, j] 1表示存在 i 指向 j 的超链接。from scipy.sparse import coo_matrix, csr_matrix data np.ones(len(src), dtypenp.float64) A_coo coo_matrix((data, (src, dst)), shape(N, N)) A A_coo.tocsr() out_degree np.diff(A.indptr).astype(np.float64) print(fnnz: {A.nnz}) print(fmax out degree: {out_degree.max()}) print(fdangling nodes: {(out_degree 0).sum()})coo_matrix本质是三元组存储三个数组分别记录非零元素的值、行号、列号构建起来最快适合从边列表直接组装。但它不支持高效的矩阵乘法所以必须调用tocsr()转成 CSR 格式。CSR 用indptr、indices、data三个数组压缩存储按行访问极快是后续所有稀疏矩阵乘法的基石。out_degree的计算方式值得说清楚CSR 的indptr数组长度为 N1第 i 行非零元素个数等于indptr[i1] - indptr[i]所以np.diff(A.indptr)直接得到每个节点的出度。这段代码跑完后打印出的dangling nodes数量通常有几千到几万个它们是后续迭代里概率泄漏的源头避坑章节会专门展开。内存占用可以提前估算indptr是 87.5 万加一长度的整型数组约 7MBindices是 510 万个整型约 20MBdata是 510 万个浮点数约 40MB整体稀疏矩阵约 70MB加上src、dst的 40MB预处理阶段峰值控制在 120MB 以内单机完全无压力。3. 方法一幂迭代法——向量矩阵乘的收敛艺术3.1 PageRank 方程与幂迭代原理从转移矩阵到不动点PageRank 的核心思想是一个网页的权重由所有指向它的网页权重按出度均摊后累加而来。写成矩阵形式就是R (1 - α) / N · e α · Mᵀ · R其中 R 是长度为 N 的 PageRank 向量e 是全 1 向量α 是阻尼系数通常取 0.85M 是行归一化后的转移矩阵。Mᵀ 的第 j 列描述了所有指向节点 j 的贡献来源矩阵向量乘Mᵀ R就是一次完整的权重传播。幂迭代的原理很简单从任意初始向量 R₀ 出发反复执行 Rₖ₊₁ (1-α)/N · e α · Mᵀ · Rₖ理论上当 k 趋近无穷时 R 收敛到唯一不动点。原因在于 Mᵀ 的谱半径小于 1乘以 α 后更小保证迭代的压缩性。但工程实现里有个绕不开的问题如果某个节点出度为 0它的行在 M 中是全零转置后对应列也是全零这部分权重在Mᵀ R中会完全丢失。标准解法是假设 dangling 节点会以 1/N 的概率跳转到所有节点也就是把丢失的权重均匀补回每个节点。3.2 稀疏矩阵版幂迭代代码dangling teleport 的补全逻辑from scipy.sparse import spdiags def pagerank_power(A, alpha0.85, tol1e-8, max_iter200): N A.shape[0] out_degree np.diff(A.indptr).astype(np.float64) dangling out_degree 0 inv_d np.where(dangling, 0.0, 1.0 / out_degree) Dinv spdiags(inv_d, 0, N, N) M (Dinv A).tocsr() # 行归一化转移矩阵dangling 行为全零 Mt M.T.tocsr() r np.full(N, 1.0 / N) teleport (1.0 - alpha) / N for it in range(max_iter): r_new alpha * (Mt r) teleport r_new alpha * r[dangling].sum() / N err np.abs(r_new - r).sum() r r_new if err tol: print(fpower iteration converged: {it 1} iters) return r print(fwarning: reach max_iter{max_iter}) return r这段代码里Dinv是出度倒数构成的对角稀疏矩阵Dinv A的效果是把每一行除以该行出度完成归一化。dangling布尔数组标记出度为 0 的节点r[dangling].sum()统计当前迭代中这些节点持有的权重总和再除以 N 均摊回所有节点这就是 teleport 补全。参数上有三个地方值得推敲。第一是 α取 0.85 是经验值表示用户以 85% 概率沿超链接跳转、15% 概率随机打开新页面α 越接近 1收敛越慢但排名对链接结构的刻画越敏锐。第二是 tol我习惯取1e-8对应 L1 范数误差也就是两次迭代间所有节点权重变化绝对值之和的阈值。第三是max_iter对 web-Google 这个规模α0.85 时通常几十次迭代就能到1e-8设 200 是留足余量。需要特别注意的是误差范数选择。有人用 L2 范数或最大分量误差但 PageRank 的权重分布极度不均匀少数高权重节点会掩盖大量低权重节点的变化用 L1 能更真实反映整体变化量。我在实践中固定用np.abs(r_new - r).sum()三种方法统一用这个判据方便横向对比收敛速度。4. 方法二与方法三Gauss-Seidel 加速与稀疏线性求解器4.1 为什么幂迭代在大图上收敛慢谱隙与逐点更新幂迭代的收敛速度取决于 Mᵀ 的第二大特征值绝对值α0.85 时理论收敛率接近 0.85。虽然每次迭代矩阵乘是 O(nnz) 的但迭代次数动辄上百轮510 万条边的规模下每轮乘法的常数开销不小。更关键的是幂迭代是同步更新所有节点用上一轮的旧值计算新值信息传播速度受限。Gauss-Seidel 的改进在于逐点更新计算第 j 个节点时引用它入边来源节点中已经算好的本轮新值而不是全部等待上一轮旧值。这样信息在同一轮内就能沿着链接传播得更远通常能显著减少迭代轮数。代价是更新顺序影响收敛路径且依赖对矩阵列的随机访问需要用 CSC 格式配合。4.2 Gauss-Seidel 的 Python 实现CSC 格式的逐列访问def pagerank_gs(A, alpha0.85, tol1e-8, max_iter200): N A.shape[0] out_degree np.diff(A.indptr).astype(np.float64) dangling_idx np.where(out_degree 0)[0] inv_d np.where(out_degree 0, 0.0, 1.0 / out_degree) Dinv spdiags(inv_d, 0, N, N) M (Dinv A).tocsr() Mcsc M.T.tocsc() # 转成 CSC按列取得入边贡献 r np.full(N, 1.0 / N) teleport (1.0 - alpha) / N for it in range(max_iter): r_old r.copy() dangling_sum r[dangling_idx].sum() base teleport alpha * dangling_sum / N for j in range(N): st, en Mcsc.indptr[j], Mcsc.indptr[j 1] if st en: r[j] base else: contrib Mcsc.data[st:en] r[Mcsc.indices[st:en]] r[j] alpha * contrib base err np.abs(r - r_old).sum() if err tol: print(fGS converged: {it 1} iters) return r return r这里把M.T转成 CSC 是最关键的一步。CSC 与 CSR 互为转置Mcsc.indptr[j]到Mcsc.indptr[j1]恰好是 Mᵀ 第 j 列的所有非零元素也就是所有指向 j 的入边来源及其权重。Mcsc.data[st:en] r[Mcsc.indices[st:en]]用一次点积完成入边权重求和因为 GS 是逐点更新r里来源节点既有本轮已更新的新值也有尚未更新的旧值这正是 Gauss-Seidel 与幂迭代的本质区别。dangling_sum必须在每轮循环开始前用当前r算一次并固定下来这是因为 dangling 节点的权重补全是一个全局均摊操作如果每更新一个节点就重新求和会破坏 GS 的逐点推进结构收敛行为反而不可控。这一点也是我在调参时踩过的坑后面避坑章节再细说。4.3 方法三把 PageRank 写成线性方程组交给 BiCGSTAB第三种思路是把 PageRank 从不动点迭代改写成线性方程组。将转移矩阵 M 的 dangling 行补成均匀分布 1/N得到随机矩阵 M此时 PageRank 是以下方程的解(I - α · Mᵀ) · R (1-α) / N · e直接构造 M 会遇到困难dangling 行全改为非零会让稀疏矩阵变成稠密填充。但 M M s · eᵀ / N其中 s 是 dangling 指示向量这是一个秩 1 修正因此可以在矩阵向量乘法中实现无需真正填充。BiCGSTAB 只需要提供 matvec 函数就能运行。from scipy.sparse.linalg import LinearOperator, bicgstab def pagerank_solve(A, alpha0.85, tol1e-10, max_iter500): N A.shape[0] out_degree np.diff(A.indptr).astype(np.float64) dangling out_degree 0 inv_d np.where(dangling, 0.0, 1.0 / out_degree) Dinv spdiags(inv_d, 0, N, N) M (Dinv A).tocsr() Mt M.T.tocsr() def matvec(v): out v - alpha * (Mt v) out - alpha * (v[dangling].sum() / N) return out Lop LinearOperator(shape(N, N), matvecmatvec, dtypenp.float64) b np.full(N, (1.0 - alpha) / N) x, info bicgstab(Lop, b, atoltol, maxitermax_iter) if info ! 0: print(fbicgstab did not converge: info{info}) return xmatvec对应的是 (I - αMᵀ) 乘以向量 v 的结果拆成三步先保留 v减去 αMᵀv 得到基础部分再减去 α 乘以 v 在 dangling 节点上的权重和均摊值这个均摊值就是秩 1 修正项。LinearOperator接收这个函数后表现为一个隐式矩阵BiCGSTAB 全程只调用matvec避免了显式构造 87 万阶满矩阵。bicgstab的atol我取1e-10而不是1e-8因为迭代求解器的误差度量与固定点迭代不同Krylov 子空间方法对残差的要求更严格。max_iter500是安全上限实际对 web-Google 数据通常在几十轮内收敛。有一点必须提醒把 dangling 修正直接加在右侧常数向量 b 上是错的因为修正项 α · (dangling_sum) / N 依赖未知解本身只能放进矩阵或者用迭代内更新处理否则解出来权重总和对不上。5. 单机 PageRank 避坑数据、参数与收敛的常见翻车点5.1 迭代次数爆炸收敛阈值不是越小越好现象把tol设成1e-12后幂迭代跑了两百多轮还没触发收敛日志里误差一直在缓慢下降单次实验耗时翻了好几倍。原因PageRank 的收敛速率由 α 和矩阵谱性质决定α0.85 时每轮误差大约压缩到原来的 0.85 倍从1e-8压到1e-12需要额外几十轮。更隐蔽的是浮点累加误差会让误差降到1e-13附近后不再下降甚至抖动阈值设得太小只会白耗算力。解决除非做学术对比实验否则tol1e-8足够。判断标准不应该是误差绝对值而是 top 排名是否稳定——可以每十轮对比一次当前 top 100 节点集合连续两轮完全不变就视为收敛这样省下的迭代次数非常可观。5.2 排名总和对不上dangling 节点把概率泄漏了现象三种方法跑完sum(r)有的接近 0.85有的接近 1.05没有一个稳定在 1.0 附近。排名结果和公开参考值差得很远。原因幂迭代和 GS 版本如果忘记补alpha * r[dangling].sum() / N每轮迭代都会丢失 dangling 节点持有的权重求解器版本如果错误地把 dangling 修正塞进右侧常数向量线性方程组本身就不等价于 PageRank 方程。解决无论哪种方法跑完第一步永远是np.abs(r.sum() - 1.0)检查。如果小于1e-6说明权重守恒没问题如果不满足优先检查 dangling 处理逻辑而不是怀疑浮点误差。补全项里分母是 N 而不是非 dangling 节点数这是最容易写错的地方。5.3 转置操作吃内存CSR 转 CSR 的隐形副本现象执行A.T.tocsr()时程序占用的内存瞬间上涨了 400 到 500MB机器接近卡死有一两次直接 OOM。原因A.T会先返回一个转置视图但调用tocsr()时 scipy 需要重新排列indices和data中间会生成多份临时数组。510 万条边在转置过程中峰值内存大约是稳态的三倍87 万阶矩阵放大效应明显。解决尽量只在预处理阶段做一次转置并复用结果不要在每个函数里反复T.tocsr()。另外可以先A.T.tocsc()再继续使用CSC 格式下按列访问本就自然省去一次显式 CSR 转换。5.4 GS 方法比幂迭代还慢Python 层循环的代价现象第一次跑pagerank_gs87 万节点的 for 循环一层套一层单次迭代耗时超过两分钟两个小时后还没收敛果断放弃。原因Gauss-Seidel 的逐点更新天然是串行的理论上能减少迭代轮数但纯 Python 循环每轮要处理 87 万个节点、500 万个非零元素解释器开销远大于稀疏矩阵乘法本身。幂迭代虽然迭代轮数多但核心计算是 C 语言实现的Mt r单轮耗时极短。解决如果非要用 GS优先考虑 numba/Cython 对这个循环做加速或者在数据规模较小十万节点以下的场景使用。对 web-Google 这个量级我更推荐把 GS 当作理解串行更新的教学工具工程落地上幂迭代和 BiCGSTAB 更实际。5.5 和 NetworkX 的结果对不上先对齐条件再对比现象用networkx.pagerank在同一个图上跑得到的 top 节点排名和手写实现差了十多名一看数值就以为代码有 bug。原因NetworkX 的默认参数和实现细节不同包括默认tol1e-6、max_iter100且它对 dangling 节点的处理方式是统一指向所有节点看起来正确但和手写补全的数学形式不完全一致。解决对比前强制统一四项α 相同、收敛阈值相同、dangling 处理方式相同、最大迭代次数相同。另外 NetworkX 内部存储的是 Python dict87 万节点上性能很差只建议用小规模子图做交叉验证不要用它对全量数据跑排名。6. 跑完别急着交差三招验证 PageRank 结果6.1 检查权重总和与 top 节点分布任何一种方法跑完第一件事是打印权重总和与 top 10 节点。以幂迭代结果为例我习惯一次性输出三项内容def sanity_check(r, A, top_k10): total r.sum() print(fsum(r) {total:.10f}) out_degree np.diff(A.indptr).astype(np.float64) top_idx np.argsort(r)[::-1][:top_k] for idx in top_idx: print(fnode {idx}: rank{r[idx]:.8f}, out_degree{int(out_degree[idx])})sum(r)偏离 1.0 超过1e-6就该回去查 dangling 补全top 节点里如果出现大量出度为零的节点说明 teleport 补全过度把低质量节点抬了上来。这个检查五秒内能完成强烈建议每次跑完都做。6.2 用子图交叉验证三种方法的一致性全量数据上很难手算验证单个节点的 PageRank 值但可以抽一个 5000 到 10000 节点的连通子图分别用三种方法求解再用 Spearman 排名相关系数对比整体排名是否一致。子图小三种方法都秒级完成。def rank_overlap(r1, r2, k500): top1 set(np.argsort(r1)[::-1][:k]) top2 set(np.argsort(r2)[::-1][:k]) return len(top1 top2) / k如果重叠率低于 0.9优先检查三种方法对 dangling 节点的处理是否等价。幂迭代和 GS 的补全项都是alpha * r[dangling].sum() / NBiCGSTAB 是隐式低秩修正数学上三者等价但数值行为会因迭代顺序和解算器的截面不同而略有差异重叠率在 0.95 以上基本可以放心。6.3 收敛曲线作为黑匣子的透视镜最后一招是把每轮误差画出来比任何日志都直观。幂迭代的误差曲线如果出现平台期说明 α 太接近 1 或者图中存在强连通集团GS 的误差曲线如果剧烈震荡说明 dangling_sum 在轮内被错误地更新了多次BiCGSTAB 的残差曲线如果先降后升多半是max_iter设得太小被提前截断。# 在迭代函数里收集 err_seq跑完直接 plot # err_seq.append(err) # import matplotlib.pyplot as plt # plt.semilogy(err_seq); plt.xlabel(iter); plt.ylabel(L1 err)从那以后我每次跑 PageRank 都强制自己走一遍这三步先查总和再抽子图做交叉验证最后扫一眼收敛曲线。尤其是收敛曲线它能把“算得慢”和“算错了”区分开省去大量无意义的调参。希望这篇笔记能帮你少踩几个我当年踩过的坑把时间留给真正有价值的调优和实验上。本文还有配套的精品资源点击获取