简介以谷歌网页数据集为样例这套资源面向大数据与图算法学习者提供PageRank计算的三种完整实现。三种方法分别包含基于Python稀疏矩阵手写迭代求解、直接调用networkx库内置pagerank函数以及在networkx基础上自行实现pagerank逻辑的定制版本覆盖从底层算法到上层接口的完整链路可帮助读者理解稀疏矩阵存储、幂迭代收敛和库函数的封装思想。资源包为RAR压缩格式共4个文件以3个Python脚本和1个数据压缩包为主整体大小约19.38MB。Python脚本可直接运行复现单机环境的PageRank计算压缩包内为待分析的网页链接数据便于二次扩展实验。已有824人学习资源特别适合搜索引擎排序、复杂网络分析及课程设计场景通过三种代码对比可快速掌握不同实现方式的性能差异与输出结果对初学者和进阶开发者均有参考价值。1. 三种方法对 web-Google.txt 做 PageRank先回答为什么必须用稀疏矩阵web-Google.txt 这份公开的 Google 网页链接图数据集大约有 87 万个网页节点、510 万条超链接边。如果按稠密矩阵去算 PageRank87 万乘 87 万的 double 数组需要约 6 TB 内存单机根本没戏但同一个图用稀疏矩阵存储只有 70MB 左右——这就是单机计算 PageRank 的唯一出路。本文用 Python 的 scipy.sparse 分别走幂迭代、特征向量、稀疏线性方程组三条路线把这份数据在普通笔记本上跑出 PageRank 值并对比三者的实现难度、收敛速度和坑。适合正在做链接分析、图算法实验或搜索引擎数据处理的开发者照着代码改一改就能复用到自己的边列表数据上。2. 准备数据解析 web-Google.txt、压缩编号并构建稀疏转移矩阵2.1 web-Google.txt 的真实格式edge list、方向与规模这份文件的每一行是两个整数前一个是出链网页的 ID后一个是它链接到的目标网页 ID行尾没有表头个别版本开头可能带#注释行。方向语义对 PageRank 很关键第 i 行u v表示网页 u 有一条指向 v 的链接那么 u 的权重会投给 v所以转移矩阵的每一行代表一个源节点每一列代表被链接的目标节点。文件规模方面常见统计是约 87.5 万个不同节点、510 万条边。这里的陷阱是节点 ID 并不是连续的 0 到 87 万原始编号最大值可能接近 92 万中间有很多空洞。直接用max(node_id) 1当矩阵维度等于多算了几万个不存在的网页而且这些空洞行还会被误判成悬挂节点把后面所有 PR 值都稀释掉。所以第一步永远是压缩编号。2.2 压缩编号把 91 万 ID 空间压缩到 87.5 万个有效索引读取文件我用np.loadtxt510 万行大约 80MB一次性读进内存没问题。为了压缩编号先np.unique拿到所有出现过的节点再做一张id - 0..N-1的映射表。这里不建议用 Python dict 做逐条映射87 万个 key 的哈希开销不小用数组索引更快。import numpy as np import scipy.sparse as sp # 读取边列表automatic 跳过 # 注释行 edges np.loadtxt(web-Google.txt, dtypenp.int32, comments#) src edges[:, 0].astype(np.int64) dst edges[:, 1].astype(np.int64) # 合并所有节点并压缩编号uniq[i] 就是新编号 i 对应的原始 ID uniq np.unique(np.concatenate([src, dst])) N len(uniq) max_id max(int(src.max()), int(dst.max())) # 用数组做 ID 映射比 dict 快 id2idx np.full(max_id 1, -1, dtypenp.int64) id2idx[uniq] np.arange(N) row id2idx[src] # 源节点新编号 col id2idx[dst] # 目标节点新编号id2idx长度是原始最大 ID 加 1不到 100 万占内存可以忽略。uniq是排序后的节点集合id2idx[uniq] np.arange(N)这一行把所有出现过的节点映射到连续的 0..N-1。注意src、dst显式转成int64避免后面和索引数组运算时发生整数溢出文件里的 ID 本身用int32读就够但索引运算统一成 64 位更安全。2.3 用 COO 组装稀疏矩阵并完成行归一化与悬挂节点标记有了row和col直接构造 COO 矩阵再转 CSR。CSR 的行切片和矩阵乘向量都高效适合后面的幂迭代和线性方程组求解。归一化时我不用np.diff(P.indptr)算行非零个数而是用P.sum(axis1)取每行真实权重这样即使数据里有重复边也不会归一化错。data np.ones(len(row), dtypenp.float64) P sp.coo_matrix((data, (row, col)), shape(N, N)).tocsr() # 行权重 每行的实际总和能兼容重复边 row_sum np.asarray(P.sum(axis1)).ravel() # 把每行的非零元素除以该行总和得到行随机矩阵 indptr P.indptr weights np.repeat(row_sum, np.diff(indptr)) P.data P.data / weights # 悬挂节点出度为 0 的行 dangling np.asarray(row_sum 0).reshape(-1)np.repeat(row_sum, np.diff(indptr))的含义是每个非零元素对应它所在行把该行的row_sum按非零个数重复对齐一次性完成整行归一化不需要 python 循环。dangling是一个布尔数组后面三种方法都要用它来把悬挂节点的权重均匀回灌到所有节点。到这一步稀疏转移矩阵 P 就绪约 70MB 内存主节点数据准备完毕。2.4 三个必调参数damping、tol、max_iterPageRank 的三个核心参数会在三种方法里反复出现提前说清楚各自的取值习惯。参数常见取值作用调参提示ddamping factor0.85控制随机跳转概率标准值 0.85调大接近 1 时收敛变慢悬挂节点风险放大tol1e-8 ~ 1e-10迭代停止的残差阈值幂迭代一般 1e-8 够用特征向量法 1e-9 起步max_iter幂迭代 100~200eigs 300防止死循环收敛慢先看 d 和悬挂节点不要盲目加迭代d 取 0.85 是搜索场景验证过的经验值它保证矩阵满足不可约非周期的随机性质主特征值唯一。tol 太小没有意义因为 PageRank 的实用精度到 1e-8 已经能保证 top N 排名稳定eigs 对 tol 更敏感太严会直接报不收敛。max_iter 在幂迭代里主要起保险作用正常几十轮就到 tol 了。3. 方法一与方法二幂迭代和特征向量法两条迭代路线的实现与选型3.1 幂迭代的最小实现P.T 乘向量 悬挂节点均匀回灌幂迭代就是把不动点方程r d * (P.T r) (1-d)/N反复迭代。P 是行随机矩阵P.T r的意思是每个节点收集所有指向它的源节点的投票。悬挂节点那一行全零如果不处理这部分权重就凭空消失所以要把r[dangling].sum()均匀回灌给全部 N 个节点。def pagerank_power(P, d0.85, tol1e-8, max_iter200): N P.shape[0] dangling (np.diff(P.indptr) 0) r np.full(N, 1.0 / N) for it in range(max_iter): r_next d * (P.T r) # 悬挂节点回灌把 lost 的权重均分给所有节点 r_next d * r[dangling].sum() / N # 随机跳转项保证全图连通数学上不可约 r_next (1 - d) / N err np.abs(r_next - r).sum() r r_next if err tol: break return r, it 1这段代码的收敛判据是 L1 残差也就是两轮迭代向量差的绝对值之和。用 L1 而不是 L2是因为 PR 值总和恒为 1L1 残差能直接反映有多少比例的权重还在流动。P.T r是稀疏矩阵乘稠密向量scipy 内部用 CSR 的转置逻辑高效执行一轮大约几十毫秒几十轮收敛这一条路线在三种方法里速度最快、最稳。参数上d不建议超过 0.9否则特征值间隔变小迭代轮数会明显增加tol从 1e-8 开始试如果只需要 top 100 排名1e-7 都能接受。返回的r自动满足求和为 1不需要额外归一化。3.2 特征向量法用 LinearOperator 表达 PageRank 矩阵避免构造补全矩阵第二种方法从特征分解视角出发PageRank 的稳态向量其实是矩阵A d*M^T (1-d)/N * 1*1^T的主特征向量其中 M 是补全了悬挂行后的行随机矩阵。理论上直接求 A 的最大特征值对应的特征向量即可。但这里有个工程陷阱A 是 87 万乘 87 万的稠密补全悬挂列要填均匀分布显式构造这个矩阵会瞬间打爆内存。正确做法是用LinearOperator只定义矩阵乘向量的行为不真正存放 A。from scipy.sparse.linalg import eigs, LinearOperator def matvec_A(x): y d * (P.T x) # 悬挂节点列的均匀重分配每个悬挂列等价于全 1/N y d * x[dangling].sum() / N # 随机跳转项全 1 矩阵外积的作用 y (1 - d) * x.sum() / N return y A_op LinearOperator((N, N), matvecmatvec_A, dtypenp.float64) # 求模最大的特征值对应特征向量 v0 np.random.default_rng(42).random(N) vals, vecs eigs(A_op, k1, whichLM, v0v0, maxiter300, tol1e-9) pr_eig vecs[:, 0].real pr_eig pr_eig / pr_eig.sum()matvec_A把一次矩阵乘向量拆成了三部分普通链接投票、悬挂节点均匀回灌、随机跳转。eigs的whichLM表示找模最大的特征值PageRank 矩阵的谱半径是 1其他特征值模不超过 d所以 LM 一定能抓到主特征值。这里不要用sigma1.0的 shift-invert 模式它会在内部做稀疏 LU 分解87 万维的 LU 填充量足以让单机卡死。v0必须给一个非零随机向量全零初值会让 ARPACK 直接失败固定随机种子便于复现结果。返回的特征值是复数取实部后归一化即可因为主特征值在实数轴上。3.3 何时选幂迭代何时选特征向量法特征向量法的数学表达很优雅但工程上它是三条路线里最娇贵的一条。ARPACK 对非对称矩阵的收敛行为不算友好87 万维下每步要做正交化实际耗时比幂迭代高一到两个数量级。它的价值在于验证当你怀疑幂迭代的悬挂节点处理写错了时用特征向量法算一次结果做交叉检查比肉眼盯排名可靠得多。对比项幂迭代法特征向量法实现难度低10 行核心逻辑中需要理解 LinearOperator每轮开销一次稀疏乘向量同样乘向量 ARPACK 正交化调参敏感度低d/tol 即可高v0/maxiter/tol 都要照顾典型耗时几十轮秒级到分钟级分钟级起视图规模而定适用场景生产环境首选结果验证、教学演示我的习惯是线上跑就用幂迭代它不挑数据、不挑参数特征向量法只在小图或抽样图上做交叉验证。如果你在做论文复现特征向量法能直接对应到线性代数教材的表述写进文档里更易读但别拿它在 87 万节点的完整图上反复调参那是给自己找麻烦。4. 方法三把 PageRank 转成稀疏线性方程组用 bicgstab 单机求解4.1 从不动点方程推导 (I - dM^T)r b第三种方法是把 PageRank 从迭代问题变成线性方程组。原来自不动点方程r d*M^T*r (1-d)/N * 1把含 r 的项移到左边就得到(I - d*M^T) * r (1-d)/N * 1。这里的 M 是补全悬挂行后的行随机矩阵和特征向量法里用的 M 是同一个对象。右边列向量 b 的每个分量都是常数(1-d)/N。这个形式的好处是解线性方程组得到的是精确解在数值误差范围内不受迭代初值影响也不存在收敛了几十轮但不知道够不够准的问题。代价是系数矩阵K I - d*M^T是一个 87 万维的非对称稀疏矩阵而且因为悬挂列的补全K 同样不容易直接以 CSR 形式高效构造——和特征向量法一样的工程门槛。解法上直接高斯消元不现实迭代求解器才是正路。4.2 求解bicgstab 配合 LinearOperator不显式构造系数矩阵bicgstab是稳定双共轭梯度法专门处理非对称稀疏系统。它只需要矩阵乘向量的能力不需要矩阵本身存量正好用 LinearOperator 承接 4.1 里的 K 算子。注意不能用cg因为 K 非对称gmres也能跑但 bicgstab 在这类图上通常更省内存。from scipy.sparse.linalg import bicgstab, LinearOperator def matvec_K(x): y x - d * (P.T x) # 悬挂列M^T 会从这些列均匀分配权重所以是减一个常数 y - d * x[dangling].sum() / N return y K_op LinearOperator((N, N), matvecmatvec_K, dtypenp.float64) b np.full(N, (1 - d) / N) # 初值用均匀分布bicgstab 对初值不敏感 x0 np.full(N, 1.0 / N) pr_lin, info bicgstab(K_op, b, x0x0, rtol1e-8, maxiter1000) if info ! 0: print(fbicgstab 未收敛info{info}需要调整预条件或放宽 rtol) else: pr_lin pr_lin / pr_lin.sum()matvec_K的核心是x - d*(P.Tx)然后再减去悬挂节点贡献的均匀项。注意第二行y - d * x[dangling].sum() / N是一个标量减到向量上Python 会自动广播到每个分量因为悬挂列对每个目标节点的贡献相等。rtol1e-8控制相对残差收敛maxiter1000是保险丝正常几百步内应收敛。info是求解器返回的标志0表示收敛非零值要排查是矩阵算子写错还是系统本身病态。解出来的pr_lin理论上总和为 1但迭代误差会破坏这个性质所以统一归一化后再用。4.3 spsolve 直接法为什么在这条数据上要谨慎很多人看到线性方程组第一反应是scipy.sparse.linalg.spsolve也就是用 SuperLU 做稀疏 LU 分解。在小图上这没错但 87 万维的 K 矩阵进行 LU 分解时填充元fill-in会把原本几十 MB 的稀疏矩阵撑大几十倍单机 16GB 内存很可能直接 OOM而且分解时间以分钟甚至小时计。这不是代码写错是直接法在超大规模稀疏问题上固有的空间爆炸。所以方法三的落地我推荐bicgstab它每步只有一次K_op矩阵乘向量内存峰值稳定在几个向量加一个稀疏 P 上。如果实测迭代不收敛优先检查悬挂节点是否正确标记再考虑用spilu做不完全 LU 预条件——但预条件本身也有内存成本在这份数据上是最后手段。三种方法对比来看幂迭代最快、最稳bicgstab 提供了「精确解」的数学保证适合做结果校对spsolve 只适合节点数在 10 万级以下的小图。5. 四个高频排查与避坑编号空洞、悬挂节点、eigs 不收敛与内存峰值5.1 坑 1节点编号不连续矩阵维度虚高近 5%现象跑出来的 PR 值整体偏小排名靠后的节点数值接近 0而且幂迭代收敛轮数比预期多。原因很直接直接把max(node_id)当成 N 建矩阵原始编号最大值约 92 万有效节点只有 87.5 万凭空多出 4.5 万个全零行。这些行在dangling判断里全部是真悬挂节点均匀回灌的目标变多每个节点的 PR 都被稀释。解决压缩编号。对我来说最有效的做法是 2.2 里的uniq id2idx数组映射一行代码拿到真实 N之后所有算子都基于压缩后的索引计算。顺带说一句压缩编号后np.unique的结果还能当 ID 反查表用最后输出 Top N 时直接uniq[top_idx]就能还原原始网页 ID。5.2 坑 2悬挂节点把 PR 值吸走尾部排名整体失真现象第一轮迭代后r.sum()明显小于 1多轮迭代后尾部几千个节点的 PR 值全部趋近相同的极小值。原因悬挂节点的行在 P 里是全零P.T r不会分发这些节点的权重如果不回灌这些权重就永久消失了。对于有大量低出度页面的真实图这会让尾部排名彻底失真——虽然 Top 10 变化不大但中间和尾部的相对序会被打乱。解决在幂迭代里补上r_next d * r[dangling].sum() / N在特征向量法里补上y d * x[dangling].sum() / N在 bicgstab 的 K 算子中把这个贡献作为均匀减项写入。判断自己有没有踩坑最简单的办法是打印第一轮迭代后的r.sum()如果比 1 小超过 1e-3悬挂节点处理就有问题。5.3 坑 3eigs 返回的特征值不是 1或直接抛 ARPACK error现象eigs抛出RuntimeError: ARPACK error或者返回的特征值实部明显偏离 1比如 0.8 附近。原因通常是三个初值v0给成了零向量或常量向量maxiter太小矩阵数据里混入了 NaN 或 inf。我当时排查了很久最后发现是np.loadtxt读文件时没有加comments#把注释行当数据读进去导致 P.data 里有脏值。解决固定随机种子生成v0maxiter300tol1e-9并在构造 P 之后马上执行np.isfinite(P.data).all()做一次数据健康检查。还要强调一个常见误用不要用sigma1.0去做 shift-invert那会在内部触发稀疏 LU 分解87 万维的分解直接让程序卡死在构造阶段连迭代的机会都没有。5.4 坑 4内存峰值来自 toarray()、LIL 列切片和过大的直接解现象程序跑着跑着内存占用突然暴涨甚至MemoryError。原因基本逃不出三个在某个环节调了P.toarray()或np.asarray(P)把 70MB 稀疏矩阵展开成 6TB 稠密数组用P.tolil()[:, dangling]改悬挂列LIL 的列切片会触发大量 Python 对象分配或者用spsolve直接解 K 矩阵LU 填充失控。解决全程保持 CSR 格式需要按列操作时就转一次P.T.tocsr()缓存起来矩阵乘向量一律交给 scipy 的稀疏算子线性系统用bicgstab而不是spsolve。想快速估算内存可以用P.data.nbytes P.indices.nbytes P.indptr.nbytes这个值在 80MB 左右是正常的如果哪个变量单独占了几个 GB一定是踩了上面某个坑。6. 验证三种方法的结果并输出 Top N迁移到更大图数据前的两个小技巧6.1 一致性验证三份 PR 向量两两做 L1 距离三种方法跑完后先别急着看排名做一次数值一致性检查。幂迭代和 bicgstab 的残差阈值都设在 1e-8 附近理论上两份 PR 向量的 L1 距离应该在 1e-6 量级特征向量法由于 ARPACK 的迭代精度误差会稍大但也不该超过 1e-4。如果偏差超标一定是悬挂节点处理在某个方法里写漏了。def l1_dist(r1, r2): r1 r1 / r1.sum() r2 r2 / r2.sum() return float(np.abs(r1 - r2).sum()) print(power vs eig:, l1_dist(pr_power, pr_eig)) print(power vs lin:, l1_dist(pr_power, pr_lin))这里做归一化再比是为了消除各方法对总和数值的微小差异。L1 距离在 1e-5 以内说明三套逻辑彼此印证可以放心使用。最后输出 Top N 时用np.argsort(pr_power)[::-1][:10]拿到压缩编号再映射回uniq[top_idx]还原原始网页 ID。注意 web-Google.txt 本身只有节点 ID 没有 URL想展示成网址需要自己准备一份 ID 到 URL 的映射表。6.2 迁移到更大图数据时的两个小技巧第一个技巧是缓存转置。P.T x这个算子在三套方法里被反复调用但每次P.T都会重新生成 CSR 对象。我一般在数据准备阶段直接算好Pt P.T.tocsr()后面所有 matvec 都用Pt x能省掉重复构建转置的时间。第二个技巧是按方法降级数据集换成更大规模时选型顺序是幂迭代 bicgstab eigs前者对内存和参数的容错性最好后者最娇贵。我第一版跑这套管线时N 直接取了原始最大 ID结果内存瞬间爆掉后来压缩编号、缓存转置才在一台普通笔记本上稳定复现。希望帮到你。本文还有配套的精品资源点击获取