目录

1·1 搬沙子的问题：什么是最优传输 Abstract

这一节回答全书最基本的问题：最优传输到底在求什么？我们跟着原书摘要里 Monge 的铲土工人，把一个工程小问题一步步抽象成“比较两个概率分布”的数学问题，再用几个小演示提前感受它为什么在今天如此有用。Abstract

Monge 的铲土工人 A worker with a shovel

原书开篇借用了法国数学家 Gaspard Monge（加斯帕尔·蒙日，1746–1818）的说法：工地上有一大堆沙子，一位手持铁铲的工人要用全部这些沙子堆出一个形状事先规定好的目标沙堆——比如一座巨大的沙堡。她当然希望总力气最小，力气可以用“每一铲沙搬运的距离之和”或“花掉的总时间”来衡量。

下面的演示把工地画成侧视图：彩色方块是现在的沙堆，每个方块代表一铲沙；虚线轮廓是要堆成的沙堡。两者体积相同（都是 42 铲），而且沙堡的左塔恰好压在沙堆的右坡上——真实工地上挖方和填方区域往往也是交错的。

同样是把沙堆变成沙堡，搬法不同，总代价可以差很多：在“距离²”下，反着搬比最优方案多花 50% 以上，随手乱搬通常也要多花两三成。最优方案堆好的沙堡上，颜色仍然按“蓝→绿→黄→红”从左到右排列——搬运路线互不交叉。再切换局部代价：用“距离”时各方案差距几乎消失；用“√距离”时，“按顺序配对”反而不再最优。局部代价怎么定，决定了哪种搬法最省力。

从沙堆到概率分布 Comparing two probability distributions

数学家把这个问题换一种说法：两堆沙就是两个概率分布(probability distribution)——只要把总量都规定为 1。“沙子在哪里、有多少”对应“概率质量在哪里、有多少”；把每一列的铲数除以总铲数 42，就得到一个非负、总和为 1 的向量，也就是一个直方图(histogram)。“两堆沙体积相同”对应“两个分布的总质量都是 1”。

接着，原书用了两个关键的词：

• 局部代价(local cost)：把一粒沙从位置 $x$ 搬到位置 $y$ 要花多少，记作 $c(x,y)$。它由“地面”空间的几何决定，所以书中叫它地面代价(ground cost)。
• 全局代价(global cost)：一种搬法（运输方案transport plan）规定“哪里的沙搬去哪里、搬多少”，把所有局部代价按搬运量加起来，就是这种搬法的总代价。

“把第一堆变成第二堆”的方式有无穷多种；数学家关心的是其中最便宜的那一种：它有什么性质？怎样高效地算出来？下面用一个只有 3 个位置的小例子，把这几步走一遍。

亲手试一试：3 个仓库，3 家商店 A tiny transport problem

“搬沙子”也可以是“运货”。原书第 2 章用“仓库与工厂”的例子解释最优传输 Rem. 2.10，这正是 Hitchcock 1941 年研究的运输问题。下面是一个缩小版：3 个仓库分别存有 40、35、25 吨货，3 家商店分别需要 35、30、35 吨；每条路线的运费单价（元/吨）约等于两地距离（公里）取整。你来当调度员：在表格里填运量，看总运费能压到多低。

“西北角法”（335 元，从表格左上角开始依次填满）和“先填最便宜路线”的贪心法（390 元）都给出可行方案，但都不是最优；最优是 285 元。贪心法先把两条单价为 2 的路线填满，结果把 W1 剩下的 35 吨货逼上了单价 7 的最贵路线——每一步都省 ≠ 整体最省，这就是需要“优化算法”的原因。另外注意：3×3 的方案表里只有 4 个格子能自由填，其余 5 个被行和、列和“锁死”（见习题 1.4）。

把这个小游戏翻译成本书的记号（1·3 节会系统整理），你会发现后面各章反复出现的符号，其实你已经全部“用过”了：

最小代价本身就是一种“距离” A distance between distributions

原书摘要说：这个最小代价不仅定义了分布之间的一种距离，还在“所有概率分布组成的空间”上带来了丰富的几何结构；而且这种结构是自然的(canonical)——它借用了分布所在的地面空间(ground space)的几何性质。什么意思？看一个最简单的实验：把一个峰一点点向右平移，看不同的“距离”怎么反应。

OT 距离 $\mathrm{W}_1$ 与平移量成正比：峰搬得越远，距离越大。逐格比较的 $L^1$、$L^2$ 距离只看“同一个格子里两个直方图差多少”，两峰一旦不再重叠就到达上限、再也不变——它们分不清“隔壁”和“天涯”；KL 散度更极端，峰稍微一挪就变成无穷大。这就是“OT 继承了地面几何”的含义；第 8 章会系统比较这几类距离与散度。

分布空间里的几何：插值与重心 Interpolation and barycenters

有了距离，还能在分布之间“走路”。原书举例说：当地面空间是欧氏空间时，插值(interpolation)、重心(barycenter)、凸性(convexity)、函数的梯度(gradients of functions)这些概念，都能自然地推广到装备了 OT 几何的分布空间。下面比较两种从 $\mathbf{a}$ 走到 $\mathbf{b}$ 的方式。

线性混合是“淡出 + 淡入”：中途总有两个峰，一个在变矮、一个在变高。OT 插值是“把沙搬过去”：中途只有一个峰在滑动，宽度也逐渐变成目标的宽度；目标选“两个峰”时，它会平滑地一分为二。两种方式的平均位置（均值）其实一样，差别全在形状上。$t$ 可以看作重心权重 $\lambda=(1-t,\,t)\in\Sigma_2$：OT 插值恰好就是 $\mathbf{a}$、$\mathbf{b}$ 按这组权重的 Wasserstein 重心。

这本书关心什么 What this survey covers

摘要最后一段交代了本书的定位：以数值方法为主线，只讲那些能指导算法设计的理论性质；重点介绍最近这一波让 OT 在数据科学中流行起来的高效算法；给近几年提出的各种推广留出显著篇幅，并把它们与统计推断、核方法、信息论中的相关方法联系起来。书中所有插图都能用配套网站（optimaltransport.github.io）上的代码复现，那里还有幻灯片等资源。

1·2 历史脉络与本书视角 Introduction

这一节回答三个问题：最优传输为什么会在两百多年里被一再重新发现？它和日常生活有什么关系？这本书站在什么视角、应该怎么读？Introduction

最短路径原理：从“一个”到“一群” The shortest path principle

原书引言从一个朴素的原则说起：最短路径原理(shortest path principle)主导着生活与科学中的大多数决策——一件商品、一个人或一个比特的信息要从某处送到另一处，就应该尽量少费力气。在平面上这意味着走直线，在更复杂的度量空间里意味着走测地线(geodesic)。

最优传输把这个直觉推广到“同时搬运很多件物品（甚至是连续分布的物品）”的情形：从一种配置(configuration)搬成另一种配置。正如老师们都深有体会的——安排一群学生走到指定的位置（每人到了要恰好站成规定的队形），比带一个学生过去复杂得多。以“群体”或“分布”为单位来思考，需要更高级的数学形式，这最早见于 Monge 1781 年的开创性工作。

这个问题和日常生活联系紧密：无论运人、运货还是传输信息，很少只搬一件东西。物流、生产计划、网络路由这些重要的经济问题，都是在“搬运分布”。难怪最优传输的几篇奠基性文献——Tolstoi（1930）、Hitchcock（1941）、Kantorovich（1942）——都源于实际需要。

一再被重新发现的理论 A rich history

原书说，OT 在许多场合、以不同形式被（重新）发现，因而有着丰富的历史。下面的时间轴按原书摘要与引言的线索整理，每个节点都注明它在本书哪里再次出现。

注意两段长长的“空白”：Monge 之后一百五十年无人问津，直到物流与经济学的实际需要让它复活；1949 年有了算法之后，又过了约四十年，Brenier 等分析学家才发现它与凸性、偏微分方程、统计的深刻联系。每一次“复活”都来自一个新领域的新需求——最近一次来自数据科学，靠的是能处理大规模问题的快速近似算法。

今天用在哪里 Applications

摘要列举了近年 OT 大显身手的三个领域。它们的共同点是：要比较或变换的对象本身就是“分布”——一张图的颜色、一个形状上的点、一批数据样本。

本书的视角：计算优先 A computational viewpoint

引言接着说：到了 1980 年代以后，数学家们（尤其是 Brenier 1991 年的工作）发现这一理论是研究的沃土，与凸性、偏微分方程、统计都有深刻联系；世纪之交，计算机、图像乃至更广泛的数据科学研究者意识到，OT 提供了研究分布的强大工具——他们手里的数据恰好常以特征包(bags-of-features)或描述子(descriptors)的分布形式出现。

已有好几本 OT 专著，理论部分本身就值得几百页。既然 OT 已逐步成为一种应用工具，两位作者（Gabriel Peyré 与 Marco Cuturi）选择用计算的视角来平衡这些文献，重心放在数据科学（特别是图像科学与机器学习）上：先概述支撑 OT 实际有效性的主要理论洞见，再花更多篇幅讲如何把这些洞见变成快速的计算方案。

三层叙述：直方图、离散测度、一般测度 Three layers of exposition

这是读这本书最重要的一条“使用说明”。原书把同一套内容分三层来写：

• ① 正文：第 2、3、4、9、10 章的主体只研究最优传输在概率向量（离散直方图）空间上诱导的几何。
• ② 浅灰框：面向更进阶的读者，在同样几章里给出针对离散测度(discrete measures)的更一般的叙述。离散测度不仅由概率权重决定，还由这些权重所在的位置决定；位置通常取在连续的度量空间里——这为刻画随机现象提供了“第二个重要的自由度”。
• ③ 深灰框：最技术性的一层，处理任意测度——不必是离散的，特别地可以关于某个基准测度有密度。这是大多数经典 OT 教材的默认设定，但在实际应用中通常没那么重要。第 5–8 章讨论连续测度与离散测度之间的相互作用，因此面向数学更强的读者。

“等宽切分”时，位置等距、形状信息全在权重里；“随机抽样”时，第 ① 层的直方图变成一排一样高的柱子（每个都是 1/n），形状信息全在位置里——只看直方图完全看不出沙堆长什么样。这就是原书说位置是“第二个重要自由度”的意思。n 越大，第 ② 层越接近第 ③ 层的连续沙堆。

同一个“最优传输”，在三层里写出来长这样（切换选项卡对照）：

全书路线图 Roadmap of the book

结合原书引言的说明（第 2、3、4、9、10 章正文讲直方图；第 5–8 章讲连续与离散的相互作用）和各章开头的介绍，可以把全书画成下面这张图。

第 2 章是所有章节的地基；第 3 章（精确算法）和第 4 章（熵正则化）是两大计算工具，后面许多章都会用到，尤其是第 4 章。第 5–8 章各自深入一个“连续 ↔ 离散”的专题，可以按兴趣选读。

1·3 记号速查 Notation

这一节把原书开头的记号表Notation整理成“能查、能玩、能读”的三件工具：先记住两条排版规则，再用可筛选的速查表查任何符号，用记号实验室亲手算几个常用运算，最后逐个符号读懂全书最核心的公式 (2.11)。以后读到任何一章卡在符号上，都可以回到这里。

两条排版规则 Typographical conventions

原书的记号看似繁多，其实大多成对出现：同一个概念，在“直方图层”和“一般测度层”各有一个写法。记住两条规则就能猜出八九成：

1. 粗体 = 离散（向量与矩阵）：$\mathbf{a},\mathbf{b}$（直方图）、$\mathbf{C}$（代价矩阵）、$\mathbf{P}$（耦合矩阵）、$\mathbf{f},\mathbf{g}$（对偶变量）、$\mathbf{u},\mathbf{v},\mathbf{K}$（Sinkhorn）。
2. 希腊字母与花体 = 一般测度：测度用 $\alpha,\beta,\pi$；它们的集合与数值用花体 $\mathcal{U},\mathcal{R},\mathcal{L},\mathcal{W}$，空间用 $\mathcal{X},\mathcal{Y}$；对应的直方图版本用粗体或正体 $\mathbf{U},\mathbf{R},\mathrm{L},\mathrm{W}$。

可筛选的记号速查表 Notation table

下表按原书记号表的顺序列出全部符号，并补充了中文含义、原文、小例子和它主要出现的章节。可以按类别筛选，也可以直接搜索中文、英文或符号名（如 simplex、耦合、kernel）。

初读时选“直方图层”就够用：剩下的大约二十个符号覆盖了第 3、4 章算法所需的全部记号。测度层的符号在灰框里出现，读到时再回来查。

记号实验室：⊗、⊕、⊙ 与 Gibbs 核 Notation lab

光看定义不如动手算一算。下面的实验室把几个最容易混淆的运算做成“表格游戏”：左边一列是第一个向量（对应矩阵的行），上面一行是第二个向量（对应矩阵的列），中间是运算结果。

$\mathbf{a}\otimes\mathbf{b}$ 的每一行都与 $\mathbf{b}$ 成比例、每一列都与 $\mathbf{a}$ 成比例（“秩为 1”）；$\mathbf{f}\oplus\mathbf{g}$ 则是每行加同一个数、每列加同一个数。在 Gibbs 核里，把 ε 调小，远离“对角带”的格子迅速变成 0——Sinkhorn 就是在这张核上反复缩放行与列（$\operatorname{diag}(\mathbf{u})\mathbf{K}\operatorname{diag}(\mathbf{v})$），直到行和、列和分别等于 $\mathbf{a}$、$\mathbf{b}$（第 4 章）。

逐个符号读懂一个公式 Reading formula (2.11)

最后，用刚学的记号读懂全书最核心的定义——第 2 章的 Kantorovich 问题 (2.11)。它在 1·1 节里其实已经出现过好几次了，只是没写成公式。点「下一步」，每一步高亮一个部分。

本章小结

自测

下面的题目考的是理解而不是记忆；答错了可以点“重做”，并回到对应的演示再玩一玩。

习题

原书没有习题，以下为本教材根据本章内容自拟，建议先独立完成再展开提示与解答。

2.1 直方图与测度 Histograms and Measures

这一节回答：最优传输（OT）要比较的「东西」到底是什么？答案有两种写法——只有权重的直方图，和带位置的离散测度。把它们分清楚，后面的公式就不会乱。

本书把直方图(histogram)和概率向量(probability vector)当作同一个东西：$n$ 个非负数，加起来等于 1。所有这样的向量组成概率单纯形(probability simplex)：

原书说，这篇综述的很大一部分就是在研究「OT 在单纯形上诱导出的几何」——也就是：单纯形里两个点（两个直方图）之间，用 OT 量出来的「距离」长什么样。

三角形的顶点就是「只有一个 Dirac」的测度，边上是有一个权重为 0 的直方图。拖动右边的竖线时，直方图 $a$ 完全不变，但测度 $\alpha$ 变了（质心在移动）。OT 会同时用到两者：权重告诉它「有多少」，位置（通过代价 $C_{ij}=c(x_i,y_j)$）告诉它「搬多远」。

2.2 指派问题与 Monge 问题 Assignment and Monge Problem

最朴素的搬运方案是「一对一」：每个源点整体送到某一个目标点。这一节看这种确定性方案——指派问题和它的推广 Monge 问题——能走多远，又在哪里碰壁。

最优指派 Optimal Assignment

设有 $n$ 个源点、$m=n$ 个目标点，代价矩阵 $C$ 的元素 $C_{i,j}$ 表示「把第 $i$ 个源点配给第 $j$ 个目标点」的代价。最优指派问题(optimal assignment problem)在所有排列 $\sigma\in\mathrm{Perm}(n)$ 中找平均代价最小的一个：

这里 $\sigma(i)$ 就是「第 $i$ 个源点去哪儿」；排列保证每个目标恰好收到一个源点。最直接的办法是把所有排列试一遍——可惜排列有 $n!$ 个：

原书（引 Dantzig, 1983）指出：$n$ 只有 70 时，排列就超过 $10^{100}$ 个（比可观测宇宙中的原子数还多得多）。所以只有存在高效算法时，这个问题才真正「可解」——例如第 3 章的拍卖算法 §3.7。

把一个点慢慢拖过去，最优配对会在某一刻突然「整体换人」——指派是离散的选择，最优值却连续变化。选距离代价 $\|x-y\|$ 时，最优连线从不交叉（两条交叉的线段总能「解开」而不增加总长，这是三角不等式）；选距离平方时偶尔会交叉。点「正方形」看 Rem. 2.3：两种指派代价相同，最优解不唯一。读数最后一行预告了 Prop. 2.1：允许拆分质量的线性规划（§2.3）在这里并不能做得更好。

离散测度之间的 Monge 问题 Monge Problem between Discrete Measures

指派问题只能处理「点数相同、每个点质量相同」的情形。Monge 在 1781 年的原始问题更一般：两边的点带不同的质量。

推前算子 Push-forward Operator

记号 $T_\sharp\alpha$ 读作「$\alpha$ 在 $T$ 下的推前(push-forward)」。对离散测度，它的意思非常朴素：每个质点跟着 $T$ 搬家，质量不变 Rem. 2.5：

若几个 $x_i$ 被送到同一个 $y_j$，它们的质量就在 $y_j$ 处累加——这正是 (2.4) 的右边。所以「$T_\sharp\alpha=\beta$」只是在说：搬完之后，每个位置上的质量恰好等于 $\beta$ 要求的量。

2.3 Kantorovich 松弛 Kantorovich Relaxation

Monge 的「整堆搬运」太死板。Kantorovich（1942）的关键一步是：允许把一个点的质量拆开，送往多个目标。这一松，问题就从棘手的组合优化变成了一个线性规划——永远有解，而且可以高效求解。

原书先总结了指派 / Monge 表述的三个毛病：① 排列只能比较均匀、同样大小的直方图；② 对一般权重，Monge 的推前表述可能没有满足质量守恒 (2.4) 的可行解；③ 指派问题是组合问题，Monge 问题 (2.9) 的可行集是非凸的——两者按原始形式都很难解。

Kantorovich 放弃了「确定性」：源点 $x_i$ 的质量可以分派到若干位置。这种概率性的运输(probabilistic transport)就是今天常说的质量拆分(mass splitting)。描述它的对象不再是排列 $\sigma$ 或映射 $T$，而是一个耦合矩阵(coupling matrix) $P\in\mathbb{R}_+^{n\times m}$：$P_{i,j}$ = 从 $i$ 流向 $j$ 的质量。

三种耦合的行和、列和完全一样（都在 $U(a,b)$ 里），代价却差很多。独立耦合是一片「云」，每个源点都往所有目标送一点，代价比最优耦合高得多；反单调耦合把最小的 $x$ 送到最大的 $y$，是最贵的。最优耦合则是一条从左上到右下的细「阶梯」——小的 $x$ 送到小的 $y$，这就是 §2.6 要讲的一维单调重排，也是原书 Fig. 2.6 左图里耦合集中在 Monge 映射图像上的原因。注意最优耦合的非零格子很少：它是多面体 $U(a,b)$ 的一个顶点（第 3 章 §3.4）。

Monge 映射 vs Kantorovich 耦合 Monge Maps vs Kantorovich Couplings

同一对点云，两种搬法放在一起比一比：Monge 只能「整堆倒」，Kantorovich 可以「分着倒」。

「图 2.2 右」这组数据里，Monge 映射存在，而且 Kantorovich 的最优解恰好不拆分——两个最优值相同。点「反过来」或「一个点拆给两个」：Monge 映射根本不存在，Kantorovich 却照样给出答案，橙色虚线圈标出被拆分的源点。随便单击几个点改质量，你会发现 Monge 映射大多数时候都不存在——而一旦存在，永远有 Kantorovich ≤ Monge（把 Monge 映射写成耦合 $P_{i,T(i)}=a_i$，它就是 $U(a,b)$ 中的一个元素）。

置换矩阵也是耦合 Permutation Matrices as Couplings

指派问题其实是 Kantorovich 问题的一个「子问题」。对排列 $\sigma$，定义对应的（缩放过的）置换矩阵(permutation matrix)：

于是 $\langle C,P_\sigma\rangle=\frac1n\sum_i C_{i,\sigma_i}$——指派问题 (2.2) 就是把 Kantorovich 问题 (2.11) 的耦合限制为置换矩阵。置换矩阵只是 Birkhoff 多面体(Birkhoff polytope) $U(\mathbb{1}_n/n,\mathbb{1}_n/n)$ 的一小部分（例如 $\mathbb{1}_n\mathbb{1}_n^\top/n^2$ 可行但不是置换矩阵），在更大的集合上取最小只会更小：

下面的命题说：其实两边相等——对指派问题，松弛是「紧」的（tight）。

代价 $\langle C,P(t)\rangle$ 是 $t$ 的一次函数，斜率为 $C_{11}-C_{12}-C_{21}+C_{22}$：斜率为正，最优在左端点；为负，在右端点——总之在线段的端点，即多面体的顶点。点「a = b = (½, ½)」：两个顶点恰好是 ½ × 单位矩阵和 ½ × 交换矩阵——两个置换矩阵（Prop. 2.1）。选「正方形」：斜率为 0，整条线段都是最优解，这就是 Rem. 2.3 的「最优解不唯一」。

原书 Fig. 2.4 对比了两种情形：左边是两团同样点数的经验分布，最优耦合是一个完美匹配（Prop. 2.1）；右边两团带不同权重的点云无法一对一匹配，只能用一般的 Kantorovich 耦合——上面「Monge vs Kantorovich」演示里的「随机质量」就是这种情况。

应用一瞥 Rem. 2.12：最优运输方案本身（耦合 $P$ 或 Monge 映射 $T$）在数据科学中用途很广，尤其是图像处理——对比度均衡、纹理合成、图像匹配与融合、医学影像、形状配准、数字水印；在天体物理中用于重建早期宇宙，在音乐转录、经济学中的匹配数据分析里也有应用；OT 计算出的映射还可以用来做采样和贝叶斯推断。

2.4 最优传输的度量性质 Metric Properties of Optimal Transport

OT 不只是一个优化问题：当代价来自一个距离时，最优值本身就是直方图之间的一个距离。换句话说，OT 提供了一种标准的方法，把「点与点之间的距离」提升（lift）为「分布与分布之间的距离」。

先看直方图的情形。借用 Rubner 等人（2000）的说法，代价矩阵 $C$ 称为地面距离(ground metric)：它是固定的，描述「柱与柱之间的替换代价」，被所有要比较的直方图共用。

三角不等式的证明是本章最漂亮的构造之一。下面的步骤播放器用一个 5 格的小例子把「接力运输」画出来（位置 $0,1,2,3,4$，$D_{ij}=|i-j|$，$p=2$）。

应用 Rem. 2.16：OT 能自动把「柱之间的地面距离」提升为「直方图之间的距离」，这使它成为计算机视觉与机器学习中比较直方图的首选工具——例如把局部特征（图像描述子）汇总成「特征袋」（bag of features）直方图来做检索、聚类、分类。Rubner 等人（2000）的地标性论文开启了检索与聚类方面的应用（推土机距离(earth mover's distance, EMD)），近年又被用于词袋（bag-of-words）的降维与文本分类，以及训练输出词袋的多类分类器的损失函数。

几何直觉与弱收敛 Geometric Intuition and Weak Convergence

$\mathrm{W}_p$ 最重要的性质是：它是一个弱距离(weak distance) Rem. 2.18——它能比较支撑互不重叠的奇异分布（比如离散分布），并且能量化两个分布之间的空间位移。「经典」的距离或散度做不到这一点：$L^2$ 范数只能用于有密度的测度；离散的 $\ell^2$ 范数要求所有点落在一个事先固定的网格上。

最极端的例子是两个 Dirac：$\mathcal{U}(\delta_x,\delta_y)$ 里只有一个耦合 $\delta_{(x,y)}$（全部质量从 $x$ 搬到 $y$），于是对任意 $p$，

当 $x\to y$ 时 $\mathcal{W}_p(\delta_x,\delta_y)\to0$——而逐点比较的量对两个不同位置的 Dirac 永远给出「完全不同」。

$\mathrm{W}_1$ 和 $\mathrm{W}_2$ 都恰好等于 $\tau$：一条直线，平移多少就「看见」多少。逐格比较的 L² 距离和 0-1 代价的 OT（= 总变差，§2.6）一旦两个峰不再重叠就饱和了——平移 0.2 和平移 0.6 在它们眼里一样远。KL 散度更糟：对方块和尖峰，只要有一点平移就是 $\infty$。选「尖峰」看 Rem. 2.18：两个 Dirac 之间 $\mathrm{W}_p=|x-y|$，随 $\tau\to0$ 连续地趋于 0。

平移与均值 Translations

在欧氏空间 $\mathbb{R}^d$ 上取 $c(x,y)=\|x-y\|^2$ 时，$\mathcal{W}_2$ 可以把「平移」干净地分离出来 Rem. 2.19。记 $T_\tau:x\mapsto x-\tau$ 为平移，$m_\alpha=\int x\,\mathrm{d}\alpha(x)$ 为 $\alpha$ 的均值，则

特别地，把两个分布都挪到均值为 0 的位置（$\tilde\alpha=T_{m_\alpha\sharp}\alpha$，$\tilde\beta$ 同理），得到一个很好用的分解：

2.5 对偶问题 Dual Problem

每个线性规划都有一个「影子问题」——对偶问题。Kantorovich 问题是带约束的凸（线性）最小化，它的对偶是带约束的凹最大化。对 OT 来说，对偶变量有非常具体的含义：价格。这一节先写出对偶，再用原书的「承运商定价」故事把它讲透。

承运商的定价故事 A Pricing Story Rem. 2.21

回到 Rem. 2.10 的仓库与工厂。运营者本来要自己解线性规划 (2.11)，按 $P^\star$ 运货，付给运输公司 $\langle P^\star,C\rangle$。

• 外包物流(outsourcing logistics)：运营者没有能力解线性规划，于是把任务外包给一家承运商。承运商把物流拆成「取货」和「送货」：在仓库 $i$ 每取一单位收 $f_i$（不管送到哪），在工厂 $j$ 每送一单位收 $g_j$（不管从哪来）。仓库 $i$ 恰有 $a_i$ 单位、工厂 $j$ 恰需 $b_j$ 单位，所以总报价是 $\langle f,a\rangle+\langle g,b\rangle$。
• 定价(setting prices)：价格可以是负的！把所有 $f_i$ 减去同一个数、所有 $g_j$ 加上同一个数，总价不变（因为总取货量 = 总送货量）。承运商可以假装「倒贴钱帮你取货」，再在送货时加倍收回来。承运商当然想把 $f,g$ 定得越高越好。
• 检查价格(checking prices)：运营者没有竞争报价可比，又算不起 $\mathrm{L}_C(a,b)$，怎么判断报价是否合理？有个便宜得多的办法：按这套价格，从 $i$ 运一单位到 $j$ 实际花 $f_i+g_j$；而运输公司的公道价是 $C_{i,j}$。只要有某一对 $(i,j)$ 满足 $f_i+g_j\gt C_{i,j}$，承运商就在这条路线上多收了钱，运营者应当拒绝。
• 最优价格 = 对偶问题(optimal prices as a dual problem)：若 $n\times m$ 个不等式 $f_i+g_j\le C_{i,j}$ 全部成立，运营者可以放心接受：对任何运输方案 $P$（包括最优的 $P^\star$），利用行和、列和约束，

也就是说，自己动手一定不会比承运商的报价更便宜（这就是弱对偶）。于是承运商要在约束 $f_i+g_j\le C_{i,j}$ 下最大化 $\langle f,a\rangle+\langle g,b\rangle$——正是问题 (2.20)。§3.1 会证明（强对偶）：承运商能拿到的最高报价恰好等于运营者自己运的最低成本 $\mathrm{L}_C(a,b)$。

从全 0 开始逐个提价：报价上升，直到某些路线变成绿色（紧）——再提就会变红被拒绝。所能达到的最高报价恰好是 $\mathrm{L}_C(a,b)$（强对偶）。点「最优定价」后注意：粗框路线（$P^\star$ 真正在用的）全是绿色——这就是互补松弛 (2.23)。点「f 全减 1、g 全加 1」：每个 $f_i+g_j$ 和总报价都不变，所以最优价格不唯一，价格可以为负。地图上每个点旁边标的就是它的「价格」，原书 Fig. 2.8 用颜色画的正是这件事。

互补松弛：最优方案只走「紧」路线 Complementary Slackness

由拉格朗日函数 (2.22) 的原始-对偶最优性条件，可以定位最优运输方案的支撑（亦见 §3.3）：

Brenier 定理：什么时候最优耦合就是一个映射 Brenier's Theorem Rem. 2.24

Kantorovich 松弛允许拆分质量；可是在很重要的一类情形下，最优耦合其实不拆分——它就是一个 Monge 映射。这就是 Brenier 定理（1991；同一时期还有 Cuesta–Matrán、Rüschendorf–Rachev 的类似结果，以及 Knott–Smith 的先驱工作）。

2.6 特殊情形 Special Cases

一般情况下，计算 OT 需要专门的数值算法（第 3、4、7 章）。但有几类问题可以直接写出答案：0-1 代价、一维分布、高斯分布。它们既是很好的直觉来源，也是检验算法的「标准答案」。

0-1 代价与总变差 Binary Cost and Total Variation

推广到一般测度 Rem. 2.27：代价 $c(x,y)$ 在 $x=y$ 时为 0、否则为 1（Kronecker 代价）时，OT 距离就是两个测度之间的总变差(total variation)距离（见 Example 8.2；在 $\mathrm{TV}(\alpha,\beta)=\sup_A|\alpha(A)-\beta(A)|$ 的约定下，它对离散测度恰为 $\tfrac12\|a-b\|_1$）。这样的「地面距离」没有任何几何：远近一视同仁，这正是 §2.4 平移演示中那条会饱和的红色曲线。

无论怎么拖，0-1 代价的 OT 都精确等于 ½‖a − b‖₁（不是 ‖a − b‖₁）。对比最后两个按钮：「δ₁ → δ₂」与「δ₁ → δ₅」的 0-1 代价 OT 都是 1——它分不出「搬到隔壁」和「搬到最远处」；而用 $|i-j|$ 作地面距离的 $\mathrm{W}_1$ 分别是 1 和 4。OT 的几何感完全来自地面距离。

一维情形：排序就是最优传输 1-D Case

在实数轴上，OT 有一个极其简单的答案：把两边各自排好序，按顺序配对。

对一般的（带权重的）一维测度，把「排序」说精确就要用到累积分布函数和它的「反函数」。

这就是原书 Fig. 2.11 的做法：用累积分布函数与分位数函数计算一维 OT。紫色连线永远不交叉——最优耦合是单调的。左右两张图的阴影面积相等：一个是「横着」积分（$x$ 方向比较 $\mathcal{C}_\alpha$ 与 $\mathcal{C}_\beta$），一个是「竖着」积分（$r$ 方向比较分位数）——同一块区域的两种切法。读数中的第三个数是通用线性规划求得的精确值，三者一致。选「等权、等数目」时，分位数函数都是阶梯宽度相同的台阶，(2.36) 就退化为排序后逐对相减的 (2.33)。

直方图均衡化与位移插值 Histogram Equalization and Displacement Interpolation

Rem. 2.29 一维 OT 可以做直方图均衡化(histogram equalization)，用于规范灰度图像的色调。把两幅同样大小（$n\times m$ 像素）图像的灰度值（0 = 黑，1 = 白）按固定顺序排成向量 $\bar x$、$\bar y$；各自排序得到置换 $\sigma_1,\sigma_2$，则 $\sigma=\sigma_2\circ\sigma_1^{-1}$ 就是两个灰度分布之间的最优指派。把第 $i$ 个像素的灰度换成 $\bar y_{\sigma(i)}$，就得到一幅灰度分布与 $\bar y$ 完全相同的「均衡化」图像。更进一步，$t\mapsto(1-t)\bar x_i+t\,\bar y_{\sigma(i)}$ 在原图与均衡化图像之间插值，其像素分布正是两者之间的位移插值(displacement interpolation)（定义见 (7.7)；原书 Fig. 2.10 展示了 $t=0,0.25,\dots,1$ 的图像与直方图）。

线性混合只是把 $\alpha$ 调暗、把 $\beta$ 调亮，中间时刻出现「两处都有一半」的奇怪分布；位移插值则让质量真的移动过去，形状平滑变形。两者的均值相同（都是 $(1-t)m_\alpha+tm_\beta$），差别全在形状。选「两个高斯」：位移插值在每个时刻仍是高斯，且标准差线性变化 $\sigma_t=(1-t)\sigma_\alpha+t\sigma_\beta$（原书 Fig. 2.14）；线性混合则是两个峰的叠加。选「直方图均衡化」：这正是 Fig. 2.10 下排直方图的变化过程。

高斯分布之间的距离 Distance between Gaussians

$\mathcal{W}_2^2$ 永远分成两块：均值之间距离的平方 + 协方差之间的 Bures 距离平方——这正是 Rem. 2.19「形状差异 + 均值位移²」的具体例子。点「只平移」：两个椭圆形状相同，$\mathcal{B}=0$，$A=\mathrm{I}$，映射是纯平移。点「对角协方差」：两个椭圆的轴都与坐标轴平行，Bures 距离退化为 Hellinger 距离 $\|\sqrt r-\sqrt s\|_2$（读数中两者相等）。无论怎么拖，$A$ 总是对称正定的——所以 $T$ 是凸函数 $\psi$ 的梯度，由 Brenier 定理它是最优的。插值过程中的每个中间分布也都是高斯。

对角与一维的特例。若 $\Sigma_\alpha=\mathrm{diag}(r_i)$、$\Sigma_\beta=\mathrm{diag}(s_i)$，Bures 距离就是 Hellinger 距离(Hellinger distance) $\mathcal{B}(\Sigma_\alpha,\Sigma_\beta)=\|\sqrt r-\sqrt s\|_2$。特别地，一维高斯 $\mathcal{N}(m,\sigma^2)$ 之间 $\mathcal{W}_2^2=(m_\alpha-m_\beta)^2+(\sigma_\alpha-\sigma_\beta)^2$——$\mathcal{W}_2$ 就是 $(m,\sigma)$ 平面上的欧氏距离，位移插值沿这个平面上的直线段走（下图）。可以证明 $\mathcal{B}$ 是协方差矩阵之间的距离，并且 $\mathcal{B}^2$ 对两个参数都是凸的。

本章小结

自测

前面各节已穿插了 8 道小题，下面再来几道综合题。

习题

原书没有习题，以下为本教材自拟，覆盖计算、概念辨析与简短证明。

function w1d(x, y, p){
  const xs = x.slice().sort((u, v) => u - v), ys = y.slice().sort((u, v) => u - v);
  let s = 0; for (let i = 0; i < xs.length; i++) s += Math.pow(Math.abs(xs[i] - ys[i]), p);
  return Math.pow(s / xs.length, 1 / p);
}
// 验证：const n = x.length, a = new Array(n).fill(1 / n);
// Math.pow(Book.ot.emd(a, a, Book.ot.cost(x, y, p)).cost, 1 / p) 应与 w1d(x, y, p) 相同

第 2 章把最优传输写成了 Kantorovich 问题 (2.11) 和它的对偶问题 (2.20)。这一章回答一个非常实际的问题：怎样把它精确地算出来？答案来自组合优化与线性规划（linear programming）——这两个领域的诞生本身就和「运输问题」纠缠在一起。

3.1 Kantorovich 线性规划 The Kantorovich Linear Programs

这一节要解决的问题：把 OT 问题「翻译」成线性规划的标准语言。一旦翻译完成，线性规划几十年积累的理论（对偶、顶点、单纯形法）就都能直接拿来用。

回忆第 2 章的原始问题：在所有耦合（coupling）$P\in\mathrm{U}(a,b)$ 中，找总运费最小的那个。

目标是线性的（$\sum C_{ij}P_{ij}$），约束也全是线性的（行和 = $a$、列和 = $b$、非负）——这正是线性规划(linear program, LP)。

标准形：把矩阵拉直成向量 Standard Form

线性规划的标准形(standard form)是：线性目标 $c^\top p$，等式约束 $Ap=q$，变量非负 $p\ge 0$。要把 (3.1) 写成这个样子，只需把 $n\times m$ 的矩阵 $P$ 按列拉直成长度为 $nm$ 的向量 $p$（第 $i+n(j-1)$ 个分量就是 $P_{ij}$），代价矩阵 $C$ 也同样拉直成 $c$。于是

这里 $A$ 是一个 $(n+m)\times nm$ 的 0/1 矩阵：前 $n$ 行各负责一个「行和」约束，后 $m$ 行各负责一个「列和」约束。原书用 Kronecker 积把它写成一行公式，但真正要记住的只有一句话：

$A$ 的每一列都是「一个蓝点 + 一个红点」。把所有蓝色行（行和约束）加起来、和把所有红色行（列和约束）加起来，得到的都是全 1 的行向量——所以这 $n+m$ 个约束线性相关，秩恰好是 $n+m-1$。

对偶问题与弱对偶 The Dual Program and Weak Duality

按照线性规划的对偶理论，(3.2) 的对偶问题是

把 $h$ 拆成两段 $h=(f,g)$，$f\in\mathbb R^n$ 对应行约束、$g\in\mathbb R^m$ 对应列约束。因为 $A$ 的第 $(i,j)$ 列只在第 $i$ 个和第 $n+j$ 个位置是 1，所以 $(A^\top h)_{(i,j)}=f_i+g_j$。约束 $A^\top h\le c$ 就是逐格的 $f_i+g_j\le C_{ij}$，目标就是 $\langle f,a\rangle+\langle g,b\rangle$——正是第 2 章的对偶问题 (2.20)（原书此处印作 (2.4)）。

3.2 C-变换 C-Transforms

这一节要解决的问题：对偶问题有 $n+m$ 个变量 $(f,g)$，它们之间互相牵制。能不能先固定一半，把另一半「一步到位」地求出来？可以——这就是 C-变换(C-transform)。它把对偶问题化成只关于 $f$ 的问题，是第 5 章半离散最优传输和本章拍卖算法的核心工具。

从第 2 章的对偶形式出发（$\mathrm R(C)$ 是满足 $f_i+g_j\le C_{ij}$ 的所有 $(f,g)$）：

于是对偶问题可以化成一个只关于 $f$、没有约束的问题：

因为 $f^C$ 的每个分量都是若干个关于 $f$ 的仿射函数取最小，(3.5) 的目标是分段仿射的凹函数（piecewise affine concave）——很好优化的一类函数，但它不可微，这也是第 4 章要引入熵正则的原因之一。

对称地，固定 $g$ 时最好的 $f$ 是 $g^{\bar C}$。那么从某个 $f$ 出发，交替做 C 与 C̄ 变换，目标会一路上升：

人们自然希望这样一直严格上升直到最优。可惜不行——很快就会停在一个「平台」上：

(iii) 的意思是：做完「C、C̄」一轮之后，再怎么交替都原地不动了。下面的演示用一维的点和平方代价 $C_{ij}=(x_i-y_j)^2$ 把这件事画出来：

• 对每个源点 $x_i$，画一条抛物线 $y\mapsto (y-x_i)^2-f_i$（顶点在 $(x_i,-f_i)$）；
• 每个目标点 $y_j$ 上的对偶变量画成红点 $(y_j, g_j)$；
• 约束 $f_i+g_j\le C_{ij}$ ⇔ 红点不高于任何一条蓝抛物线，即不高于它们的下包络（lower envelope）。

于是 C-变换 = 把每个红点移到下包络上（从可行的位置出发，就是往上顶，直到碰到包络）；C̄-变换 = 把每条抛物线上下平移，直到恰好碰到某个红点（从可行的位置出发，就是往下压）。

从 $f=g=0$ 出发：C-变换后每个红点都顶到了包络上（右图每列出现 0）；C̄-变换后每条抛物线都压到了某个红点上（每行也出现 0）。之后再按哪个按钮都不动了——目标值卡在最优值下面很远的地方。「放入最优」后，接触点恰好就是最优运输的路线：要到达最优，必须让一组抛物线和红点协同移动，而不是一次只动一边。

3.3 互补松弛 Complementary Slackness

这一节要解决的问题：原始问题和对偶问题可以分别求解，得到最优的 $P^\star$ 和 $(f^\star,g^\star)$。它们之间有什么关系？答案是一个极其好用的「最优性证书」——后面的网络单纯形、对偶上升、拍卖算法，停机判据全都是它。

反过来也成立。先给「互补」一个名字：

选「最优耦合」时，所有圆盘都落在绿框（松弛 = 0）里，原始代价与对偶目标完全相等——这就是可行 + 互补的最优证书。换成西北角解或独立耦合，圆盘开始落在松弛为正的格子上，「差距」一栏恰好等于 Σ 圆盘质量 × 格子松弛：每一份放错地方的质量都按它所在格子的松弛「罚款」。

3.4 运输多面体的顶点 Vertices of the Transportation Polytope

这一节要解决的问题：线性规划的最优值总能在可行集的「角」上取到，所以只需要在 $\mathrm U(a,b)$ 的角上找。那么这些角长什么样？答案非常漂亮：角就是支撑图里没有环的耦合——稀疏、呈树状，至多 $n+m-1$ 个非零元。

线性规划的基本定理（Bertsimas & Tsitsiklis 1997，Thm. 2.7）：可行集非空且有界时，最小值一定在某个顶点上取到。$\mathrm U(a,b)$ 是有界的运输多面体(transportation polytope)，所以找最优 $P$ 只需在它的顶点中找。

3.4.1 所有顶点的支撑都是树 Tree Structure of the Support of All Vertices of U(a, b)

把运输问题画成二部图(bipartite graph)：左边 $n$ 个源点 $1,\dots,n$，右边 $m$ 个汇点 $1',\dots,m'$（加撇只是为了区分两边），每对 $(i,j')$ 之间一条边，边上的代价是 $C_{ij}$。一个运输方案就是这个图上的一个流：从源点 $i$ 流出 $a_i$，流入汇点 $j'$ 的总量为 $b_j$。

记 $S(P)=\{(i,j'):P_{ij}\gt 0\}$ 为 $P$ 的支撑(support)，$G(P)$ 为只保留这些边的图。

最小的非平凡例子是 2×3：只有 $(2-1)(3-1)=2$ 个自由度，所以整个 $\mathrm U(a,b)$ 可以画在平面上。取 $x=P_{11}$、$y=P_{12}$ 作坐标，其余 4 个元素都由边缘约束算出；每个元素 $\ge 0$ 给出一个半平面，六个半平面围成一个（至多）六边形。

在多边形内部，6 个元素全为正，二部图（5 个节点、6 条边）一定有环；在边上，一个元素为 0，还剩 5 条边、仍然有一个环，沿环扰动恰好是沿着这条边的方向；只有在顶点，两个元素为 0，剩下 4 = n+m−1 条边构成一棵树，再也没有可以挪动的方向。代价是 $(x,y)$ 的线性函数，等高线是平行直线，最低的那条总是碰在某个顶点上。

3.4.2 西北角法则 The North-West Corner Rule

怎样快速得到一个顶点？西北角法则(north-west corner rule)只需不超过 $n+m$ 步，它常被用来初始化在原始问题上工作的算法（比如下一节的网络单纯形）。

跟着原书的例子 $a=(0.2,0.5,0.3)$、$b=(0.5,0.1,0.4)$ 走一遍（左：表格与剩余量；右：支撑图逐步长成一棵树）：

所有「先重排、再做西北角、再排回来」得到的方案组成集合 $\mathcal N(a,b)=\{\mathrm{NW}_{\sigma^{-1}\sigma'^{-1}}(a_\sigma,b_{\sigma'})\}$。它们都至多有 $n+m-1$ 个非零元；不同的 $(\sigma,\sigma')$ 可能给出同一张表（指数多个）；并且 $\mathcal N(a,b)$ 只是顶点集合的一个子集（通常是真子集，Brualdi 2006 Cor. 8.1.4）。

3.5 网络单纯形法的直观描述 A Heuristic Description of the Network Simplex

这一节要解决的问题：怎样从一个顶点出发，一步步走到最优顶点？把前两节拼起来，思路就出来了。设 $P$ 是一个支撑无环的可行解（一个顶点）。由 Prop. 3.3，只要找到一对与 $P$ 互补、并且可行的对偶 $(f,g)$，就证明了 $P$ 最优。网络单纯形法(network simplex)建立在两条简单的原则上：

1. 对每个顶点 $P$，都能配出一对与它互补的 $(f,g)$——沿着树「传播」一下就行；
2. 如果这对 $(f,g)$ 可行，停机，$P$ 最优；否则就修改 $P$（保持可行），让它的互补对偶「更接近可行」。

3.5.1 求与 P 互补的对偶变量 Obtaining a Dual Pair Complementary to P

互补要求：对支撑中的每条边 $(i,j')\in S(P)$，都有 $f_i+g_j=C_{ij}$。设 $S(P)=\{(i_1,j_1'),\dots,(i_s,j_s')\}$，这是 $s$ 个方程、$n+m$ 个未知数的线性方程组：

因为 $s\le n+m-1\lt n+m$，方程组总是欠定的。多出的自由度有两个来源：一是 Rem. 3.1 说的约束冗余（整体上 $f+t$、$g-t$ 永远是另一组解）；二是当 $G(P)$ 不连通、是由几棵树组成的森林时，每棵树各自多出一个自由度——未定的对偶变量个数恰好等于连通分支的个数。例如原书图 3.3 中有 $5+6=11$ 个对偶变量，却只有 8 条边、8 个方程，分成 3 棵树。

解法非常简单：在每棵树里任选一个根，令它的对偶变量为 0，然后沿树边传播（广度或深度优先都行）。每走过一条边 $(i,j')$，它的一端已知，另一端就由 $g_j=C_{ij}-f_i$ 或 $f_i=C_{ij}-g_j$ 唯一确定。

3.5.2 网络单纯形的更新 Network Simplex Update

对每个格子算一下约化代价(reduced cost)$\ r_{ij}=C_{ij}-f_i-g_j$（这是单纯形法的标准术语，原书没有起名字）。树上的边 $r_{ij}=0$；若所有 $r_{ij}\ge 0$，$(f,g)$ 可行，由 Prop. 3.3 $P$ 已经最优。否则找一个「违规」的格子 $(i,j')$：$f_i+g_j\gt C_{ij}$，即 $r_{ij}\lt 0$。

先令图 $G$ 等于 $G(P)$，再把违规的边 $(i,j')$ 加进去。会出现两种情况：

• (a) $G$ 仍是森林——$(i,j')$ 连接了两棵原本分开的树。这时 $P$ 完全不变，只是两棵树合成一棵、少了一个自由度，重新按 §3.5.1 计算 $(f,g)$ 即可。这种更新叫退化更新(degenerate update)：边进入了 $G$，但 $P_{ij}$ 仍是 0（此后 $G(P)\subseteq G$，$G$ 里可以有流量为 0 的边）。
• (b) $G$ 出现了一个环。这时要从环上删掉一条边使 $G$ 恢复成森林，同时修改 $P$ 保持可行、并且 $G(P)\subseteq G$。做法和 Prop. 3.4 的证明如出一辙：从新边出发给环定向，$i_1\to j_1'\to i_2\to j_2'\to\cdots\to j_l'\to i_{l+1}=i_1$（其中 $i_1=i$、$j_1=j$），「$i\to j'$」方向的边标 +，「$j'\to i$」方向的边标 −，然后沿环推流：

$\theta$ 取「能推多少推多少」：受限于环上 − 边里最小的流量，$\theta=\min_k P_{i_{k+1},j_k}$。取到最小值的那条 − 边流量降为 0，把它从 $G$ 中删掉（出基(leaving edge)），而 $(i,j')$ 留下（入基(entering edge)）。最后按 §3.5.1 重新计算 $(f,g)$。

3.5.3 原始解的改进 Improvement of the Primal Solution

虽然更新的动机是「让对偶更接近可行」，但它同时保证了原始代价只降不升。情形 (a) 中 $P$ 不变；情形 (b) 中：

第二个等号是因为：除新边外，环上所有边都在旧的 $G$ 里，满足 $C=f_i+g_j$；把这些 $f_i+g_j$ 代进交错和，每个节点的对偶变量恰好一加一减互相抵消，只剩下新边的 $C_{ij}-f_i-g_j$。由于 $r_{ij}\lt 0$，只要 $\theta\gt 0$ 代价就严格下降；$\theta=0$ 时（例如某条 − 边的流量本来就是 0）只是图 $G$ 变了而 $P$ 不变，这是单纯形法中常见的退化换基。（原书此处把条件印成了 $C_{ij}\lt f_i-g_j$，应为 $C_{ij}\lt f_i+g_j$。）

例 1 的西北角解支撑只有 3 棵树（边缘分布的部分和在 0.2 与 0.75 处重合），前两次加边都是情形 (a)：把树连起来，$P$ 不变。之后每次换基都沿一个环推流，代价按 $\theta\times r_{ij}$ 精确下降；当所有约化代价 $\ge 0$ 时停机，最终代价与 Book.ot.emd 独立算出的最优值一致。随机问题里偶尔会看到 $\theta=0$ 的退化换基：图变了，代价没变。

3.6 对偶上升法 Dual Ascent Methods

这一节要解决的问题：能不能从对偶这一侧出发？网络单纯形始终保持 $P$ 可行，一步步让对偶变得可行；对偶上升法(dual ascent methods)反其道而行——始终保持 $(f,g)$ 可行，每次给一组 $f_i$ 加价、一组 $g_j$ 降价，让对偶目标严格上升，直到能找到一个与之互补的可行 $P$ 为止。

这类方法比网络单纯形还早几十年：可以追溯到 Borchardt 与 Jacobi（1865），以及后来的 König 和 Egerváry（Kuhn 1955 有回顾）。其中最有名的是匈牙利算法(Hungarian algorithm)，但它只处理 $a=b=\mathbf 1_n/n$ 的特殊情形（指派问题）。下面介绍一般的对偶上升框架（原书的叙述主要依据 Bertsimas & Tsitsiklis 1997 §7.7），需要用到最大流(maximum flow)问题。

由互补松弛，最优的 $P$ 只能用平衡边运货。所以对偶上升的每一步都在问：只走平衡边，能不能把所有质量运完？能，就找到了互补的可行 $P$，结束；不能，就调价，让新的边变成平衡边。调价的方式由下面两个命题给出。

用 Prop. 3.6 构造的 $S,S'$ 作方向、Prop. 3.5 允许的步长作 $\varepsilon$，就得到经典的原始—对偶法(primal-dual method)；把它用在指派问题上，就是匈牙利算法。

每一轮：先只沿平衡边求最大流（流量从上一轮继续增加，已有的流不会失效）；运不完时，标号集合 $S$ 的总供给严格大于 $S'$ 的总需求——「卡住」的正是这批仓库。于是给 $S$ 加价、给 $S'$ 降价，对偶目标上升 $\varepsilon\,(a(S)-b(S'))$，并出现一条新的平衡边。有时流量不增加、只是标号集合变大，但目标总在上升；当流量达到 1 时，得到的 $P$ 与网络单纯形的最优代价一致。

对偶上升与「对偶网络单纯形」有相似之处，但至少有两点重要区别：

• 单纯形类方法的迭代点总是可行集（对偶时是 $\mathrm R(C)$）的顶点；对偶上升没有这个限制，迭代点可以在可行集的内部。
• 对偶网络单纯形构造的 $P$ 始终满足边缘约束，只有收敛时才满足非负；对偶上升构造的 $P$（即最大流）始终非负，但只有收敛时才满足边缘约束。

3.7 拍卖算法 Auction Algorithm

这一节要解决的问题：对于指派问题（$n$ 个人、$n$ 件物品，一一配对，使总代价 $\sum_i C_{i\sigma_i}$ 最小），有没有比前面更简单、更容易并行的算法？拍卖算法(auction algorithm)由 Bertsekas（1981）提出，后与 Eckstein（1988）一起改进，有很自然的经济学解释（Bertsekas 1992）。它也可以推广到任意边缘分布，但这里只讲指派问题。

指派问题中的互补松弛 Complementary Slackness

指派问题的顶点都是置换矩阵 $P_\sigma$，互补松弛因此变得格外简单。设 $\sigma^\star$ 与 $(f^\star,g^\star)$ 为原始、对偶最优解，则 $f^\star_i+g^\star_{\sigma^\star_i}=C_{i,\sigma^\star_i}$。又由 §3.2 的 C-变换原理，可以取 $f^\star=g^{\star\bar C}$，于是

反过来，只要某个 $g$ 和置换 $\sigma$ 满足

则 $\sigma$ 是最优指派，$(g^{\bar C},g)$ 是最优对偶。

部分指派与 ε-互补松弛 Partial Assignments and ε-Complementary Slackness

算法维护一个三元组 $(S,\xi,g)$：$S\subset\llbracket n\rrbracket$ 是已经分到物品的人，$\xi$ 是从 $S$ 到物品的单射（部分指派），$g$ 是对偶向量（负的价格）。每次迭代后保持三条性质：

• (a) ε-互补松弛（ε-CS）：对所有 $i\in S$，$\ C_{i,\xi_i}-g_{\xi_i}\le\varepsilon+\min_j\ C_{ij}-g_j$ ——每个已分配的人都「差不多满意」，离最好的选择至多差 $\varepsilon$；
• (b) $S$ 的大小不会减少；
• (c) 有某个 $g_j$ 至少下降 $\varepsilon$（某件物品的价格至少上涨 $\varepsilon$）。

拍卖更新 Auction Algorithm Updates

挑一个还没分到物品的人 $i\notin S$，找出他眼中最好和次好的物品：

$j^1_i\in\operatorname{argmin}_j\ C_{ij}-g_j,\qquad j^2_i\in\operatorname{argmin}_{j\ne j^1_i}\ C_{ij}-g_j .$

1. 出价（更新 $g$）：把 $g_{j^1_i}$ 减去「次好与最好的差 + ε」：
   $$g_{j^1_i}\ \leftarrow\ g_{j^1_i}-\underbrace{\Big[(C_{i,j^2_i}-g_{j^2_i})-(C_{i,j^1_i}-g_{j^1_i})+\varepsilon\Big]}_{\ge\,\varepsilon\,\gt\,0}\ =\ C_{i,j^1_i}-(C_{i,j^2_i}-g_{j^2_i})-\varepsilon .$$(3.9)
2. 换主人（更新 $S$ 与 $\xi$）：若有 $i'\in S$ 原本拿着 $j^1_i$，把他踢出：$S\leftarrow S\setminus\{i'\}$；然后令 $\xi_i=j^1_i$，$S\leftarrow S\cup\{i\}$。

先用一个 3×3 的小例子手算一遍：代价 $C=\left[\begin{smallmatrix}0 & 1 & 5 \\ 1 & 1 & 5 \\ 4 & 1 & 3\end{smallmatrix}\right]$，$\varepsilon=0.5$。人 1、人 2 都觉得物品 1′、2′ 便宜，人 3 最喜欢 2′——三个人抢两件好物品，会打一场小小的「价格战」。下面用价格 $p_j=-g_j$ 叙述（出价就是涨价），每一步都对应式 (3.9)。

每次出价只动一件物品的价格，而且只涨不跌；被挤掉的人之后面对更高的价格，会转向别的物品。结束时总代价与最优值（勾选「显示最优指派」可对照）至多差 $n\varepsilon$：把 ε 调大，出价次数变少，但可能停在一个稍差的指派上；把 ε 调小，结果更准，但几个人可能围绕同一件物品反复小幅加价。

算法性质 Algorithmic Properties

从 $S=\varnothing$、没有任何指派、$g=\mathbf 0$ 开始，反复执行上面两步，直到 $S=\llbracket n\rrbracket$（人人都有物品）为止。性质 (b)、(c) 从 (3.9) 一眼可见；一开始 $S$ 为空，ε-CS 自然成立，而下面的命题保证它一直成立。

「拥挤」时所有人看各物品的代价都差不多，大家轮流对同几件物品小幅加价（每次只涨 ε 左右），ε 缩小 10 倍，出价次数就大幅增加——这正是 Prop. 3.8 的 $1/\varepsilon$ 因子。ε-scaling 先用大 ε 把价格「粗调」到位，出价次数少得多。右图中所有点都在虚线 $n\varepsilon$ 之下；ε 足够小时，结果直接就是最优指派。点分布均匀时，价格战不明显，普通拍卖本来就很快。

所以拍卖算法可以看成使用 C-变换机制的另一种方式——它在 C-变换中加入了「次好」的信息和一点点 ε，从而打破了 §3.2 中交替变换的僵局。下一章介绍另一种思路：正则化，即 Sinkhorn 算法，它和拍卖算法也有不少相似之处（Schmitzer 2016b）。最后值得一提：在低维欧氏空间的规则网格上，这些经典的线性规划求解器还可以和多尺度（multiscale）策略结合，获得显著的加速（Schmitzer 2016a；Oberman & Ruan 2015）。

本章小结

自测

习题

第 3 章的精确算法（网络单纯形、拍卖算法）能算出精确的最优传输，但它们的复杂度大约是 $O(n^3\log n)$，难以并行，而且得到的最优值对输入不光滑——没法方便地求导、放进机器学习模型里。本章换一个思路：不再精确求解，而是在目标里加一个熵惩罚，把问题「软化」。原书开篇列出了这样做的五个好处，它们也正是本章的线索：

4.1 熵正则化 Entropic Regularization

这一节要解决的问题：最优传输是一个线性规划，它的解稀疏（大部分元素为 0）、可能不唯一、对输入不光滑。我们给目标加上一个「鼓励摊开」的熵项，得到一个严格凸的新问题；一个参数 $\varepsilon$ 控制摊开的程度。

离散熵与正则化问题 Discrete Entropy and the Regularized Problem

两个小细节值得记住：

• 那个「$-1$」只是为了公式好看：它让 $\partial H/\partial P_{ij}=-\log P_{ij}$。对总质量为 1 的 $P$，$H(P)$ 就等于香农熵（Shannon entropy）加 1。
• $H$ 是 1-强凹(1-strongly concave) 的：它的 Hessian 是 $-\operatorname{diag}(1/P_{ij})$，而 $P_{ij}\le 1$。所以 $-H$ 是一个「碗形」的强凸函数，加到目标上会让最优解唯一。

先在最简单的情形里看熵起了什么作用。原书Fig. 4.1在概率单纯形 $\Sigma_3$（一个三角形）上最小化线性函数 $\langle c,p\rangle$：这是一个线性规划，解总落在顶点上（最稀疏的点）。加上 $-\varepsilon H(p)$ 之后，解有闭式 $p_\varepsilon\propto e^{-c/\varepsilon}$（即 softmax），它随 $\varepsilon$ 增大离开边界、走向三角形的「熵中心」$(\tfrac13,\tfrac13,\tfrac13)$。

$\varepsilon$ 很小时解贴在顶点上（稀疏）；$\varepsilon$ 变大，解沿灰色路径进入三角形内部，最终到达中心——所有分量都为正，最「模糊」。点「让两个顶点代价相同」：这时整条底边都是线性规划的最优解（不唯一），而 $\varepsilon\to 0$ 时 $p_\varepsilon$ 恰好收敛到底边中点，也就是最优解中熵最大的那一个——这正是下面 Prop. 4.1 的内容。

ε 的两个极限 Convergence with ε

换句话说：$\varepsilon$ 很小时，熵项只在「平局」时起作用（在一堆同样便宜的方案里挑最分散的）；$\varepsilon$ 很大时代价几乎被忽略，得到的是「最随意」的方案——两个边缘分布的独立乘积，相当于每个源按比例把质量撒给所有目标。

$\varepsilon$ 小时，$P_\varepsilon$ 收缩成精确 OT 解那条细细的「单调阶梯」（勾选黑点对比），代价逼近精确值；$\varepsilon$ 大时，热力图变成 $a$ 与 $b$ 的外积——每一行都是 $b$ 的缩放，代价升到 $\langle C,ab^\top\rangle$，互信息降到 0。代价和互信息此消彼长：$\varepsilon$ 就是在「省钱」与「随机」之间调节的旋钮。

KL 投影与 Gibbs 核 KL Projection and the Gibbs Kernel

(4.2) 还有一个非常有用的改写：它是某个矩阵在 KL 散度意义下的「投影」。

离散测度与一般测度 Discrete and General Measures

$\varepsilon$ 小时每个红点几乎只连向一个蓝点（精确 OT 在等质量时是一个置换，恰好 10 条线）；$\varepsilon$ 增大，每个点的质量被撒向越来越多的目标，连线变成一团「毛线」，代价上升。打开重心投影：它把每个 $x_i$ 送到它所连目标的加权平均，$\varepsilon$ 越大，所有箭头越被拉向目标点云的中心。

4.2 Sinkhorn 算法及其收敛性 Sinkhorn's Algorithm and Its Convergence

这一节要解决的问题：怎样把 $P_\varepsilon$ 真正算出来？关键发现是最优解的结构极其简单——只需给 Gibbs 核的每一行、每一列各乘一个缩放因子。于是 $nm$ 个未知数变成 $n+m$ 个，交替地让行和、列和对上即可。这就是 Sinkhorn 算法。

解的形式：缩放后的 Gibbs 核 Form of the Solution

矩阵缩放与 Sinkhorn 迭代 Regularized OT as Matrix Scaling

把 $P=\operatorname{diag}(u)K\operatorname{diag}(v)$ 代入两个边缘约束 (4.13)，由于 $\operatorname{diag}(v)\mathbb 1_m=v$，约束化简为

$\odot$ 是逐元素乘法。在数值分析里这叫矩阵缩放问题(matrix scaling problem)。两组方程各自都很好解：固定 $v$，第一组直接给出 $u=a/(Kv)$；固定 $u$，第二组给出 $v=b/(K^\top u)$。交替进行，就是 Sinkhorn 算法：

K ← exp(−C/ε);   v ← 1
重复，直到 ‖u ⊙ (K v) − a‖₁ < tol：
    u ← a ./ (K v)        # 行缩放：之后每一行的和恰好是 a
    v ← b ./ (Kᵀ u)       # 列缩放：之后每一列的和恰好是 b
输出 P = diag(u) · K · diag(v)

• 缩放因子不唯一，耦合唯一。若 $(u,v)$ 满足 (4.13)，则 $(\lambda u,v/\lambda)$ 也满足（$\lambda\gt0$）；不同的初值会得到不同的 $u,v$，但 $\operatorname{diag}(u)K\operatorname{diag}(v)$ 总是同一个最优耦合。
• 一定收敛。原因有两种讲法：交替 KL 投影（Rem. 4.8）和 Hilbert 度量下的压缩（Rem. 4.14），下面都会讲到。
• 演化的样子。耦合从 Gibbs 核 $K$（一条对角带）出发，一步步把质量「挪」向满足边缘的位置（原书 Fig. 4.5 上排）。

下面用一个 4×4 的小例子手动走几步。位置 $x=y=(0,1,2,3)$，$C_{ij}=|x_i-y_j|^2$，$\varepsilon=1$；源 $a=(0.4,0.3,0.2,0.1)$ 偏左，目标 $b=(0.1,0.2,0.3,0.4)$ 偏右，所以质量要整体向右搬。

第一次行缩放后，红色实心条与空心条完全重合（行 ✓），蓝色的却不重合；下一步反过来。随着迭代进行，两边的误差都越来越小，热力图从对角带 $K$ 逐渐「弯」成最终的耦合。右图中误差曲线在对数坐标下接近直线——线性收敛；把 $\varepsilon$ 调小再重置，直线明显变平（变慢）。

历史与复杂度 History and Complexity

Rem. 4.5 这个迭代被独立发现过很多次，名字也很多：统计学里叫 迭代比例拟合(iterative proportional fitting procedure, IPFP)，经济学里叫 RAS 方法，指派问题里叫 softassign。收敛性由 Sinkhorn（1964）证明，算法因此得名。2013 年 Cuturi 指出它能高效、可并行地近似 OT，并给出可微的损失函数，此后在机器学习、视觉、图形学中广泛使用。

Sinkhorn 的另外几种视角 Other Viewpoints (Remarks 4.8–4.11)

近端点算法 Rem. 4.9。如果想要的是 $\varepsilon=0$ 的精确解，可以反复求解 $P^{(\ell+1)}=\operatorname*{argmin}_{P\in U(a,b)}\mathrm{KL}(P|P^{(\ell)})+\tfrac1\varepsilon\langle P,C\rangle$ (4.17)。每一步都是一个 Sinkhorn 问题，只是核换成 $e^{-C/\varepsilon}\odot P^{(\ell)}$；展开后相当于用核 $e^{-C/(\varepsilon/\ell)}$——也就是让正则化参数按 $\varepsilon/\ell$ 逐渐减小。这与「$\varepsilon$ 退火」（从大 $\varepsilon$ 开始、逐步减小，用上一次的解热启动）的思路一脉相承（Kosowsky–Yuille 1994）；在低维空间里配合多尺度稀疏网格（Schmitzer 2016）非常高效。

其他正则项 Rem. 4.10。把 $-H(P)$ 换成任何严格凸的 $R(P)$ 也可以，例如二次正则 $R(P)=\sum_{ij}P_{ij}^2$ (4.18)，或 Tsallis 熵。区别在于：熵里的 $\log$ 自动把 $P$ 挡在正象限内，不需要再加 $P\ge0$；二次正则却必须显式加上非负约束，约束集不再是仿射的，交替 Bregman 投影不再收敛，要改用 Dykstra 算法（Benamou 等 2015）。

重心投影 Rem. 4.11。耦合 $P$ 不是映射；若想要一个「把 $x_i$ 送到哪里」的映射，可以取它所连目标的加权平均（重心投影(barycentric projection)）：

它只在 $\alpha$ 的支撑点上有定义。若 $P$ 是置换矩阵，$T$ 就是 Monge 映射；$\varepsilon\to0$ 时（若 Monge 映射唯一）它收敛到 Monge 映射，$\varepsilon\gt0$ 时给出它的光滑近似。一般测度的版本见 (4.20)。应用包括图像的颜色迁移（原书 Fig. 9.6）。在上面 §4.1 的点云演示里勾选「重心投影」就能看到它。

为什么收敛：Hilbert 射影度量 Convergence via the Hilbert Projective Metric

收敛性证明的关键观察：$u$ 与 $\lambda u$ 给出同一个耦合，所以不该比较向量本身，而该比较它们的方向（从原点出发的射线）。Hilbert 射影度量正是这样一种「射线之间的距离」。

Rem. 4.13 典型应用是 Perron–Frobenius 定理的定量证明：若 $K\gt0$ 且 $K^\top\mathbb 1=\mathbb 1$（列随机矩阵），它把概率单纯形 $\Sigma_n$ 映入自身且是严格压缩，因此存在唯一的不变分布 $p^\star=Kp^\star$，并且 $d_{\mathcal H}(K^\ell p_0,p^\star)\le\lambda(K)^\ell d_{\mathcal H}(p_0,p^\star)$——线性收敛（原书 Fig. 4.8）。

实用含义：(4.23) 说明边缘误差 $\|P^{(\ell)}\mathbb 1_m-a\|_1$、$\|P^{(\ell)\top}\mathbb 1_n-b\|_1$ 是可靠的停止准则；由 (4.21)，对偶变量 $f^{(\ell)}=\varepsilon\log u^{(\ell)}$ 在 $\|\cdot\|_{\mathrm{var}}$ 下也线性收敛。（Linial 等 1998 还证明了每次迭代都会增大缩放后矩阵的积和式 permanent。）

每条曲线在对数坐标下都近似直线：线性收敛。$\varepsilon$ 每缩小约 3 倍，直线就明显变平。表中「理论」一列的 $1-\lambda(K)^2$ 小到可以忽略（虚线几乎是水平的），而实测的每步收缩要快得多——这就是 Rem. 4.15 说的「全局界很悲观」。切到「随机」代价：同样的 $\varepsilon$ 慢得多，小 $\varepsilon$ 在 600 步内几乎不动。把 $\omega$ 调到 1.5～1.6：本来收敛慢的小 $\varepsilon$ 曲线，迭代次数减少到约 1/3～1/5；而本来就快的大 $\varepsilon$ 曲线反而会变慢——最合适的 $\omega$ 取决于收敛速度本身。$\omega$ 再大（小 $\varepsilon$ 时 1.7 左右起）就会振荡甚至发散，所以实用的外推算法都要配合安全措施。

4.3 加速 Sinkhorn 迭代 Speeding Up Sinkhorn's Iterations

这一节要解决的问题：Sinkhorn 每一步的瓶颈是乘以 $K$ 和 $K^\top$，朴素实现要 $O(nm)$ 次运算。能不能更快？原书给了四个方向：同时算很多对（把矩阵–向量乘法变成矩阵–矩阵乘法）、利用网格结构（乘以 $K$ = 卷积）、在曲面上用热方程近似、以及外推加速收敛。

批量并行与 GPU Parallel and GPU-Friendly Computation

Rem. 4.16 若要同时计算 $N$ 对直方图 $(a_1,b_1),\dots,(a_N,b_N)$ 的正则化距离，把它们排成矩阵 $A=[a_1,\dots,a_N]$（$n\times N$）、$B=[b_1,\dots,b_N]$（$m\times N$），$N$ 个 Sinkhorn 可以齐步走：

除法仍是逐元素的。一次矩阵–矩阵乘法 $KV$ 就完成了 $N$ 次矩阵–向量乘法——这正是 GPU 最擅长的运算（BLAS-3）。收敛后 $N$ 个正则化代价也能用几次矩阵运算一并算出。这是 Sinkhorn 能在大规模直方图比较中流行的根本原因。

可分离核：网格上的卷积 Separable Kernels on Grids

Rem. 4.17 这是最重要的加速。设每个下标都是 $d$ 维网格上的一个多重下标 $i=(i_1,\dots,i_d)$，并且代价沿各坐标可加：$C_{ij}=\sum_k C^k_{i_kj_k}$。那么 Gibbs 核就可分离为乘积：

于是计算 $Ku$ 时根本不必构造 $n\times n$ 的大矩阵（$n=\prod_k n_k$ 随维数指数增长），只要沿每个坐标轴依次作用一维的 $K^k$。当 $n=m$ 时，复杂度从 $O(n^2)$ 降到 $O(n^{1+1/d})$。

最典型的例子：$[0,1]^d$ 上的规则网格，代价 $\|x-y\|_q^q=\sum_k|x_k-y_k|^q$，一维核 $K^k_{rs}=\exp(-|\tfrac{r-s}{n_k}|^q/\varepsilon)$ 就是一个卷积。$d=2$、$q=2$ 时它是高斯模糊：把 $u$ 排成 $n_1\times n_2$ 的图像 $U$，则 $Ku$ 对应 $K^1UK^2$——先对每一列、再对每一行做一维高斯模糊，只需 $n_1^2n_2+n_1n_2^2$ 次运算，而不是 $n_1^2n_2^2$。正是这一点让 Sinkhorn 能处理 $200^3=8\times10^6$ 个格点的三维形状（Solomon 等 2015；Bonneel 等 2016）。

上排三幅图说明：「乘以 $K$」= 用宽度约 $\sqrt{\varepsilon/2}$ 的高斯核模糊图像，而二维高斯模糊可以拆成两次一维模糊。32×32 网格已有 1024 个点，稠密的 $K$ 有一百万个元素；可分离卷积每次只需 $2N^3$ 次乘法，快十几倍，网格越大优势越大（64×64 时稠密矩阵要 134 MB，这里干脆不构造）。$\varepsilon$ 变大时，箭头变短、所有箭头都被拉向目标的中心——重心投影被「平均」掉了；右下图中一个点的质量被撒得更开。

近似卷积与热核 Approximated Convolutions and the Heat Kernel

Rem. 4.18 更一般地，只要核是平移不变的（$K_{ij}=k_{i-j}$，均匀网格上总是如此），$Kv=k\star v$ 就是卷积，可以用近线性时间的算法近似：快速傅里叶变换（假设周期边界；但 $\varepsilon$ 小时数值不稳定，常常不可接受），或用一串递归滤波器近似高斯卷积（Deriche 1993；比较见 Getreuer 2013）。

Rem. 4.19 在曲面（流形）上，测地距离的核 $e^{-d_{\mathcal M}^p/\varepsilon}$ 不是卷积，但 Varadhan 公式说它与热扩散密切相关：

• $p=2$：$\sqrt{-4t\log\mathcal H_t(x,y)}=d_{\mathcal M}(x,y)+o(t)$ (4.29)，$\mathcal H_t$ 是热方程 $\partial_tg=\Delta_{\mathcal M}g$ 的热核。所以乘以 $K$ ≈ 让热量扩散 $t=\varepsilon/4$ 的时间：用离散 Laplace 矩阵 $L$（如余切 Laplacian），做 $R$ 步隐式 Euler，即乘以 $(\mathrm{Id}-\tfrac tRL)^{-R}$。
• $p=1$：$-\tfrac{\sqrt t}{2}\log\mathcal P_t(x,y)=d_{\mathcal M}(x,y)+o(t)$，$\mathcal P_t=(\mathrm{Id}-t\Delta_{\mathcal M})^{-1}$ (4.28)：取 $\varepsilon=\sqrt t/2$，乘以 $K$ 换成解一个线性方程组 $(\mathrm{Id}-tL)w=v$。

这些线性系统每次迭代都要解，但矩阵不变，可以预先做一次稀疏 Cholesky 分解；二维曲面上每次 Sinkhorn 迭代因此只需近线性的代价（Crane 等 2013 的「热方法求测地线」；Solomon 等 2015 首先与 Sinkhorn 结合）。

外推加速 Extrapolation Acceleration

Rem. 4.20 Sinkhorn 是一个不动点迭代（Rem. 4.15），所以可以套用经典的线性或非线性外推。最简单的是过松弛(over-relaxation)，类似逐次超松弛法（SOR）：在对数域里把每次更新写成 $f\leftarrow(1-\omega)f+\omega\,\hat f$（$\hat f$ 是普通 Sinkhorn 的新值，$1\lt\omega\lt2$）。它把局部线性速度从 $O((1-\kappa)^\ell)$ 提升到 $O((1-\sqrt\kappa)^\ell)$——$\kappa$ 越小（$\varepsilon$ 越小）收益越大。回到上面 §4.2 的收敛速度演示，把 $\omega$ 调到 1.6 就能亲眼看到：$\varepsilon=0.003$ 时每轮收缩因子从约 0.95 变成约 0.80，接近 $1-\sqrt\kappa$（这里 $\kappa\approx0.045$，$\sqrt\kappa\approx0.21$），迭代次数少了约 4/5。但 $\omega$ 不能太大，否则会振荡发散；本来就收敛很快（$\kappa$ 大）时，外推也没有好处。

4.4 稳定性与对数域计算 Stability and Log-Domain Computations

这一节要解决的问题：当 $\varepsilon$ 相对于代价 $C$ 很小时，朴素的 Sinkhorn 会数值溢出或下溢（Rem. 4.7）。办法是改存缩放因子的对数 $f=\varepsilon\log u$、$g=\varepsilon\log v$。它们恰好就是熵正则问题对偶的变量：从对偶的角度看，Sinkhorn 是块坐标上升，每一步是一个 soft-min，可以用经典的 log-sum-exp 技巧稳定地算出来。

熵正则化的对偶问题 The Entropic Dual

Sinkhorn = 对偶上的块坐标上升 Sinkhorn as Block Coordinate Ascent

Rem. 4.21 记 (4.30) 的目标为 $Q(f,g)$。它对 $f$、$g$ 的梯度是

梯度 = 目标边缘 − 当前边缘。固定 $g$、令 $\nabla_fQ=0$，可以精确解出最优的 $f$；再固定 $f$ 解 $g$：

这正是 (4.15) 取对数：$(f^{(\ell)},g^{(\ell)})=\varepsilon(\log u^{(\ell)},\log v^{(\ell)})$。所以每半步 Sinkhorn 都在一组变量上把对偶目标精确最大化，对偶值单调上升。

soft-min 视角 Soft-Min Rewriting

Rem. 4.22 用 soft-min 改写 (4.35)–(4.36)：

对比第 3 章的 c-变换 $f_i=\min_j(C_{ij}-g_j)$：Sinkhorn 就是在交替做「软 c-变换」（也叫熵 c-变换，见 §5.3 的 (5.11)），再加上一个 $\varepsilon\log a_i$ 的修正。原书把对矩阵每行、每列取 soft-min 的操作记为 $\mathrm{Min}^{\mathrm{row}}_\varepsilon$、$\mathrm{Min}^{\mathrm{col}}_\varepsilon$，于是 $f^{(\ell+1)}=\mathrm{Min}^{\mathrm{row}}_\varepsilon(C-\mathbb 1_ng^{(\ell)\top})+\varepsilon\log a$ (4.40–4.41)。

$\varepsilon$ 大时 soft-min 明显低于 min，而且权重分散在所有分量上；$\varepsilon$ 小时权重几乎全压在最小的分量上，soft-min 贴近 min——所以它是 min 的光滑版本。把平移 $s$ 调大（比如 100 以上）并把 $\varepsilon$ 调小：$e^{-z_i/\varepsilon}$ 全部下溢为 0，朴素公式得到 $-\varepsilon\log 0=+\infty$；稳定公式先减去最小值，指数里的数都 $\le0$ 且最大的恰为 0，求和至少是 1，永远不会出问题。

对数域 Sinkhorn Log-Domain Sinkhorn

$u,v$ 的数量级范围大约与 $1/\varepsilon$ 成正比：$\varepsilon=10^{-2}$ 时还在 $10^{\pm 30}$ 以内，$\varepsilon=10^{-3}$ 时已达 $10^{\pm 130}$，$\varepsilon\approx3\times10^{-4}$ 时撞上 $10^{\pm308}$ 的墙，出现 $+\infty$、$0/0$（非数），朴素迭代崩溃。与此同时 $K$ 里越来越多的元素下溢为 0（$\varepsilon\lt1/745$ 时开始出现）。而对数域里的 $f,g$ 始终是普通大小的数（$|f|$ 不超过代价的量级），算法照常运行——只是 $\varepsilon$ 很小时收敛本身变慢（300 轮后误差仍较大，§4.2）。

4.5 最优传输代价的正则化近似 Regularized Approximations of the Optimal Transport Cost

这一节要解决的问题：我们真正关心的往往是精确的 OT 代价 $\mathrm L_C(a,b)$（比如 Wasserstein 距离），而 Sinkhorn 给出的是一堆近似量。它们各自是上界还是下界？差多少？停在有限步时还靠得住吗？另外，熵正则代价本身光滑、凸、梯度现成，是很好的损失函数。

对偶解给出下界 Dual Potentials Are Feasible

这是熵正则化在机器学习中大受欢迎的原因之一：跑完 Sinkhorn，梯度就免费得到了——$f^\star$ 告诉你「把一点质量从 $a$ 的哪个位置挪走，代价下降得最多」。第 9 章会把它当作损失函数来拟合模型。

Sinkhorn 散度：夹住真实的 OT 代价 Sinkhorn Divergences

有限步停止时 Finite Sinkhorn Divergences

上面的关系只在 Sinkhorn 收敛后才成立，而实践中总是跑有限的 $L$ 步就停。好消息是对偶那一侧在任何时候都靠得住：

• 原始迭代不可行 Rem. 4.26：中途的 $P^{(\ell)}$ 边缘不对，$\langle C,P^{(\ell)}\rangle$ 不是上界——前几步它甚至远小于真值（每个源的质量只在附近打转，还没被「推」到目标去）。要得到可靠的上界，先用 Rem. 4.6 的取整把它变成可行耦合，这一步同样可以批量并行（Lacombe 等 2018）。
• 非凸但可微 Rem. 4.27：$\mathfrak D^{(L)}_C(a,b)$ 一般不是 $(a,b)$ 的凸函数（数值上很容易验证），但它是可微的，可以对 $C$、$a$、$b$ 做自动微分（§9.1.3）。

红线（原始 $\mathfrak P^\varepsilon$）永远在黑线（精确值）之上，蓝线（对偶 $\mathfrak D^\varepsilon$）永远在其下，间隙 = $\varepsilon\times$香农熵，随 $\varepsilon$ 减小而收窄。注意 $t=0$ 时 $a=b$，精确代价是 0，而 $\mathfrak P^\varepsilon(a,a)\gt0$：熵正则化带有偏差；去偏后的曲线在 $t=0$ 处回到 0，并且紧贴精确值。打开「有限 L 步」：蓝色虚线 $\mathfrak D^{(L)}$ 始终在下方（下界不失效）；紫色点线 $\langle C,P^{(L)}\rangle$ 在 $L$ 很小（比如 1）时明显低于真值——它不是上界；$L$ 增大后它才逐渐爬回红线。再把 $\varepsilon$ 调大：对偶下界变得非常松，而原始上界离真值并不远。

4.6 广义 Sinkhorn Generalized Sinkhorn

这一节要解决的问题：OT 的许多变体——边缘质量可以不守恒的非平衡 OT、Wasserstein 重心、梯度流——都有同样的结构：熵项，加上「只依赖于两个边缘的函数」。能否用同一个 Sinkhorn 框架求解？可以：只需把两步缩放里的「$a/(Kv)$」换成一个邻近算子。

这套迭代的价值在于 $\operatorname{Prox}^{\mathrm{KL}}_F$ 常常有闭式或很容易算。特别地，若 $F(s)=\sum_iF_i(s_i)$ 可分离，邻近算子就是逐个坐标的标量优化问题。对某些 $F,G$ 可以证明线性收敛，也能推广到一般测度，并同样可以在对数域中稳定地实现（Chizat 等 2018）。

$\lambda$ 大时解接近平衡 OT：右下那团质量几乎全部被长途运到中间（热力图下方、第二团对应的那几行出现一块），运输代价高。$\lambda$ 小时，把它运过去不如「违约」划算：它大部分被销毁（左侧红色实心条在下方几乎消失），中间的目标也收不满（蓝色实心条低于空心条），总质量小于 1。把 $b$ 的质量调得与 $a$ 不同：平衡 OT 此时根本无解，而非平衡 OT 照样给出方案，总质量落在两者之间。（$\lambda$ 再小到与 $\varepsilon$ 相当时，熵项本身会开始「凭空」增加质量，所以滑块下限取在 $\lambda\approx0.02$。）

(4.49) 还能进一步推广到多于两个函数、多于一个耦合的情形（Chizat 等 2018）：多边缘问题的 Sinkhorn（§10.1 的 (10.2)）、熵正则的 Wasserstein 重心（§9.2 的 (9.15) 与迭代 (9.18)）、梯度流（§9.3 的 (9.26)）都是它的特例。第 10 章的非平衡 OT 迭代 (10.9) 正是上面那个例子。

本章小结

自测

习题

前几章处理的大多是"离散对离散"：两堆点之间求一个耦合矩阵 $\mathbf P$。本章换一个非常常见、而且几何味道很浓的场景：一边是连续分布 $\alpha$（有密度，比如正方形上的均匀分布、一张图片的亮度），另一边是离散分布 $\beta=\sum_{j=1}^m \mathbf b_j\delta_{y_j}$（m 个"站点"，第 j 个需要质量 $\mathbf b_j$）。这就是半离散最优传输(semidiscrete optimal transport)。原书第 2 章的 Fig. 2.5 已经预告过这三种情形：

原书开篇一段话点出了本章的三条主线：

• 几何：低维时，从连续密度到离散测度的最优传输是一个分片常数映射——每个站点 $y_j$ 的"原像"是一块（或几块不相连的）区域。代价为欧氏距离平方时，这些区域就是计算几何里的 Laguerre 单元(Laguerre cells)：把 Voronoi 单元的边界按常数平移而得。于是可以借用计算几何的快速算法。
• 随机优化：高维时没法把区域精确切出来，但半离散问题可以看成一个随机规划(stochastic programming)问题，只靠从 $\alpha$ 采样来求解；再加一点熵正则效果更好——这把第 4 章的熵正则推广到了新场景。
• 工具：一切都建立在 c-变换(c-transform) 之上。第 3 章 §3.2 只对代价矩阵讲过它，这里推广到一般的代价函数。c-变换是凸分析中 Legendre 变换(Legendre transform) 的推广，在 OT 的理论和算法中都处于枢纽地位。

5.1 c-变换与 c̄-变换 c-Transform and c̄-Transform

这一节要解决的问题：对偶问题有两个未知函数 f、g，还被"对所有 (x, y) 都成立"的约束捆在一起，不好直接优化。c-变换告诉我们：只要知道其中一个，另一个的最佳选择可以直接写出来。于是两个未知函数减成一个，约束也消失了——这是整章的出发点。

回顾：对偶问题 The Dual Problem (2.24)

回忆第 2 章的对偶问题 (2.24)：

原书把约束用指示函数 $\iota_{\mathcal R(c)}(f,g)$（满足约束时为 0，否则为 $-\infty$，记号见 (4.50)）写进目标，记整个目标为 $\mathcal E(f,g)$。借用 Rem. 2.21 的"物流外包"比喻：f(x) 是在 x 收货的单价，g(y) 是在 y 送货的单价，约束说任何一对 (x, y) 的"收 + 送"都不能贵过真实运费 c(x, y)；在这个前提下，报价方希望总收入越高越好。

定义：c-变换与 c̄-变换 Definition (5.1)–(5.2)

一句话直觉：收货价 f 定好之后，y 处的送货价最多能定多高？必须对每一个可能的来源 x 都满足 $g(y)\le c(x,y)-f(x)$，所以最高只能取这些上限中最小的那个。

注意：部分最大化只在 $\beta$（或 $\alpha$）的支撑上决定了最优值，而定义 (5.1) 却给出了定义在整个 $\mathcal Y$（或 $\mathcal X$）上的函数——这是把对偶解自然地延拓到全空间的一种方式。

几何图像：一排"碗"的下包络 Geometric Picture

当 $\mathcal X=\mathbb R^d$、$c(x,y)=\|x-y\|^p$ 时，$f^c$ 叫做 $-f$ 与 $\|\cdot\|^p$ 的下确界卷积(inf-convolution)，这个公式也常被称为 Hopf–Lax 公式(Hopf–Lax formula)（它给出 Hamilton–Jacobi 方程的解）。下面的演示中，f 只在 m 个点上有值（离散的 $\alpha$），包络由 m 只碗构成——这正是下一节半离散情形要用到的；如果 f 是连续函数，就是无穷多只碗、碗底沿曲线 $-f$ 滑动。

初始时有一只碗（第 2 个站点）整个悬在包络上方——它的单元是空的，完全"没生意"。勾选"c̄c 往返"：每只碗被往下压，直到刚好碰到包络（空心圆是新碗底）；包络本身一点没变，但对偶目标 E 变大了。点"用 g^c̄c 替换 g"后再做一次往返：什么都不动了——这就是"只能改进一次"。再试试 p = 1/2，并把第 3 个碗底往下拖到 −0.3 左右：碗变成尖顶的"钉子"，第 3 个站点的单元被别的站点切成好几段不相连的区间。

c-凹函数与 Legendre 变换 c-Concavity

映射 $(f,g)\mapsto(g^{\bar c},f^c)$ 把一对对偶势替换成"更好"的一对（对偶目标 $\mathcal E$ 不降）。能写成 $f^c$ 或 $g^{\bar c}$ 形式的函数称为 c-凹函数(c-concave function) 或 c̄-凹函数。特例 $c(x,y)=\langle x,y\rangle$、$\mathcal X=\mathcal Y=\mathbb R^d$ 时，$f^c(y)=\inf_x\langle x,y\rangle-f(x)$ 是一族线性函数的下包络，这时 c-凹就是通常意义下的凹函数；而且 $f^c(y)=-(-f)^*(-y)$，其中 $\varphi^*(y)=\sup_x\langle x,y\rangle-\varphi(x)$ 是 Legendre 变换——除了几处正负号，c-变换就是 Legendre 变换。所以说 c-变换是 Legendre 变换的推广。

只能改进一次 One Improvement Only

既然 c-变换能改进对偶势，能不能反复做 $f\to f^c\to f^{c\bar c}\to\cdots$ 一路改进到最优？不能。把 Prop. 3.1（矩阵情形）自然推广到连续情形，有

这条性质说明：用这种方式"改进"对偶势只能改进一次。换句话说，对 f、g 做交替最大化不会收敛到最优——它马上就进入循环。这是约束（不光滑的耦合）把两组变量捆在一起时的典型现象。与之形成鲜明对比的是第 4 章的熵正则：正则化后对偶目标 (4.30) 是光滑的，交替最大化恰好就是 Sinkhorn 迭代 (4.43)(4.44)，而且会收敛。把 Sinkhorn 写成对偶变量的形式，它做的正是熵光滑版的 c-变换：把 min 换成 soft-min（见 5.3 节）。

只剩一个势函数：半对偶形式 Single-Potential Formulation (5.4)–(5.5)

利用 (5.3)，可以把 (2.24) 改写成只关于一个势函数的无约束凸规划（对 f 是凹最大化）：

约束去哪儿了？它被 c-变换"自动满足"了：$(f,f^c)$ 永远是可行的一对。由于映射 $(f,g)\mapsto(g^{\bar c},f^c)$ 可以迭代，还可以额外要求 f 是 c̄-凹的、g 是 c-凹的而不损失最优值。这一点在理论上很重要：c-凹函数自带正则性（例如 $c=d^p$ 时是 Lipschitz 的），借此可以证明对偶问题 (2.24) 的解存在（第 2 章 §2.5 提到的"c-变换的机制"就是指这个）。

5.2 半离散形式 Semidiscrete Formulation

这一节要解决的问题：当 $\beta=\sum_j\mathbf b_j\delta_{y_j}$ 只有 m 个点时，(5.5) 里的未知函数 g 其实只是 m 个数。于是一个关于连续分布的无穷维问题，变成了关于向量 $\mathbf g\in\mathbb R^m$ 的有限维凹最大化。本节讲清它的几何（Laguerre 单元）和解法（梯度上升、牛顿法）。（若 $\alpha$ 是离散的而 $\beta$ 连续，交换两者的角色即可。）

离散的 c̄-变换 Discrete c̄-Transform (5.6)

这正是上一节演示里的"m 只碗的下包络"。原书 Fig. 5.1 上半部分画的就是一维的例子（代价 $|x-y|^p$，p 从 1/2 到 2，红点在 $(y_j,-\mathbf g_j)$），下半部分是二维的例子：背景颜色是 $\mathbf g^{\bar c}$ 的数值，粗黑线围出的区域就是下面要定义的 Laguerre 单元。下面的演示 2 可以复现它。

有限维的对偶问题 A Finite-Dimensional Dual (5.7)

关键的一步：把离散 c̄-变换代入半对偶 (5.5)，得到一个只有 m 个未知数的优化问题：

$\mathcal E$ 是凹函数：对每个固定的 x，$\mathbf g\mapsto\mathbf g^{\bar c}(x)$ 是 m 个仿射函数的最小值（凹），对 x 积分、再加上线性项，仍然是凹的。所以这是一个没有约束、没有局部极大陷阱的问题。

Laguerre 单元 Laguerre Cells

一句话直觉：$\mathbf g^{\bar c}(x)$ 里的 min 在哪个 j 处取到，x 就"归属"于哪个站点。把归属于同一站点的点收集起来，就是它的"势力范围"。

有了单元，目标函数可以方便地改写为（在每个单元内，min 就是那一项）：

$\mathbf g$ 全为 0 时是 Voronoi 图：边界是相邻站点连线的垂直平分线。提高某个 $\mathbf g_j$，它的单元就向外扩张、边界平行地推向邻居（p = 2 时边界始终是直线，单元是凸多边形）；权重足够大时，一个站点甚至能"吞掉"别的站点所在的位置——站点不一定在自己的单元里。给所有 $\mathbf g_j$ 同时加同一个数，单元不变。换成 p = 1 或 1/2，边界变成曲线。打开等高线可以看到 $\mathbf g^{\bar c}$ 的"碗"形地形，单元边界正是两只碗的交线。

梯度：目标质量 − 单元质量 Gradient of E

这个公式非常好读：

• 最优条件：$\nabla\mathcal E(\mathbf g)=0$ ⟺ 每个单元的质量恰好等于目标 $\alpha(\mathbb L_j(\mathbf g))=\mathbf b_j$。
• 梯度上升 $\mathbf g\leftarrow\mathbf g+\tau\,(\mathbf b-\alpha(\mathbb L(\mathbf g)))$：客流不足的站点加补贴（单元扩张），爆满的站点减补贴（单元收缩）——就像市场上的价格调节。
• 梯度各分量之和为 $\sum_j\mathbf b_j-\sum_j\alpha(\mathbb L_j)=1-1=0$，所以梯度上升自动保持 $\sum_j\mathbf g_j$ 不变。

最优传输映射是分片常数的 Piecewise Constant Transport Map

求出最优的 $\mathbf g$ 后，从 $\alpha$ 到 $\beta$ 的最优传输映射 T 就是：把 $\mathbb L_j(\mathbf g)$ 里的每个点都送到 $y_j$。它是分片常数的——整个单元被"压缩"到一个点上。对 $c=\|x-y\|^2$，这与第 2 章的 Brenier 定理（Rem. 2.24）一致：$T=\nabla\varphi$，其中 $\varphi(x)=\max_j\,\langle x,y_j\rangle+\tfrac12(\mathbf g_j-\|y_j\|^2)$ 是一个分片线性的凸函数，它在 $\mathbb L_j$ 上的梯度恰好是 $y_j$。

梯度上升：单元像"呼吸"一样一点点调整，误差曲线在对数坐标下近似一条直线（线性收敛）；均匀分布下把 τ 调到 0.5 以上就会来回震荡（步长超过了 2/曲率）。阻尼牛顿通常 3–6 步就把误差压到 $10^{-10}$ 以下（曲线几乎垂直下落）。换成高斯团块的 $\alpha$：密度高的地方边界稍一移动，质量就变化很多（曲率大），τ = 0.3 在均匀分布下能收敛，在这里却会来回震荡；而离团块远的站点必须把单元伸得很长才能"够到"足够的质量——曲率忽大忽小，正是二阶方法的用武之地。勾选"显示传输映射"：每个单元的质量（箭头从单元的 α-重心出发）整体运往自己的站点。

幂图与凸包 Power Diagrams and Convex Hulls

当 $c(x,y)=\|x-y\|^2$ 时，Laguerre 剖分也叫幂图(power diagram)。单元 j 与 k 的分界满足 $\|x-y_j\|^2-\mathbf g_j=\|x-y_k\|^2-\mathbf g_k$，展开后二次项抵消：

这是一张垂直于连线 $y_jy_k$ 的超平面（二维时是直线）。$\mathbf g_j$ 每增加 δ，这条边界就朝 $y_k$ 平移 $\delta/(2\|y_j-y_k\|)$。所以每个单元都是若干半空间的交——凸多面体，可以用计算几何算法高效求出（见 [Aurenhammer, 1987]）。

最常用的算法基于一个"抬升"技巧：因为 $\|x-y_j\|^2-\mathbf g_j=\|x\|^2-h_j(x)$，其中 $h_j(x)=2\langle x,y_j\rangle-\|y_j\|^2+\mathbf g_j$ 是 x 的仿射函数，所以"碗的下包络"减去公共的 $\|x\|^2$ 之后，就变成了"一族超平面的上包络"。上包络是 $\mathbb R^{d+1}$ 中的凸多面体曲面，它的每个面竖直投影下来恰好是一个 Laguerre 单元。这个上包络与抬升点集 $\big(y_j,\ \|y_j\|^2-\mathbf g_j\big)_{j=1}^m\subset\mathbb R^{d+1}$ 的（下）凸包互为对偶，所以幂图可以通过计算凸包得到。求凸包的算法很多，例如 Chan [1996] 的算法在二维和三维的复杂度是 $O(m\log Q)$，Q 是凸包的顶点数。

5.3 熵正则的半离散形式 Entropic Semidiscrete Formulation

这一节要解决的问题：把第 4 章的熵正则搬到半离散问题上。效果一句话：c-变换里的 min 换成 soft-min。于是 Laguerre 单元的"硬边界"变成渐变的"软边界"，目标函数变得处处光滑——可以放心地用牛顿法、L-BFGS，下一节还能用随机梯度。

熵正则的对偶与软 c-变换 Smoothed c-Transform (5.9)

一般测度之间熵正则问题 (4.9) 的对偶，是一个光滑、无约束的优化问题（参见 (4.45)）：

其中 $(f\oplus g)(x,y)\overset{\text{def.}}{=}f(x)+g(y)$。和 (2.24) 对比：硬约束 $f\oplus g\le c$ 被换成了一个指数罚项——违反约束不再是"不可能"，而是代价随违反量指数增长。和不带正则时 (5.1) 一样，可以对 f 或 g 之一显式地求最大，得到熵光滑的 c-变换：

Sinkhorn = 交替计算软 c-变换 Sinkhorn as Alternate Smoothed c-Transforms (5.11)

当 $\alpha$ 也是离散的，可以类似地定义 $\mathbf f$ 的软变换。第 4 章用 soft-min 改写的 Sinkhorn 迭代 (4.40)(4.41)，恰好就是交替地计算熵光滑的 c-变换：

把它和 5.1 节放在一起看，就是本章最重要的对比之一：$\varepsilon=0$ 时交替做 c-变换往返一次就停住；$\varepsilon\gt0$ 时交替做软 c-变换会一步步收敛到最优——这就是 Sinkhorn。下面的演示把两者放在一起跑。

红线（硬）在第 1 步跳一下就再也不动了——这就是 $f^{c\bar cc}=f^c$：往返一次之后势函数已经是 c-凹的，交替最大化原地打转，通常停在精确值的一半到九成之间（取决于点的位置）。蓝线（软）则一步步爬升：ε 越小，最终越接近 100%，但爬得越慢；ε 太大时 soft-min 太"软"，得到的 $\mathbf f$ 离最优较远。左图的连线是 Sinkhorn 在当前迭代的耦合 $\mathbf P$：一开始很模糊，逐渐集中到少数几条，接近精确的最优耦合。

光滑的半离散目标 Smooth Semidual (5.12)–(5.13)

与其求解 (5.9)，不如直接求解有限维问题

它在不带正则的情形 $\varepsilon=0$ 下依然有效：此时 $\mathbf g^{\bar c,0}=\mathbf g^{\bar c}$ 就是 (5.6) 的 c̄-变换，(5.12) 就是 (5.8)。它的梯度是

$\chi^\varepsilon_j$ 是 Laguerre 单元指示函数 $\chi^0_j=\mathbb 1_{\mathbb L_j(\mathbf g)}$ 的光滑版本（(5.13) 在 $\varepsilon=0$ 时仍然成立，就退回 5.2 节的梯度）。函数族 $(\chi^\varepsilon_j)_j$ 是一个单位分解(partition of unity)：$\chi^\varepsilon_j\ge0$ 且 $\sum_j\chi^\varepsilon_j=1$。换句话说，点 x 不再"整个"属于某个站点，而是按比例 $\chi^\varepsilon_j(x)$ 分给各个站点——从"硬分配"变成了"软分配"。

原书 Fig. 5.3 展示了这一切：上图是一维的软变换 $\mathbf g^{\bar c,\varepsilon}$（代价 $|x-y|$，ε 从 0 到 0.3），下图是二维的情形（代价 $\|x-y\|$），颜色表示软分配 $\chi^\varepsilon_j$，黑色曲线是 $\mathbf g^{\bar c,\varepsilon}$ 的等高线。下面两个演示分别对应上、下两行。

ε 很小时软包络几乎贴着硬包络，彩色带是一块块"纯色"——就是硬的 Laguerre 区间。ε 增大，包络在两碗交界处的尖角被磨圆、向下沉，而碗底附近变化较小（ε 很大时整条曲线都会明显下沉）；彩色带在交界处出现渐变：交界附近的点被按比例分给两个站点。无论 ε 多大，每一列彩色带的总高度都是 1（单位分解）。

ε = 0 时是纯色的硬单元；ε 增大，边界晕开成宽度约与 ε 成正比的渐变带，但"最可能的站点"的分界线（虚线，argmax 边界）依然是直线——软分配就是一个线性的多类 logistic 回归分类器（见下面的 Rem. 5.1）。取消"保持最优 g"再拖 ε：g 不变时，ε 越大各单元的软质量越偏离目标，说明最优的 g 本身随 ε 而变。

二阶方法与 logistic 回归 Second-Order Methods and Logistic Regression

5.4 随机优化方法 Stochastic Optimization Methods

这一节要解决的问题：在高维空间里，没法把 Laguerre 单元精确地切出来再对它们积分。但半离散目标 (5.8) 和它的光滑版本 (5.12) 有一个好性质：它们都是关于 $\alpha$ 的期望。所以只要能从 $\alpha$ 里抽样本，就能用随机梯度法求解——完全不需要把 $\alpha$ 离散化。

目标是一个期望 The Objective as an Expectation

其中 X 是按 $\alpha$ 分布的随机向量。对单个样本 x 求梯度：

读法："目标比例 − 这个样本的分配比例"。$\varepsilon=0$ 时 $\chi^0(x)$ 是一个独热向量 $\mathbf e_{j^*}$（x 落在哪个单元，哪一位就是 1），于是随机梯度就是 $\mathbf b-\mathbf e_{j^*}$。它的期望恰好是 (5.13) 的真梯度 $\mathbf b-\int\chi^\varepsilon\mathrm d\alpha$——一个无偏估计。

这样做的好处（Genevay 等 [2016] 提出）：可以得到有收敛性保证的算法，而不必对 $\alpha$ 做任何人为的离散化（既不用 Dirac 和来近似 $\alpha$，也不用求积公式算积分）。$\alpha$ 只被当作一个能产生独立样本的黑盒——在统计与机器学习的高维应用中，这是最自然的计算设定（比如 $\alpha$ 是一个庞大的数据集，或者一个生成模型）。这类方法已被 Staib 等 [2017b] 推广到 Wasserstein 重心的计算（§9.2）。

随机梯度法 Stochastic Gradient Descent (5.14)–(5.15)

从 $\mathbf g^{(0)}=\mathbf 0$ 出发，随机梯度法(stochastic gradient descent, SGD；这里用作最大化方法)在第 ℓ 步按 $\alpha$ 独立地抽取一个样本 $x_\ell$（与之前、之后的样本都独立），然后更新

步长 $\tau_\ell$ 必须足够快地衰减到 0，才能抵消"用样本梯度代替真梯度"带来的噪声。典型的选择是

其中 $\ell_0$ 大致表示"热身阶段"（warmup）的迭代次数：前 $\ell_0$ 步步长基本不变，之后按 $1/\ell$ 衰减。可以证明

其中 $\mathbf g^\star$ 是 (5.12) 的解，期望是对独立抽样 $(x_\ell)_\ell$ 取的。原书 Fig. 5.4 展示了一个简单的二维例子：$\alpha$ 是 $[0,1]^2$ 上的均匀分布、$\varepsilon=0$，几条不同颜色的曲线是不同随机样本序列下的能量间隙，下方的小图是 Laguerre 单元随迭代的演变。

带平均的随机梯度 SGD with Averaging (SGA)

SGD 慢，是因为步长 $\tau_\ell$ 衰减得太快。改进办法是对过去的迭代做平均：等价于在辅助变量 $\tilde{\mathbf g}^{(\ell)}$ 上跑一个"经典"的 SGD，

并以平均值 $\mathbf g^{(\ell)}$ 作为输出。这就是带平均的随机梯度法(stochastic gradient descent with averaging, SGA)。不必存下所有迭代，只要维护一个滑动平均：

注意这里的步长衰减慢得多，按 $\ell^{-1/2}$ 衰减。Bach [2014] 证明 SGA 比 SGD 收敛更快（涉及的常数更小），因为与 SGD 不同，SGA 能自适应目标函数局部的强凸性（对最大化问题是强凹性）。

开始时步长较大，样本噪声甚至会让间隙先上升一段；随后 SGD（红）的步长按 $1/\ell$ 迅速变小，下降明显放缓。SGA 的原始迭代（浅蓝虚线）步长衰减得慢，一直在最优点附近抖动、间隙降不下去；但它的平均（蓝）又稳又快，默认参数下通常比 SGD 低 5–10 倍。每条曲线只是一次随机运行，所以锯齿很多；多按几次"重置"可以比较不同的样本序列（灰线）。灰色虚线只标出 $1/\sqrt\ell$ 的斜率（高度随意取）：理论只保证期望间隙至少按这个速率下降，实际曲线后期往往下降得更快。再试试调参数：τ₀ = 0.02、ℓ₀ = 10 时 SGD 的步长过早变小，g 被"冻结"在离最优较远处；τ₀ = 0.3、ℓ₀ = 1000 时它长期迈大步、噪声很大。SGA 对这两个参数则稳健得多——这正是"自适应"的好处。

两边都连续怎么办 Continuous–Continuous Problems

本章小结

术语卡片 Glossary Cards

自测

习题

本章导读：当代价就是距离 Introduction

这一节要回答：为什么要专门拿出一章讨论「代价 = 距离」这种情形？它和前几章常用的平方代价有什么本质区别？

本章只讨论一种地面代价(ground cost)：两点之间搬一单位质量的代价恰好等于一个距离 $d(x,y)$，而不是它的平方。对应的最优传输值记作 $\mathcal{W}_1(\alpha,\beta)$，它就是第 2 章 $p$-Wasserstein 距离 (2.17) 在 $p=1$ 时的特例：

也就是「每一单位质量付的钱 = 它走过的路程」，总账 = 质量 × 距离之和。

和平方代价（第 7 章会详细研究）相比，距离代价有得有失。下表把原书导言里的要点并排放在一起：

先用一个经典的小例子直观感受第一行——为什么 W₁ 的最优方案常常不唯一。

绿色的 W₁ 代价是一条水平线：「整体平移」「只搬一本」以及两者之间的任何混合都一样便宜（代价都是 1），最优方案有无穷多个。紫色的平方代价却随 $t$ 上升，只有 $t=0$ 的平移最优。书越多，「只搬一本」在平方代价下越贵（代价 $=n$），在 W₁ 下却始终是 1。

6.1 度量空间上的 W₁ W1 on Metric Spaces

这一节要回答：当代价 $c=d$ 是一个距离时，Kantorovich 对偶问题会变成什么样？答案出奇地干净：两个对偶势合并成一个 $f$，唯一的约束是「$f$ 的坡度不超过 1」。由此还能看出 W₁ 是一个范数，并在离散与一维情形得到非常简单的线性规划。

本节始终假设 $d$ 是空间 $\mathcal{X}=\mathcal{Y}$ 上的一个距离，求解地面代价为 $c(x,y)=d(x,y)$ 的最优传输问题。

Lipschitz 常数 Lipschitz Constant

命题 6.1：c-变换只留下 1-Lipschitz 函数 c-Transform When the Cost Is a Distance

回忆第 5 章的 c-变换(c-transform) (5.1)：$f^c(y)=\inf_{x\in\mathcal{X}}\,c(x,y)-f(x)$。在对偶问题里，固定 $f$ 之后最优的第二个势就是 $f^c$，于是对偶问题可以只用一个势来写（单势形式 (5.4)）。当代价是距离时，c-变换有一个特别漂亮的刻画：

选「平缓」时，所有帐篷都藏在 $f$ 下面，绿线与蓝线重合，下半部红线与蓝色点线重合：$f^c=-f$。选「有一段陡坡」，陡坡顶上的帐篷探出来，把左侧垫成坡度 1 的斜面（绿色阴影），$f^c$ 也随之低于 $-f$。点「把 $f$ 换成包络」后再看：包络已经是 1-Lipschitz 的，再做 c-变换就不再变化——这就是「只需改进一次」。

这个命题让对偶问题再简化一步。从单势形式 (5.4) 出发，可以把势对 $(g,g^c)$ 再「改进」一次换成 $(g^c,(g^c)^c)$。命题 6.1 说 $g^c$ 是 1-Lipschitz 的，于是 $(g^c)^c=-g^c$：势对变成 $(g^c,-g^c)$——也就等价于任意一对 $(f,-f)$，其中 $\operatorname{Lip}(f)\le1$。代入对偶目标 $\int f\,\mathrm{d}\alpha+\int g\,\mathrm{d}\beta$，就得到下面的 W₁ 公式。

Kantorovich–Rubinstein 对偶 Kantorovich–Rubinstein Duality

把势对 $(f,-f)$ 代入对偶问题，得到 W₁ 的另一种表达——Kantorovich–Rubinstein 对偶(Kantorovich–Rubinstein duality)：

一维离散情形可以亲手玩一玩这个「最大化游戏」。下面的点 $z_1\lt\dots\lt z_8$ 上，差 $\alpha-\beta$ 的质量为 $m_k$（正 = $\alpha$ 多出来，负 = $\beta$ 多出来），你的任务是选势 $f_k$ 让 $\sum_k f_km_k$ 最大——这正是 6.1 节末尾的离散公式 (6.2)。

想让 $\sum f_km_k$ 变大，就要把蓝柱处的 $f$ 抬高、红柱处压低，但坡度卡在 1。点「显示最优 $f$」：凡是有净质量越过的间隙（有绿色箭头），$f$ 的斜率都恰好是 $\mp1$（橙色），而且沿箭头方向下降；没有质量越过的间隙，斜率可以随意——所以最优 $f$ 一般不唯一，但最优值 W₁ 唯一，并与一维公式、Book.ot.emd1d 的结果一致。

W₁ 是一个范数 W1 Is a Norm

仔细看 (6.1) 的右边：$\alpha$ 与 $\beta$ 只以差 $\alpha-\beta$ 的形式出现。因此

并且这个表达式对任何不必为正的测度（带符号测度(signed measure)）都有效，只要 $\int_{\mathcal{X}}\alpha=\int_{\mathcal{X}}\beta$。这个范数常被称为 Kantorovich–Rubinstein 范数(Kantorovich–Rubinstein norm) [1958]。几个立刻可以读出来的推论：

• 共同质量不影响距离：对任意 $\gamma$，$\mathcal{W}_1(\alpha+\gamma,\beta+\gamma)=\mathcal{W}_1(\alpha,\beta)$。
• 只有「多出来的部分」需要搬：$\mathcal{W}_1(\alpha,\beta)=\mathcal{W}_1\big((\alpha-\beta)_+,(\alpha-\beta)_-\big)$，其中 $(\alpha-\beta)_+$、$(\alpha-\beta)_-$ 是差的正部与负部——两边共有的质量原地不动即可。
• 齐次与三角不等式：$\mathcal{W}_1(\lambda\alpha,\lambda\beta)=\lambda\,\mathcal{W}_1(\alpha,\beta)$（$\lambda\ge0$）；三角不等式直接来自范数的三角不等式。

无论把 $\gamma$ 放在哪里、放多少，W₁ 纹丝不动——最下排的差 $\alpha-\beta$ 也纹丝不动。平方代价却会变小：$\gamma$ 成了「中转站」，一次长途被拆成两段短途（$2^2=4\gt1^2+1^2$）。切换到「直接搬」：它对 W₁ 仍然最优（代价相同），对平方代价则不是。

离散测度与一维情形 Discrete Measures and the 1-D Case

对离散测度（第 2 章 (2.1) 那样的 Dirac 之和），把差写成 $\alpha-\beta=\sum_k\mathbf{m}_k\delta_{z_k}$，其中 $z_k\in\mathcal{X}$、$\sum_k\mathbf{m}_k=0$（$\mathbf{m}_k\gt0$ 表示 $\alpha$ 在该点多出来，$\lt0$ 表示 $\beta$ 多出来）。优化 (6.1) 就变成有限维问题：

原书称它为「带二次锥约束的有限维凸规划」，可以用内点法求解，或像后文处理类似问题那样用近端方法(proximal methods)。（注：(6.2) 的每个约束 $|\mathbf{f}_k-\mathbf{f}_\ell|\le d(z_k,z_\ell)$ 其实就是两个线性不等式；二阶锥约束 $\|\nabla f\|_2\le1$ 是在欧氏情形 (6.3) 离散化之后才出现的。）

(6.2) 有 $K^2$ 个约束。但在实轴上（$\mathcal{X}=\mathbb{R}$，$d(x,y)=|x-y|$），把点排好序 $z_1\le z_2\le\cdots$，只需保留相邻点之间的约束：

这是一个线性规划。原因很简单：相邻约束一段段串起来，由三角不等式自动推出任意两点间的约束 $|\mathbf{f}_k-\mathbf{f}_\ell|\le\sum_{j=k}^{\ell-1}(z_{j+1}-z_j)=z_\ell-z_k$。（6.3 节在图上还会用同一个技巧：只需要「边」上的约束。）

一维还有闭式解，就是第 2 章的 (2.37)：记 $C_\alpha(x)=\int_{-\infty}^x\mathrm{d}\alpha$ 为累积分布函数(cumulative distribution function)，则

阴影面积与 Book.ot.emd1d（单调重排，$p=1$）算出的最优代价始终相等。选「整体平移」：两条曲线形状相同、错开 $\Delta$，面积 = $\Delta$，就是平移距离。注意箭头：在蓝线高于红线的地方质量向右流，反之向左——同一个一维问题里，质量可以在不同位置朝不同方向流。

注记 6.1：$0\lt p\le1$ 的 W_p Wp with 0 < p ≤ 1

6.2 欧氏空间上的 W₁ W1 on Euclidean Spaces

这一节要解决：在 $\mathbb{R}^d$ 中，「任意两点之间坡度 $\le1$」要对无穷多对点逐一检查，能不能换成逐点的条件？可以——换成「每一点的梯度长度 $\le1$」。再对它取对偶，就得到一个关于流场的问题：Beckmann 形式。它把「谁搬给谁」这种全局配对问题，变成了「每个地方流量多大、往哪流」这种局部问题。

局部化：梯度不超过 1 Gradient Constraint

在欧氏空间 $\mathcal{X}=\mathcal{Y}=\mathbb{R}^d$ 上取 $c(x,y)=\|x-y\|$，(6.1) 中全局的 Lipschitz 约束可以改写成对 $f$ 的梯度的一致界：

这里 $\|\nabla f\|_\infty\le1$ 的意思是：在每一点 $x$，梯度的欧氏长度 $\|\nabla f(x)\|_2\le1$——站在任何地方，最陡方向上的坡度也不超过 1。

Beckmann 形式：流与散度 Beckmann Formulation

对 (6.3) 取对偶，得到一个在固定散度约束下的优化问题：

它常被称为 Beckmann 形式(Beckmann formulation) [Beckmann, 1952]。

绿色箭头总是垂直穿过等高线，而且沿箭头方向每走 0.5 就恰好跨过一条等高线：流沿最陡的下坡走，坡度恰为 1（互补松弛）。把检测圈只套住一个蓝点：净流出 = 1/4；套住一蓝一红：净流出 = 0；放在空白处：流进多少就流出多少，净通量为 0——这就是 $\operatorname{div}(s)=\alpha-\beta$。读数里原始值 $\sum P\|x-y\|$ 与对偶值始终相等。（这里画的势是 $f=g^c$ 形式，恰好处处坡度为 1；约束只要求 $\le1$，真正必须「绷紧」的只是有流经过的地方。）

数值求解与近似 Numerical Methods and Approximations

用有限元等方法恰当离散之后，(6.3) 与 (6.4) 都是非光滑凸优化问题。可以直接调用现成的内点法二次锥规划求解器，但正如 §7.3 所提倡的，大规模问题更适合简单而贴合结构的一阶方法，例如 Douglas–Rachford (DR) 迭代 (7.14) 或与之密切相关的交替方向乘子法（ADMM）。在均匀网格上，投影到散度约束 $\{\operatorname{div}(s)=\alpha-\beta\}$ 可以借助快速傅里叶变换（FFT）高效完成。Solomon 等 [2014a] 详细介绍了这些方法及其在三角网格上的应用；Li 等 [2018a]、Ryu 等 [2017b,a] 用原始-对偶分裂得到类似的方案。

另一条路是放松对 $f$ 的 Lipschitz 约束以得到近似：例如约束 $f$ 的小波系数，可以得到显式公式 [Shirdhonkar & Jacobs, 2008]；或者只考虑用「整流线性」$\max(0,\cdot)$ 激活、并截断（clip）权重的多层神经网络来参数化 $f$ [Arjovsky et al., 2017]。

6.3 图上的 W₁ W1 on a Graph

这一节要解决：当空间本身是一张图（道路网、像素网格、三角网格），距离是最短路长度时，能不能不去构造 $n\times n$ 的耦合矩阵，只在边上「局部地」计算 W₁？可以——把梯度与散度搬到图上，(6.3)(6.4) 就变成一对线性规划，其中 (6.6) 正是运筹学里经典的最小费用流问题。

测地距离：图上的最短路 Geodesic Distance on a Graph

公式 (6.3)(6.4) 可以推广到 $\mathcal{X}$ 是测地空间(geodesic space)的情形，即 $c(x,y)=d(x,y)$ 是测地距离（两点之间最短曲线的长度）；$\mathcal{X}$ 为黎曼流形时的理论分析见 Feldman & McCann [2002]。当 $\mathcal{X}=\llbracket1,n\rrbracket$ 是一个离散集合，配有无向边 $(i,j)\in\mathcal{E}\subset\mathcal{X}^2$、每条边标有权重（长度）$\mathbf{w}_{i,j}$ 时，就得到一个重要情形：$\mathcal{X}$ 是一张图，配备测地距离（最短路距离）(geodesic distance / shortest path metric)：

其中 $i\to j$ 表示 $i_1=i$、$i_K=j$，即路径从 $i$ 出发、到 $j$ 结束。

考虑两个向量 $(\mathbf{a},\mathbf{b})\in(\mathbb{R}^n)^2$，它们定义图上的（带符号的）离散测度，满足 $\sum_i\mathbf{a}_i=\sum_i\mathbf{b}_i$（这些权重不必为正）。目标是计算以最短路距离为地面度量的 $\mathrm{W}_1(\mathbf{a},\mathbf{b})$（即 (2.17) 取 $p=1$、$C=\mathbf{D}$）。关键在于：不要去算完整的 $n\times n$ 耦合 $\mathbf{P}$，而是借助图的连通结构，只用局部算子——梯度与散度这两个微分算子的离散版本。

图上的梯度与散度 Gradient and Divergence on a Graph

离散的对偶 Kantorovich 势 $\mathbf{f}\in\mathbb{R}^n$ 是定义在所有顶点上的向量。梯度算子(gradient operator) $\nabla:\mathbb{R}^n\to\mathbb{R}^{\mathcal{E}}$ 定义为

流(flow) $\mathbf{s}=(\mathbf{s}_{i,j})_{i,j}$ 定义在边上：每条无向边看成方向相反的两条有向边，各带一个非负流量，$\mathbf{s}_{i,j}$ 是从 $i$ 流向 $j$ 的量。散度算子(divergence operator) $\operatorname{div}:\mathbb{R}^{\mathcal{E}}\to\mathbb{R}^n$ 是梯度的伴随，它把流映成顶点上的向量：

最小费用流 Min-Cost Flow

有了这两个局部算子，(6.3) 在图上变成

与之相伴的对偶问题（对应 (6.4)）是

• (6.5) 只要边约束就够了：和一维「只需相邻约束」一样，沿最短路把边上的约束一段段串起来，就得到任意两点之间 $|\mathbf{f}_i-\mathbf{f}_j|\le\mathbf{D}_{i,j}$——恰是 (6.2) 以 $\mathbf{D}$ 为距离时的全部 $n^2$ 个约束。
• (6.6) 是图上的 Beckmann 形式：找一个非负流，使每个点「净流出 = 供给 − 需求」，并让「流量 × 边长」的总和最小。它是一个线性规划，更确切地说是最小费用流(min-cost flow)问题的一个实例。针对它已有非常高效的专用单纯形求解器（例如 [Ling & Okada, 2007]）。(6.6) 正是所谓的 Beckmann 形式 [Beckmann, 1952]，它还被推广用来定义和研究交通拥堵模型（例如 [Carlier et al., 2008]）。

四种视角，同一个小例子 Four Views of One Example

下面用一张 9 个点的小图，把本章的各种视角放在一起对照：传输（谁搬给谁）、流（每条路上走多少）、势（每个点的海拔），以及把它们连起来的互补松弛。

看懂了这个例子，就可以自己动手改：换供需、改边长，观察流如何绕路、势如何变化，以及结果如何与「构造完整耦合」的做法对上。

三个数永远相等：图上的最小费用流 $\sum\mathbf{w}\mathbf{s}$（只用边变量）、对偶值 $\sum\mathbf{f}(\mathbf{a}-\mathbf{b})$、以及 Book.ot.emd 在最短路距离矩阵 $\mathbf{D}$ 上解出的完整耦合的代价。把 D–G 这座「桥」的长度调大，流会改走 D–F–H 或 E–G 等路线；把 A、H 各加一份再清空，结果只取决于净供给 $\mathbf{a}-\mathbf{b}$——W₁ 是范数。

原书图 6.1 Figure 6.1