图算法是内存受限的
不规则的内存访问模式、低算术强度使图算法天然受限于内存系统。优化局部性和减少随机访问是提升性能的关键。
CS267 · 第 22 课
图是表示实体间关系的通用抽象。从道路网络到基因组组装、从互联网拓扑到社交网络，图算法是高性能数据分析的核心。本课深入探讨如何为并行计算机设计高效的图算法，并分析五种经典问题的并行化策略。
图是实体（顶点）及实体间关系（边）的抽象表示。形式化地，定义图 \(G = (V, E)\)，其中 \(V\) 是顶点集合，\(E\) 是边集合。我们用以下符号刻画图的规模：
| 符号 | 含义 |
|---|---|
| \(n = |V|\) | 顶点数量 |
| \(m = |E|\) | 边数量 |
| \(D\) | 直径（任意两顶点间最大跳数） |
边可以是有向或无向的，可以带权（有权重）或不带权。边甚至可以携带属性（语义图）。由边序列 \(\langle u_1,u_2\rangle, \langle u_2,u_3\rangle, \ldots, \langle u_{n-1},u_n\rangle\) 构成的序列称为从 \(u_1\) 到 \(u_n\) 的路径，其长度等于路径上所有边权重之和。
注意：图算法的并行化面临独特挑战：局部性差（不规则内存访问模式）、算术强度低（每条边仅少量计算）、负载不均衡（scale-free 图的度分布极度不均）。这些特性使图算法天然受限于内存系统而非计算单元。
图抽象在科学计算和数据分析中无处不在，以下几类应用尤为突出：
点对点最短路径查询是导航系统的核心。Bast 等人（2007, Science）提出的 Transit Nodes 方法将朴素算法的 15 秒查询时间缩短至 10 微秒，提升超过 6 个数量级。
万维网天然表示为有向图（页面为顶点，超链接为边）。图算法支撑：
互联网拓扑（路由器网络）也自然建模为图。
图抽象对分析复杂数据集极为有用。数据来源包括模拟、实验设备、互联网和传感器网络。挑战包括数据规模、异构性、不确定性和数据质量。
Isomap（Tenenbaum 等, 2000, Science）利用图上测地距离（geodesic distances）保持数据的内在几何结构，实现非线性降维。所需工具包括 k-近邻图、全源最短路径（APSP）和 top-k 特征值计算。
科学计算中产生的图规模可能与要求解的稀疏线性系统或执行的模拟同量级。例如，二分图匹配可用于将矩阵变换为稠密对角或块三角形式；图划分用于并行模拟中的动态负载均衡。
PRAM（Parallel Random Access Machine）是 1980 年代提出的理想并行共享内存系统模型，是当时设计并行图算法的主要理论框架。
| 优点 | 缺点 |
|---|---|
| 语义简洁清晰，易于推理 | 不现实——通信模型过于强大 |
| 绝大多数理论并行算法基于 PRAM 设计 | 忽略通信开销 |
| 独立于通信网络拓扑 | 处理器完全同步，无局部内存 |
| 大 O 记号常具有误导性（忽略常数因子） |
尽管 PRAM 不切实际，但它为并行算法设计提供了理论基础。现代实践者在设计分布式内存图算法时，仍借鉴 PRAM 的分析框架，但需额外考虑通信开销和负载均衡。
CSR 是缓存高效的邻接表表示，由三个数组构成：
顶点: 1 2 3 4
row ptr: 1 3 7 12 19 (进入邻接数组的索引)
col id: 1 12 | 3 26 2 19 2 14 | 4 7 12 26 19 14 7
weight: ...对应权重...在分布式内存系统中，图数据需要划分到各处理器：
| 策略 | 描述 | 通信模式 |
|---|---|---|
| 全复制 | 每个处理器存储完整图 | 无通信，但内存受限 |
| 1D 分区 | 每个处理器存储 \(n/p\) 个顶点及其所有邻接出边 | 需要 all-to-all 通信交换邻接 |
| 2D 棋盘分区 | \(p_r \times p_c\) 处理器网格，每个处理器持有子矩阵（边块） | 行/列内通信，减少通信量 |
1D 分区中，可通过图划分算法（递归优化 conductance）或随机打乱顶点 ID 来均衡边数/处理器。
2D 棋盘分区将所需的 frontier 段缩小为 \(1/p_c\)（通常 \(1/\sqrt{p}\)），无需依赖网络内快速归约。这对 BFS 等遍历算法的通信复杂度有显著改善。
BFS 是图算法的基石，是（二分）匹配、最大流、（强）连通分量、介数中心性等算法的基础构建块。
procedure BFS(s):
queue Q ← {s}
dist[s] ← 0
while Q not empty:
v ← Q.dequeue()
for each neighbor w of v:
if dist[w] == ∞:
dist[w] ← dist[v] + 1
Q.enqueue(w)
内存需求（机器字数）：
稀疏图表示: m + n
访问顶点栈: n
距离数组: nDFS 的教训：深度优先搜索本质上是串行的（Reif, 1985），其跨度（span）为 \(O(n)\)。因此并行图算法通常基于 BFS 而非 DFS。
两种主要策略适用于不同直径的图：
算法步骤：
1. 确定当前 frontier 邻接的所有者 [计算]
2. 通过 all-to-all 交换邻接 [通信]
3. 更新未访问顶点的距离/父节点 [计算]算法步骤：
1. 收集处理器列中的顶点 [通信]
2. 确定当前 frontier 邻接的所有者 [计算]
3. 在处理器行内交换邻接 [通信]
4. 更新未访问顶点的距离/父节点 [计算]利用半环代数（min-plus 或 select-min），BFS 可表达为矩阵-向量乘法。将标量操作替换为：乘法→选择（select），加法→取最小值。
核心操作: x = A^T ⊗ x (在 min-plus 半环上)
其中 A^T 是邻接矩阵转置，x 是 frontier 向量
每轮迭代将 frontier 扩展一层经典 BFS 是自顶向下（Top-Down）的：对 frontier 中每个顶点的所有邻居尝试认领为子节点。但在 frontier 峰值时（如 Youtube 社交网络），几乎所有边的检查都"失败"——因为目标顶点已被访问。
Top-Down: for all v in frontier
attempt to parent all neighbors(v)
Bottom-Up: for all v in unvisited
find any parent (neighbor(v) in frontier)结合两者优势：frontier 小时用 Top-Down，frontier 大时切换为 Bottom-Up。切换启发式：当已发现顶点的多数时切换（详见 Beamer 等, SC 2012）。
2D 分区减少了 frontier 段（\(1/\sqrt{p}\)），但 bottom-up 内层循环是串行的（找到一个父节点后希望提前退出）。解决方案：
| 技术 | 描述 |
|---|---|
| 时间划分 | 每个顶点最多被一个处理器处理（避免重复检查） |
| 脉动旋转（Systolic Rotation） | 将完成信息发送给下一个处理器，使其知道要跳过哪些顶点 |
在 ORNL Titan（Cray XK6, Gemini 互连）上，对 160 亿顶点、2560 亿边的 Kronecker (Graph500) 图，方向优化 + 2D 分解实现了约 8 倍的加速。第一个量子跃迁来自 2D 算法的采用，第二个来自 bottom-up 搜索的引入。
在基因组组装中，需要遍历 de Bruijn 图寻找 UU contigs（UU 节点链）。每个处理器 Pi 从分布式哈希表中随机选取 k-mer 作为种子，通过连续查找向两端延伸，直到到达 UU contig 的端点。
多个处理器可能从同一 UU contig 选取初始种子，需要通过基于状态机的轻量级同步方案协作拼接子 contigs，避免冗余工作。
低直径图（如 R-MAT 合成网络）与长窄图（如基因组 de Bruijn 图）需要不同的遍历策略。实践者必须根据图的拓扑特征选择算法变体。
两种经典串行算法：
| 算法 | 类型 | 适用条件 | 串行复杂度 |
|---|---|---|---|
| Bellman-Ford | label correcting | 任意图（可处理负权） | \(O(nm)\) |
| Dijkstra | label setting | 非负权重 | \(O(m + n \log n)\) |
目前无已知 PRAM 算法能在亚线性时间和 \(O(m + n \log n)\) 工作量内解决 SSSP。Ullman-Yannakakis 提出了随机化方法。
Δ-stepping 是"近似的 Dijkstra 桶实现"——一种 label correcting 算法，也可从未定顶点松弛边。顶点按宽度为 Δ 的桶排序，每个桶可并行处理。
基本松弛操作: Relax(e(u,v)):
d(v) = min { d(v), d(u) + w(u, v) }
Δ 的选择:
Δ < min w(e) → 退化为 Dijkstra
Δ > max w(e) → 退化为 Bellman-Ford
并行阶段 (while bucket not empty):
i) 检查轻边 (w < Δ)
ii) 构建请求集合 R
iii) 清空当前桶
iv) 记录被删除的顶点 S
v) 松弛 R 中的请求对
松弛 S 中的重请求对
转到下一个桶考虑 7 个顶点的有向图，边权重归一化到 [0,1]，设 Δ = 0.1：
对于随机权重 [0,1] 的图，Δ-stepping 的复杂度为 \(O(n + m + D \times L)\)，其中 \(L\) 是从源点到任意节点的最大距离。
在高直径图（如网格图、道路网络）中，Δ-stepping 需要过多的阶段（phases）。Level-synchronous BFS 面临同样的困境。实际中，Chakaravarthy 等人（IPDPS 2014）将 Δ-stepping 与方向优化 BFS 结合，在超算规模图上实现了可扩展的 SSSP。
Floyd-Warshall 是经典的动态规划算法，使用三重嵌套循环。然而，之前被忽视的递归版本比三重嵌套循环更易并行化。
给定距离矩阵 D 分块为 [A B; C D]:
A* = solve(A) // 递归求解
B = B × A* // min-plus 矩阵乘
C = A* × C
D = D + C × B // min-plus: + 是 min, × 是加
D* = solve(D) // 递归求解
B = B × D*
C = D* × C
A = A + B × C
其中: + → min × → +这个分块递归形式天然适合并行化：四个子块操作可独立执行。
Buluç, Gilbert & Budak（2010）证明，将 Floyd-Warshall 移植到 GPU 并使用正确的基元（矩阵乘法），相比朴素递归实现可获得 480 倍加速。关键洞察：矩阵乘法是高度优化的 GPU 基元，min-plus 代数下的矩阵乘法可以复用这些优化。
Solomonik, Buluç & Demmel（2013）提出 2.5D 矩阵乘法方法用于 APSP：
2.5D 矩阵乘法:
复制因子 c（顶点数）
W_bc-2.5D(n, p) = O(n²/√(cp))
S_bc-2.5D(p) = O(c log²(p)) (latency)
对任意内存大小都最优！APSP 可归约为 min-plus 矩阵乘法半环上的矩阵幂运算。这一联系使得 APSP 能直接受益于矩阵乘法在 GPU 和分布式系统上的优化，是代数图算法（GraphBLAS）的核心思想。
独立集是图中互不相邻的顶点集合。最大（maximal）独立集无法再添加顶点而保持独立（局部最优）；最大（maximum）独立集是所有独立集中最大的（全局最优，NP-hard）。
1. S = 空集
2. for 顶点 v = 1 to n:
3. if (v 没有邻居在 S 中):
4. 将 v 加入 S
work ~ O(n), 但 span ~ O(n) —— 完全串行1. S = 空集; C = V
2. while C 不为空:
3. 为 C 中每个 v 标记随机数 r(v)
4. for all v in C (并行):
5. if r(v) < min(r(v 的邻居)):
6. 将 v 从 C 移入 S
7. 将 v 的邻居从 C 中移除
定理: 该算法"极大概率"在 O(log n) 轮内完成
work ~ O(n log n), span ~ O(log n)Luby 算法通过随机标记打破对称性，每轮以常数比例消除顶点。随机化使串行跨度 \(O(n)\) 降为对数跨度 \(O(\log n)\)，是并行算法设计的经典范例。GraphBLAS 提供了该算法的变体实现。
有向图中，强连通分量是极大子图，其中任意两顶点互相可达。SCC 分解可通过对称置换转化为块三角形式。
基于深度优先搜索，时间复杂度 \(O(m+n)\)。但 DFS 本质上是串行的（Reif, 1985），难以直接并行化。
核心思想：选取 pivot 顶点 \(u\)，利用前向/后向可达集将图分为三个独立子图。
定义:
FW(u): 从顶点 u 可达的顶点集合
BW(u): 能到达顶点 u 的顶点集合
引理: FW(u) ∩ BW(u) 是唯一的 SCC
对于其他任意 SCC s，要么:
(a) s ⊂ FW(u) \ BW(u) (从 u 可达但不可回)
(b) s ⊂ BW(u) \ FW(u) (可到达 u 但 u 不可达)
(c) s ⊂ V \ (FW(u) ∪ BW(u)) (互相不可达)
算法:
1. 随机选取 pivot u
2. 计算 FW(u) 和 BW(u)
3. 将图分为三个子图 (a)(b)(c)
4. 递归处理每个子图Hong, Rodia & Olukotun（SC 2013）观察到：真实世界图常有巨大的 SCC，串行 FW/BW 的跨度 \(O(N)\) 成为瓶颈。解决方案：
真实世界图（如 Web 图、社交网络）通常有一个巨大的"巨型 SCC"包含大部分顶点。这使得朴素的递归并行 SCC 算法效率低下，必须结合并行 BFS 和弱连通分量分解。
匹配是图中没有公共端点的边集。基数指匹配中边的数量。最大基数匹配是边数最多的匹配。
增广路径是起点和终点都是未匹配顶点的交错路径（匹配边/非匹配边交替）。翻转增广路径上的边状态，可使匹配基数增加 1。
1. 初始匹配（任意）
2. 从未匹配顶点 x3 搜索增广路径（交错 BFS）
3. 翻转增广路径上的边
4. 重复直到无增广路径同时从多个未匹配顶点搜索顶点不相交的增广路径，并行翻转。通过 BFS 森林（alternating BFS forest）实现。
Azad, Buluç & Pothen（IPDPS 2015）提出的创新：当从某个顶点出发的增广路径搜索失败时，不销毁整棵树，而是将其"嫁接"到相邻的搜索森林上，避免重复搜索。
算法步骤:
1. 并行方向优化 BFS (Beamer et al. SC 2012)
- frontier 大时使用 bottom-up BFS
- 维护 visited 数组
- 顶点不相交路径：每轮迭代中顶点只访问一次
- 使用线程安全的原子操作
2. 增广路径顶点不相交 → 可并行翻转
3. 未匹配顶点尝试"附着"到活跃顶点
- 无需同步树嫁接算法相比之前的 Pothen-Fan（Azad 等 IPDPS 2012）和 Push-Relabel（Langguth 等 2014）方法有显著提升，尤其在真实世界二分图上。关键在于避免重复搜索和高效的并行 BFS 基元。
| 问题 | 串行复杂度 | 并行复杂度 | 关键挑战 |
|---|---|---|---|
| BFS | \(O(m+n)\) | \(O(D)\) 步（每步并行） | frontier 大小变化、负载均衡 |
| SSSP (Δ-stepping) | \(O(m + n \log n)\) | \(O(n + m + D \times L)\) | 高直径图阶段过多 |
| APSP (Floyd-Warshall) | \(O(n^3)\) | \(O(n^3/p)\) + 通信 | 通信量、内存分布 |
| MIS (Luby) | \(O(n)\) | \(O(\log n)\) 轮 | 随机化、竞争条件 |
| SCC | \(O(m+n)\) | 最坏 \(O(n)\)，实际亚线性 | 巨型 SCC、递归深度 |
| Max Matching | \(O(nm)\) | 多轮增广 | 路径不相交、搜索森林 |
Q1：CSR 格式存储稀疏图时，不包括以下哪个数组？
Q2：Direction-optimizing BFS 在什么条件下从 Top-Down 切换到 Bottom-Up？
Q3：Δ-stepping 算法中，当 Δ 小于最小边权重时，退化为哪个算法？
Q4：Luby 随机化 MIS 算法的期望并行轮数（span）是多少？
Q5：Fleischer-Hendrickson-Pınar 并行 SCC 算法中，选取 pivot 顶点 u 后，图被分为几个独立子图？
不规则的内存访问模式、低算术强度使图算法天然受限于内存系统。优化局部性和减少随机访问是提升性能的关键。
根据问题规模动态调整算法方向（top-down ↔ bottom-up）是图算法性能优化的通用策略，不仅适用于 BFS，也适用于 SSSP。
将图算法表达为矩阵运算（min-plus 半环、GraphBLAS）能复用高度优化的矩阵乘法基元，在 GPU 上实现数量级的加速。
图的拓扑特征（直径、度分布、SCC 结构）深刻影响算法性能。实践者必须根据具体问题选择合适的算法变体。