元信息：Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, Clifford Stein《Introduction to Algorithms》（第 3 版，MIT Press 2009）。负责 PDF 第 667–883 页。页码换算：原书页码 = PDF 页码 − 21（例如 PDF p.667 = 原书 p.646）。本块从第 24 章引言中段（原书 p.646）开始，到第 29 章 29.2 节“多商品流”小节开头（原书 p.862）结束。记号约定：$G=(V,E)$ 为图，$w$ 为边权函数，$\delta(u,v)$ 为最短路径权重，$v.d$ 为最短路径估计，$v.\pi$ 为前驱。复杂度中的 $V,E$ 分别指 $|V|,|E|$。

第 24 章单源最短路径（Single-Source Shortest Paths）（接上一块）

第 24 章引言（后半）（PDF p.667–671）

（接上一块：本章引言前半部分——最短路径问题定义、路径权重 $w(p)=\sum w(v_{i-1},v_i)$、$\delta(u,v)$ 定义、问题变体、最优子结构、负权边——在上一块。）

负权边示例（图 24.1）：图中顶点 e、f 构成从 s 可达的负权环，所以 $\delta(s,e)=\delta(s,f)=-\infty$；g 从一个 $-\infty$ 顶点可达，也是 $-\infty$；h、i、j 从 s 不可达，$\delta=+\infty$，尽管它们自己在一个负权环上。Dijkstra 算法要求所有边权非负；Bellman-Ford 允许负权边，只要从源点不可达负权环就能给出正确答案，并能检测出可达负权环。

最短路径不含环（Cycles）：

不能含负权环（否则无最短路，权重为 $-\infty$）。
不能含正权环：若路径 $p=\langle v_0,\dots,v_k\rangle$ 上有环 $c=\langle v_i,\dots,v_j\rangle$（$v_i=v_j$，$w(c)>0$），删掉环得 $p'$，$w(p')=w(p)-w(c)<w(p)$，矛盾。
0 权环可以删去而不改变权重。
因此不失一般性，最短路径都是简单路径（simple path），至多含 $|V|$ 个不同顶点、$|V|-1$ 条边。这一事实是 Bellman-Ford 做 $|V|-1$ 轮的根据。

最短路径的表示：每个顶点维护前驱 $v.\pi$（另一顶点或 NIL）。沿前驱链回溯即得从 s 到 v 的最短路径（用 22.2 节的 PRINT-PATH）。前驱子图（predecessor subgraph） $G_\pi=(V_\pi,E_\pi)$：

V_\pi=\{v\in V: v.\pi\ne \text{NIL}\}\cup\{s\},\qquad E_\pi=\{(v.\pi,v)\in E: v\in V_\pi-\{s\}\}.

最短路径树（shortest-paths tree）：设 G 无从 s 可达的负权环。以 s 为根的最短路径树是有向子图 $G'=(V',E')$，$V'\subseteq V,E'\subseteq E$，满足：(1) $V'$ 是 G 中从 s 可达的顶点集；(2) $G'$ 是以 s 为根的有根树；(3) 对每个 $v\in V'$，$G'$ 中从 s 到 v 的唯一简单路径就是 G 中 s 到 v 的一条最短路径。它与 BFS 树类似，只是以边权而非边数度量。最短路径和最短路径树都不一定唯一（图 24.2 给出同一图、同一根的两棵不同最短路径树）。

松弛（relaxation）——本章所有算法的核心操作。每个顶点维护 $v.d$，是 $\delta(s,v)$ 的上界，称为最短路径估计（shortest-path estimate）。

def initialize_single_source(G, s):     # Θ(V)
    for v in G.V:
        v.d = INF
        v.pi = None
    s.d = 0

def relax(u, v, w):                      # O(1)
    if v.d > u.d + w(u, v):
        v.d = u.d + w(u, v)
        v.pi = u

图 24.3 例：$w(u,v)=2$。(a) $u.d=5, v.d=9$，松弛后 $v.d=7$；(b) $u.d=5,v.d=6$，$6\le 7$，不变。

名称来历（脚注）：“松弛”实为收紧上界，是历史用语——它可看作对约束 $v.d\le u.d+w(u,v)$ 的放松：若该约束已满足，就没有“压力”。

各算法的区别只在于松弛哪些边、多少次、什么顺序：Dijkstra 与 DAG 算法对每条边恰好松弛一次；Bellman-Ford 对每条边松弛 $|V|-1$ 次。

最短路径与松弛的性质（证明在 24.5 节；后五条假设先调用 INITIALIZE-SINGLE-SOURCE，之后只通过松弛改变 d 与 π）：

三角不等式（Triangle inequality，引理 24.10）：对任意边 $(u,v)\in E$，$\delta(s,v)\le\delta(s,u)+w(u,v)$。
上界性质（Upper-bound property，引理 24.11）：始终有 $v.d\ge\delta(s,v)$；一旦 $v.d$ 达到 $\delta(s,v)$ 就不再改变。
无路径性质（No-path property，推论 24.12）：若 s 到 v 无路径，则始终 $v.d=\delta(s,v)=\infty$。
收敛性质（Convergence property，引理 24.14）：若 $s\leadsto u\to v$ 是最短路径，且在松弛边 $(u,v)$ 之前的任意时刻 $u.d=\delta(s,u)$，则此后始终 $v.d=\delta(s,v)$。
路径松弛性质（Path-relaxation property，引理 24.15）：若 $p=\langle v_0,\dots,v_k\rangle$ 是 $s=v_0$ 到 $v_k$ 的最短路径，且按 $(v_0,v_1),(v_1,v_2),\dots,(v_{k-1},v_k)$ 的顺序松弛过这些边，则 $v_k.d=\delta(s,v_k)$。无论中间穿插了多少其他松弛都成立。
前驱子图性质（Predecessor-subgraph property，引理 24.17）：一旦对所有 v 有 $v.d=\delta(s,v)$，前驱子图就是以 s 为根的最短路径树。

章节安排：24.1 Bellman-Ford（一般情形，可检测负环）；24.2 DAG 上的线性时间算法；24.3 Dijkstra（更快，需非负权）；24.4 用 Bellman-Ford 解线性规划的特例（差分约束）；24.5 证明上述性质。

无穷大算术约定：对实数 $a\ne-\infty$，$a+\infty=\infty+a=\infty$；对 $a\ne\infty$，$a+(-\infty)=(-\infty)+a=-\infty$。图均以邻接表存储，边权随边存放，遍历时 O(1) 取权。

24.1 Bellman-Ford 算法（PDF p.672–676）

功能：一般情形（边权可负）的单源最短路径。返回布尔值：若存在从源点可达的负权环则返回 FALSE（无解）；否则返回 TRUE 并给出最短路径及权重。

def bellman_ford(G, w, s):
    initialize_single_source(G, s)           # Θ(V)
    for i in range(1, len(G.V)):             # |V|-1 轮
        for (u, v) in G.E:                   # 每轮 Θ(E)
            relax(u, v, w)
    for (u, v) in G.E:                       # O(E) 负环检测
        if v.d > u.d + w(u, v):
            return False
    return True

复杂度：时间 $O(VE)$（初始化 $\Theta(V)$，$|V|-1$ 轮每轮 $\Theta(E)$，检测 $O(E)$）；空间 $O(V)$（d、π 数组）。

例（图 24.4）：5 个顶点 s,t,x,y,z；每轮按 $(t,x),(t,y),(t,z),(x,t),(y,x),(y,z),(z,x),(z,s),(s,t),(s,y)$ 的顺序松弛。边权：$s\to t=6, s\to y=7, t\to x=5, t\to y=8, t\to z=-4, x\to t=-2, y\to x=-3, y\to z=9, z\to x=7, z\to s=2$。最终 $d(s,t,x,y,z)=(0,2,4,7,-2)$，返回 TRUE。四轮后各值依次为：第 1 轮 $t=6,y=7$；第 2 轮 $x=4,z=2$；第 3 轮 $t=2$；第 4 轮 $z=-2$。

引理 24.2：若 G 无从 s 可达的负权环，则 $|V|-1$ 轮循环后，对所有从 s 可达的 v，$v.d=\delta(s,v)$。证明：取 s 到 v 的最短路径 $p=\langle v_0,\dots,v_k\rangle$，它是简单路径，故 $k\le|V|-1$。第 i 轮松弛了所有边，其中包括 $(v_{i-1},v_i)$，由路径松弛性质得 $v.d=\delta(s,v)$。

推论 24.3：在同样假设下，s 到 v 有路径 ⇔ Bellman-Ford 结束时 $v.d<\infty$（证明为习题 24.1-2）。

定理 24.4（Bellman-Ford 正确性）：若 G 无从 s 可达的负权环，算法返回 TRUE，所有 $v.d=\delta(s,v)$，前驱子图是最短路径树；若有，返回 FALSE。证明要点：

无负环时：可达顶点由引理 24.2，不可达顶点由无路径性质，得 $v.d=\delta(s,v)$；再由前驱子图性质得最短路径树。对每条边，$v.d=\delta(s,v)\le\delta(s,u)+w(u,v)=u.d+w(u,v)$（三角不等式），检测不会触发，返回 TRUE。
有可达负环 $c=\langle v_0,\dots,v_k\rangle$，$v_0=v_k$，$\sum_{i=1}^k w(v_{i-1},v_i)<0$（式 24.1）。反设返回 TRUE，则 $v_i.d\le v_{i-1}.d+w(v_{i-1},v_i)$。沿环求和：$\sum v_i.d\le\sum v_{i-1}.d+\sum w(v_{i-1},v_i)$。由于 $v_0=v_k$，两个 d 之和相等，且由推论 24.3 都有限，消去得 $0\le\sum w(v_{i-1},v_i)$，与 (24.1) 矛盾。

直观理解：第 i 轮后，所有“边数 ≤ i 的最短路径”都已求对，相当于按“允许的边数”做动态规划。负环检测的含义是：若 $|V|-1$ 轮后仍能松弛，说明存在更长（边数 ≥ |V|）的更短路径，只能来自负环。

常见误区：Bellman-Ford 只能检测从源点可达的负环；若想检测全图负环，可加超级源点 $v_0$ 连到所有顶点（权 0），见 24.4 节和 25.3 节 Johnson 算法。

习题 24.1：

24.1-1 手工运行（以 z 为源；再把 $w(z,x)$ 改为 4 以 s 为源，此时出现负环）。
24.1-2 证明推论 24.3。
24.1-3 设 m 为所有顶点“最短路径中最少边数”的最大值，修改算法使其在 m+1 轮后终止（提示：某轮无任何 d 改变即提前退出）。
24.1-4 修改算法，使负环可达的所有顶点 $v.d=-\infty$。
24.1-5★ 求每个 v 的 $\delta^*(v)=\min_{u}\delta(u,v)$，要求 $O(VE)$。
24.1-6★ 存在负环时，列出一个负环的顶点并证明正确性（从最后一轮被松弛的顶点沿 π 回溯 |V| 步后必在环上）。

24.2 有向无环图中的单源最短路径（PDF p.676–679）

思想：对带权 DAG（directed acyclic graph） 按拓扑序松弛各顶点的出边，$\Theta(V+E)$ 时间求出单源最短路径。DAG 中即使有负权边也不存在负环，最短路径总有定义。

def dag_shortest_paths(G, w, s):
    order = topological_sort(G)        # Θ(V+E)，见 22.4 节
    initialize_single_source(G, s)     # Θ(V)
    for u in order:                    # 每个顶点一次
        for v in G.Adj[u]:             # 每条边恰好松弛一次（聚合分析）
            relax(u, v, w)

复杂度：时间 $\Theta(V+E)$，即邻接表规模的线性时间；空间 $O(V)$。

例（图 24.5）：顶点按拓扑序 r,s,t,x,y,z，源 s。边：$r\to s=5,r\to t=3,s\to t=2,s\to x=6,t\to x=7,t\to y=4,t\to z=2,x\to y=-1,x\to z=1,y\to z=-2$。最终 $d(r,s,t,x,y,z)=(\infty,0,2,6,5,3)$。

定理 24.5：若 G 无环、源点为 s，则 DAG-SHORTEST-PATHS 结束时所有 $v.d=\delta(s,v)$，且前驱子图是最短路径树。证明：不可达顶点由无路径性质；可达顶点取最短路径 $\langle v_0,\dots,v_k\rangle$，拓扑序保证依次松弛 $(v_0,v_1),\dots,(v_{k-1},v_k)$，由路径松弛性质成立；再用前驱子图性质。

应用：PERT 图中的关键路径（critical path）（PERT = program evaluation and review technique）。边代表任务，边权代表耗时；若边 $(u,v)$ 进入 v、$(v,x)$ 离开 v，则任务 $(u,v)$ 须先于 $(v,x)$。关键路径是 DAG 中的最长路径，其权重是完成全部任务所需总时间的下界。求法：(1) 把边权取负后运行 DAG-SHORTEST-PATHS；或 (2) 初始化时把 $\infty$ 换成 $-\infty$，RELAX 中把“>”换成“<”。注意：一般图的最长简单路径是 NP 难问题，只有 DAG 上才能这样线性求解。

习题 24.2：

24.2-1 以 r 为源运行。
24.2-2 证明只处理拓扑序前 $|V|-1$ 个顶点仍正确（最后一个顶点没有出边）。
24.2-3 顶点带权（任务在顶点上、边表示先后约束）时，线性时间求最长路径。
24.2-4 计数 DAG 中路径总数（按拓扑序做 DP）。

24.3 Dijkstra 算法（PDF p.679–685）

适用条件：所有边权非负，$w(u,v)\ge0$。好的实现比 Bellman-Ford 快。

思想：维护集合 S，其中顶点的最终最短路径权重已确定。反复从 $V-S$ 中选 d 值最小的顶点 u 加入 S，并松弛 u 的所有出边。用以 d 为键的最小优先队列 Q。这是贪心策略（greedy strategy）。

def dijkstra(G, w, s):
    initialize_single_source(G, s)
    S = set()
    Q = MinPriorityQueue(G.V, key=lambda v: v.d)   # 隐含 |V| 次 INSERT
    while Q:                                        # 恰好 |V| 次
        u = Q.extract_min()
        S.add(u)
        for v in G.Adj[u]:                          # 总计 |E| 次
            relax(u, v, w)                          # 隐含 DECREASE-KEY

循环不变式：每次 while 迭代开始时 $Q=V-S$。顶点只在第 3 行入队，每个顶点恰好出队并加入 S 一次，故 while 恰好执行 $|V|$ 次。第一次取出的是 s。

例（图 24.6）：顶点 s,t,x,y,z；边 $s\to t=10,s\to y=5,t\to x=1,t\to y=2,x\to z=4,y\to t=3,y\to x=9,y\to z=2,z\to s=7,z\to x=6$。依次取出 s(0)、y(5)、z(7)、t(8)、x(9)。中间过程：取 s 后 $t=10,y=5$；取 y 后 $t=8,x=14,z=7$；取 z 后 $x=13$；取 t 后 $x=9$。最终 $d=(s0,t8,x9,y5,z7)$。

定理 24.6（Dijkstra 正确性）：在非负权有向图上，Dijkstra 结束时对所有 u 有 $u.d=\delta(s,u)$。证明（循环不变式：每次迭代开始时，S 中所有 v 满足 $v.d=\delta(s,v)$）：

初始化：$S=\emptyset$，平凡成立。
保持：反设 u 是第一个加入 S 时 $u.d\ne\delta(s,u)$ 的顶点。$u\ne s$（s 第一个加入且 $s.d=0=\delta(s,s)$），因此此时 $S\ne\emptyset$。s 到 u 必有路径（否则由无路径性质 $u.d=\delta=\infty$），取最短路径 p。设 y 是 p 上第一个属于 $V-S$ 的顶点，x 为其在 p 上的前驱（$x\in S$），分解 $p: s\overset{p_1}{\leadsto}x\to y\overset{p_2}{\leadsto}u$（图 24.7）。
- 断言 $y.d=\delta(s,y)$：x 加入 S 时 $x.d=\delta(s,x)$（u 是第一个出错者），此时松弛了 $(x,y)$，由收敛性质得证。
- 因 y 在最短路径上位于 u 之前且边权非负（特别是 $p_2$ 上），$\delta(s,y)\le\delta(s,u)$，于是 $y.d=\delta(s,y)\le\delta(s,u)\le u.d$（上界性质）。（式 24.2）
- 但 u 与 y 都在 $V-S$ 中而算法选了 u，故 $u.d\le y.d$。两个不等式取等，$y.d=\delta(s,y)=\delta(s,u)=u.d$，与假设矛盾。
终止：$Q=\emptyset$，结合 $Q=V-S$ 得 $S=V$。

推论 24.7：结束时前驱子图是以 s 为根的最短路径树（定理 24.6 + 前驱子图性质）。

复杂度分析：INSERT 与 EXTRACT-MIN 各 $|V|$ 次；每个顶点只进 S 一次，所以每条邻接边只检查一次，DECREASE-KEY 至多 $|E|$ 次（聚合分析）。

优先队列实现	INSERT	EXTRACT-MIN	DECREASE-KEY	总时间
以顶点编号为下标的数组	O(1)	O(V)	O(1)	$O(V^2+E)=O(V^2)$
二叉最小堆（需维护顶点与堆元素的互相句柄）	建堆总 O(V)	O(lg V)	O(lg V)	$O((V+E)\lg V)$，全部可达时为 $O(E\lg V)$
斐波那契堆（第 19 章）	O(1)	摊还 O(lg V)	摊还 O(1)	$O(V\lg V+E)$

二叉堆在 $E=o(V^2/\lg V)$（稀疏图）时优于数组实现。历史上，斐波那契堆的提出正是因为 Dijkstra 中 DECREASE-KEY 调用远多于 EXTRACT-MIN，把前者降到 $o(\lg V)$ 摊还而不增加后者代价，就能渐近加速。空间 $O(V)$（不计图本身）。

与其他算法的关系：像 BFS——S 对应 BFS 中的黑色顶点，都已得到最终距离；像 Prim 算法——都用最小优先队列找集合外“最轻”的顶点，加入集合后调整其余顶点的键。

常见误区：存在负权边时 Dijkstra 可能出错（证明中 $\delta(s,y)\le\delta(s,u)$ 依赖非负权）；但如果只有从源点出发的边为负、且无负环，Dijkstra 仍正确（习题 24.3-10）。

习题 24.3：

24.3-1 在图 24.2 上分别以 s、z 为源手算。
24.3-2 构造负权边使 Dijkstra 出错的例子，并说明证明在哪一步失效。
24.3-3 把循环条件改成 $|Q|>1$（只循环 $|V|-1$ 次）是否正确（正确，最后一个顶点无需再松弛出边）。
24.3-4 $O(V+E)$ 时间验证给定的 d、π 是否对应某棵最短路径树。
24.3-5 构造反例说明 Dijkstra 并不总是按路径顺序松弛最短路径上的边。
24.3-6 可靠性 $r(u,v)\in[0,1]$ 独立，求最可靠路径（取 $-\log r$ 作权后跑 Dijkstra，或直接对乘积做“最大化”版本）。
24.3-7 权值为 $\{1..W\}$ 的整数时把每条边拆成单位边，证明 BFS 染黑顺序与 Dijkstra 出队顺序相同。
24.3-8 权值为 $\{0..W\}$ 的整数，$O(WV+E)$ 实现（桶队列 / Dial 算法）。
24.3-9 改进到 $O((V+E)\lg W)$（任一时刻 $V-S$ 中不同的估计值至多 W+1 个）。
24.3-10 只有源点出边可负时，Dijkstra 仍正确。

24.4 差分约束与最短路径（Difference constraints and shortest paths）（PDF p.685–691）

背景：线性规划（linear programming）。给定 $m\times n$ 矩阵 A、m 维向量 b、n 维向量 c，求 n 维向量 x，使目标函数 $\sum_{i=1}^n c_ix_i$ 最大，满足 $Ax\le b$。单纯形法（第 29 章）最坏情况下不是多项式时间，但存在多项式时间的线性规划算法。了解线性规划建模的两个理由：(1) 能把问题写成多项式规模的 LP，就立刻有多项式算法；(2) 许多 LP 特例有更快的专用算法，例如单对最短路径（习题 24.4-4）和最大流（习题 26.1-5）。有时不关心目标函数，只求可行解（feasible solution），即满足 $Ax\le b$ 的任意 x，或判定不存在——本节就是这样的可行性问题。

差分约束系统（system of difference constraints）：A 的每一行恰有一个 1 和一个 −1，其余为 0。于是 $Ax\le b$ 是 m 个形如

x_j-x_i\le b_k\quad(1\le i,j\le n,\ i\ne j,\ 1\le k\le m)

的约束。

例：5 个未知数、8 个约束：

x_1-x_2\le0,\ x_1-x_5\le-1,\ x_2-x_5\le1,\ x_3-x_1\le5,\ x_4-x_1\le4,\ x_4-x_3\le-1,\ x_5-x_3\le-3,\ x_5-x_4\le-3\quad(24.3\text{–}24.10)

一个解是 $x=(-5,-3,0,-1,-4)$；另一个是 $x'=(0,2,5,4,1)$，每个分量都大 5。这并非巧合：

引理 24.8：若 x 是差分约束系统 $Ax\le b$ 的解，d 为任意常数，则 $x+d=(x_1+d,\dots,x_n+d)$ 也是解。证明：$(x_j+d)-(x_i+d)=x_j-x_i$。

应用举例：$x_i$ 表示事件发生时间，约束表示两事件之间至少/至多间隔多久。例如胶水在 $x_1$ 时刻涂上，需要 2 小时凝固后才能在 $x_2$ 安装零件：$x_2\ge x_1+2$，即 $x_1-x_2\le-2$。若要求零件在涂胶之后、但不晚于胶水凝固一半时安装：$x_2\ge x_1$ 且 $x_2\le x_1+1$，即 $x_1-x_2\le0$、$x_2-x_1\le1$。

约束图（constraint graph）：把 $m\times n$ 矩阵 A 视为一个 n 顶点 m 条边的图的关联矩阵（incidence matrix）的转置。正式定义：$G=(V,E)$，

V=\{v_0,v_1,\dots,v_n\},\qquad E=\{(v_i,v_j): x_j-x_i\le b_k\text{ 是约束}\}\cup\{(v_0,v_1),\dots,(v_0,v_n)\}.

约束 $x_j-x_i\le b_k$ 对应边 $(v_i,v_j)$，权 $w(v_i,v_j)=b_k$；从附加源点 $v_0$ 出发的边权均为 0。$v_0$ 保证有一个顶点能到达所有其他顶点。图 24.8 是上例的约束图，各顶点的 $\delta(v_0,v_i)$ 恰为 $(-5,-3,0,-1,-4)$。

定理 24.9：设 G 是差分约束系统 $Ax\le b$ 的约束图。若 G 无负权环，则

x=(\delta(v_0,v_1),\delta(v_0,v_2),\dots,\delta(v_0,v_n))\qquad(24.11)

是可行解；若 G 有负权环，则系统无可行解。证明：

无负环：对任意边 $(v_i,v_j)$，三角不等式给出 $\delta(v_0,v_j)\le\delta(v_0,v_i)+w(v_i,v_j)$，即 $x_j-x_i\le w(v_i,v_j)=b_k$，满足对应约束。
有负环 $c=\langle v_1,\dots,v_k\rangle$，$v_1=v_k$（$v_0$ 无入边，不可能在环上）。它对应约束 $x_2-x_1\le w(v_1,v_2),\dots,x_k-x_{k-1}\le w(v_{k-1},v_k)$。若有解，把这 k−1 个不等式相加，左边每个未知数一加一减，和为 0（因 $x_1=x_k$），右边为 $w(c)$，得 $0\le w(c)<0$，矛盾。

求解方法：在约束图上以 $v_0$ 为源点跑 Bellman-Ford。因为 $v_0$ 到所有顶点都有边，任何负环都从 $v_0$ 可达。返回 TRUE 时最短路径权重即可行解（上例 $x=(-5,-3,0,-1,-4)$，由引理 24.8，$(d-5,d-3,d,d-1,d-4)$ 也是可行解）；返回 FALSE 时无可行解。

复杂度：m 个约束、n 个未知数，约束图有 n+1 个顶点、n+m 条边，Bellman-Ford 用时 $O((n+1)(n+m))=O(n^2+nm)$。习题 24.4-5 要求改成 $O(nm)$（即使 $m\ll n$）。

直观：差分约束 $x_j\le x_i+b_k$ 与松弛后的状态 $v.d\le u.d+w(u,v)$ 完全同构，最短路径距离是“最紧”的一组满足所有约束的值。由习题 24.4-8/24.4-9，Bellman-Ford 给出的解在 $x_i\le0$ 约束下使 $\sum x_i$ 最大，并使 $\max x_i-\min x_i$ 最小。

习题 24.4：

24.4-1、24.4-2 两个具体的差分约束系统，求可行解或判定无解。
24.4-3 约束图中从 $v_0$ 出发的最短路径权重能否为正（不能，至多为 0，因为有权 0 的直达边）。
24.4-4 把单对最短路径写成线性规划。
24.4-5 修改 Bellman-Ford 使差分约束求解为 $O(nm)$。
24.4-6 加入等式约束 $x_i=x_j+b_k$（拆成两个不等式）。
24.4-7 不用附加顶点 $v_0$ 的 Bellman-Ford 式解法（初始化所有 d 为 0）。
24.4-8★ 证明 Bellman-Ford 在 $Ax\le b,\ x_i\le0$ 下最大化 $\sum x_i$。
24.4-9★ 证明它最小化 $\max\{x_i\}-\min\{x_i\}$，并说明在施工排期中的用处（总工期最短）。
24.4-10 加入单变量约束 $x_i\le b_k$ 或 $-x_i\le b_k$（借助 $v_0$ 表示常数 0）。
24.4-11 b 为实数、要求 x 为整数（把 $b_k$ 向下取整）。
24.4-12★ 只要求部分变量为整数。

24.5 最短路径性质的证明（PDF p.692–698）

本节逐一证明章首的六条性质。

引理 24.10（三角不等式）：对任意边 $(u,v)\in E$，$\delta(s,v)\le\delta(s,u)+w(u,v)$。证明：s 到 v 的最短路径权重不超过“先走 s 到 u 的最短路径再走边 (u,v)”这条特定路径的权重。（无最短路径即 $\pm\infty$ 的情形留作习题 24.5-3。）它推广了 BFS 中的引理 22.1。

引理 24.11（上界性质）：初始化后，对所有 v 有 $v.d\ge\delta(s,v)$，任何松弛序列都保持该不变式；且一旦 $v.d=\delta(s,v)$ 就不再改变。证明（对松弛步数归纳）：基础——初始化后 $v.d=\infty$，$s.d=0\ge\delta(s,s)$（若 s 在负环上 $\delta(s,s)=-\infty$，否则为 0）。归纳——松弛 $(u,v)$ 只可能改变 $v.d$，若改变则 $v.d=u.d+w(u,v)\ge\delta(s,u)+w(u,v)\ge\delta(s,v)$。达到下界后不能再减（刚证的不变式），也不能增（松弛从不增大 d）。

推论 24.12（无路径性质）：若 s 到 v 无路径，则始终 $v.d=\delta(s,v)=\infty$。证明：$\infty=\delta(s,v)\le v.d$。

引理 24.13：执行 RELAX(u,v,w) 之后立即有 $v.d\le u.d+w(u,v)$。证明：若之前 $v.d>u.d+w$，则被赋值为等号；否则两者都不变，不等式本已成立。

引理 24.14（收敛性质）：若 $s\leadsto u\to v$ 是最短路径，且在调用 RELAX(u,v,w) 之前某时刻 $u.d=\delta(s,u)$，则调用后始终 $v.d=\delta(s,v)$。证明：由上界性质 $u.d=\delta(s,u)$ 一直保持；松弛后 $v.d\le u.d+w(u,v)=\delta(s,u)+w(u,v)=\delta(s,v)$（引理 24.13 与最短路径的最优子结构引理 24.1）；又 $v.d\ge\delta(s,v)$，故相等并保持。

引理 24.15（路径松弛性质）：对最短路径 $p=\langle v_0,\dots,v_k\rangle$，若松弛序列中按顺序包含 $(v_0,v_1),\dots,(v_{k-1},v_k)$，则此后 $v_k.d=\delta(s,v_k)$，无论中间穿插什么其他松弛。证明：对 i 归纳，第 i 条边松弛后 $v_i.d=\delta(s,v_i)$；基础 $v_0.d=s.d=0=\delta(s,s)$ 且不再变；归纳步用收敛性质。

引理 24.16：设 G 无从 s 可达的负权环。初始化后前驱子图 $G_\pi$ 是以 s 为根的有根树，且任何松弛序列都保持这一点。证明要点：

$G_\pi$ 无环：反设某次松弛产生环 $c=\langle v_0,\dots,v_k\rangle$，$v_k=v_0$，$v_i.\pi=v_{i-1}$，不妨设是松弛 $(v_{k-1},v_k)$ 造成的。环上每个顶点都有非 NIL 前驱，赋予前驱时 d 有限，由上界性质其最短路径权重有限，因而从 s 可达。松弛前，对 $i=1..k-1$，$v_i.d$ 最后一次更新为 $v_{i-1}.d+w(v_{i-1},v_i)$，此后 $v_{i-1}.d$ 只可能减小，于是 $v_i.d\ge v_{i-1}.d+w(v_{i-1},v_i)$（24.12）；而 $v_k.\pi$ 被改变说明 $v_k.d>v_{k-1}.d+w(v_{k-1},v_k)$（严格）。求和后两边 d 之和相同，得 $0>\sum w(v_{i-1},v_i)$，即 c 是可达负环，矛盾。
从 s 到 $V_\pi$ 每个顶点有路径（习题 24.5-6 归纳证明）。
路径唯一：若有两条简单路径 $s\leadsto u\leadsto x\to z\leadsto v$ 与 $s\leadsto u\leadsto y\to z\leadsto v$（$x\ne y$），则 $z.\pi=x$ 且 $z.\pi=y$，矛盾（图 24.9）。由附录 B.5-2，这三点说明 $G_\pi$ 是以 s 为根的有根树。

引理 24.17（前驱子图性质）：设 G 无从 s 可达的负环，初始化后经若干松弛使所有 $v.d=\delta(s,v)$，则 $G_\pi$ 是以 s 为根的最短路径树。证明三条定义：

$V_\pi$ 恰为可达顶点：$\delta(s,v)$ 有限 ⇔ v 可达；$v\ne s$ 的 d 有限 ⇔ $v.\pi\ne$ NIL。
有根树由引理 24.16 给出。
树上路径 $p=\langle v_0=s,\dots,v_k=v\rangle$ 是最短路径：对每个 i，$v_i.d=\delta(s,v_i)$ 且 $v_i.d\ge v_{i-1}.d+w(v_{i-1},v_i)$，故 $w(v_{i-1},v_i)\le\delta(s,v_i)-\delta(s,v_{i-1})$；求和后望远镜相消：$w(p)\le\delta(s,v_k)-\delta(s,v_0)=\delta(s,v_k)$。而 $\delta$ 是下界，所以 $w(p)=\delta(s,v)$。

习题 24.5：

24.5-1 为图 24.2 再给出两棵最短路径树。
24.5-2 构造图使每条边都既属于某棵最短路径树又不属于另一棵。
24.5-3 把引理 24.10 的证明补全到 $\pm\infty$ 情形。
24.5-4 证明：若松弛使 $s.\pi$ 变为非 NIL，则 G 含负环。
24.5-5 非负权图中，若允许 $v.\pi$ 取任意最短路径上的前驱，可能出现 $G_\pi$ 有环（0 权环），说明这种赋值不能由松弛产生。
24.5-6 证明 $V_\pi$ 中每个顶点从 s 可达且松弛保持该性质。
24.5-7 证明存在恰好 $|V|-1$ 次松弛的序列使所有 d 达到 δ（沿最短路径树的 BFS 序）。
24.5-8 有可达负环时，构造无限的松弛序列，每次松弛都改变某个估计值。

第 24 章思考题（Problems）（PDF p.699–703）

24-1 Yen 对 Bellman-Ford 的改进：给顶点任意线性序 $v_1,\dots,v_{|V|}$，把边分为 $E_f=\{(v_i,v_j):i<j\}$ 与 $E_b=\{(v_i,v_j):i>j\}$。(a) 证明 $G_f$、$G_b$ 都是 DAG，拓扑序分别为 $\langle v_1,\dots,v_{|V|}\rangle$ 和逆序。每轮先按 $v_1..v_{|V|}$ 顺序松弛 $E_f$ 出边，再按逆序松弛 $E_b$ 出边。(b) 证明无可达负环时只需 $\lceil|V|/2\rceil$ 轮。(c) 渐近时间是否改善（否，仍是 $O(VE)$，只改善常数）。
24-2 嵌套盒子（Nesting boxes）：d 维盒子 $(x_1..x_d)$ 嵌入 $(y_1..y_d)$ 当存在排列 π 使 $x_{\pi(i)}<y_i$ 对所有 i 成立。(a) 证明传递性；(b) 高效判断嵌套（各自排序后逐维比较）；(c) 求最长嵌套序列（建 DAG 求最长路径，$O(n^2d+nd\lg d)$）。
24-3 套利（Arbitrage）：利用汇率差异把一单位货币换回多于一单位。例：1 美元换 49 卢比，1 卢比换 2 日元，1 日元换 0.0107 美元，则 $49\times2\times0.0107=1.0486$，获利 4.86%。给定 n 种货币和汇率表 $R[i,j]$。(a) 判断是否存在货币序列 $\langle c_{i_1},\dots,c_{i_k}\rangle$ 使 $R[i_1,i_2]\cdot R[i_2,i_3]\cdots R[i_{k-1},i_k]\cdot R[i_k,i_1]>1$，分析时间；(b) 若存在则输出该序列。（标准解法：取边权 $w(i,j)=-\ln R[i,j]$，乘积 >1 ⇔ 权和 <0，即负权环；用 Bellman-Ford，$O(n^3)$；输出环用前驱回溯。）
24-4 Gabow 缩放算法（scaling algorithm）：非负整数边权，$W=\max w$，目标 $O(E\lg W)$。设 $k=\lceil\lg(W+1)\rceil$，$w_i(u,v)=\lfloor w(u,v)/2^{k-i}\rfloor$（取最高 i 位）。例：k=5，$w=25=\langle11001\rangle$，则 $w_3=\langle110\rangle=6$；$w=4=\langle00100\rangle$，$w_3=1$。依次计算 $\delta_1,\delta_2,\dots,\delta_k=\delta$。(a) 若所有 $\delta(s,v)\le|E|$，可 $O(E)$ 求出（桶式 Dijkstra）；(b) $O(E)$ 求 $\delta_1$；(c) 证明 $w_i=2w_{i-1}$ 或 $2w_{i-1}+1$，并且 $2\delta_{i-1}(s,v)\le\delta_i(s,v)\le2\delta_{i-1}(s,v)+|V|-1$；(d) 定义重新赋权 $\hat w_i(u,v)=w_i(u,v)+2\delta_{i-1}(s,u)-2\delta_{i-1}(s,v)$，证明它是非负整数；(e) 证明 $\delta_i(s,v)=\hat\delta_i(s,v)+2\delta_{i-1}(s,v)$ 且 $\hat\delta_i(s,v)\le|E|$；(f) 由此每步 $O(E)$，总计 $O(E\lg W)$。
24-5 Karp 最小平均权环算法：环 $c=\langle e_1..e_k\rangle$ 的平均权重 $\mu(c)=\frac1k\sum w(e_i)$，$\mu^*=\min_c\mu(c)$。设所有顶点从 s 可达，$\delta_k(s,v)$ 为恰好 k 条边的最短路径权重（不存在则 ∞）。(a) 若 $\mu^*=0$，则无负环且 $\delta(s,v)=\min_{0\le k\le n-1}\delta_k(s,v)$；(b) 若 $\mu^*=0$，则 $\max_{0\le k\le n-1}\frac{\delta_n(s,v)-\delta_k(s,v)}{n-k}\ge0$；(c) 0 权环上 u 到 v 沿环权 x，则 $\delta(s,v)=\delta(s,u)+x$；(d) $\mu^*=0$ 时每个最小平均环上存在顶点使上式等于 0；(e) $\mu^*=0$ 时 $\min_v\max_k\frac{\delta_n-\delta_k}{n-k}=0$；(f) 每条边加常数 t 则 $\mu^*$ 增加 t，从而一般地 $\mu^*=\min_{v\in V}\max_{0\le k\le n-1}\frac{\delta_n(s,v)-\delta_k(s,v)}{n-k};$ (g) 给出 $O(VE)$ 算法（DP 计算所有 $\delta_k$）。
24-6 双调最短路径（Bitonic shortest paths）：序列先单调增后单调减（或循环移位后如此）称为双调，如 $\langle1,4,6,8,3,-2\rangle$、$\langle9,2,-4,-10,-5\rangle$、$\langle1,2,3,4\rangle$ 是，$\langle1,3,12,4,2,10\rangle$ 不是。边权互不相同，且已知每条最短路径的边权序列是双调的，求尽可能高效的单源最短路径算法（按权排序后做一次递增、一次递减的松弛，$O(E\lg E)$）。

第 24 章注记（Chapter notes）（PDF p.703–704）

Dijkstra 算法发表于 1959 年，原文未提优先队列。Bellman-Ford 源自 Bellman 和 Ford 各自的工作；Bellman 论述了最短路径与差分约束的关系。DAG 线性时间算法由 Lawler 记述，视为民间流传。小的非负整数权时，Dijkstra 的 EXTRACT-MIN 返回值单调递增，可用更快的数据结构：Ahuja–Mehlhorn–Orlin–Tarjan 的 $O(E+V\sqrt{\lg W})$；Thorup 的 $O(E\lg\lg V)$；Raman 的 $O(E+V\min\{(\lg V)^{1/3+\epsilon},(\lg W)^{1/4+\epsilon}\})$（空间依赖机器字长，可用随机哈希降为线性）。无向整数权图上 Thorup 有 $O(V+E)$ 算法（不是 Dijkstra 的实现）。负权图：Gabow–Tarjan $O(\sqrt V E\lg(VW))$，Goldberg $O(\sqrt V E\lg W)$。Cherkassky、Goldberg、Radzik 做了大规模实验比较。

第 24 章小结

本章要点 单源最短路径的所有算法都建立在“初始化 + 松弛”的框架上，差别只在松弛的顺序与次数。最短路径有最优子结构，且可取为简单路径（至多 |V|−1 条边）。Bellman-Ford 适用于任意实数权，$O(VE)$，并能检测可达负环；DAG 上按拓扑序松弛一次即可，$\Theta(V+E)$，还能求最长路径（关键路径）；Dijkstra 要求非负权，用二叉堆 $O((V+E)\lg V)$、斐波那契堆 $O(V\lg V+E)$。差分约束系统 $x_j-x_i\le b_k$ 等价于约束图上的最短路径问题：无负环时 $x_i=\delta(v_0,v_i)$ 是可行解，有负环则不可行。正确性证明依靠三角不等式、上界、无路径、收敛、路径松弛、前驱子图六条性质。

与量化交易的关联

外汇/加密货币三角套利检测：思考题 24-3 本身就是量化场景。把汇率取 $-\ln R$ 作为边权，套利机会即负权环，Bellman-Ford 或 SPFA 可在 $O(n^3)$ 内检测；实际系统中还要把买卖价差、手续费并入 R，并用习题 24.1-6 的方法回溯出具体货币环路。交易所多、币种多时，常用增量版本（只对变动的边重新松弛）。
最优兑换/最优执行路径：在多交易所、多交易对之间把资产 A 换成资产 B 的“最便宜路径”，就是以 $-\ln(\text{成交率})$ 为权的单源最短路径；全部权非负时用 Dijkstra。
事件调度与系统时序约束：差分约束可用来排定交易系统任务（数据到达、因子计算、下单截止时间）之间的“至少/至多间隔”约束，并检查约束是否自相矛盾（负环）。DAG 关键路径可用于估计每日批处理流水线（行情落地 → 清洗 → 因子 → 优化 → 下单）的最短完成时间并找出瓶颈。
动态规划视角：Bellman-Ford 的“按边数递推”与最优停止、多期决策中的值迭代同构，便于理解后面的定价与执行算法。
平均权最小环（24-5）可用于在汇率图中找“单位步收益最大”的循环套利，但实际中受流动性约束，参考意义有限。

推荐习题

24.1-3（提前终止的 Bellman-Ford，工程中常用）；24.1-6（输出负环，套利实现必需）。
24.2-4（DAG 路径计数，DP 练习）。
24.3-2、24.3-10（理解 Dijkstra 对非负权的依赖）；24.3-6（概率乘积转对数和）；24.3-8（桶队列）。
24.4-1/24.4-2（差分约束建模手算）；24.4-9（最小化跨度的调度含义）。
思考题 24-3（套利，强烈推荐）；24-5（Karp 最小平均环，DP 与图论结合）。

第 25 章所有结点对的最短路径（All-Pairs Shortest Paths）

第 25 章引言（PDF p.705–707）

问题：对带权有向图 $G=(V,E)$、$w:E\to\mathbb R$，求每一对顶点 $u,v$ 间的最短路径，通常以表格形式输出（u 行 v 列为 $\delta(u,v)$），例如公路地图册的城市距离表。

朴素做法——跑 |V| 次单源算法：

非负权用 Dijkstra：数组优先队列 $O(V^3+VE)=O(V^3)$；二叉堆 $O(VE\lg V)$（稀疏图更好）；斐波那契堆 $O(V^2\lg V+VE)$。
有负权边只能用 Bellman-Ford：$O(V^2E)$，稠密图为 $O(V^4)$。本章要做得更好，并考察与矩阵乘法的关系及其代数结构。

输入表示：与单源算法用邻接表不同，本章多数算法用邻接矩阵（Johnson 算法除外）。顶点编号 $1..n$，输入 $n\times n$ 矩阵 $W=(w_{ij})$：

w_{ij}=\begin{cases}0 & i=j\\ \text{边 }(i,j)\text{ 的权} & i\ne j,(i,j)\in E\\ \infty & i\ne j,(i,j)\notin E\end{cases}\qquad(25.1)

允许负权边，暂假定无负权环。

输出：$n\times n$ 矩阵 $D=(d_{ij})$，结束时 $d_{ij}=\delta(i,j)$；以及前驱矩阵（predecessor matrix） $\Pi=(\pi_{ij})$：$i=j$ 或无路径时为 NIL，否则为 i 到 j 某条最短路径上 j 的前驱。第 i 行诱导的子图应是以 i 为根的最短路径树：$V_{\pi,i}=\{j:\pi_{ij}\ne\text{NIL}\}\cup\{i\}$，$E_{\pi,i}=\{(\pi_{ij},j):j\in V_{\pi,i}-\{i\}\}$。

def print_all_pairs_shortest_path(Pi, i, j):
    if i == j:
        print(i)
    elif Pi[i][j] is None:
        print(f"no path from {i} to {j} exists")
    else:
        print_all_pairs_shortest_path(Pi, i, Pi[i][j])
        print(j)

章节安排：25.1 基于矩阵乘法的动态规划，用“重复平方”达到 $\Theta(V^3\lg V)$；25.2 Floyd-Warshall 算法 $\Theta(V^3)$，以及传递闭包；25.3 Johnson 算法 $O(V^2\lg V+VE)$，适合大规模稀疏图。记号：$n=|V|$；矩阵用大写字母，元素用带下标的小写字母；迭代序号用上标括号，如 $L^{(m)}=(l_{ij}^{(m)})$；$A.rows$ 存放矩阵阶数 n。

25.1 最短路径与矩阵乘法（PDF p.707–714）

按第 15 章动态规划步骤：刻画最优解结构 → 递归定义最优值 → 自底向上计算（构造最优解留作习题）。

最优解结构：由引理 24.1，最短路径的子路径也是最短路径。若 i 到 j 的最短路径 p 至多含 m 条边（无负环时 m 有限）：$i=j$ 时权 0；否则分解为 $i\overset{p'}{\leadsto}k\to j$，$p'$ 至多 m−1 条边且是 i 到 k 的最短路径，故 $\delta(i,j)=\delta(i,k)+w_{kj}$。

递归式：令 $l_{ij}^{(m)}$ 为 i 到 j 至多含 m 条边的路径的最小权重。

l_{ij}^{(0)}=\begin{cases}0&i=j\\\infty&i\ne j\end{cases},\qquad l_{ij}^{(m)}=\min\Big(l_{ij}^{(m-1)},\min_{1\le k\le n}\{l_{ik}^{(m-1)}+w_{kj}\}\Big)=\min_{1\le k\le n}\{l_{ik}^{(m-1)}+w_{kj}\}\quad(25.2)

第二个等号因为 $w_{jj}=0$。无负环时，最短路径是简单路径，至多 n−1 条边，因此

\delta(i,j)=l_{ij}^{(n-1)}=l_{ij}^{(n)}=l_{ij}^{(n+1)}=\cdots\qquad(25.3)

自底向上：$L^{(1)}=W$，依次算 $L^{(2)},\dots,L^{(n-1)}$。核心过程把已求得的最短路径再延伸一条边：

def extend_shortest_paths(L, W):            # Θ(n^3)
    n = len(L)
    Lp = [[INF] * n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            for k in range(n):
                Lp[i][j] = min(Lp[i][j], L[i][k] + W[k][j])
    return Lp

与矩阵乘法的对应：$C=A\cdot B$，$c_{ij}=\sum_k a_{ik}\cdot b_{kj}$（25.4）。在 (25.2) 中做替换 $l^{(m-1)}\to a$，$w\to b$，$l^{(m)}\to c$，$\min\to+$，$+\to\cdot$，并把 min 的单位元 $\infty$ 换成加法单位元 0，就得到 4.2 节的 SQUARE-MATRIX-MULTIPLY。即这是 (min, +) 半环 上的矩阵乘法（热带代数）。用 $A\cdot B$ 表示 EXTEND-SHORTEST-PATHS(A,B) 的结果：

L^{(1)}=L^{(0)}\cdot W=W,\ L^{(2)}=W^2,\ \dots,\ L^{(n-1)}=W^{n-1}.

def slow_all_pairs_shortest_paths(W):       # Θ(n^4)
    n = len(W)
    L = W
    for m in range(2, n):                   # m = 2 .. n-1
        L = extend_shortest_paths(L, W)
    return L

例（图 25.1）：5 个顶点，边 $1\to2=3,1\to3=8,1\to5=-4,2\to4=1,2\to5=7,3\to2=4,4\to1=2,4\to3=-5,5\to4=6$。最终

L^{(4)}=\begin{pmatrix}0&1&-3&2&-4\\3&0&-4&1&-1\\7&4&0&5&3\\2&-1&-5&0&-2\\8&5&1&6&0\end{pmatrix},

且 $L^{(5)}=L^{(4)}\cdot W=L^{(4)}$，此后不变。

重复平方（repeated squaring）：只需 $L^{(n-1)}$；无负环时对所有 $m\ge n-1$ 有 $L^{(m)}=L^{(n-1)}$。(min,+) 矩阵乘法满足结合律（习题 25.1-4），故计算 $W,W^2,W^4,W^8,\dots,W^{2^{\lceil\lg(n-1)\rceil}}$，共 $\lceil\lg(n-1)\rceil$ 次乘积即可，因 $2^{\lceil\lg(n-1)\rceil}\ge n-1$。

def faster_all_pairs_shortest_paths(W):     # Θ(n^3 lg n)
    n = len(W)
    L, m = W, 1
    while m < n - 1:
        L = extend_shortest_paths(L, L)     # L^(2m) = (L^(m))^2
        m *= 2
    return L

最后一次迭代算出某个 $L^{(2m)}$，$n-1\le2m<2n-2$，由 (25.3) 它等于 $L^{(n-1)}$。时间 $\Theta(n^3\lg n)$，代码紧凑，常数小；按原写法保存 $\lceil\lg(n-1)\rceil$ 个矩阵空间为 $\Theta(n^2\lg n)$，只用两个矩阵可降到 $\Theta(n^2)$（习题 25.1-8）。

习题 25.1：

25.1-1 在图 25.2（6 顶点）上手算 SLOW 与 FASTER 版本。
25.1-2 为什么要求 $w_{ii}=0$（使 $l^{(m)}$ 表示“至多 m 条边”，并让 (25.2) 的两种写法等价）。
25.1-3 $L^{(0)}$（对角为 0、其余 ∞）对应普通矩阵乘法中的单位矩阵。
25.1-4 证明 (min,+) 乘法满足结合律。
25.1-5 把单源最短路径写成矩阵与向量的乘积，对应类似 Bellman-Ford 的算法。
25.1-6 由最终 L 在 $O(n^3)$ 内求前驱矩阵 Π。
25.1-7 在计算 $L^{(m)}$ 的同时计算 $\Pi^{(m)}$。
25.1-8 只用两个 $n\times n$ 矩阵，空间 $\Theta(n^2)$。
25.1-9 修改 FASTER 版本以检测负环（再乘一次看对角线是否出现负值或是否变化）。
25.1-10 求边数最少的负权环的边数。

25.2 Floyd-Warshall 算法（PDF p.714–721）

适用：可有负权边，无负权环；时间 $\Theta(V^3)$。

最优解结构——按中间顶点刻画：简单路径 $p=\langle v_1,\dots,v_l\rangle$ 的**中间顶点（intermediate vertex）**是除 $v_1,v_l$ 外的顶点。考虑所有中间顶点都取自 $\{1,\dots,k\}$ 的 i→j 路径，设 p 为其中权重最小者（简单路径）：

若 k 不是 p 的中间顶点，则 p 的中间顶点都在 $\{1..k-1\}$ 中，p 也是“中间顶点限于 $\{1..k-1\}$”的最短路径。
若 k 是中间顶点，分解 $i\overset{p_1}{\leadsto}k\overset{p_2}{\leadsto}j$（图 25.3）。$p_1$、$p_2$ 的中间顶点都在 $\{1..k-1\}$ 中（k 只出现一次），且分别是相应限制下 i→k、k→j 的最短路径。

递归式：$d_{ij}^{(k)}$ 为中间顶点都在 $\{1..k\}$ 中的 i→j 最短路径权重。$k=0$ 时路径无中间顶点，至多一条边：

d_{ij}^{(k)}=\begin{cases}w_{ij}&k=0\\\min\big(d_{ij}^{(k-1)},\ d_{ik}^{(k-1)}+d_{kj}^{(k-1)}\big)&k\ge1\end{cases}\qquad(25.5)

$D^{(n)}$ 即答案：$d_{ij}^{(n)}=\delta(i,j)$。

def floyd_warshall(W):
    n = len(W)
    D = [row[:] for row in W]               # D^(0) = W
    for k in range(n):
        for i in range(n):
            for j in range(n):
                if D[i][k] + D[k][j] < D[i][j]:
                    D[i][j] = D[i][k] + D[k][j]
    return D

原书为每个 k 新建 $D^{(k)}$，空间 $\Theta(n^3)$；习题 25.2-4 证明去掉上标、原地更新（如上）仍正确，空间 $\Theta(n^2)$（原因：第 k 轮中 $d_{ik}$、$d_{kj}$ 不会被改变，因为 $d_{kk}=0$）。时间 $\Theta(n^3)$，三重循环、无复杂数据结构，常数很小，中等规模图上很实用。图 25.4 给出图 25.1 上的 $D^{(0)}..D^{(5)}$ 与 $\Pi^{(0)}..\Pi^{(5)}$，$D^{(5)}$ 即上面的 $L^{(4)}$。

构造最短路径：方法一，先算 D 再在 $O(n^3)$ 内推出 Π（习题 25.1-6）。方法二，同时计算 $\Pi^{(0)},\dots,\Pi^{(n)}$，$\pi_{ij}^{(k)}$ 是中间顶点限于 $\{1..k\}$ 的最短路径上 j 的前驱：

\pi_{ij}^{(0)}=\begin{cases}\text{NIL}&i=j\text{ 或 }w_{ij}=\infty\\ i&i\ne j\text{ 且 }w_{ij}<\infty\end{cases}\qquad(25.6)

\pi_{ij}^{(k)}=\begin{cases}\pi_{ij}^{(k-1)}&d_{ij}^{(k-1)}\le d_{ik}^{(k-1)}+d_{kj}^{(k-1)}\\ \pi_{kj}^{(k-1)}&d_{ij}^{(k-1)}>d_{ik}^{(k-1)}+d_{kj}^{(k-1)}\end{cases}\qquad(25.7)

（走经过 k 的路径时，j 的前驱沿用 k→j 那段路径中 j 的前驱。）把它并入算法及证明 $G_{\pi,i}$ 是最短路径树留作习题 25.2-3；方法三见习题 25.2-7（记录最大编号中间顶点 $\phi_{ij}^{(k)}$，类似矩阵链乘法的 s 表）。

有向图的传递闭包（transitive closure）：$G^*=(V,E^*)$，$E^*=\{(i,j):G\text{ 中有 } i\text{ 到 } j\text{ 的路径}\}$。

方法一：所有边权设为 1 跑 Floyd-Warshall，有路径则 $d_{ij}<n$，否则 $\infty$。$\Theta(n^3)$。
方法二：把 min、+ 换成逻辑或 ∨、逻辑与 ∧。$t_{ij}^{(k)}=1$ 当且仅当存在中间顶点都在 $\{1..k\}$ 中的 i→j 路径： $t_{ij}^{(0)}=\begin{cases}0&i\ne j\text{ 且 }(i,j)\notin E\\1&i=j\text{ 或 }(i,j)\in E\end{cases},\qquad t_{ij}^{(k)}=t_{ij}^{(k-1)}\vee\big(t_{ik}^{(k-1)}\wedge t_{kj}^{(k-1)}\big)\qquad(25.8)$

def transitive_closure(n, edges):           # Θ(n^3) 时间
    T = [[1 if i == j or (i, j) in edges else 0 for j in range(n)] for i in range(n)]
    for k in range(n):
        for i in range(n):
            for j in range(n):
                T[i][j] = T[i][j] or (T[i][k] and T[k][j])
    return T

同为 $\Theta(n^3)$，但单比特逻辑运算在某些机器上更快，存储量也比整数版少一个字长的倍数（可用位并行，每次处理一整行）。图 25.5 的 4 顶点例子（边 2→3,2→4,3→2,4→1,4→3）最终 $T^{(4)}$：第 1 行 (1,0,0,0)，其余三行全为 1。

习题 25.2：

25.2-1 在图 25.2 上运行并给出每轮 $D^{(k)}$。
25.2-2 用 25.1 节的方法求传递闭包。
25.2-3 并入 $\Pi^{(k)}$ 计算并严格证明 $G_{\pi,i}$ 是最短路径树。
25.2-4 原地更新版本正确，空间 $\Theta(n^2)$。
25.2-5 把 (25.7) 中的相等情形改为取 $\pi_{kj}$ 是否仍正确。
25.2-6 用输出检测负环（对角线出现负值 $d_{ii}<0$）。
25.2-7 用最大编号中间顶点 $\phi$ 重建路径。
25.2-8 $O(VE)$ 求传递闭包（从每个顶点 BFS/DFS）。
25.2-9 若 DAG 的传递闭包可在 $f(|V|,|E|)$ 内求出，则一般有向图可在 $f+O(V+E)$ 内求出（强连通分量缩点）。

25.3 稀疏图上的 Johnson 算法（PDF p.721–726）

概述：$O(V^2\lg V+VE)$，稀疏图上渐近快于重复平方和 Floyd-Warshall；返回最短路径矩阵，或报告存在负环。以 Bellman-Ford 和 Dijkstra 为子程序。核心技术是重新赋权（reweighting）：若所有边权非负，对每个顶点跑一次 Dijkstra（斐波那契堆）即 $O(V^2\lg V+VE)$；若有负权边但无负环，就先构造新的非负权 $\hat w$，要求：

对所有 u,v，p 是 w 下的最短路径 ⇔ p 是 $\hat w$ 下的最短路径；
对所有边，$\hat w(u,v)\ge0$。预处理求 $\hat w$ 只需 $O(VE)$。

引理 25.1（重新赋权不改变最短路径）：对任意 $h:V\to\mathbb R$，定义

\hat w(u,v)=w(u,v)+h(u)-h(v).\qquad(25.9)

则对任意路径 $p=\langle v_0,\dots,v_k\rangle$：$w(p)=\delta(v_0,v_k)\iff\hat w(p)=\hat\delta(v_0,v_k)$；且 G 在 w 下有负环 ⇔ 在 $\hat w$ 下有负环。证明：望远镜求和

\hat w(p)=\sum_{i=1}^k\big(w(v_{i-1},v_i)+h(v_{i-1})-h(v_i)\big)=w(p)+h(v_0)-h(v_k).\qquad(25.10)

$h(v_0)-h(v_k)$ 与路径无关，所以路径之间的大小关系不变。对环 $v_0=v_k$，$\hat w(c)=w(c)$，负环性质不变。

构造非负权：新建图 $G'=(V',E')$，$V'=V\cup\{s\}$，$E'=E\cup\{(s,v):v\in V\}$，$w(s,v)=0$。s 无入边，所以除以 s 为源的路径外，任何最短路径都不经过 s；$G'$ 无负环 ⇔ G 无负环。令 $h(v)=\delta(s,v)$，由三角不等式 $h(v)\le h(u)+w(u,v)$，于是 $\hat w(u,v)=w(u,v)+h(u)-h(v)\ge0$。

例（图 25.6，用图 25.1 的图）：$h(1..5)=(0,-1,-5,0,-4)$。重新赋权后：$\hat w(1,2)=4,\hat w(1,3)=13,\hat w(1,5)=0,\hat w(2,4)=0,\hat w(2,5)=10,\hat w(3,2)=0,\hat w(4,1)=2,\hat w(4,3)=0,\hat w(5,4)=2$，全部非负。再以每个顶点为源跑 Dijkstra，图中每个顶点标注 $\hat\delta(u,v)/\delta(u,v)$，并有 $d_{uv}=\delta(u,v)=\hat\delta(u,v)+h(v)-h(u)$。

def johnson(G, w):
    # 1. 加超级源 s，到每个顶点一条 0 权边
    Gp = G.copy(); s = Gp.add_vertex()
    for v in G.V:
        Gp.add_edge(s, v, 0)
    # 2. Bellman-Ford 求 h，并检测负环           O(VE)
    if not bellman_ford(Gp, w, s):
        raise ValueError("the input graph contains a negative-weight cycle")
    h = {v: v.d for v in Gp.V}
    # 3. 重新赋权                                O(E)
    w_hat = {(u, v): w(u, v) + h[u] - h[v] for (u, v) in Gp.E}
    # 4. 每个顶点跑一次 Dijkstra                 |V| × O(V lg V + E)
    D = {}
    for u in G.V:
        delta_hat = dijkstra(G, w_hat, u)
        for v in G.V:
            D[u, v] = delta_hat[v] + h[v] - h[u]   # 还原原权重
    return D

复杂度：斐波那契堆 $O(V^2\lg V+VE)$；二叉堆 $O(VE\lg V)$，稀疏图上仍渐近快于 Floyd-Warshall。空间 $O(V^2)$（输出矩阵）+ $O(V+E)$。

常见误区：

不能简单地把所有边权都加上 $-w^*$（最小边权）来消除负权：不同边数的路径加的量不同，会改变最短路径（习题 25.3-4）。势函数 h 之所以可行，是因为它的修正量只取决于起点和终点。
超级源点 s 必不可少：若随便取已有顶点为 s，可能有顶点不可达，h 为 ∞（习题 25.3-6）；若图强连通则可以。

习题 25.3：

25.3-1 在图 25.2 上运行 Johnson，给出 h 与 $\hat w$。
25.3-2 新增顶点 s 的作用。
25.3-3 若所有 $w\ge0$，w 与 $\hat w$ 的关系（h≡0，二者相同）。
25.3-4 指出“减去最小边权”方法的错误。
25.3-5 证明 0 权环上每条边重新赋权后为 0。
25.3-6 不加新源点时的反例，以及强连通时正确。

第 25 章思考题与注记（PDF p.726–728）

25-1 动态图的传递闭包：边逐条插入，用布尔矩阵维护传递闭包。(a) 每插入一条边 $O(V^2)$ 更新；(b) 举例说明某次插入必须 $\Omega(V^2)$；(c) 设计算法，使任意 n 次插入总时间 $\sum t_i=O(V^3)$（只对新变为可达的对做传播，每对至多由 0 变 1 一次）。
25-2 ε-稠密图上的最短路径：$|E|=\Theta(V^{1+\epsilon})$，$0<\epsilon\le1$。用 d 叉最小堆（思考题 6-2）可在不使用斐波那契堆的情况下达到同样时间。(a) d 叉堆 INSERT、EXTRACT-MIN、DECREASE-KEY 的时间关于 d、n 的表达式，取 $d=\Theta(n^\alpha)$ 时如何，与斐波那契堆摊还代价比较（INSERT/DECREASE-KEY $O(\log_d n)$，EXTRACT-MIN $O(d\log_d n)$）；(b) 非负权单源最短路径 $O(E)$（取 $d=V^\epsilon$）；(c) 非负权全源 $O(VE)$；(d) 有负权无负环时全源 $O(VE)$（结合 Johnson）。

注记：Lawler 有关于全源最短路径的良好讨论，矩阵乘法算法属民间流传。Floyd-Warshall 来自 Floyd，基于 Warshall 关于布尔矩阵传递闭包的定理。Johnson 算法见 [192]。改进：Fredman 用 $O(V^{5/2})$ 次比较，得到 $O(V^3(\lg\lg V/\lg V)^{1/3})$；Han 改进到 $O(V^3(\lg\lg V/\lg V)^{5/4})$。利用快速矩阵乘法（$O(n^\omega)$，$\omega<2.376$），Galil–Margalit、Seidel 在无向无权图上达到 $O(V^\omega p(V))$（p 为多对数）；Shoshan–Zwick 对整数权 $\{1..W\}$ 无向图 $O(WV^\omega p(VW))$。Karger–Koller–Phillips 与 McGeoch 给出依赖于“参与某条最短路径的边集 $E^*$”的界 $O(VE^*+V^2\lg V)$。Baswana–Hariharan–Sen 研究递减（删除边）算法；Demetrescu–Italiano 处理插入删除混合。Aho–Hopcroft–Ullman 定义的**闭半环（closed semiring）**给出有向图路径问题的统一代数框架，Floyd-Warshall 与传递闭包算法都是其实例。

第 25 章小结

本章要点 全源最短路径有三条路线：(1) 把“再走一条边”的递推看作 (min,+) 半环上的矩阵乘法，朴素 $\Theta(n^4)$，重复平方 $\Theta(n^3\lg n)$；(2) Floyd-Warshall 按“允许使用的中间顶点集合 $\{1..k\}$”做动态规划，$\Theta(n^3)$ 时间、$\Theta(n^2)$ 空间，代码最简，换成布尔运算即得传递闭包；(3) Johnson 算法用势函数 $h(v)=\delta(s,v)$ 重新赋权，把负权变非负，再跑 |V| 次 Dijkstra，$O(V^2\lg V+VE)$，适合稀疏图。重新赋权的核心等式 $\hat w(p)=w(p)+h(v_0)-h(v_k)$ 说明它不改变最短路径，也不改变环的权重。负环检测：Floyd-Warshall 看对角线是否为负，Johnson 由 Bellman-Ford 报告。

与量化交易的关联

多币种/多资产全对套利与最优兑换表：对 n 种货币，以 $-\ln R[i,j]$ 为权跑 Floyd-Warshall，$O(n^3)$ 一次得到任意两种货币间的最优兑换率；对角线 $d_{ii}<0$ 直接标出参与套利环的货币。几十到几百种币种时 $n^3$ 完全可接受。
(min,+) 代数视角：最优执行、多期库存控制等问题的 DP 递推常可写成 (min,+) 或 (max,+) 矩阵乘积，可用重复平方加速长期限的价值递推，也便于向量化（numpy 广播实现 EXTEND）。
传递闭包：可用于系统实现中的依赖分析（因子之间、任务之间谁依赖谁），或判断资产之间是否存在可兑换链路。
势函数重新赋权的思想在优化中常见（对偶变量、影子价格），与第 29 章对偶性呼应；在最小费用流（订单路由、资金调拨）中也用同样技巧让 Dijkstra 能处理负成本边。
对相关性网络、资产聚类等研究，本章算法本身用得不多；那类问题更常用最小生成树（第 23 章，在上一块）。

推荐习题

25.1-4（证明 (min,+) 结合律，理解半环）；25.1-8（空间优化）；25.1-9（负环检测）。
25.2-4（原地 Floyd-Warshall 正确性，实现必会）；25.2-6（负环检测）；25.2-3（前驱矩阵）。
25.3-4、25.3-6（理解重新赋权为何必须用势函数和超级源点）。
思考题 25-2（d 叉堆，工程中替代斐波那契堆的实用选择）。

第 26 章最大流（Maximum Flow）

第 26 章引言（PDF p.729–730）

把有向图看作流网络（flow network），回答物料流动问题：物料从源点（source）以稳定速率产生，流经系统到汇点（sink）被消耗。每条有向边是一条管道，有给定容量（capacity）（如每小时 200 加仑、导线 20 安培）；顶点是管道交汇点，除源汇外物料不在顶点积存——进入速率等于离开速率，称为流量守恒（flow conservation），相当于电路中的基尔霍夫电流定律。最大流问题：在不违反容量约束的前提下，求从源到汇的最大输送速率。可建模液体管道、装配线零件、电网电流、通信网络信息流等。

章节安排：26.1 形式化流网络与最大流问题；26.2 Ford-Fulkerson 方法；26.3 应用于二分图最大匹配；26.4 推送-重贴标签（push-relabel）方法，是许多最快网络流算法的基础；26.5 “前置重贴标签（relabel-to-front）”算法，$O(V^3)$，虽非最快但展示了最快算法用到的技巧，实践中也相当高效。

26.1 流网络（Flow networks）（PDF p.730–735）

定义：流网络 $G=(V,E)$ 是有向图，每条边 $(u,v)$ 有非负容量 $c(u,v)\ge0$。要求：若 $(u,v)\in E$ 则反向边 $(v,u)\notin E$（不允许反平行边，后面有绕过办法）；$(u,v)\notin E$ 时约定 $c(u,v)=0$；不允许自环。指定源点 s、汇点 t。为方便假定每个顶点都在某条 $s\leadsto v\leadsto t$ 路径上，因此图连通且 $|E|\ge|V|-1$。

**流（flow）**是函数 $f:V\times V\to\mathbb R$，满足：

容量约束（capacity constraint）：对所有 u,v，$0\le f(u,v)\le c(u,v)$。
流量守恒：对所有 $u\in V-\{s,t\}$，$\sum_{v\in V}f(v,u)=\sum_{v\in V}f(u,v)$（流入 = 流出）。 $(u,v)\notin E$ 时 $f(u,v)=0$。

流的值（value）：

|f|=\sum_{v\in V}f(s,v)-\sum_{v\in V}f(v,s)\qquad(26.1)

即流出源点的总量减去流入源点的总量（这里 $|\cdot|$ 表示流值，不是绝对值或基数）。通常没有进入源点的边，第二项为 0；保留它是因为在残存网络中流入源点的流有意义。最大流问题：给定 G、s、t，求值最大的流。

例（图 26.1，Lucky Puck 公司运输问题）：温哥华工厂（s）生产冰球，温尼伯仓库（t）存货，经埃德蒙顿（$v_1$）、卡尔加里（$v_2$）、萨斯卡通（$v_3$）、里贾纳（$v_4$）中转，每对城市间每天至多运 $c(u,v)$ 箱。容量：$s\to v_1=16,s\to v_2=13,v_1\to v_3=12,v_2\to v_1=4,v_2\to v_4=14,v_3\to v_2=9,v_3\to t=20,v_4\to v_3=7,v_4\to t=4$。图 (b) 给出一个值为 19 的流（标记 $f/c$：11/16, 8/13, 12/12, 1/4, 11/14, 4/9, 15/20, 7/7, 4/4）。公司关心的是每天能稳定发出并到达多少箱 p，不关心单个冰球花多久；稳态下中转城市不能积压，正对应流量守恒。

反平行边建模（antiparallel edges）：如果又能租用埃德蒙顿→卡尔加里每天 10 箱的运力，就同时有 $(v_1,v_2)$ 与 $(v_2,v_1)$，违反假设。处理：选其中一条如 $(v_1,v_2)$，引入新顶点 $v'$，替换为 $(v_1,v')$ 与 $(v',v_2)$，容量都等于原容量 10（图 26.2）。习题 26.1-1 证明等价。

多源多汇（multiple sources and sinks）：m 个工厂 $\{s_1..s_m\}$、n 个仓库 $\{t_1..t_n\}$。加**超级源点（supersource）s，连边 $(s,s_i)$ 容量 ∞；加超级汇点（supersink）**t，连边 $(t_i,t)$ 容量 ∞（图 26.3）。两个问题的流一一对应（习题 26.1-2）。

习题 26.1：

26.1-1 证明拆边得到等价网络。
26.1-2 把流的定义推广到多源多汇，证明与超级源汇网络等价。
26.1-3 若某顶点 u 不在任何 $s\leadsto u\leadsto t$ 路径上，则存在最大流使 u 相关的流全为 0。
26.1-4 证明流构成凸集：$\alpha f_1+(1-\alpha)f_2$ 仍是流。
26.1-5 把最大流写成线性规划。
26.1-6 Adam 教授两个孩子不愿走同一条街区（路口交叉可以），判断能否去同一所学校——建模为最大流（无向边拆成两条有向边，容量 1，看最大流是否 ≥2）。
26.1-7 顶点容量 $l(v)$：把每个顶点拆成入点和出点，中间连容量 $l(v)$ 的边；新网络 $2|V|$ 个顶点、$|E|+|V|$ 条边。

26.2 Ford-Fulkerson 方法（PDF p.735–752）

称为“方法”而非“算法”，因为它包含多种运行时间不同的实现。三个关键思想：残存网络（residual network）、增广路径（augmenting path）、割（cut），它们是最大流最小割定理（定理 26.6）的基础。

基本框架：从 $f\equiv0$ 开始，每次在残存网络 $G_f$ 中找一条增广路径，沿它增加流值；直到没有增广路径。过程中单条边上的流可能增加也可能减少——减少某些边上的流可能是为了让更多流到达汇点。

FORD-FULKERSON-METHOD(G, s, t):
    f ← 0
    while 残存网络 G_f 中存在增广路径 p:
        沿 p 增广 f
    return f

残存容量（residual capacity）：

c_f(u,v)=\begin{cases}c(u,v)-f(u,v)&(u,v)\in E\\ f(v,u)&(v,u)\in E\\0&\text{其他}\end{cases}\qquad(26.2)

由于不允许反平行边，每个有序对恰好适用一种情形。例：$c(u,v)=16,f(u,v)=11$，则还能再加 $c_f(u,v)=5$，也能把最多 11 单位“退回”，$c_f(v,u)=11$。反向的残存边表示减少原边上的流。

残存网络：$G_f=(V,E_f)$，$E_f=\{(u,v)\in V\times V:c_f(u,v)>0\}$（26.3）。$E_f$ 中的边要么是原边，要么是原边的反向，故 $|E_f|\le2|E|$。$G_f$ 除了可能同时含 $(u,v)$ 与 $(v,u)$ 外，和流网络性质相同，可在其中定义关于 $c_f$ 的流。

流的增广（augmentation）：f 是 G 中的流，$f'$ 是 $G_f$ 中的流，定义

(f\uparrow f')(u,v)=\begin{cases}f(u,v)+f'(u,v)-f'(v,u)&(u,v)\in E\\0&\text{其他}\end{cases}\qquad(26.4)

在反向残存边上推流即抵消（cancellation）：u→v 送 5 箱、v→u 送 2 箱，等价于 u→v 送 3 箱。抵消对任何最大流算法都至关重要。

引理 26.1：$f\uparrow f'$ 是 G 中的流，且 $|f\uparrow f'|=|f|+|f'|$。证明要点：

容量下界：$f'(v,u)\le c_f(v,u)=f(u,v)$，故 $(f\uparrow f')(u,v)\ge f'(u,v)\ge0$。
容量上界：$\le f(u,v)+f'(u,v)\le f(u,v)+c_f(u,v)=c(u,v)$。
守恒：f 与 $f'$ 各自守恒，展开求和后流入等于流出。
流值：令 $V_1=\{v:(s,v)\in E\}$，$V_2=\{v:(v,s)\in E\}$，二者不交。代入展开、重组后把求和范围扩展到 V（多出来的项为 0），得 $|f|+|f'|$（式 26.5–26.7）。

增广路径：$G_f$ 中从 s 到 t 的简单路径 p。其残存容量 $c_f(p)=\min\{c_f(u,v):(u,v)\in p\}$。图 26.4(b) 中增广路径的残存容量为 $c_f(v_2,v_3)=4$。

引理 26.2：定义 $f_p(u,v)=c_f(p)$（若 (u,v) 在 p 上），否则为 0（26.8）。则 $f_p$ 是 $G_f$ 中的流，值 $|f_p|=c_f(p)>0$（证明为习题 26.2-7）。

推论 26.3：$f\uparrow f_p$ 是 G 中的流，值 $|f|+|f_p|>|f|$。图 26.4(c) 是把 26.1(b) 中值 19 的流沿该路径增广 4 后的流（值 23），(d) 是其残存网络。

割（cut）：流网络的割 $(S,T)$ 是 V 的划分，$T=V-S$，$s\in S,t\in T$（与第 23 章不同，这里是有向图且要求 s、t 分属两边）。

穿过割的净流量（net flow）：$f(S,T)=\sum_{u\in S}\sum_{v\in T}f(u,v)-\sum_{u\in S}\sum_{v\in T}f(v,u)$（26.9）。
割的容量：$c(S,T)=\sum_{u\in S}\sum_{v\in T}c(u,v)$（26.10），只算 S→T 方向。
最小割（minimum cut）：容量最小的割。定义的不对称是有意的：容量只算 S→T 的边，流算 S→T 减 T→S。例（图 26.5）：$S=\{s,v_1,v_2\}$，$T=\{v_3,v_4,t\}$，净流量 $f(v_1,v_3)+f(v_2,v_4)-f(v_3,v_2)=12+11-4=19$，容量 $c(v_1,v_3)+c(v_2,v_4)=12+14=26$。

引理 26.4：对任意割，$f(S,T)=|f|$。证明：在 $|f|$ 的定义上加上 $S-\{s\}$ 中各顶点的守恒式（每个都等于 0），重组求和，把对 V 的求和拆成对 S 与对 T；S 内部两项相同而抵消，剩下 $f(S,T)$。

推论 26.5：任意流的值不超过任意割的容量：$|f|=f(S,T)\le\sum_{u\in S,v\in T}f(u,v)\le c(S,T)$。因此最大流 ≤ 最小割容量。

定理 26.6（最大流最小割定理，Max-flow min-cut theorem）：以下三条等价：

f 是 G 的最大流；
残存网络 $G_f$ 不含增广路径；
存在某个割使 $|f|=c(S,T)$。证明：(1)⇒(2) 若有增广路径，由推论 26.3 可严格增大流值，矛盾。(2)⇒(3) 令 $S=\{v:G_f\text{ 中 s 可达 }v\}$，$T=V-S$，因 t 不可达所以是割。对 $u\in S,v\in T$：若 $(u,v)\in E$ 必有 $f(u,v)=c(u,v)$（否则 $(u,v)\in E_f$，v 会在 S 中）；若 $(v,u)\in E$ 必有 $f(v,u)=0$（否则 $c_f(u,v)=f(v,u)>0$）；故 $f(S,T)=c(S,T)$，再用引理 26.4。(3)⇒(1) 由推论 26.5，$|f|\le c(S,T)$ 对所有割成立，取等即最大。

基本 Ford-Fulkerson 算法：

def ford_fulkerson(G, s, t):
    for (u, v) in G.E:
        f[u, v] = 0
    while (p := find_path_in_residual(G, f, s, t)) is not None:   # DFS 或 BFS
        cf_p = min(cf(u, v) for (u, v) in p)
        for (u, v) in p:
            if (u, v) in G.E:
                f[u, v] += cf_p        # 原边：加流
            else:
                f[v, u] -= cf_p        # 反向边：抵消
    return f

图 26.6 给出示例运行：在图 26.1(a) 网络上经过若干次增广，最终最大流值 23，最后的残存网络中 s 到 t 无路径。

分析：

增广路径选得不好，算法可能不终止，流值甚至不收敛到最大值——这只会在容量为无理数时发生（脚注）。
整数容量（有理数可缩放为整数）：设 $f^*$ 为最大流，每轮流值至少增加 1，while 循环至多 $|f^*|$ 次。维护图 $G'=(V,E')$，$E'=\{(u,v):(u,v)\in E\text{ 或 }(v,u)\in E\}$，残存网络的边是 $G'$ 中 $c_f>0$ 的边，用 DFS/BFS 找路径 $O(V+E')=O(E)$。总时间 $O(E|f^*|)$。
坏例子（图 26.7）：s→u、s→v、u→t、v→t 容量各 1,000,000，u→v 容量 1。最大流 2,000,000。若交替选 $s\to u\to v\to t$ 与 $s\to v\to u\to t$（后者利用反向边），每次只增 1，要 2,000,000 次增广。

Edmonds-Karp 算法：用 BFS 找增广路径，即在残存网络中取（以边数计）最短的 s→t 路径。时间 $O(VE^2)$。记 $\delta_f(u,v)$ 为 $G_f$ 中单位边权的最短距离。

引理 26.7：Edmonds-Karp 运行时，对所有 $v\in V-\{s,t\}$，$\delta_f(s,v)$ 随每次增广单调不减。证明：反设某次增广（f→f'）使某些距离减小，取其中 $\delta_{f'}(s,v)$ 最小的 v。设 $G_{f'}$ 中最短路径 $s\leadsto u\to v$，$\delta_{f'}(s,u)=\delta_{f'}(s,v)-1$（26.12），且 u 的距离没减小：$\delta_{f'}(s,u)\ge\delta_f(s,u)$（26.13）。

若 $(u,v)\in E_f$：$\delta_f(s,v)\le\delta_f(s,u)+1\le\delta_{f'}(s,u)+1=\delta_{f'}(s,v)$，矛盾。
所以 $(u,v)\notin E_f$ 而 $(u,v)\in E_{f'}$，说明增广增加了 v→u 的流，即增广路径（最短路径）包含边 (v,u)，故 $\delta_f(s,v)=\delta_f(s,u)-1\le\delta_{f'}(s,u)-1=\delta_{f'}(s,v)-2$，也矛盾。

定理 26.8：Edmonds-Karp 的增广总次数为 $O(VE)$。证明：若 $c_f(p)=c_f(u,v)$，称 (u,v) 是增广路径 p 上的关键边（critical edge）；增广后关键边从残存网络消失，每条增广路径至少有一条关键边。(u,v) 第一次成为关键边时 $\delta_f(s,v)=\delta_f(s,u)+1$；它要再出现，必须先有 (v,u) 出现在某条增广路径上（流 f'），此时 $\delta_{f'}(s,u)=\delta_{f'}(s,v)+1\ge\delta_f(s,v)+1=\delta_f(s,u)+2$。所以两次成为关键边之间，u 的距离至少增加 2。u 的距离从 ≥0 开始，在 u 变得不可达之前至多 $|V|-2$（路径中间顶点不含 s、u、t），故每条边至多成为关键边 $|V|/2$ 次。残存网络中可能有边的顶点对共 $O(E)$ 个，关键边总次数 $O(VE)$。 总时间：每次增广 BFS $O(E)$，共 $O(VE^2)$；空间 $O(V+E)$。后面 26.4 节的推送-重贴标签给出 $O(V^2E)$，26.5 节 $O(V^3)$。

习题 26.2：

26.2-1 证明 (26.6) 与 (26.7) 的求和相等。
26.2-2 图 26.1(b) 中割 $(\{s,v_2,v_4\},\{v_1,v_3,t\})$ 的流与容量。
26.2-3 在图 26.1(a) 上运行 Edmonds-Karp。
26.2-4 图 26.6 中最大流对应的最小割；哪条增广路径抵消了流。
26.2-5 证明超级源汇构造中，原边容量有限则流值有限。
26.2-6 每个源恰好产出 $p_i$、每个汇恰好消耗 $q_j$（$\sum p_i=\sum q_j$）——把超级源边容量设为 $p_i$，汇边设为 $q_j$，检查最大流是否饱和。
26.2-7 证明引理 26.2。
26.2-8 残存网络不允许进入 s 的边时 Ford-Fulkerson 仍正确。
26.2-9 两个流的增广 $f\uparrow f'$ 是否满足守恒与容量约束（守恒满足，容量可能违反）。
26.2-10 证明最大流可以由至多 |E| 条增广路径得到（求出最大流后分解）。
26.2-11 用至多 |V| 次最大流求无向图的边连通度（edge connectivity）（树为 1，环为 2）。
26.2-12 有进入 s 的边且 $f(v,s)=1$ 时，$O(E)$ 构造等值且 $f'(v,s)=0$ 的流。
26.2-13 在所有最小割中找边数最少的：修改容量（如 $c'=c\cdot(|E|+1)+1$）。

26.3 二分图最大匹配（Maximum bipartite matching）（PDF p.753–757）

定义：无向图 G 中的匹配（matching）是边子集 $M\subseteq E$，每个顶点至多与 M 中一条边关联。被 M 中某边关联的顶点称为已匹配（matched），否则为未匹配。**最大匹配（maximum matching）**是基数最大的匹配。二分图（bipartite graph）：$V=L\cup R$，L、R 不交，所有边在 L 与 R 之间；假定每个顶点至少有一条边。应用：L 为机器，R 为可同时执行的任务，边 (u,v) 表示机器 u 能做任务 v，最大匹配让尽可能多的机器有活干。图 26.8：(a) 基数 2 的匹配，(b) 基数 3 的最大匹配。

构造流网络 $G'=(V',E')$：$V'=V\cup\{s,t\}$，

E'=\{(s,u):u\in L\}\cup\{(u,v):(u,v)\in E,\ \text{方向 } L\to R\}\cup\{(v,t):v\in R\},

所有边容量为 1。因每个顶点至少有一条边，$|E|\ge|V|/2$，故 $|E|\le|E'|=|E|+|V|\le3|E|$，$|E'|=\Theta(E)$。

称流 f 整数值（integer-valued），若所有 $f(u,v)$ 都是整数。

引理 26.9：若 M 是 G 的匹配，则 $G'$ 中存在整数值流 f，$|f|=|M|$；反之，若 f 是 $G'$ 中整数值流，则存在匹配 M，$|M|=|f|$。证明：正向——对 $(u,v)\in M$ 令 $f(s,u)=f(u,v)=f(v,t)=1$，其余为 0；这些 $s\to u\to v\to t$ 路径除 s、t 外顶点不交，穿过割 $(L\cup\{s\},R\cup\{t\})$ 的净流量为 $|M|$。反向——令 $M=\{(u,v):u\in L,v\in R,f(u,v)>0\}$。每个 $u\in L$ 只有一条入边 (s,u)，容量 1，整数流下至多 1 单位流入、由守恒恰从一条边流出；R 侧对称。所以 M 是匹配，且穿过上述割的净流量等于 $|M|$。

定理 26.10（整数性定理，Integrality theorem）：若容量全为整数，Ford-Fulkerson 求出的最大流 f 满足 $|f|$ 为整数，且所有 $f(u,v)$ 为整数（对迭代次数归纳，习题 26.3-2）。

推论 26.11：二分图最大匹配的基数 = 对应流网络最大流的值。证明：若最大匹配 M 对应的流不是最大流，则存在更大的最大流 $f'$，由整数性定理可取整数值，对应更大的匹配，矛盾；反向同理。

算法与复杂度：构造 $G'$，跑 Ford-Fulkerson，直接从整数最大流读出匹配。任何匹配基数至多 $\min(|L|,|R|)=O(V)$，故最大流值 $O(V)$，时间 $O(VE')=O(VE)$。（更快的 Hopcroft-Karp $O(\sqrt VE)$ 见思考题 26-6。）

习题 26.3：

26.3-1 在图 26.8(c) 上手算，每次取字典序最小的增广路径。
26.3-2 证明整数性定理。
26.3-3 增广路径长度上界（$2\min(|L|,|R|)+1$）。
26.3-4★ Hall 定理：$|L|=|R|$ 时存在完美匹配（perfect matching） ⇔ 对每个 $A\subseteq L$，$|A|\le|N(A)|$，其中 $N(X)$ 为 X 的邻域。
26.3-5★ d-正则二分图（每个顶点度为 d，必有 $|L|=|R|$）存在基数 |L| 的匹配：证明对应网络的最小割容量为 |L|。

★26.4 推送-重贴标签算法（Push-relabel algorithms）（PDF p.757–769）

定位：迄今许多渐近最快的最大流算法以及最快的实际实现都基于推送-重贴标签方法，它也能高效求解最小费用流等问题。本节介绍 Goldberg 的“通用（generic）”最大流算法，简单实现 $O(V^2E)$，优于 Edmonds-Karp 的 $O(VE^2)$；26.5 节细化为 $O(V^3)$。

与 Ford-Fulkerson 的区别：更“局部”——不在整个残存网络中找增广路径，而是一次只处理一个顶点，只看它在残存网络中的邻居；执行过程中不保持流量守恒，只维护预流（preflow）：满足容量约束，以及守恒的放松形式

\sum_{v}f(v,u)-\sum_{v}f(u,v)\ge0\quad\forall u\in V-\{s\},

即流入可以超过流出。超额流（excess flow）

e(u)=\sum_{v}f(v,u)-\sum_vf(u,v)\qquad(26.14)

若 $u\in V-\{s,t\}$ 且 $e(u)>0$，称 u 溢出（overflowing）。

直观（流体类比）：边是管道，顶点是接头；每个顶点有一根通向任意大蓄水池的出流管来存放超额流；每个顶点连同蓄水池和管道坐落在一个平台上，平台高度随算法推进而升高。只向下坡推流（从高处到低处）。源点高度固定为 |V|，汇点固定为 0，其余初始为 0。开始时从源点尽量往下推，恰好把源点每条出边填满（即割 $(s,V-\{s\})$ 的容量），流进中间顶点后先积在蓄水池里，再逐步往下推。若某溢出顶点 u 所有未饱和出管都通向同高或更高的顶点，就要重贴标签（relabel）——把 u 的高度提升到“它有未饱和管道相连的最低邻居高度 + 1”，这样至少有一根出管可以继续推流。最终所有能到达汇点的流都到了（受割容量限制，不能再多），为把预流变成合法的流，算法继续抬高溢出顶点到超过源点高度 |V|，把蓄水池中多余的流送回源点。所有蓄水池清空后，预流不仅是合法的流，而且是最大流。

高度函数（height function）：$h:V\to\mathbb N$，满足 $h(s)=|V|$、$h(t)=0$，且对每条残存边 $(u,v)\in E_f$，$h(u)\le h(v)+1$。（文献中常称为“距离函数”和“距离标号”；高度与在转置图 $G^T$ 中到汇点 t 的 BFS 距离有关。）

引理 26.12：若 $h(u)>h(v)+1$，则 (u,v) 不是残存边。

PUSH(u,v)——适用条件：u 溢出、$c_f(u,v)>0$、$u.h=v.h+1$。

def push(u, v):
    delta = min(u.e, cf(u, v))
    if (u, v) in E:
        f[u, v] += delta          # 原边加流
    else:
        f[v, u] -= delta          # 反向边：减少原流
    u.e -= delta
    v.e += delta

O(1)。推送后 f 仍是预流。只允许高度差恰好为 1 的推送——由引理 26.12，高度差大于 1 的顶点间本来就没有残存边。推送后若 $c_f(u,v)=0$ 称为饱和推送（saturating push），该边从残存网络消失；否则为非饱和推送（nonsaturating push）。

引理 26.13：非饱和推送后 u 不再溢出（推送量等于 u.e）。

RELABEL(u)——适用条件：u 溢出，且对所有残存边 $(u,v)\in E_f$ 有 $u.h\le v.h$（无法向下推）。s、t 按定义不会溢出，不会被重贴标签。

def relabel(u):
    u.h = 1 + min(v.h for v in residual_neighbors(u))

最小值集合非空：u 溢出说明有某个 v 使 $f(v,u)>0$，于是 $c_f(u,v)>0$。RELABEL 给 u 的是高度函数约束允许的最大高度。

初始化预流：

def initialize_preflow(G, s):
    for v in G.V:
        v.h = 0; v.e = 0
    for (u, v) in G.E:
        f[u, v] = 0
    s.h = len(G.V)
    for v in G.Adj[s]:
        f[s, v] = c(s, v)         # 源点出边全部填满
        v.e = c(s, v)
        s.e -= c(s, v)

即 $f(u,v)=c(u,v)$ 若 $u=s$，否则 0（26.15）；$h(u)=|V|$ 若 $u=s$，否则 0（26.16）。这是高度函数：唯一满足 $h(u)>h(v)+1$ 的边是源点出边，而它们已饱和，不在残存网络中。

GENERIC-PUSH-RELABEL(G):
    INITIALIZE-PREFLOW(G, s)
    while 存在可执行的 push 或 relabel 操作:
        任选一个可执行的操作并执行

引理 26.14：对任意溢出顶点 u，要么可推送，要么可重贴标签。（若不能推送，则所有残存边满足 $h(u)<h(v)+1$ 即 $h(u)\le h(v)$。）

正确性：

引理 26.15：顶点高度从不减少；每次重贴标签至少加 1。
引理 26.16：算法始终维持 h 为高度函数。RELABEL(u) 后 u 的出残存边满足 $u.h\le v.h+1$；入残存边 (w,u) 原有 $w.h\le u.h+1$，u 升高后更满足。PUSH(u,v) 可能新增残存边 (v,u)，而 $v.h=u.h-1<u.h+1$；删除残存边只是去掉约束。
引理 26.17：若 h 是高度函数，则残存网络中不存在 s 到 t 的路径。证明：若有简单路径 $\langle v_0=s,\dots,v_k=t\rangle$，$k<|V|$，沿路径累加 $h(v_i)\le h(v_{i+1})+1$ 得 $h(s)\le h(t)+k=k<|V|$，矛盾。
定理 26.18（通用算法正确性）：若算法终止，求得的预流是最大流。循环不变式“每次 while 测试时 f 是预流”。终止时不再有溢出顶点（引理 26.14），故 f 是流；h 仍是高度函数，由引理 26.17 残存网络无 s→t 路径；由最大流最小割定理 f 是最大流。

运行时间分析（分别界定重贴标签、饱和推送、非饱和推送三类操作）：

引理 26.19：对任意溢出顶点 x，残存网络中存在 x 到 s 的简单路径。证明：设 U 为 x 在 $G_f$ 中可达的顶点集，反设 $s\notin U$。对 U 中超额求和并展开，得 $\sum_{u\in U}e(u)=\sum_{u\in U}\sum_{v\in\bar U}f(v,u)-\sum_{u\in U}\sum_{v\in\bar U}f(u,v)$；左边 > 0（$e(x)>0$，除 s 外超额均非负），所以存在 $u'\in U,v'\in\bar U$ 使 $f(v',u')>0$，则残存边 $(u',v')$ 存在，$v'$ 也可达，矛盾。（直觉：溢出的流来自源点，可以沿原路退回。）
引理 26.20：任何时刻所有顶点 $u.h\le2|V|-1$。重贴标签时 u 溢出，有到 s 的简单路径（长度 $\le|V|-1$），沿路径累加得 $u.h\le s.h+|V|-1=2|V|-1$。
推论 26.21（重贴标签次数）：每个顶点至多 $2|V|-1$ 次，总计至多 $(2|V|-1)(|V|-2)<2|V|^2$。
引理 26.22（饱和推送次数）：少于 $2|V||E|$。对每对 u,v：u→v 饱和推送时 $v.h=u.h-1$；要再次 u→v 推送，必须先 v→u 推送，需要 $v.h=u.h+1$，即 v.h 至少增加 2；反过来同理。高度不超过 $2|V|-1$，每个顶点高度增加 2 的次数少于 |V|，所以每对之间饱和推送少于 $2|V|$ 次，乘以边数得界。
引理 26.23（非饱和推送次数）：少于 $4|V|^2(|V|+|E|)$。用势函数（potential function） $\Phi=\sum_{v:e(v)>0}v.h$。初始 0。重贴标签使 Φ 增加少于 $2|V|$；饱和推送使 Φ 增加少于 $2|V|$（只可能让 v 新溢出）；非饱和推送使 u 不再溢出（减少 u.h），v 可能新溢出（增加 v.h=u.h−1），净减少至少 1。总增量 $<(2|V|)(2|V|^2)+(2|V|)(2|V||E|)=4|V|^2(|V|+|E|)$，而 $\Phi\ge0$，故非饱和推送次数不超过该数。
定理 26.24：基本操作总数 $O(V^2E)$。
推论 26.25：存在 $O(V^2E)$ 的实现（每次重贴标签 O(V)、每次推送 O(1)、O(1) 选择可执行操作，习题 26.4-2）。

习题 26.4：

26.4-1 证明初始化后 $s.e\le-|f^*|$。
26.4-2 实现：重贴标签 O(V)、推送 O(1)、选择 O(1)，总 $O(V^2E)$。
26.4-3 所有 $O(V^2)$ 次重贴标签总共只花 $O(VE)$ 时间。
26.4-4 由推送-重贴标签得到的最大流快速找最小割。
26.4-5 用推送-重贴标签求二分图最大匹配并分析。
26.4-6 容量都在 $\{1..k\}$ 时的运行时间（每条边饱和前至多支持 k 次非饱和推送）。
26.4-7 把 s.h 初始化为 |V|−2 不影响正确性和渐近性能。
26.4-8 证明：$u.h<|V|$ 时 $u.h\le\delta_f(u,t)$；$u.h\ge|V|$ 时 $u.h-|V|\le\delta_f(u,s)$。
26.4-9★ 维护上述为等式，额外总时间 $O(VE)$（精确距离标号）。
26.4-10 $|V|\ge4$ 时非饱和推送至多 $4|V|^2|E|$ 次。

★26.5 前置重贴标签算法（The relabel-to-front algorithm）（PDF p.769–781）

推送-重贴标签允许任意顺序执行基本操作；精心选择顺序并高效管理数据结构，可以突破 $O(V^2E)$。前置重贴标签算法时间 $O(V^3)$，渐近上不差于 $O(V^2E)$，在稠密网络上更好。

思路：维护顶点链表 L，从表头扫描，每次选一个溢出顶点 u 并**释放（discharge）**它——反复推送与重贴标签直到 u 无超额。一旦某顶点被重贴标签，就把它移到表头（算法名由此而来），并从它之后继续扫描。

可容许边（admissible edge）：$c_f(u,v)>0$ 且 $h(u)=h(v)+1$；否则不可容许。可容许网络 $G_{f,h}=(V,E_{f,h})$，由能推流的边组成。

引理 26.26：可容许网络无环（DAG）。证明：若有环，沿环累加 $h(v_{i-1})=h(v_i)+1$ 得 $0=k$，矛盾。

引理 26.27：u 溢出且 (u,v) 可容许，则 PUSH(u,v) 适用；推送不会产生新的可容许边（唯一可能新增的残存边 (v,u) 有 $v.h=u.h-1$，不可容许），但饱和推送会使 (u,v) 变为不可容许。

引理 26.28：u 溢出且无可容许出边，则 RELABEL(u) 适用；重贴标签后 u 至少有一条可容许出边（取到最小值的那个邻居），但没有可容许入边（若 (v,u) 可容许则重贴前 $v.h>u.h+1$，由引理 26.12 不存在该残存边；重贴标签不改变残存网络）。

邻居表（neighbor lists）：$u.N$ 是单链表，包含所有满足 $(u,v)\in E$ 或 $(v,u)\in E$ 的 v，即可能存在残存边 (u,v) 的所有顶点。$u.N.head$ 指向表头，$v.next\text{-}neighbor$ 指向下一个。每个顶点有指针 $u.current$ 指向当前考察的邻居，初始为表头；扫描顺序任意但固定。

DISCHARGE(u)：

def discharge(u):
    while u.e > 0:
        v = u.current
        if v is None:                       # 走到邻居表末尾
            relabel(u)
            u.current = u.N.head
        elif cf(u, v) > 0 and u.h == v.h + 1:   # 可容许边
            push(u, v)                      # current 不前进
        else:
            u.current = v.next_neighbor

每次迭代恰好做三件事之一：到表尾则重贴标签并重置指针；当前边可容许则推送；否则指针后移。DISCHARGE 对溢出顶点调用时，最后一个动作一定是推送（只有推送改变 u.e）。

例（图 26.9）：释放顶点 y（初始超额 19，高度 0），邻居表 ⟨s, x, z⟩，共 15 次迭代：迭代 1–3 发现无可容许边，迭代 4 指针为 NIL，重贴标签到高度 1；迭代 7 向 z 推 8 单位（饱和 (y,z)）；迭代 9 再次重贴标签；迭代 11 向 x 推 5 单位；迭代 14 第三次重贴标签；迭代 15 向 s 推回 6 单位，y 超额为 0，结束。

引理 26.29：DISCHARGE 调用 PUSH 时推送确实适用；调用 RELABEL 时重贴标签确实适用。后者的关键：一个完整的“遍”（指针从表头走到 NIL）中，每条出边都曾被判为不可容许；而在这一遍期间，推送不会制造可容许边（引理 26.27），u 自己没被重贴标签，其他被重贴标签的顶点也没有可容许入边（引理 26.28），所以遍结束时 u 的所有出边仍不可容许。

算法：

def relabel_to_front(G, s, t):
    initialize_preflow(G, s)
    L = LinkedList(v for v in G.V if v not in (s, t))   # 任意顺序
    for u in L:
        u.current = u.N.head
    u = L.head
    while u is not None:
        old_height = u.h
        discharge(u)
        if u.h > old_height:            # 被重贴标签了
            L.move_to_front(u)
        u = u.next                      # 若已移到表头，则继续处理其后继

例（图 26.10）：源点 s 初始送出 26 单位（$s.e=-26$），L = ⟨x, y, z⟩。(a) 释放 x：重贴标签到 1，推 5 给 y、推 7 给 t，移到表头（结构不变）。(b) 释放 y（即图 26.9 过程），y 被重贴标签，移到表头。(c) x 在 y 后，再次释放，把 5 单位全推给 t，未重贴标签，位置不变。(d) 释放 z：重贴标签到 1，8 单位推给 t，移到表头。(e) y、x 都无超额，DISCHARGE 立即返回，到达表尾结束；此时无溢出顶点，预流即最大流，t 收到 20 单位，$s.e=-20$。

正确性（循环不变式）：每次测试 while 条件时，L 是可容许网络的一个拓扑排序，且 L 中位于 u 之前的顶点都没有超额。

初始化：初始化后除 s 外高度全为 0，s.h = |V| ≥ 2，没有可容许边，任何顺序都是拓扑序；u 是表头，前面没有顶点。
保持：只有重贴标签能产生可容许边（引理 26.27）；重贴标签后 u 无可容许入边、可能有可容许出边（引理 26.28），把 u 移到表头即保持拓扑序。若 u 被重贴标签，下一轮的 $u'$ 前面只有 u 本身（已无超额）；若未被重贴标签，释放期间 L 始终拓扑有序，推送只把超额推往表中更靠后的顶点（或 s、t），前面的顶点不会获得超额。
终止：u 越过表尾，所有顶点超额为 0，不再有基本操作适用。它是通用算法的一个实现，因此求得最大流。

定理 26.30：RELABEL-TO-FRONT 在任意流网络上运行时间为 $O(V^3)$。证明：定义“阶段”为两次相邻重贴标签之间的时间，共 $O(V^2)$ 个阶段。每个阶段至多 |V| 次 DISCHARGE 调用（不重贴标签则下一个调用在表中更靠后，表长 < |V|；重贴标签则进入新阶段）。故 DISCHARGE 共调用 $O(V^3)$ 次，主循环自身开销 $O(V^3)$。DISCHARGE 内部：重贴标签总计 $O(VE)$（习题 26.4-3）；指针前移每次重贴标签后 $O(\deg u)$，每个顶点总共 $O(V\deg u)$，由握手引理总计 $O(VE)$；饱和推送 $O(VE)$ 次；非饱和推送会使超额变 0 并立即返回，每次调用至多一次，总计 $O(V^3)$。合计 $O(V^3+VE)=O(V^3)$。空间 $O(V+E)$。

习题 26.5：

26.5-1 在图 26.1(a) 上按给定邻居表手算。
26.5-2★ FIFO 队列维护溢出顶点的推送-重贴标签，$O(V^3)$ 实现。
26.5-3 RELABEL 改为 $u.h=u.h+1$ 仍正确，对分析的影响。
26.5-4★ 总是释放最高的溢出顶点（highest-label），可达 $O(V^3)$。
26.5-5 间隙启发式（gap heuristic）：若存在 $0<k\le|V|-1$ 没有任何顶点高度为 k，则所有高度 > k 的顶点都在某个最小割的源点一侧；把这些顶点（s 除外）高度设为 $\max(v.h,|V|+1)$ 后 h 仍是高度函数。这对推送-重贴标签实际性能至关重要。

第 26 章思考题（PDF p.781–786）

26-1 逃脱问题（Escape problem）：$n\times n$ 网格（内部顶点 4 个邻居），给定 $m\le n^2$ 个起点，问能否找到 m 条顶点不相交的路径从起点到 m 个不同的边界点（图 26.11：(a) 可以逃脱，(b) 不能）。(a) 同时有边容量与顶点容量的流网络可归约为普通最大流（拆点）；(b) 设计高效算法并分析（拆点后建网络，超级源连起点，边界点连超级汇，看最大流是否为 m）。
26-2 最小路径覆盖（Minimum path cover）：有向图的路径覆盖是一组顶点不相交的路径，每个顶点恰在一条路径上（长度可为 0）。(a) DAG 的最小路径覆盖：构造 $V'=\{x_0..x_n\}\cup\{y_0..y_n\}$，$E'=\{(x_0,x_i)\}\cup\{(y_i,y_0)\}\cup\{(x_i,y_j):(i,j)\in E\}$，跑最大流，最小覆盖数 = $n-$ 最大流；(b) 对有环图是否有效（否，会得到环）。
26-3 算法咨询公司（Algorithmic consulting）：n 个子领域 $A_k$，雇专家费用 $c_k$；m 个项目 $J_i$，需要子领域集合 $R_i$，收入 $p_i$；专家可同时做多个项目。目标：最大化净收入。网络：$s\to A_k$ 容量 $c_k$，$J_i\to t$ 容量 $p_i$，$A_k\to J_i$（若 $A_k\in R_i$）容量 ∞。(a) 有限容量割中若 $J_i\in T$ 则其所需 $A_k\in T$；(b) 最大净收入 $=\sum p_i-$ 最小割容量；(c) 给出接受哪些项目、雇哪些专家的算法并用 m、n、$r=\sum|R_i|$ 表示时间。（这是经典的项目选择/闭包问题。）
26-4 更新最大流：整数容量、已知最大流。(a) 某边容量加 1，$O(V+E)$ 更新（找一次增广路径）；(b) 某边容量减 1，$O(V+E)$ 更新（若该边原本满载，先沿路径退回 1 单位，再尝试增广）。
26-5 缩放最大流（Maximum flow by scaling）：$C=\max c(u,v)$。(a) 最小割容量至多 $C|E|$；(b) $O(E)$ 找残存容量 ≥ K 的增广路径（只保留 $c_f\ge K$ 的边做 BFS）；算法：K 从 $2^{\lfloor\lg C\rfloor}$ 开始，每阶段反复增广容量 ≥ K 的路径，然后 K 减半；(c) 证明返回最大流；(d) 每次检查外层循环条件时残存网络最小割容量 ≤ $2K|E|$；(e) 每个 K 内层循环 $O(E)$ 次；(f) 总时间 $O(E^2\lg C)$。
26-6 Hopcroft-Karp 二分图匹配算法，$O(\sqrt VE)$。关于匹配 M 的增广路径：从 L 中未匹配顶点出发、到 R 中未匹配顶点结束、边交替属于 $E-M$ 与 M 的简单路径（与流网络的增广路径相关但不同）。对称差 $A\oplus B=(A-B)\cup(B-A)$。(a) $M\oplus P$ 是匹配且大小加 1；k 条顶点不相交增广路径同时翻转则加 k。算法：反复找一个极大的顶点不相交最短增广路径集，全部翻转，直到没有增广路径。(b) 两个匹配的对称差中每个顶点度 ≤ 2，是不相交的简单路径与环的并，边交替；若 $|M|\le|M^*|$，则 $M\oplus M^*$ 至少含 $|M^*|-|M|$ 条关于 M 的顶点不相交增广路径。(c)(d) 设最短增广路径长 l，翻转一个极大集后新的最短增广路径长 > l。(e) 若最短增广路径长为 l，最大匹配至多 $|M|+|V|/(l+1)$。(f) repeat 循环至多 $2\sqrt{|V|}$ 次（$\sqrt{|V|}$ 次迭代后最短路径长 ≥ $\sqrt{|V|}$，剩余至多 $\sqrt{|V|}$ 次增长）。(g) $O(E)$ 找极大集（BFS 分层 + DFS），总时间 $O(\sqrt VE)$。

第 26 章注记（PDF p.786–787）

参考书：Ahuja–Magnanti–Orlin、Even、Lawler、Papadimitriou–Steiglitz、Tarjan；Goldberg–Tardos–Tarjan 综述；Schrijver 的历史回顾。Ford-Fulkerson 方法由 Ford 与 Fulkerson 提出，开创了网络流的形式化研究（含最大流与二分匹配）。Edmonds–Karp 与 Dinic 独立证明 BFS 增广是多项式的；Dinic 还提出阻塞流（blocking flow）。Karzanov 首创预流；推送-重贴标签由 Goldberg、Goldberg–Tarjan 提出，后者给出基于队列的 $O(V^3)$ 算法和基于动态树的 $O(VE\lg(V^2/E+2))$ 算法。之后有缩放（Ahuja–Orlin、Ahuja–Orlin–Tarjan）、最高标号优先（Cheriyan–Maheshwari）、随机排列邻居表及其去随机化，King–Rao–Tarjan 达 $O(VE\log_{E/(V\lg V)}V)$。迄今渐近最快的是 Goldberg–Rao 的 $O(\min(V^{2/3},E^{1/2})E\lg(V^2/E+2)\lg C)$，基于阻塞流，并给高容量边赋长度 0、低容量边赋长度 1，使最短路径倾向高容量，减少迭代。实践中推送-重贴标签优于增广路径或线性规划方法；Cherkassky–Goldberg 指出两个关键启发式：定期在残存网络上做 BFS 以获得精确高度（全局重贴标签），以及间隙启发式；最佳变体是释放最高的溢出顶点。二分图最大匹配最好的算法是 Hopcroft–Karp $O(\sqrt VE)$。匹配问题参考 Lovász–Plummer。

第 26 章小结

本章要点 流网络由容量约束和流量守恒定义，流值等于穿过任意割的净流量，且不超过任意割的容量。最大流最小割定理把三件事等同起来：f 是最大流；残存网络无增广路径；存在割使流值等于割容量。残存网络中的反向边允许“抵消”已有的流，这是所有最大流算法的关键。Ford-Fulkerson 在整数容量下为 $O(E|f^*|)$，用 BFS 增广（Edmonds-Karp）后为 $O(VE^2)$，关键在于残存网络中的 BFS 距离单调不减、每条边成为关键边至多 $|V|/2$ 次。整数性定理保证整数容量下得到整数流，从而二分图最大匹配可归约为最大流（$O(VE)$）。推送-重贴标签方法维护预流与高度函数，只向下坡推送，通用版本 $O(V^2E)$，前置重贴标签版本 $O(V^3)$；分析依靠高度上界 $2|V|-1$ 和势函数。建模技巧：反平行边拆点、多源多汇加超级源汇、顶点容量拆点、项目选择问题化为最小割。

与量化交易的关联

最小割/闭包问题与组合选择：思考题 26-3 的“项目选择”模型可直接用于有依赖关系的选择问题，例如选择一组策略或因子时，某些策略依赖共同的数据源/基础设施（有固定成本），求净收益最大的子集，用一次最小割精确求解。
资金调拨与清算网络：在多个账户、交易所、托管行之间转移资金或证券，各通道有额度上限，最大流给出单位时间最大可调拨量，最小割指出瓶颈通道；多源多汇构造对应多个资金来源与多个需求方。
订单撮合与分配：二分图匹配可用于把一批委托分配到不同券商通道/账户（每个通道容量有限），或在内部交叉撮合（crossing）中匹配买卖单；带容量的版本就是流问题。更一般的“带成本”版本（最小费用流）用于最优订单路由，本章的残存网络和推送-重贴标签是其基础。
运输问题与线性规划：最大流是线性规划的特例（习题 26.1-5），理解其对偶（最小割）有助于理解第 29 章 LP 对偶与组合优化中的影子价格。
风险传染分析：银行间敞口网络中，最小割可用于识别“切断后能隔离违约传染”的最小敞口集合，属于系统性风险研究中的工具。
对日常因子研究、回测本身，最大流没有直接用处；它主要出现在执行、清算与资源分配这类有网络结构的环节。

推荐习题

26.1-5（最大流写成 LP）、26.1-7（顶点容量拆点）——建模基本功。
26.2-4、26.2-11（最小割与边连通度）；26.2-10（流分解）。
26.3-4（Hall 定理）。
26.4-4（由推送-重贴标签结果求最小割）；26.5-5（间隙启发式，实现必用）。
思考题 26-3（项目选择/最小割建模，强烈推荐）；26-5（容量缩放）；26-6（Hopcroft-Karp）。

第七部分算法问题选编（Part VII Selected Topics）

第七部分导言（PDF p.789–792）

本部分精选若干扩展前文的专题：有的引入新的计算模型（电路、并行计算机），有的涉及专门领域（计算几何、数论），最后两章讨论高效算法设计的已知局限及应对技术。各章内容：

第 27 章多线程算法：基于动态多线程的并行计算模型，以工作量（work）和跨度（span）度量并行性，介绍矩阵乘法与归并排序的多线程算法。
第 28 章矩阵运算：用 LU 与 LUP 分解以 $O(n^3)$ 时间通过高斯消元解线性方程组；证明矩阵求逆与矩阵乘法一样快；无精确解时求最小二乘近似解。
第 29 章线性规划：在有限资源和相互竞争的约束下最大化/最小化目标；讲如何建模与求解，方法是最古老的单纯形算法——最坏情况非多项式，但实践中相当高效、广泛使用。
第 30 章多项式与快速傅里叶变换（FFT）：$O(n\lg n)$ 时间乘两个 n 次多项式，包括高效实现与并行电路。
第 31 章数论算法：欧几里得 gcd、模线性方程、模幂；RSA 公钥密码系统（加密与数字签名）；Miller-Rabin 随机素性测试；Pollard rho 因数分解启发式。
第 32 章字符串匹配：朴素法、Rabin-Karp、有限自动机、Knuth-Morris-Pratt。
第 33 章计算几何：基本原语、扫描线判定线段相交、凸包（Graham 扫描、Jarvis 步进）、平面最近点对。
第 34 章 NP 完全性：判定 NP 完全的技术；哈密顿回路、布尔可满足性、子集和、旅行商问题等经典 NP 完全问题。
第 35 章近似算法：顶点覆盖（无权/加权）、MAX-3-CNF、旅行商、集合覆盖、子集和。

第 27 章多线程算法（Multithreaded Algorithms）

第 27 章引言（PDF p.793–795）

本书多数算法是串行算法（serial algorithms），适合单处理器一次执行一条指令。本章扩展到并行算法（parallel algorithms），在允许多条指令并发执行的多处理器上运行，重点是便于设计、分析且实践中可高效实现的**动态多线程（dynamic multithreading）**模型。

并行计算机：从廉价的单芯片多核（chip multiprocessor，多个“核”共享内存），到由普通 PC 通过专用网络互连的集群（cluster），再到使用定制架构和网络的超级计算机。串行计算早已统一于 RAM 模型，但并行计算没有广泛公认的单一模型，主要因为厂商架构不统一：共享内存（shared memory）——每个处理器可直接访问任意内存位置；分布式内存（distributed memory）——各处理器内存私有，须显式发消息访问他人内存。多核普及后，每台新笔记本和台式机都是共享内存并行机，趋势倾向共享内存，本章采用此模型。

静态线程（static threading）：软件抽象出共享内存的“虚拟处理器”即线程，各有程序计数器，可独立执行；操作系统负责装载和切换。创建与销毁线程较慢，所以线程通常在整个计算期间存在，故称“静态”。直接用静态线程编程困难且易错，尤其是负载均衡——动态地把工作均匀分给各线程需要复杂的通信协议和调度器。于是出现了并发平台（concurrency platforms）：一层负责协调、调度、管理并行资源的软件，可以是运行时库，也可以是带编译器和运行时支持的完整并行语言。

动态多线程：程序员只需指定应用中的逻辑并行性，不必处理通信协议与负载均衡；平台的**调度器（scheduler）**自动负载均衡。几乎所有动态多线程环境都支持两个特性：嵌套并行（nested parallelism）——可以“派生（spawn）”子程序，调用者继续执行而子程序同时计算；并行循环（parallel loops）——迭代可并发执行的 for 循环。

该模型的优点：

是串行编程模型的简单扩展：伪代码只加三个并发关键字 parallel、spawn、sync；删去它们得到同一问题的串行伪代码，称为多线程算法的串行化（serialization）。
基于“工作量”与“跨度”可以在理论上干净地量化并行性。
许多嵌套并行算法自然来自分治法，并可像串行分治一样用递归式分析。
贴合并行计算实践的发展：Cilk、Cilk++、OpenMP、Task Parallel Library、Threading Building Blocks 等都支持某种动态多线程。

章节：27.1 模型与工作量、跨度、并行度指标；27.2 多线程矩阵乘法；27.3 多线程归并排序。

27.1 动态多线程基础（The basics of dynamic multithreading）（PDF p.795–813）

引例：斐波那契数。$F_0=0,F_1=1,F_i=F_{i-1}+F_{i-2}$。

def fib(n):                 # 串行
    if n <= 1:
        return n
    x = fib(n - 1)
    y = fib(n - 2)
    return x + y

图 27.1 是 FIB(6) 的递归树：FIB(5) 和 FIB(6) 都调用 FIB(4)，大量重复计算（未做备忘）。$T(n)=T(n-1)+T(n-2)+\Theta(1)$，用代入法（假设 $T(n)\le aF_n-b$，b 足够大以吸收常数）得

T(n)=\Theta(\phi^n),\qquad\phi=(1+\sqrt5)/2\ \text{（黄金比）}\qquad(27.1)

指数级，是很差的算法（思考题 31-3 有快得多的方法），但适合说明多线程分析的概念：两次递归调用相互独立，可以并行。

P-FIB(n):
    if n <= 1: return n
    x = spawn P-FIB(n-1)     # 子任务可与父任务并行
    y = P-FIB(n-2)           # 父任务继续执行
    sync                     # 等待所有派生的子任务完成
    return x + y

串行化：删去 spawn、sync（以及 parallel）后得到串行算法，求解同一问题。
嵌套并行：spawn 放在过程调用前。执行 spawn 的**父（parent）可以与被派生的子（child）**并行，而不是像串行调用那样等待。递归下去形成巨大的并行子计算树。
spawn 表示“可以”并行而非“必须”并行；关键字表达的是逻辑并行性，哪些子计算真正并发由运行时调度器决定。
sync：过程必须等待所有已派生的子任务完成才能继续，才能安全使用它们的返回值（否则可能在 x 算出前就相加）。每个过程返回前隐式执行 sync。

多线程执行模型——计算 DAG（computation dag） $G=(V,E)$：顶点是指令，边 $(u,v)$ 表示 u 必须先于 v 执行。把不含并行控制（spawn、sync、从 spawn 返回）的指令链合并为一个链（strand）。若一个链有两个后继，其中之一必是被派生的；有多个前驱说明在 sync 处汇合。若 DAG 中有 u 到 v 的有向路径，u、v （逻辑上）串行，否则**（逻辑上）并行**。

边的分类（图 27.2，P-FIB(4)）：延续边（continuation edge） $(u,u')$ 连接同一过程实例内相继的链（水平向右）；派生边（spawn edge） $(u,v)$（向下）；调用边（call edge）（向下）——派生与调用的区别是派生还产生一条从 u 到其后继 $u'$ 的延续边，表示 $u'$ 可与 v 同时执行；返回边（return edge） $(u,x)$，x 是调用过程中下一个 sync 之后的链（向上）。计算从单一初始链开始、到单一最终链结束。

理想并行计算机（ideal parallel computer）：一组处理器 + **顺序一致（sequentially consistent）**的共享内存——实际可能同时有许多读写，但结果如同每一步只有一个处理器执行一条指令，且全局线性顺序保持每个处理器自身的指令顺序。对动态多线程计算，相当于指令被交错成一个与计算 DAG 偏序一致的线性序；不同运行的顺序可以不同。性能假设：各处理器计算能力相同，忽略调度开销（对并行度充分的算法，实践中调度开销通常很小）。

性能指标：

工作量（work） $T_1$：在单处理器上执行整个计算的总时间，即各链时间之和；单位时间链时等于 DAG 顶点数。
跨度（span） $T_\infty$：DAG 中最长路径上链的执行时间之和，单位时间链时等于**关键路径（critical path）**上的顶点数（可用 24.2 节在 $\Theta(V+E)$ 内求）。图 27.2 中工作量 17、跨度 8。
$T_P$：P 个处理器上的运行时间；$T_1$ 即工作量，$T_\infty$ 即无限多处理器时的时间。

两条下界：

\textbf{工作量定律（work law）：}\ T_P\ge T_1/P\qquad(27.2)

（P 个处理器每步至多做 P 单位工作。）

\textbf{跨度定律（span law）：}\ T_P\ge T_\infty\qquad(27.3)

（无限处理器的机器可以模拟 P 处理器的机器。）

加速比（speedup） $T_1/T_P\le P$。$T_1/T_P=\Theta(P)$ 称线性加速，$=P$ 称完美线性加速。
并行度（parallelism） $T_1/T_\infty$。三种理解：关键路径上每一步平均可并行执行的工作量；任意多处理器下可能达到的最大加速比；能否达到完美线性加速的界限——一旦 $P>T_1/T_\infty$，由跨度定律 $T_1/T_P\le T_1/T_\infty<P$；若 $P\gg T_1/T_\infty$，加速比远小于 P。例：P-FIB(4) 并行度 17/8 = 2.125，再多处理器也难以超过约 2 倍加速。
（并行）松弛度（slackness） $(T_1/T_\infty)/P=T_1/(PT_\infty)$：并行度超出处理器数的倍数。小于 1 不可能完美线性加速；大于 1 时每处理器工作量成为主要约束，松弛度越大，好的调度器越接近完美线性加速。

调度（scheduling）：编程模型不指定哪个链在哪个处理器上执行，由平台调度器把动态展开的计算映射到处理器（实际上映射到静态线程，再由操作系统调度，此处可忽略这一层）。调度器必须在线（on-line）工作（事先不知道何时派生、何时完成），好的调度器还是分布式的。可证明良好的在线分布式调度器存在，但分析复杂；本节分析在线集中式的贪心调度器（greedy scheduler）：每个时间步分配尽可能多的链。若就绪链 ≥ P，称完全步（complete step），任选 P 个执行；否则为不完全步（incomplete step），所有就绪链都执行。

定理 27.1：P 处理器理想并行计算机上，贪心调度器执行工作量 $T_1$、跨度 $T_\infty$ 的计算，时间

T_P\le T_1/P+T_\infty.\qquad(27.4)

证明：完全步每步做 P 单位工作，若完全步多于 $\lfloor T_1/P\rfloor$ 个，则总工作 $\ge P(\lfloor T_1/P\rfloor+1)=T_1-(T_1\bmod P)+P>T_1$，矛盾。不完全步（设每链单位时间，长链可拆成单位链）：最长路径必从入度 0 的顶点开始，不完全步执行了未执行子图 $G'$ 中全部入度 0 的链，故剩余子图的最长路径长度减 1；因此不完全步至多 $T_\infty$ 个。

推论 27.2：贪心调度在最优的 2 倍以内：$T_P^*\ge\max(T_1/P,T_\infty)$，故 $T_P\le T_1/P+T_\infty\le2\max(T_1/P,T_\infty)\le2T_P^*$。

推论 27.3：若 $P\ll T_1/T_\infty$，则 $T_P\approx T_1/P$，加速比约为 P。经验法则：松弛度至少 10 通常就够了——此时跨度项不到每处理器工作量项的 10%。只有 10 或 100 个处理器时，并行度 1,000,000 与 10,000 没什么区别；有时牺牲极端的并行度可以换来其他方面更好的算法（思考题 27-2）。

分析多线程算法：工作量就是串行化的运行时间。跨度的组合规则（图 27.3）：两个子计算串联时工作量相加、跨度相加；并联时工作量相加、跨度取最大值。P-FIB：$T_1(n)=\Theta(\phi^n)$；跨度 $T_\infty(n)=\max(T_\infty(n-1),T_\infty(n-2))+\Theta(1)=T_\infty(n-1)+\Theta(1)=\Theta(n)$；并行度 $\Theta(\phi^n/n)$，随 n 急剧增长，中等 n 即可在最大的并行机上接近完美线性加速。

并行循环（parallel for）：例如矩阵-向量乘法 $y_i=\sum_j a_{ij}x_j$：

MAT-VEC(A, x):
    n = A.rows; y = new vector(n)
    parallel for i = 1 to n: y[i] = 0
    parallel for i = 1 to n:
        for j = 1 to n:
            y[i] = y[i] + a[i][j] * x[j]
    return y

编译器用嵌套并行把 parallel for 实现为分治子程序：

MAT-VEC-MAIN-LOOP(A, x, y, n, i, i'):
    if i == i':
        for j = 1 to n: y[i] += a[i][j] * x[j]
    else:
        mid = floor((i + i') / 2)
        spawn MAT-VEC-MAIN-LOOP(A, x, y, n, i, mid)
        MAT-VEC-MAIN-LOOP(A, x, y, n, mid+1, i')
        sync

形成以单次迭代为叶的二叉树（图 27.4）。

工作量：串行化为 $\Theta(n^2)$。递归派生的开销：满二叉树内部结点比叶子少 1，每个内部结点常数工作，可摊到迭代上，只增加常数因子。实际平台常把若干次迭代合并到一个叶子（粗化（coarsening））以减少开销，代价是降低并行度，但只要松弛度足够就不影响近乎完美的加速。
跨度：n 次迭代、第 i 次迭代跨度为 $iter_\infty(i)$ 的并行循环， $T_\infty(n)=\Theta(\lg n)+\max_{1\le i\le n}iter_\infty(i).$ MAT-VEC 初始化循环跨度 $\Theta(\lg n)$，主循环每次外层迭代含 n 次串行内层迭代，跨度 $\Theta(n)$；总跨度 $\Theta(n)$，并行度 $\Theta(n^2)/\Theta(n)=\Theta(n)$。（习题 27.1-6 要求并行度 $\Theta(n^2/\lg n)$。）

竞争条件（race conditions）：多线程算法若对同一输入无论如何调度都行为相同，则是确定性的（deterministic），否则非确定性的。本应确定的算法常因**确定性竞争（determinacy race）**而失败：两条逻辑上并行的指令访问同一内存位置，且至少一条是写。著名竞争 bug：Therac-25 放疗机（3 人死亡多人受伤）、2003 年北美大停电（逾 5000 万人断电）；这类 bug 极难发现，实验室测试几天不出错，现场却偶发崩溃。

RACE-EXAMPLE():
    x = 0
    parallel for i = 1 to 2:
        x = x + 1
    print x

串行化总打印 2，但并行时可能打印 1：x = x+1 不是原子操作，分为“读入寄存器 → 寄存器加 1 → 写回”。若处理器 1 读 x=0 并加到 r1=1，处理器 2 也读 x=0、加 1、写回 x=1，处理器 1 再写回 1，就丢失了一次更新（图 27.5，执行序 1,2,3,4,5,6,7,8）。顺序 ⟨1,2,3,7,4,5,6,8⟩ 或 ⟨1,4,5,6,2,3,7,8⟩ 则正确。多数交错是对的，只有少数交错出错，所以很难测出。处理办法有互斥锁等同步手段；本章的约定是让并行的链相互独立：parallel for 的各次迭代独立；spawn 与对应 sync 之间，子任务代码与父任务代码（包括其他子任务）独立。注意派生子任务的参数在父任务中、派生之前求值，因此参数求值与派生后对参数的访问是串行的。

反例 MAT-VEC-WRONG：把内层循环也改为 parallel for 以得到 $\Theta(\lg n)$ 跨度，但所有 j 并发更新 $y_i$，产生竞争，结果错误。有竞争的代码有时也正确（例如两个线程写入同一个值），但本书一般视有竞争的代码为非法。

国际象棋程序的教训（★Socrates）：原型在 32 处理器机器上开发，最终跑在 512 处理器超级计算机上。某项“优化”使 32 处理器上的基准从 $T_{32}=65$ 秒降到 $T'_{32}=40$ 秒。但原版 $T_1=2048,T_\infty=1$，优化版 $T'_1=1024,T'_\infty=8$。用近似 $T_P\approx T_1/P+T_\infty$：$T_{32}=2048/32+1=65$，$T'_{32}=1024/32+8=40$；而 $T_{512}=2048/512+1=5$，$T'_{512}=1024/512+8=10$——在 512 处理器上“优化”版反而慢一倍，因为跨度 8 在 512 处理器时成为主导项。开发者据此放弃了该优化。教训：用工作量和跨度外推性能，比用实测时间更可靠。

习题 27.1：

27.1-1 把第 4 行 P-FIB(n−2) 也改为 spawn，对渐近工作量、跨度、并行度的影响（无影响）。
27.1-2 画 P-FIB(5) 的 DAG，求工作量、跨度、并行度，给出 3 处理器的贪心调度。
27.1-3 证明更强的界 $T_P\le(T_1-T_\infty)/P+T_\infty$（27.5）。
27.1-4 构造 DAG，使同一处理器数下两次贪心调度的时间相差近 2 倍。
27.1-5 Karan 教授声称 $T_4=80,T_{10}=42,T_{64}=10$ 秒，用工作量定律、跨度定律与 (27.5) 证明不可能。
27.1-6 矩阵-向量乘法达到 $\Theta(n^2/\lg n)$ 并行度且工作量 $\Theta(n^2)$（内层用并行归约）。
27.1-7 原地转置 P-TRANSPOSE（两层 parallel for 交换 $a_{ij},a_{ji}$）的工作量、跨度、并行度（$\Theta(n^2)$、$\Theta(\lg n)$）。
27.1-8 内层改为普通 for 后的分析（跨度 $\Theta(n)$）。
27.1-9 两版象棋程序在多少处理器时同样快（$2048/P+1=1024/P+8$，P≈146）。

27.2 多线程矩阵乘法（Multithreaded matrix multiplication）（PDF p.813–818）

并行化三重循环：

P-SQUARE-MATRIX-MULTIPLY(A, B):
    n = A.rows; C = new n×n matrix
    parallel for i = 1 to n:
        parallel for j = 1 to n:
            c[i][j] = 0
            for k = 1 to n:
                c[i][j] += a[i][k] * b[k][j]
    return C

工作量 $T_1=\Theta(n^3)$（串行化即 SQUARE-MATRIX-MULTIPLY）；跨度 $\Theta(\lg n)+\Theta(\lg n)+\Theta(n)=\Theta(n)$；并行度 $\Theta(n^2)$。习题 27.2-3 要求把内层也并行化达到 $\Theta(n^3/\lg n)$——不能直接用 parallel for，否则对 $c_{ij}$ 有竞争。

分治多线程矩阵乘法：把 A、B、C 各分成四个 $n/2\times n/2$ 子矩阵，

\begin{pmatrix}C_{11}&C_{12}\\C_{21}&C_{22}\end{pmatrix}=\begin{pmatrix}A_{11}B_{11}&A_{11}B_{12}\\A_{21}B_{11}&A_{21}B_{12}\end{pmatrix}+\begin{pmatrix}A_{12}B_{21}&A_{12}B_{22}\\A_{22}B_{21}&A_{22}B_{22}\end{pmatrix}\qquad(27.6)

8 次 $n/2$ 阶乘法 + 1 次 n 阶加法。

P-MATRIX-MULTIPLY-RECURSIVE(C, A, B):      # 输出矩阵作为参数，避免多余分配
    n = A.rows
    if n == 1: c11 = a11 * b11
    else:
        T = new n×n matrix
        分块 A, B, C, T
        spawn P-MMR(C11, A11, B11); spawn P-MMR(C12, A11, B12)
        spawn P-MMR(C21, A21, B11); spawn P-MMR(C22, A21, B12)
        spawn P-MMR(T11, A12, B21); spawn P-MMR(T12, A12, B22)
        spawn P-MMR(T21, A22, B21);       P-MMR(T22, A22, B22)   # 最后一个在主链中执行
        sync
        parallel for i = 1 to n:
            parallel for j = 1 to n:
                c[i][j] += t[i][j]

工作量：$M_1(n)=8M_1(n/2)+\Theta(n^2)=\Theta(n^3)$（主定理情形 1），与三重循环相同。
跨度：分块 $\Theta(1)$，被加法循环的 $\Theta(\lg n)$ 主导；8 个并行递归规模相同，取其一： $M_\infty(n)=M_\infty(n/2)+\Theta(\lg n)\qquad(27.7)$ 不属于主定理任何情形，但满足习题 4.6-2 的条件，解为 $\Theta(\lg^2n)$。
并行度 $\Theta(n^3/\lg^2n)$，非常高。缺点是需要临时矩阵 T（思考题 27-2 去掉 T，跨度变为 $\Theta(n)$）。

多线程 Strassen 方法：

按 (27.6) 分块，下标计算 $\Theta(1)$ 工作量和跨度。
构造 10 个 $n/2\times n/2$ 矩阵 $S_1..S_{10}$（两块之和或差），双重 parallel for，$\Theta(n^2)$ 工作量、$\Theta(\lg n)$ 跨度。
递归派生计算 7 个乘积 $P_1..P_7$。
由 $P_i$ 加减组合得到 $C_{11},C_{12},C_{21},C_{22}$，同样 $\Theta(n^2)$ 工作量、$\Theta(\lg n)$ 跨度。工作量 $\Theta(n^{\lg7})$（串行化即原算法）；跨度同样满足 (27.7)，为 $\Theta(\lg^2n)$；并行度 $\Theta(n^{\lg7}/\lg^2n)$，略低于 P-MATRIX-MULTIPLY-RECURSIVE。

习题 27.2：

27.2-1、27.2-2 画 2×2 情形的计算 DAG 并分析。
27.2-3 工作量 $\Theta(n^3)$、跨度 $\Theta(\lg n)$ 的矩阵乘法。
27.2-4 $p\times q$ 乘 $q\times r$，即使某维为 1 也高度并行。
27.2-5 分治原地转置（不用 parallel for）。
27.2-6 多线程 Floyd-Warshall（k 循环串行，i、j 两层并行，工作量 $\Theta(n^3)$、跨度 $\Theta(n\lg n)$）。

27.3 多线程归并排序（Multithreaded merge sort）（PDF p.818–826）

朴素版本：

MERGE-SORT'(A, p, r):
    if p < r:
        q = floor((p + r) / 2)
        spawn MERGE-SORT'(A, p, q)
        MERGE-SORT'(A, q+1, r)
        sync
        MERGE(A, p, q, r)          # 串行合并，Θ(n) 工作量与跨度

工作量 $MS'_1(n)=2MS'_1(n/2)+\Theta(n)=\Theta(n\lg n)$；跨度 $MS'_\infty(n)=MS'_\infty(n/2)+\Theta(n)=\Theta(n)$；并行度只有 $\Theta(\lg n)$——排 1000 万个元素也许在几个处理器上能线性加速，但无法扩展到几百个处理器。瓶颈是串行的 MERGE。

多线程合并思路（图 27.6）：把 $T[p_1..r_1]$（长 $n_1$）与 $T[p_2..r_2]$（长 $n_2$）合并到 $A[p_3..r_3]$（$n_3=n_1+n_2$），不妨设 $n_1\ge n_2$。

取较长子数组的中位元素 $x=T[q_1]$，$q_1=\lfloor(p_1+r_1)/2\rfloor$。
二分查找 $q_2$，使 x 插入 $T[q_2-1]$ 与 $T[q_2]$ 之间后第二个子数组仍有序。
$q_3=p_3+(q_1-p_1)+(q_2-p_2)$（x 之前的元素个数），$A[q_3]=x$。
递归并行地合并 $T[p_1..q_1-1]$ 与 $T[p_2..q_2-1]$ 到 $A[p_3..q_3-1]$，以及 $T[q_1+1..r_1]$ 与 $T[q_2..r_2]$ 到 $A[q_3+1..r_3]$。基本情形 $n_1=n_2=0$；由于 $n_1\ge n_2$，只需检查 $n_1=0$。

def binary_search(x, T, p, r):
    # 返回 p..r+1 中第一个满足 x <= T[q] 的位置；空区间返回 p
    low, high = p, max(p, r + 1)
    while low < high:
        mid = (low + high) // 2
        if x <= T[mid]:
            high = mid
        else:
            low = mid + 1
    return high                       # Θ(lg n) 工作量与跨度

P-MERGE(T, p1, r1, p2, r2, A, p3):
    n1 = r1 - p1 + 1; n2 = r2 - p2 + 1
    if n1 < n2: 交换 (p1,r1,n1) 与 (p2,r2,n2)     # 保证 n1 >= n2
    if n1 == 0: return                            # 两个都为空
    q1 = floor((p1 + r1) / 2)
    q2 = BINARY-SEARCH(T[q1], T, p2, r2)
    q3 = p3 + (q1 - p1) + (q2 - p2)
    A[q3] = T[q1]
    spawn P-MERGE(T, p1, q1-1, p2, q2-1, A, p3)
    P-MERGE(T, q1+1, r1, q2, r2, A, q3+1)
    sync

与 MERGE 不同：两个子数组不必相邻，输出写到另一个数组 A；$r_3=p_3+(r_1-p_1)+(r_2-p_2)+1$ 不作为输入。第二个子数组为空时也能正确处理：每次递归把较长子数组的中位元素放入输出，直到它也为空。

P-MERGE 分析（$n=n_1+n_2$）：

跨度：最坏情况下任一递归调用至多涉及 $3n/4$ 个元素。因 $n_2\le n_1$，$n_2\le n/2$；最坏时一个递归合并 $\lfloor n_1/2\rfloor$ 个与全部 $n_2$ 个：$\lfloor n_1/2\rfloor+n_2\le n_1/2+n_2/2+n_2/2=n/2+n_2/2\le3n/4$。加上二分查找： $PM_\infty(n)=PM_\infty(3n/4)+\Theta(\lg n)\qquad(27.8)$ 由习题 4.6-2，$PM_\infty(n)=\Theta(\lg^2n)$。
工作量：每个元素都要复制，$\Omega(n)$；上界：两个递归共处理至多 $n-1$ 个元素，单个至多 $3n/4$： $PM_1(n)=PM_1(\alpha n)+PM_1((1-\alpha)n)+O(\lg n),\quad 1/4\le\alpha\le3/4\qquad(27.9)$ （α 每层可变）。代入法设 $PM_1(n)\le c_1n-c_2\lg n$，得 $\le c_1n-c_2\lg n-(c_2(\lg n+\lg(\alpha(1-\alpha)))-\Theta(\lg n))\le c_1n-c_2\lg n$（$c_2$ 足够大）。故 $PM_1(n)=\Theta(n)$。
并行度 $\Theta(n/\lg^2n)$。

多线程归并排序：

P-MERGE-SORT(A, p, r, B, s):         # 把 A[p..r] 排序后写入 B[s..s+r-p]
    n = r - p + 1
    if n == 1: B[s] = A[p]
    else:
        T = new array[1..n]
        q = floor((p + r) / 2); q' = q - p + 1
        spawn P-MERGE-SORT(A, p, q, T, 1)
        P-MERGE-SORT(A, q+1, r, T, q'+1)
        sync
        P-MERGE(T, 1, q', q'+1, n, B, s)

工作量：$PMS_1(n)=2PMS_1(n/2)+\Theta(n)=\Theta(n\lg n)$（与串行相同）。
跨度：$PMS_\infty(n)=PMS_\infty(n/2)+\Theta(\lg^2n)$（27.10），由习题 4.6-2 得 $\Theta(\lg^3n)$。
并行度：$\Theta(n\lg n)/\Theta(\lg^3n)=\Theta(n/\lg^2n)$，远好于 MERGE-SORT' 的 $\Theta(\lg n)$。
空间：每层递归分配临时数组，总计 $O(n\lg n)$（可优化）。实践中应粗化基本情形，规模足够小时改用串行排序（如快速排序），牺牲一点并行度以降低常数。

习题 27.3：

27.3-1 粗化 P-MERGE 的基本情形。
27.3-2 改为求两个有序子数组全体的中位数（用习题 9.3-8）来划分，给出并分析。
27.3-3 围绕主元并行划分数组（可用辅助数组、多遍扫描，借助并行前缀和）。
27.3-4 多线程 RECURSIVE-FFT。
27.3-5★ 多线程 RANDOMIZED-SELECT（用 27.3-3 的划分）。
27.3-6★ 多线程 SELECT（9.3 节）。

第 27 章思考题与注记（PDF p.826–833）

27-1 用嵌套并行实现并行循环：SUM-ARRAYS 对 $C[i]=A[i]+B[i]$ 做 parallel for。(a) 仿照 MAT-VEC-MAIN-LOOP 用 spawn/sync 改写并分析并行度；另一种实现 SUM-ARRAYS' 按 grain-size（粒度） 把数组切成 $r=\lceil n/\text{grain-size}\rceil$ 段，串行 for 循环逐段 spawn ADD-SUBARRAY；(b) grain-size = 1 时的并行度；(c) 用 n 和 grain-size 表示跨度，求使并行度最大的粒度（跨度 $\Theta(n/g+g)$，取 $g=\sqrt n$，并行度 $\Theta(\sqrt n)$）。
27-2 矩阵乘法中节省临时空间：P-MATRIX-MULTIPLY-RECURSIVE 并行度极高，1000×1000 时约 $1000^3/10^2=10^7$，远超实际处理器数，而临时矩阵 T 影响常数。(a) 计算 $C=C+AB$，并行初始化 C，在恰当位置插入 sync（先并行做 4 个乘积累加到 C，sync，再做另外 4 个），去掉 T，跨度变为 $\Theta(n)$；(b) 写出并解工作量与跨度递归式（$\Theta(n^3)$、$\Theta(n)$）；(c) 1000×1000 时并行度约 $10^6$，仍然足够。
27-3 多线程矩阵算法：多线程化 LU-DECOMPOSITION、LUP-DECOMPOSITION、LUP-SOLVE，以及基于式 (28.13) 的对称正定矩阵求逆，分析工作量、跨度、并行度。
27-4 多线程归约与前缀计算：⊗ 为结合运算，⊗-归约（reduction） $y=x[1]\otimes\cdots\otimes x[n]$。(a) P-REDUCE：$\Theta(n)$ 工作量、$\Theta(\lg n)$ 跨度（分治）。⊗-前缀计算（prefix computation / scan）：$y[i]=x[1]\otimes\cdots\otimes x[i]$；串行 SCAN 的循环有依赖，直接 parallel for 会产生竞争。(b) P-SCAN-1：每个位置独立调用 P-REDUCE，工作量 $\Theta(n^2)$、跨度 $\Theta(\lg n)$；(c) P-SCAN-2：递归求左右两半前缀，再用 parallel for 把左半最后值并入右半，工作量 $\Theta(n\lg n)$、跨度 $\Theta(\lg^2n)$；(d)(e) P-SCAN-3 两遍法：上行（P-SCAN-UP，在 t 中记录各子区间的归约，返回 t[k] ⊗ right）与下行（P-SCAN-DOWN，左半传 v，右半传 v ⊗ t[k]），不变式为传入的 $v=x[1]\otimes\cdots\otimes x[i-1]$；工作量 $\Theta(n)$、跨度 $\Theta(\lg n)$、并行度 $\Theta(n/\lg n)$。
27-5 简单模板（stencil）计算的多线程化：$A[i,j]$ 只依赖于左上方已算的元素（如 LCS），填一个元素 $\Theta(1)$。(a) 2×2 分块：先 $A_{11}$，再并行 $A_{12}$、$A_{21}$，最后 $A_{22}$；工作量 $\Theta(n^2)$，跨度 $S(n)=3S(n/2)+\Theta(1)=\Theta(n^{\lg3})$，并行度 $\Theta(n^{2-\lg3})\approx n^{0.415}$。(b) 3×3 分块（按反对角线波前并行），跨度 $5S(n/3)$，即 $\Theta(n^{\log_35})$。(c) 一般 b×b 分块：跨度 $\Theta(n^{\log_b(2b-1)})$，并行度 $o(n)$。(d) 设计并行度 $\Theta(n/\lg n)$ 的算法，并论证问题内在并行度为 $\Theta(n)$（按反对角线逐条并行）。
27-6 随机化多线程算法：(a) 工作量定律、跨度定律、贪心界改为期望形式；(b) 例：1% 情况 $T_1=10^4,T_{10000}=1$，99% 情况 $T_1=T_{10000}=10^9$，说明加速比应定义为 $E[T_1]/E[T_P]$ 而非 $E[T_1/T_P]$；(c) 并行度定义为 $E[T_1]/E[T_\infty]$；(d)(e) 用嵌套并行多线程化 RANDOMIZED-QUICKSORT（不并行化划分）并分析（期望跨度 $\Theta(n)$，并行度 $\Theta(\lg n)$）。

注记：前几版讲过排序网络与 PRAM 模型；**数据并行（data-parallel）模型以向量、矩阵运算为原语。Graham 与 Brent 证明存在达到定理 27.1 界的调度器；Eager–Zahorjan–Lazowska 证明任何贪心调度器都达到该界，并提出用工作量与跨度分析并行算法；Blelloch 基于工作量与“深度（depth）”建立数据并行编程模型。Blumofe–Leiserson 提出基于随机工作窃取（work-stealing）**的分布式调度，达到 $E[T_P]\le T_1/P+O(T_\infty)$。伪代码与模型深受 MIT Cilk 项目和 Cilk++ 影响；多线程归并排序受 Akl 启发；顺序一致性概念源自 Lamport。

第 27 章小结

本章要点 动态多线程只用 spawn、sync、parallel 三个关键字表达逻辑并行，删掉它们即串行算法。计算被建模为链组成的 DAG：工作量 $T_1$ 是总时间，跨度 $T_\infty$ 是关键路径长度，并行度 $T_1/T_\infty$ 是可能的最大加速比。工作量定律 $T_P\ge T_1/P$、跨度定律 $T_P\ge T_\infty$ 给出下界，贪心调度达到 $T_P\le T_1/P+T_\infty$，在最优的 2 倍以内；松弛度约 10 以上即接近线性加速。分析时串联跨度相加、并联取最大；并行循环的跨度为 $\Theta(\lg n)$ 加单次迭代最大跨度。必须避免确定性竞争。代表算法：矩阵乘法（三重循环并行度 $\Theta(n^2)$，分治 $\Theta(n^3/\lg^2n)$，Strassen $\Theta(n^{\lg7}/\lg^2n)$），归并排序（串行合并时只有 $\Theta(\lg n)$，用二分 + 分治的 P-MERGE 后为 $\Theta(n/\lg^2n)$）。

与量化交易的关联

回测与参数扫描的并行化：大量回测、蒙特卡洛模拟、参数网格搜索天然是 parallel for，各次迭代独立；用工作量/跨度估算能否在给定核数上线性加速，避免像国际象棋程序那样被跨度（例如串行的数据加载、最终汇总）卡住。
竞争条件：多线程撮合引擎、风控计数器、持仓与资金的并发更新是竞争 bug 的高发区，本章的 RACE-EXAMPLE 就是“两个线程同时加仓导致持仓记错”的原型。实盘系统应使用原子操作、锁或单线程事件循环（很多交易系统正是因此采用单写者设计）。
并行前缀和（思考题 27-4）：累计收益、累计成交量、VWAP、滚动统计等都是前缀计算，GPU/向量化库（如 cumsum、CUDA scan）正是用两遍法实现的；理解它有助于在大规模 tick 数据上高效计算。
并行矩阵运算：协方差矩阵估计、因子暴露矩阵乘法、组合优化中的大矩阵运算都依赖并行 BLAS；本章解释了分块矩阵乘法为何并行度高。
任务调度 DAG：日终批处理（行情清洗 → 因子 → 模型 → 优化）可看作计算 DAG，跨度即关键路径，决定加多少机器还有没有用。
性能外推：用 $T_P\approx T_1/P+T_\infty$ 估算扩容收益，比直接测速更可靠。

推荐习题

27.1-5（用三条定律识破不可能的测量结果）；27.1-9（象棋程序交叉点）；27.1-6（并行归约降低跨度）。
27.2-3（$\Theta(\lg n)$ 跨度矩阵乘法）；27.2-6（并行 Floyd-Warshall）。
27.3-3（并行划分，需前缀和）。
思考题 27-1（粒度选择，工程中最实用）；27-4（并行前缀和，强烈推荐）；27-5（模板计算/动态规划的波前并行）。

第 28 章矩阵运算（Matrix Operations）

第 28 章引言（PDF p.834）

矩阵运算是科学计算的核心。本章讨论矩阵乘法与解线性方程组（矩阵基础见附录 D）。28.1 用 LUP 分解解线性方程组；28.2 讨论矩阵乘法与求逆的紧密关系；28.3 讨论对称正定矩阵，并用它求超定线性方程组的最小二乘解。数值稳定性（numerical stability）：浮点精度有限，舍入误差可能在计算中被放大导致错误结果，称为数值不稳定。本章偶尔提及但不重点讨论，详见 Golub & Van Loan。

28.1 求解线性方程组（Solving systems of linear equations）（PDF p.834–848）

问题：n 个方程 n 个未知数

\sum_{j=1}^na_{ij}x_j=b_i\ (i=1..n)\quad\Longleftrightarrow\quad Ax=b.\qquad(28.1)(28.2)

若 A 非奇异，则 $x=A^{-1}b$（28.3）是唯一解（若 $Ax=Ax'=b$，则 $x=A^{-1}Ax=A^{-1}Ax'=x'$）。本节只讨论 A 非奇异（等价地，秩为 n）。若方程数少于未知数或秩 < n，为**欠定（underdetermined）系统，通常无穷多解，也可能因矛盾而无解；方程数多于未知数为超定（overdetermined）**系统，可能无解（28.3 节求近似解）。

先求 $A^{-1}$ 再乘 b 的方法数值不稳定；LUP 分解既数值稳定，实践中也更快。

LUP 分解概述：找 $n\times n$ 矩阵 L、U、P 使

PA=LU,\qquad(28.4)

其中 L 是单位下三角矩阵（unit lower-triangular）（对角线全为 1），U 是上三角矩阵，P 是置换矩阵（permutation matrix）。每个非奇异矩阵都有 LUP 分解。有了它，$Ax=b$ 两边左乘 P 得 $PAx=Pb$（即重排方程），即 $LUx=Pb$。令 $y=Ux$：

**前代（forward substitution）**解下三角系统 $Ly=Pb$（28.5）；
**回代（back substitution）**解上三角系统 $Ux=y$（28.6）。验证：$A=P^{-1}LU$（28.7），$Ax=P^{-1}LUx=P^{-1}Ly=P^{-1}Pb=b$。

前代与回代：用数组 $\pi[1..n]$ 紧凑表示 P：$P_{i,\pi[i]}=1$，于是 PA 第 i 行第 j 列为 $a_{\pi[i],j}$，Pb 第 i 个元素为 $b_{\pi[i]}$。由于 L 单位下三角：

y_i=b_{\pi[i]}-\sum_{j=1}^{i-1}l_{ij}y_j;\qquad x_i=\Big(y_i-\sum_{j=i+1}^nu_{ij}x_j\Big)\Big/u_{ii}.

def lup_solve(L, U, pi, b):          # Θ(n^2) 时间，O(n) 额外空间
    n = len(L)
    y = [0.0] * n; x = [0.0] * n
    for i in range(n):                                    # 前代
        y[i] = b[pi[i]] - sum(L[i][j] * y[j] for j in range(i))
    for i in reversed(range(n)):                          # 回代
        x[i] = (y[i] - sum(U[i][j] * x[j] for j in range(i + 1, n))) / U[i][i]
    return x

例：

A=\begin{pmatrix}1&2&0\\3&4&4\\5&6&3\end{pmatrix},\ b=\begin{pmatrix}3\\7\\8\end{pmatrix};\quad L=\begin{pmatrix}1&0&0\\0.2&1&0\\0.6&0.5&1\end{pmatrix},\ U=\begin{pmatrix}5&6&3\\0&0.8&-0.6\\0&0&2.5\end{pmatrix},\ P=\begin{pmatrix}0&0&1\\1&0&0\\0&1&0\end{pmatrix}.

前代解 $Ly=Pb=(8,3,7)^T$ 得 $y=(8,1.4,1.5)^T$；回代解 $Ux=y$ 得 $x=(-1.4,2.2,0.6)^T$。

LU 分解（无置换，$P=I_n$）——高斯消元（Gaussian elimination）：从其余方程中减去第一个方程的倍数以消去第一个变量，再用第二个方程消去后续方程中的第二个变量……直到剩下上三角形式，即 U；L 由消元所用的行乘数组成。递归表述：$n=1$ 时 $L=I_1,U=A$；$n>1$ 时分块

A=\begin{pmatrix}a_{11}&w^T\\v&A'\end{pmatrix}=\begin{pmatrix}1&0\\v/a_{11}&I_{n-1}\end{pmatrix}\begin{pmatrix}a_{11}&w^T\\0&A'-vw^T/a_{11}\end{pmatrix},\qquad(28.8)

其中 $v=(a_{21},\dots,a_{n1})^T$ 为列向量，$w^T=(a_{12},\dots,a_{1n})$ 为行向量，$vw^T/a_{11}$ 是外积除以 $a_{11}$。

A'-vw^T/a_{11}\qquad(28.9)

称为 A 关于 $a_{11}$ 的舒尔补（Schur complement）。

断言：A 非奇异则舒尔补非奇异。否则舒尔补行秩 < n−1，(28.8) 右侧第二个矩阵下面 n−1 行（第一列全为 0）行秩 < n−1，整个矩阵行秩 < n，进而 A 秩 < n，与非奇异矛盾。
递归求 $A'-vw^T/a_{11}=L'U'$，则 $A=\begin{pmatrix}1&0\\v/a_{11}&L'\end{pmatrix}\begin{pmatrix}a_{11}&w^T\\0&U'\end{pmatrix}=LU.$
除数称为主元（pivots），位于 U 的对角线上。若 $a_{11}=0$ 或某步舒尔补左上角为 0，此法因除以 0 而失败。引入 P 正是为了避免除以 0（或除以很小的数而引起数值不稳定），这叫选主元（pivoting）。对称正定矩阵无需选主元，LU 分解总能进行（28.3 节证明）。

def lu_decomposition(A):              # Θ(n^3) 时间；尾递归改写为循环
    n = len(A)
    L = [[1.0 if i == j else 0.0 for j in range(n)] for i in range(n)]
    U = [[0.0] * n for _ in range(n)]
    for k in range(n):
        U[k][k] = A[k][k]                       # 主元
        for i in range(k + 1, n):
            L[i][k] = A[i][k] / A[k][k]         # v_i / 主元
            U[k][i] = A[k][i]                   # w_i
        for i in range(k + 1, n):               # 舒尔补写回 A
            for j in range(k + 1, n):
                A[i][j] -= L[i][k] * U[k][j]
    return L, U

标准优化：把 L（$i>j$ 部分）和 U（$i\le j$ 部分）原地存放在 A 中（把 l、u 的引用都换成 a 即可）。

例（图 28.1）：

\begin{pmatrix}2&3&1&5\\6&13&5&19\\2&19&10&23\\4&10&11&31\end{pmatrix}=\begin{pmatrix}1&0&0&0\\3&1&0&0\\1&4&1&0\\2&1&7&1\end{pmatrix}\begin{pmatrix}2&3&1&5\\0&4&2&4\\0&0&1&2\\0&0&0&3\end{pmatrix}.

主元依次为 2、4、1、3。

LUP 分解：即使 A 非奇异，也要避免除以很小的值，所以选绝对值最大的元素作主元（部分选主元，partial pivoting）。第一列不可能全为 0（否则行列式为 0，A 奇异）。把第一列绝对值最大的 $a_{k1}$ 所在行与第 1 行交换，相当于左乘置换矩阵 Q：

QA=\begin{pmatrix}a_{k1}&w^T\\v&A'\end{pmatrix}=\begin{pmatrix}1&0\\v/a_{k1}&I_{n-1}\end{pmatrix}\begin{pmatrix}a_{k1}&w^T\\0&A'-vw^T/a_{k1}\end{pmatrix}

（v 中 $a_{11}$ 代替了 $a_{k1}$，$w^T=(a_{k2},\dots,a_{kn})$）。舒尔补非奇异，递归得 $P'(A'-vw^T/a_{k1})=L'U'$。令 $P=\begin{pmatrix}1&0\\0&P'\end{pmatrix}Q$（置换矩阵之积仍为置换矩阵），则

PA=\begin{pmatrix}1&0\\P'v/a_{k1}&L'\end{pmatrix}\begin{pmatrix}a_{k1}&w^T\\0&U'\end{pmatrix}=LU.

与 LU 不同，列向量 $v/a_{k1}$ 也要乘 $P'$——所以实现时整行交换。

def lup_decomposition(A):             # Θ(n^3)，原地；返回置换数组 pi
    n = len(A)
    pi = list(range(n))
    for k in range(n):
        p, kp = 0.0, None
        for i in range(k, n):                    # 找第 k 列绝对值最大的元素
            if abs(A[i][k]) > p:
                p, kp = abs(A[i][k]), i
        if p == 0:
            raise ValueError("singular matrix")
        pi[k], pi[kp] = pi[kp], pi[k]
        A[k], A[kp] = A[kp], A[k]                # 交换整行
        for i in range(k + 1, n):
            A[i][k] /= A[k][k]                   # L 的元素
            for j in range(k + 1, n):
                A[i][j] -= A[i][k] * A[k][j]     # 舒尔补
    return pi                                    # 结束时 a_ij = l_ij (i>j)，u_ij (i<=j)

时间 $\Theta(n^3)$，与 LU 分解相同——选主元最多只多花常数因子。

例（图 28.2）：

\begin{pmatrix}0&0&1&0\\1&0&0&0\\0&0&0&1\\0&1&0&0\end{pmatrix}\begin{pmatrix}2&0&2&0.6\\3&3&4&-2\\5&5&4&2\\-1&-2&3.4&-1\end{pmatrix}=\begin{pmatrix}1&0&0&0\\0.4&1&0&0\\-0.2&0.5&1&0\\0.6&0&0.4&1\end{pmatrix}\begin{pmatrix}5&5&4&2\\0&-2&0.4&-0.2\\0&0&4&-0.5\\0&0&0&-3\end{pmatrix}.

第一步第一列最大元为第 3 行的 5，交换第 1、3 行；随后各步类推，第四步无变化。

习题 28.1：

28.1-1 用前代解一个单位下三角系统。
28.1-2 求一个 3×3 矩阵的 LU 分解。
28.1-3 用 LUP 分解解一个 3×3 系统。
28.1-4 对角矩阵的 LUP 分解。
28.1-5 置换矩阵的 LUP 分解并证明唯一。
28.1-6 对所有 n ≥ 1 存在有 LU 分解的奇异矩阵。
28.1-7 LU/LUP 分解中 k = n 的那次外循环是否必要。

28.2 矩阵求逆（Inverting matrices）（PDF p.848–853）

实践中一般不用逆矩阵解方程（LUP 更稳定），但有时确实需要逆。本节用 LUP 分解求逆，并证明矩阵乘法与矩阵求逆在渐近意义上同样难——因此可用 Strassen 算法求逆（Strassen 原论文的动机正是证明线性方程组能比常规方法更快求解）。

由 LUP 分解求逆：有了 $PA=LU$，每解一个 $Ax=b$ 只需 $\Theta(n^2)$；k 个只差右端项的方程组共 $\Theta(kn^2)$。$AX=I_n$（28.10）可视为 n 个方程组 $AX_i=e_i$（$X_i$ 为 X 的第 i 列，$e_i$ 为单位向量）。每列 $\Theta(n^2)$，共 $\Theta(n^3)$；加上分解 $\Theta(n^3)$，求逆总计 $\Theta(n^3)$。

定理 28.1（乘法不比求逆难）：若能在 $I(n)$ 时间内求 $n\times n$ 矩阵的逆，$I(n)=\Omega(n^2)$ 且满足正则性条件 $I(3n)=O(I(n))$，则可在 $O(I(n))$ 时间内计算两个 $n\times n$ 矩阵之积。证明：构造 $3n\times3n$ 矩阵

D=\begin{pmatrix}I_n&A&0\\0&I_n&B\\0&0&I_n\end{pmatrix},\qquad D^{-1}=\begin{pmatrix}I_n&-A&AB\\0&I_n&-B\\0&0&I_n\end{pmatrix},

AB 就是 $D^{-1}$ 右上角的 $n\times n$ 子块。构造 D 用 $\Theta(n^2)=O(I(n))$，求逆 $O(I(3n))=O(I(n))$。$I(n)=\Theta(n^c\lg^dn)$（$c>0,d\ge0$）都满足正则性。

定理 28.2（求逆不比乘法难）：若能在 $M(n)$ 时间内乘两个 $n\times n$ 实矩阵，$M(n)=\Omega(n^2)$，且满足 $M(n+k)=O(M(n))$（$0\le k\le n$）以及 $M(n/2)\le cM(n)$（某常数 $c<1/2$），则任何实非奇异矩阵可在 $O(M(n))$ 内求逆。证明：

可设 n 是 2 的幂：$\begin{pmatrix}A&0\\0&I_k\end{pmatrix}^{-1}=\begin{pmatrix}A^{-1}&0\\0&I_k\end{pmatrix}$，补到下一个 2 的幂，第一个正则条件保证只增加常数因子。
先设 A 对称正定，分块 $A=\begin{pmatrix}B&C^T\\C&D\end{pmatrix},\quad A^{-1}=\begin{pmatrix}R&T\\U&V\end{pmatrix}.\qquad(28.11)$ 令 $S=D-CB^{-1}C^T$（28.12）为 A 关于 B 的舒尔补，则 $A^{-1}=\begin{pmatrix}B^{-1}+B^{-1}C^TS^{-1}CB^{-1}&-B^{-1}C^TS^{-1}\\-S^{-1}CB^{-1}&S^{-1}\end{pmatrix}.\qquad(28.13)$ 由引理 28.4、28.5，B 与 S 对称正定，因而可逆（引理 28.3），且逆也对称。计算步骤（全为 $n/2$ 阶）：(1) 取出 B、C、$C^T$、D；(2) 递归求 $B^{-1}$；(3) $W=CB^{-1}$，其转置 $W^T=B^{-1}C^T$；(4) $X=WC^T=CB^{-1}C^T$，$S=D-X$；(5) 递归求 $S^{-1}$，令 $V=S^{-1}$；(6) $Y=S^{-1}W=S^{-1}CB^{-1}$，$Y^T=B^{-1}C^TS^{-1}$，令 $T=-Y^T$，$U=-Y$；(7) $Z=W^TY=B^{-1}C^TS^{-1}CB^{-1}$，$R=B^{-1}+Z$。共 2 次递归求逆、4 次乘法、$O(n^2)$ 其他开销： $I(n)\le2I(n/2)+4M(n/2)+O(n^2)=2I(n/2)+\Theta(M(n))=O(M(n))$ （第二个正则条件给出 $4M(n/2)<2M(n)$，并使主定理情形 3 适用）。
一般非奇异 A：$A^TA$ 对称正定，且 $A^{-1}=(A^TA)^{-1}A^T$（因为 $((A^TA)^{-1}A^T)A=I_n$ 且逆唯一）。先算 $A^TA$，用上法求逆，再乘 $A^T$，三步各 $O(M(n))$。

推论用法：解 $Ax=b$ 也可两边乘 $A^T$ 得 $(A^TA)x=A^Tb$，对对称正定的 $A^TA$ 做不选主元的 LU 分解。理论上正确，但实践中 LUP 分解更好：运算量少一个常数因子，数值性质也更好（$A^TA$ 会使条件数平方）。

习题 28.2：

28.2-1 矩阵乘法与矩阵平方难度相同。
28.2-2 $M(n)$ 时间的乘法蕴含 $O(M(n))$ 的 LUP 分解。
28.2-3 矩阵乘法与求行列式难度相同。
28.2-4 布尔矩阵乘法 $M(n)$ 蕴含 $O(M(n)\lg n)$ 的传递闭包，反之传递闭包 $T(n)$ 蕴含 $O(T(n))$ 的布尔矩阵乘法。
28.2-5 基于定理 28.2 的求逆在模 2 整数域上是否可行（不行，“正定”在 GF(2) 中无意义，$A^TA$ 可能奇异）。
28.2-6★ 推广到复矩阵：用共轭转置 $A^*$ 与 Hermite 矩阵。

28.3 对称正定矩阵与最小二乘逼近（PDF p.853–861）

对称正定（symmetric positive-definite）矩阵：$A=A^T$ 且对所有非零 x，$x^TAx>0$。它们非奇异，可做 LU 分解而不必担心除以 0。

引理 28.3：正定矩阵非奇异。（若奇异，存在 $x\ne0$ 使 $Ax=0$，则 $x^TAx=0$。）

第 k 个顺序主子矩阵（leading submatrix） $A_k$：前 k 行与前 k 列的交。

引理 28.4：对称正定矩阵的每个顺序主子矩阵都对称正定。证明：分块 $A=\begin{pmatrix}A_k&B^T\\B&C\end{pmatrix}$（28.14），若有 $x_k\ne0$ 使 $x_k^TA_kx_k\le0$，取 $x=(x_k^T,0)^T$，则 $x^TAx=x_k^TA_kx_k\le0$，矛盾。

一般舒尔补：A 关于 $A_k$ 的舒尔补

S=C-BA_k^{-1}B^T,\qquad(28.15)

k = 1 时与 (28.9) 一致。

引理 28.5（舒尔补引理）：A 对称正定，则 S 对称正定。证明：对称性显然。把 x 分为与 $A_k$、C 对应的 y、z，“配方”得

x^TAx=(y+A_k^{-1}B^Tz)^TA_k(y+A_k^{-1}B^Tz)+z^T(C-BA_k^{-1}B^T)z.\qquad(28.16)

对任意 $z\ne0$ 取 $y=-A_k^{-1}B^Tz$，第一项为 0，得 $z^TSz=x^TAx>0$。

推论 28.6：对称正定矩阵的 LU 分解不会除以 0。更强地，每个主元都严格为正：第一个主元 $a_{11}=e_1^TAe_1>0$；LU 第一步得到关于 $A_1=(a_{11})$ 的舒尔补，由引理 28.5 归纳即得。

最小二乘逼近（least-squares approximation）：给定 m 个带测量误差的数据点 $(x_1,y_1),\dots,(x_m,y_m)$，求函数 F 使逼近误差 $\eta_i=F(x_i)-y_i$（28.17）尽量小。设 F 是基函数的线性组合 $F(x)=\sum_{j=1}^nc_jf_j(x)$，n 与 $f_j$ 由问题背景决定；常用 $f_j(x)=x^{j-1}$，即 n−1 次多项式 $F(x)=c_1+c_2x+\cdots+c_nx^{n-1}$。

取 n = m 可精确穿过每个点，但高次 F 会“拟合噪声”，对新的 x 预测效果差（即过拟合）。通常取 n 远小于 m，期望抓住数据的主要模式而不过分关注噪声。选 n 的理论超出本书范围。此时得到超定方程组。
记 $A=(a_{ij})$，$a_{ij}=f_j(x_i)$（$m\times n$），系数向量 c，则 Ac 是预测值向量，$\eta=Ac-y$ 是误差向量。最小化 $\|\eta\|^2=\|Ac-y\|^2=\sum_{i=1}^m\Big(\sum_{j=1}^na_{ij}c_j-y_i\Big)^2,$ 对每个 $c_k$ 求导令其为 0： $\frac{d\|\eta\|^2}{dc_k}=\sum_{i=1}^m2\Big(\sum_{j=1}^na_{ij}c_j-y_i\Big)a_{ik}=0\qquad(28.18)$ 即 $(Ac-y)^TA=0$，等价于 $A^T(Ac-y)=0$，得正规方程（normal equation） $A^TAc=A^Ty.\qquad(28.19)$ $A^TA$ 对称；若 A 列满秩，则 $A^TA$ 正定，可逆， $c=\big((A^TA)^{-1}A^T\big)y=A^+y,\qquad(28.20)$ $A^+=(A^TA)^{-1}A^T$ 称为 A 的伪逆（pseudoinverse），把逆矩阵推广到非方阵（对比精确解 $A^{-1}b$）。

例（图 28.3）：5 个点 $(-1,2),(1,1),(2,1),(3,0),(5,3)$，拟合二次多项式 $F(x)=c_1+c_2x+c_3x^2$。

A=\begin{pmatrix}1&-1&1\\1&1&1\\1&2&4\\1&3&9\\1&5&25\end{pmatrix},\quad A^+=\begin{pmatrix}0.500&0.300&0.200&0.100&-0.100\\-0.388&0.093&0.190&0.193&-0.088\\0.060&-0.036&-0.048&-0.036&0.060\end{pmatrix},

$c=A^+y=(1.200,-0.757,0.214)^T$，即 $F(x)=1.200-0.757x+0.214x^2$ 是最小二乘意义下最接近的二次曲线。 实际做法：先算 $A^Ty$，再对 $A^TA$ 做 LU 分解后前代、回代求 c。A 满秩时 $A^TA$ 对称正定，必非奇异。复杂度：形成 $A^TA$ 需 $O(mn^2)$，分解 $O(n^3)$。

习题 28.3：

28.3-1 对称正定矩阵对角元均为正。
28.3-2 2×2 对称正定矩阵 $\begin{pmatrix}a&b\\b&c\end{pmatrix}$ 用配方证明 $ac-b^2>0$。
28.3-3 对称正定矩阵的最大元素在对角线上。
28.3-4 每个顺序主子矩阵的行列式为正。
28.3-5 第 k 个主元等于 $\det(A_k)/\det(A_{k-1})$（$\det A_0=1$）。
28.3-6 用 $F(x)=c_1+c_2x\lg x+c_3e^x$ 对 $(1,1),(2,1),(3,3),(4,8)$ 做最小二乘拟合。
28.3-7 证明伪逆满足四个 Moore-Penrose 条件：$AA^+A=A$，$A^+AA^+=A^+$，$(AA^+)^T=AA^+$，$(A^+A)^T=A^+A$。

第 28 章思考题与注记（PDF p.861–863）

28-1 三对角线性方程组（Tridiagonal systems）：以 5×5 矩阵（对角线 1,2,2,2,2，次对角线全 −1）为例：(a) 求 LU 分解；(b) 用前代回代解 $Ax=(1,1,1,1,1)^T$；(c) 求逆；(d) 对任意对称正定三对角矩阵，用 LU 分解 $O(n)$ 解 $Ax=b$，并说明任何先求 $A^{-1}$ 的方法最坏情况下渐近更慢（逆矩阵一般是稠密的，$\Omega(n^2)$）；(e) 对任意非奇异三对角矩阵，用 LUP 分解 $O(n)$ 求解。（即 Thomas 算法。）
28-2 样条（Splines）：用三次样条（cubic spline）插值 n+1 个点 $(x_i,y_i)$，$x_0<\cdots<x_n$；曲线由 n 段三次多项式 $f_i(x)=a_i+b_ix+c_ix^2+d_ix^3$ 组成，$x_i\le x\le x_{i+1}$ 时 $f(x)=f_i(x-x_i)$，拼接点称为节点（knots），先设 $x_i=i$。连续性 $f_i(0)=y_i,f_i(1)=y_{i+1}$；一阶导连续 $f'_i(1)=f'_{i+1}(0)$。(a) 已知各节点一阶导 $D_i$ 时，用 $y_i,y_{i+1},D_i,D_{i+1}$ 表示 $a_i,b_i,c_i,d_i$，$O(n)$ 求出 4n 个系数；要求二阶导也连续，且端点 $f''(x_0)=f''(x_n)=0$（自然三次样条，natural cubic spline）。(b) 证明 $D_{i-1}+4D_i+D_{i+1}=3(y_{i+1}-y_{i-1})$（$i=1..n-1$，28.21）；(c) $2D_0+D_1=3(y_1-y_0)$（28.22），$D_{n-1}+2D_n=3(y_n-y_{n-1})$（28.23）；(d) 写成关于 $D=(D_0..D_n)$ 的矩阵方程——系数矩阵对称、三对角、严格对角占优（正定）；(e) 由 28-1 可 $O(n)$ 插值；(f) 推广到不等距节点。

注记：推荐 George–Liu、Golub–Van Loan、Press 等《Numerical Recipes》、Strang。Golub–Van Loan 指出 $\det(A)$ 不是衡量稳定性的好指标，建议用 $\|A\|_\infty\|A^{-1}\|_\infty$（$\|A\|_\infty=\max_i\sum_j|a_{ij}|$，即条件数），并讨论不求 $A^{-1}$ 而估计它的方法。高斯消元是最早的系统解线性方程组方法之一，通常归功于 Gauss（1777–1855）。Strassen 证明 $n\times n$ 矩阵可在 $O(n^{\lg7})$ 时间求逆；Winograd 首先证明乘法不比求逆难，反方向由 Aho–Hopcroft–Ullman 证明。另一重要分解是奇异值分解（SVD）：$A=Q_1\Sigma Q_2^T$，Σ 为 $m\times n$ 且只有对角线非零，$Q_1$（$m\times m$）、$Q_2$（$n\times n$）的列两两标准正交（内积为 0 且范数为 1）。

第 28 章小结

本章要点 解 $Ax=b$ 不应先求逆，而应做 LUP 分解 $PA=LU$（$\Theta(n^3)$），再用前代和回代（各 $\Theta(n^2)$）求解；分解一次可以复用到多个右端项。LU 分解本质是高斯消元，每步把问题化为对舒尔补的分解；选绝对值最大的主元（部分选主元）既避免除以 0 又提高数值稳定性。求逆可由 n 次 LUP-SOLVE 完成，总 $\Theta(n^3)$；理论上矩阵求逆与矩阵乘法同样难（定理 28.1、28.2），可借 Strassen 加速。对称正定矩阵的顺序主子矩阵和舒尔补仍对称正定，LU 分解时主元全为正，无需选主元。最小二乘拟合归结为正规方程 $A^TAc=A^Ty$，解为伪逆 $c=A^+y$。

与量化交易的关联

线性回归与因子模型：OLS 估计正是本章最小二乘，正规方程 $X^TX\beta=X^Ty$ 是截面回归（Fama-MacBeth）、时序回归求因子暴露、Barra 风格模型因子收益估计的核心计算。工程上应避免显式求 $(X^TX)^{-1}$，用 LU/Cholesky 或 QR 求解；$X^TX$ 会把条件数平方，因子高度共线时更应使用 QR 或 SVD（注记中的条件数与 SVD 正与此相关）。
协方差矩阵与对称正定性：风险模型中的协方差矩阵必须对称正定（或半正定）。舒尔补引理对应条件协方差 $\Sigma_{22}-\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}$（已知一组资产收益后另一组的条件协方差），是多元正态条件分布、对冲比率与最小方差对冲的公式基础；分块求逆公式 (28.13) 也用于增量更新协方差逆。
组合优化：均值-方差最优解 $w\propto\Sigma^{-1}\mu$ 实际应通过解线性系统 $\Sigma w=\mu$ 获得；Cholesky 分解（对称正定矩阵的 LU 特例）还用于生成相关正态随机数做蒙特卡洛模拟。
过拟合：书中“高次多项式拟合噪声”的警告就是因子挖掘和策略参数调优中的过拟合问题。
样条（28-2）：用于利率期限结构（收益率曲线）插值、波动率曲面平滑，自然三次样条只需解三对角系统，$O(n)$ 完成。
三对角求解（28-1）：有限差分法求解 Black-Scholes 偏微分方程（隐式格式、Crank-Nicolson）每一步都要解三对角系统，Thomas 算法使每步 $O(n)$。

推荐习题

28.1-3（手算 LUP 分解与求解）；28.1-7。
28.2-3（行列式与乘法等价）。
28.3-2、28.3-5（配方法与主元的行列式解释）；28.3-6（非多项式基函数的最小二乘）；28.3-7（伪逆的 Moore-Penrose 条件）。
思考题 28-1（三对角系统，PDE 定价必备）；28-2（自然三次样条，曲线插值必备）。

第 29 章线性规划（Linear Programming）（续见下一块）

第 29 章引言（PDF p.864–871）

许多问题是在有限资源和相互竞争的约束下最大化或最小化某个目标。若目标是变量的线性函数，约束是关于变量的线性等式或不等式，就是**线性规划（linear programming, LP）**问题。

政治问题示例：选区有城市、郊区、农村三类地区，登记选民分别 100,000、200,000、50,000，希望在每类地区至少赢得半数选票，即 50,000、100,000、25,000 张。四个议题：修路、枪支管制、农业补贴、专用于公共交通的汽油税。每在某议题上花 1000 美元广告，赢得（负数为失去）的选票数（千张）如下（图 29.1）：

政策	城市	郊区	农村
修路	−2	5	3
枪支管制	8	2	−5
农业补贴	0	0	10
汽油税	10	0	−2

试错方案：修路 20、枪支 0、补贴 4、汽油税 9（千美元），则城市 $20(-2)+0+4(0)+9(10)=50$，郊区 $20(5)=100$，农村 $20(3)+0+4(10)+9(-2)=82$（甚至超过农村选民总数），花费 33 千美元。是否能更省？设 $x_1..x_4$ 为四项广告支出（千美元），得线性规划：

\begin{aligned}\min\ & x_1+x_2+x_3+x_4 &(29.6)\\ \text{s.t.}\ & -2x_1+8x_2+0x_3+10x_4\ge50 &(29.7)\\ &5x_1+2x_2+0x_3+0x_4\ge100 &(29.8)\\ &3x_1-5x_2+10x_3-2x_4\ge25 &(29.9)\\ &x_1,x_2,x_3,x_4\ge0 &(29.10)\end{aligned}

（负面广告存在，但没有负成本的广告，故非负。）

一般线性规划：线性函数 $f(x_1..x_n)=\sum_ja_jx_j$；$f=b$ 为线性等式，$f\le b$、$f\ge b$ 为线性不等式，统称线性约束（linear constraints）；不允许严格不等式。线性规划问题：在有限个线性约束下最小化（最小化线性规划）或最大化（最大化线性规划）一个线性函数。已有多项式时间算法，但本章研究最古老的单纯形算法（simplex algorithm）——最坏情况非多项式，但相当高效、应用广泛。

两种规范形式：标准型（standard form）——在线性不等式约束下最大化线性函数；松弛型（slack form）——在线性等式约束下最大化。通常用标准型表达 LP，描述单纯形算法细节时用松弛型更方便。

二维示例：

\max\ x_1+x_2\quad\text{s.t.}\ 4x_1-x_2\le8,\ 2x_1+x_2\le10,\ 5x_1-2x_2\ge-2,\ x_1,x_2\ge0.\qquad(29.11\text{–}29.15)

满足所有约束的取值称为可行解（feasible solution）；可行解集合在平面上构成凸区域（convex region）（区域内任意两点连线上的点都在区域内），称为可行域（feasible region）；要最大化的函数称为目标函数（objective function），在某点的取值为目标值（objective value），取最大目标值的点为最优解（optimal solution）。可行域通常含无穷多个点，不能逐点计算。二维可用图解：$x_1+x_2=z$ 是斜率 −1 的直线，画出 z = 0、4、8 的直线（图 29.2(b)）；可行域有界，存在使直线与可行域相交的最大 z。本例最优解为 $x_1=2,x_2=6$，目标值 8。

最优解在顶点：使直线仍与可行域相交的最大 z 必在边界上，交集要么是单个顶点（唯一最优解），要么是一条线段（线段上各点目标值相同，端点即顶点也是最优解）。三维时每个约束对应一个半空间，交集为可行域，等目标值集合是平面；若目标系数全非负且原点可行，沿目标函数法向远离原点时目标值增加。n 维时每个约束定义 n 维空间的半空间，它们的交称为单纯形（simplex），目标函数是超平面，由于凸性，最优解仍在单纯形的某个顶点上。

单纯形算法思路：从单纯形的某个顶点出发，每次沿一条边移到目标值不更小（通常更大）的相邻顶点，直到达到局部最大（所有相邻顶点目标值都更小）。因为可行域凸、目标线性，局部最优就是全局最优（29.4 节用对偶性（duality）证明）。29.3 节用代数而非几何描述：把 LP 写成松弛型，用一些变量（基本变量，basic variables）表示另一些变量（非基本变量，nonbasic variables）；从一个顶点到另一个顶点就是让一个基本变量变为非基本、一个非基本变量变为基本，称为转轴（pivot），代数上就是把 LP 改写成等价的另一个松弛型。还需处理：无可行解的 LP、无有限最优解的 LP、原点不可行的 LP。

应用：运筹学教材中遍布 LP 例子，也是商学院标准工具。例如：航空公司在 FAA 约束（连续工作小时数上限、每个机组每月只飞一种机型等）下用最少机组人员安排所有航班；石油公司在有限预算下选择钻井地点使期望出油量最大。LP 也可用来建模和求解图论与组合问题：24.4 节的差分约束就是特例；29.2 节把若干图和网络流问题写成 LP；35.4 节用 LP 求另一图问题的近似解。

线性规划算法：单纯形算法精心实现时通常很快，但对精心构造的输入可能需要指数时间。第一个多项式时间算法是椭球算法（ellipsoid algorithm），实践中很慢。第二类多项式算法是内点法（interior-point methods）：不同于单纯形法沿可行域外表面移动、每步维持一个顶点可行解，内点法穿过可行域内部，中间解可行但不一定是顶点，最终解是顶点；对大规模输入，内点法与单纯形法一样快，有时更快。若再要求所有变量取整数，就是整数线性规划（integer linear program），仅判断是否存在可行解就是 NP 难的（习题 34.5-3），没有已知多项式算法；而一般 LP 可在多项式时间内求解。记号：变量 $x=(x_1..x_n)$ 的某个具体取值记为 $\bar x=(\bar x_1..\bar x_n)$。

29.1 标准型与松弛型（Standard and slack forms）（PDF p.871–879）

标准型：给定实数 $c_1..c_n$、$b_1..b_m$、$a_{ij}$，求 $x_1..x_n$：

\max\sum_{j=1}^nc_jx_j\quad(29.16)\qquad\text{s.t.}\ \sum_{j=1}^na_{ij}x_j\le b_i\ (i=1..m)\quad(29.17),\qquad x_j\ge0\ (j=1..n)\quad(29.18)

(29.16) 为目标函数，n+m 个不等式为约束，其中 (29.18) 为非负约束（nonnegativity constraints）——一般 LP 不必有，但标准型要求有。紧凑形式：

\max\ c^Tx\quad\text{s.t.}\ Ax\le b,\ x\ge0\qquad(29.19\text{–}29.21)

标准型可用三元组 $(A,b,c)$ 表示（A 为 $m\times n$）。术语：满足所有约束的 $\bar x$ 是可行解，否则是不可行解（infeasible solution）；目标值 $c^T\bar x$；目标值最大的可行解为最优解，其目标值为最优目标值。没有可行解的 LP 称为不可行的（infeasible），否则可行的；有可行解但无有限最优目标值的称为无界的（unbounded）。可行域无界时最优目标值仍可能有限（习题 29.1-9）。

转化为标准型：LP 不是标准型的四种原因及处理：

最小化目标 → 目标系数取负变为最大化（可行解集相同，目标值互为相反数）。两个 LP 等价的定义：对一个的每个目标值为 z 的可行解，另一个有目标值为 z（最小化对最大化时为 −z）的可行解，反之亦然（不要求一一对应）。
变量无非负约束 → 用 $x_j'-x_j''$ 替换 $x_j$，加 $x_j',x_j''\ge0$；目标与约束中的 $c_jx_j$、$a_{ij}x_j$ 相应替换。新解 $\hat x$ 对应原解 $\bar x_j=\hat x_j'-\hat x_j''$；原解对应 $\hat x_j'=\bar x_j,\hat x_j''=0$（若 $\bar x_j\ge0$）或 $\hat x_j''=-\bar x_j,\hat x_j'=0$（若 $\bar x_j<0$）。
等式约束 $f=b$ → 拆成 $f\le b$ 与 $f\ge b$。
≥ 约束 → 两边乘 −1：$\sum_ja_{ij}x_j\ge b_i\iff\sum_j-a_{ij}x_j\le-b_i$。

完整例子：

\min\ -2x_1+3x_2\quad\text{s.t.}\ x_1+x_2=7,\ x_1-2x_2\le4,\ x_1\ge0.

取负目标得 $\max\ 2x_1-3x_2$；$x_2$ 无非负约束，替换为 $x_2'-x_2''$（29.22）；等式拆成 ≤ 与 ≥（29.23），再把 ≥ 取负；重命名 $x_2'\to x_2$、$x_2''\to x_3$，得标准型

\max\ 2x_1-3x_2+3x_3\quad\text{s.t.}\ x_1+x_2-x_3\le7,\ -x_1-x_2+x_3\le-7,\ x_1-2x_2+2x_3\le4,\ x_1,x_2,x_3\ge0.\qquad(29.24\text{–}29.28)

转化为松弛型：为用单纯形法高效求解，希望除非负约束外都是等式。对不等式 $\sum_ja_{ij}x_j\le b_i$（29.29），引入新变量 s：

s=b_i-\sum_ja_{ij}x_j,\quad s\ge0.\qquad(29.30)(29.31)

s 称为松弛变量（slack variable），度量左右两边的差（松弛量）。从标准型转换时用 $x_{n+i}$ 表示第 i 个约束的松弛变量：$x_{n+i}=b_i-\sum_ja_{ij}x_j$，$x_{n+i}\ge0$（29.32）。上例引入 $x_4,x_5,x_6$：

\max\ 2x_1-3x_2+3x_3;\quad x_4=7-x_1-x_2+x_3,\ x_5=-7+x_1+x_2-x_3,\ x_6=4-x_1+2x_2-2x_3;\ \text{全部变量}\ge0.\qquad(29.33\text{–}29.37)

每个等式左边一个变量、右边其余变量；各等式右边变量集合相同，且只有它们出现在目标函数中。左边的称为基本变量，右边的称为非基本变量。省略“maximize”“subject to”和显式非负约束，并用 z 表示目标值，得到松弛型：

z=2x_1-3x_2+3x_3,\quad x_4=7-x_1-x_2+x_3,\quad x_5=-7+x_1+x_2-x_3,\quad x_6=4-x_1+2x_2-2x_3.\qquad(29.38\text{–}29.41)

松弛型的紧凑记号：N 为非基本变量下标集，B 为基本变量下标集，$|N|=n$，$|B|=m$，$N\cup B=\{1..n+m\}$；方程用 B 中下标索引，右边变量用 N 中下标索引。再加目标函数中的可选常数项 ν（便于读出目标值）。松弛型用六元组 $(N,B,A,b,c,\nu)$ 表示：

z=\nu+\sum_{j\in N}c_jx_j,\qquad x_i=b_i-\sum_{j\in N}a_{ij}x_j\ \ (i\in B),\qquad(29.42)(29.43)

所有变量非负。注意因为是减去 $\sum a_{ij}x_j$，$a_{ij}$ 是松弛型中“看起来”的系数的相反数。

例：松弛型

z=28-\tfrac{x_3}{6}-\tfrac{x_5}{6}-\tfrac{2x_6}{3},\quad x_1=8+\tfrac{x_3}{6}+\tfrac{x_5}{6}-\tfrac{x_6}{3},\quad x_2=4-\tfrac{8x_3}{3}-\tfrac{2x_5}{3}+\tfrac{x_6}{3},\quad x_4=18-\tfrac{x_3}{2}+\tfrac{x_5}{2}

中 $B=\{1,2,4\}$，$N=\{3,5,6\}$，

A=\begin{pmatrix}a_{13}&a_{15}&a_{16}\\a_{23}&a_{25}&a_{26}\\a_{43}&a_{45}&a_{46}\end{pmatrix}=\begin{pmatrix}-1/6&-1/6&1/3\\8/3&2/3&-1/3\\1/2&-1/2&0\end{pmatrix},\ b=\begin{pmatrix}8\\4\\18\end{pmatrix},\ c=(-1/6,-1/6,-2/3)^T,\ \nu=28.

下标不必连续，取决于 B 和 N。例如 $x_1$ 方程中有 $+x_3/6$，而 $a_{13}=-1/6$。

习题 29.1：

29.1-1 写出 (29.24)–(29.28) 的 n、m、A、b、c。
29.1-2 给出该 LP 的三个可行解及目标值。
29.1-3 写出 (29.38)–(29.41) 的 N、B、A、b、c、ν。
29.1-4 把一个最小化、含等式和无符号变量的 LP 化为标准型。
29.1-5 把一个 LP 化为松弛型并指出基本与非基本变量。
29.1-6 证明某 LP 不可行（$x_1+x_2\le2$ 与 $-2x_1-2x_2\le-10$ 矛盾）。
29.1-7 证明某 LP 无界。
29.1-8 一般 LP（n 变量、m 约束）化为标准型后变量数与约束数的上界（至多 2n 个变量、2m 个不等式约束，不计非负约束）。
29.1-9 举例：可行域无界但最优目标值有限。

29.2 把问题表述为线性规划（Formulating problems as linear programs）（PDF p.880–883）（续见下一块）

能识别何时可以把问题写成 LP 很重要：一旦写成多项式规模的 LP，就可用椭球法或内点法在多项式时间内求解，也可交给现成 LP 软件包。本节例子：已学过的单源最短路径和最大流，然后是最小费用流，最后是多商品流（其唯一已知的多项式算法基于 LP）。记号：LP 中用下标而非属性，如 $d_v$ 代替 $v.d$，$f_{uv}$ 代替 $(u,v).f$；输入量仍记为 $w(u,v)$、$c(u,v)$。

最短路径（单对，单源推广为习题 29.2-3）：给定带权有向图、源 s、终点 t，求 $d_t$。Bellman-Ford 终止时对每条边有 $d_v\le d_u+w(u,v)$，且 $d_s=0$。于是

\max\ d_t\quad\text{s.t.}\ d_v\le d_u+w(u,v)\ \ \forall(u,v)\in E,\quad d_s=0.\qquad(29.44\text{–}29.46)

为什么是最大化：若最小化，令所有 $\bar d_v=0$ 就是最优解，却没有解决最短路问题。最短路问题的解令每个 $\bar d_v=\min_{u:(u,v)\in E}\{\bar d_u+w(u,v)\}$，即不超过集合 $\{\bar d_u+w(u,v)\}$ 中所有值的最大值；要在这些约束下使 s→t 最短路径上所有顶点的 $d_v$ 最大，最大化 $d_t$ 即可达到。规模：|V| 个变量，|E|+1 个约束。（这与 24.4 节差分约束是同一结构。）

最大流：

\max\ \sum_{v}f_{sv}-\sum_vf_{vs}\quad\text{s.t.}\ f_{uv}\le c(u,v),\ \ \sum_vf_{vu}=\sum_vf_{uv}\ (u\in V-\{s,t\}),\ \ f_{uv}\ge0.\qquad(29.47\text{–}29.50)

$|V|^2$ 个变量（每对顶点一个流量），$2|V|^2+|V|-2$ 个约束。为了记号方便对非边也设流量与容量为 0；更高效的写法只需 $O(V+E)$ 个约束（习题 29.2-5）。

最小费用流（minimum-cost flow）：用 LP 求解已有高效专用算法的问题（Dijkstra、推送-重贴标签）通常不如专用算法高效，LP 真正的威力在于求解新问题（如开篇的政治问题），以及我们还不知道高效算法的问题变体。最小费用流是最大流的推广：每条边除容量 $c(u,v)$ 外还有实数费用 $a(u,v)$，在边上送 $f_{uv}$ 单位流的费用为 $a(u,v)f_{uv}$；给定需求 d，要从 s 送 d 单位流到 t 且总费用 $\sum_{(u,v)\in E}a(u,v)f_{uv}$ 最小。 例（图 29.3）：边（容量 c，费用 a）：$s\to x$（5,2），$s\to y$（2,5），$x\to y$（1,3），$x\to t$（2,7），$y\to t$（4,1），要送 4 单位。最优解：$f_{sx}=2,f_{sy}=2,f_{xy}=1,f_{xt}=1,f_{yt}=3$，总费用 $2\cdot2+5\cdot2+3\cdot1+7\cdot1+1\cdot3=27$。存在专门的多项式时间算法（超出本书范围），但可写成 LP：

\min\sum_{(u,v)\in E}a(u,v)f_{uv}\quad\text{s.t.}\ f_{uv}\le c(u,v);\ \sum_vf_{vu}-\sum_vf_{uv}=0\ (u\ne s,t);\ \sum_vf_{sv}-\sum_vf_{vs}=d;\ f_{uv}\ge0.\qquad(29.51)(29.52)

多商品流（multicommodity flow）（本块在此中断，续见下一块）：Lucky Puck 公司扩展产品线，除冰球外还运送球棍和头盔，每种产品有自己的工厂和仓库（球棍从温哥华运往萨斯卡通，头盔从埃德蒙顿运往里贾纳），但共享同一运输网络、容量不变。这是多商品流问题的实例：给定有向图，每条边容量 $c(u,v)\ge0$（非边容量为 0，无反平行边）……

第 24 章 单源最短路径（Single-Source Shortest Paths）（接上一块）

第 24 章引言（后半）（PDF p.667–671）

24.1 Bellman-Ford 算法（PDF p.672–676）

24.2 有向无环图中的单源最短路径（PDF p.676–679）

24.3 Dijkstra 算法（PDF p.679–685）

24.4 差分约束与最短路径（Difference constraints and shortest paths）（PDF p.685–691）

24.5 最短路径性质的证明（PDF p.692–698）

第 24 章思考题（Problems）（PDF p.699–703）

第 24 章注记（Chapter notes）（PDF p.703–704）

第 24 章小结

第 25 章 所有结点对的最短路径（All-Pairs Shortest Paths）

第 25 章引言（PDF p.705–707）

25.1 最短路径与矩阵乘法（PDF p.707–714）

25.2 Floyd-Warshall 算法（PDF p.714–721）

25.3 稀疏图上的 Johnson 算法（PDF p.721–726）

第 25 章思考题与注记（PDF p.726–728）

第 25 章小结

第 26 章 最大流（Maximum Flow）

第 26 章引言（PDF p.729–730）

26.1 流网络（Flow networks）（PDF p.730–735）

26.2 Ford-Fulkerson 方法（PDF p.735–752）

26.3 二分图最大匹配（Maximum bipartite matching）（PDF p.753–757）

★26.4 推送-重贴标签算法（Push-relabel algorithms）（PDF p.757–769）

★26.5 前置重贴标签算法（The relabel-to-front algorithm）（PDF p.769–781）

第 26 章思考题（PDF p.781–786）

第 26 章注记（PDF p.786–787）

第 26 章小结

第七部分 算法问题选编（Part VII Selected Topics）

第七部分导言（PDF p.789–792）

第 27 章 多线程算法（Multithreaded Algorithms）

第 27 章引言（PDF p.793–795）

27.1 动态多线程基础（The basics of dynamic multithreading）（PDF p.795–813）

27.2 多线程矩阵乘法（Multithreaded matrix multiplication）（PDF p.813–818）

27.3 多线程归并排序（Multithreaded merge sort）（PDF p.818–826）

第 27 章思考题与注记（PDF p.826–833）

第 27 章小结

第 28 章 矩阵运算（Matrix Operations）

第 28 章引言（PDF p.834）

28.1 求解线性方程组（Solving systems of linear equations）（PDF p.834–848）

28.2 矩阵求逆（Inverting matrices）（PDF p.848–853）

28.3 对称正定矩阵与最小二乘逼近（PDF p.853–861）

第 28 章思考题与注记（PDF p.861–863）

第 28 章小结

第 29 章 线性规划（Linear Programming）（续见下一块）

第 29 章引言（PDF p.864–871）

29.1 标准型与松弛型（Standard and slack forms）（PDF p.871–879）

29.2 把问题表述为线性规划（Formulating problems as linear programs）（PDF p.880–883）（续见下一块）

第 24 章单源最短路径（Single-Source Shortest Paths）（接上一块）

第 25 章所有结点对的最短路径（All-Pairs Shortest Paths）

第 26 章最大流（Maximum Flow）

第七部分算法问题选编（Part VII Selected Topics）

第 27 章多线程算法（Multithreaded Algorithms）

第 28 章矩阵运算（Matrix Operations）

第 29 章线性规划（Linear Programming）（续见下一块）