浅谈二分图带权匹配——从 OI 到计算机图形学
前置知识：二分图匹配，费用流，模拟费用流，Primal-Dual，KM 算法
相关问题：差分约束，线性规划，Voronoi Diagram，Power Diagram，Optimal Transport
0. 前言
本文的起点是我在科研过程中遇到的一些算法问题。经过一些调研和思考之后，我竟然对自己在遥远的 OI 在役时期没有掌握好的一些算法有了新的理解。我感觉这些想法还算有趣，因此写下这篇文章与大家分享。
想不到在退役多年以后，我竟然还能发一篇 UOJ 博客，来分享一些我对 OI 问题的理解和思考；本文也是我的第一篇 UOJ 博客。
由于过去我从来都不是一位站在“科技前沿”的选手，退役之后更是好多年没有接触过 OI 界的新进展了，文章中难免有一些过时或理解不当的内容，欢迎大家批评指正，或在评论区友好交流。
1. 从原始问题说起
数据集中有 \(n\) 个散点与 \(m\) 个站点，每个站点 \(s_i\) 有一个容量限制 \(C_i\)（我们不妨设 \(\sum C_i=n\)）；将散点 \(p_i\) 分配给站点 \(s_i\) 的代价为 \(f(i,j)\)。现需要求出一组分配方案，最小化总代价；我们不妨设 \(n\gg m\)。
显然这是一个二分图带权多重匹配问题；熟悉 OI 相关算法的朋友们一定能在一瞬间想到一个费用流建模：散点在左，站点在右；散点与源点 \(S\) 连容量为 \(1\) 的边，站点与汇点 \(T\) 连容量为 \(C_j\) 的边，散点、站点之间连容量为 \(1\)、费用为 \(f(i,j)\) 的边。跑费用流。
这样做的复杂度太高了，我们不能接受；优化策略是从增广路的形态入手。
众所周知，匹配问题中的增广路形态是 \(S\rightarrow p_1\rightarrow s_1\rightarrow p_2\rightarrow s_2\rightarrow\cdots\rightarrow p_k\rightarrow s_k\rightarrow T\)；它的实际意义是，我们选出一个未分配的点 \(p_1\)，将其分配给站点 \(s_1\)，然后选出一个原本分配给 \(s_1\) 的散点 \(p_2\)，将其重新分配给 \(s_2\)，……，依此类推。
此时有一个性质是，如果我们需要选出一个原本分配给 \(s_x\) 的散点 \(p_i\)，将其重新分配给站点 \(s_y\)，那么这个散点需要“最小化重分配带来的代价变化”。也就是说，我们需要最小化 \(f(i,y)-f(i,x)\)。当 \(s_x,s_y\) 均确定时，实际上 \(p_i\) 的选取也是 固定 的。
于是我们根本没有必要在原始的二分图上做最短路。考虑到 \(m\ll n\)，我们只需要考虑“失配链”依次经过了哪些 站点，确定了这个之后，链上的散点也就随之固定了。
具体地，我们用 \(m\) 个站点来构造一张有向图，图上 \(s_x\) 与 \(s_y\) 之间的边权 \(w(x,y)=\min(f(i,y)-f(i,x))\)，其中 \(i\) 是一个先前被分配给了 \(s_x\) 的散点。
增广过程可以直接在这个图上完成，建一个超级源点跑 Bellman-Ford 即可。
至于重分配代价最小的散点怎么找，我们维护 \(m^2\) 个优先队列即可。这个过程大概就是优先队列 \(Q(x,y)\) 里存的是所有的 \(\Delta(x,y)=f(i,y)-f(i,x)\)。
每次增广之后我们要更新一下这些优先队列，这部分的复杂度显然是 \(O(m^2\log n)\)。
于是我们把这个问题做到了 \(O(n(m^3+m^2\log n))\)。
然而我们不太满足于这个复杂度——众所周知最小费用最大流可以用 Primal Dual 来优化。我们给每个站点赋势能 \(\varphi(i)\)，将增广过程中的 Bellman-ford 修改成 Dijkstra，就得到了 \(O(nm^2\log n)\) 的复杂度，此时的复杂度瓶颈在于 \(Q(x,y)\) 的更新。
空间复杂度 \(O(nm)\)，瓶颈在于优先队列的存储。
2. 一些简单的变式
2.1 当 \(m\) 不再远小于 \(n\)
当 \(m\) 与 \(n\) 同阶时，上述算法的复杂度退化至 \(O(n^3\log n)\)；发现此时 \(Q(x,y)\) 反而成了负优化。
我们直接把优先队列这个优化删掉即可，增广过程中“求重分配代价最小的点”这一步直接暴力完成。
事实上也可以把所有“优化”都去掉，只留一个 Primal-Dual，复杂度不变——于是我们就得到了一种 \(O(n^3)\) 的最大权完美匹配算法。
我们说，其实这种算法与经典的 KM 算法是等价的。KM 算法通过更新“顶标”来增广，此处“顶标”与 Primal-Dual 中用到的势能函数等价。
非常熟悉 Primal-Dual 与 KM 算法的朋友其实不难观察到这一点；下面我再给出一个非常直观的理解方式。
KM 算法的顶标满足如下条件：\(lx(u)+ly(v)\le w(u,v)\)（此处按照 OI-wiki 上的方式，记 \(lx(u)\) 表示左部点 \(u\) 的顶标，\(ly(v)\) 表示右部点 \(v\) 的顶标）。
对于匹配边 \((u,v)\)，我们有：
对于非匹配边 \((u',v)\)，有：
将 \((1)\) 中的 \(ly(v)\) 带入 \((2)\) 中并移项，我们得到：
如果将匹配的左部点视为站点，右侧点视为散点，那么 \((3)\) 式左侧正是我们在 \(Q(x,y)\) 中维护的“转移代价” \(\Delta(x,y)\)，也就是增广过程中构建的压缩图上 \(u'\)、\(u\) 两点的权值差。
于是我们可以立刻意识到，若将此处的顶标 \(l\) 直接作为增广过程中的势能函数，它同样满足 Primal-Dual 的“边权非负”性质！当然，不等号的方向是反的，因为这里我们用 KM 求解的是最大权匹配，而原问题是最小权匹配。
当然这种解释非常简陋，但我们着实可以由此认识到，费用流与 KM 殊途同归，它们本质上是在做同一件事情。
2.2 当多重匹配变成单一匹配
再考虑这样一种变式：所有站点的容量均为 \(1\)（或均为常数级别之类的，下文以容量均为 \(1\) 为例），我们不再需要将所有散点都分配给站点，只需要从 \(n\) 个散点中挑选 \(m\) 个进行匹配；\(m\ll n\) 的前提仍然成立。
此时我们发现一条非常有用的性质——对于每个站点来说，只有权重前 \(m\) 小的边是有用的。这个其实非常显然，此处省略证明。
这启示我们，实际上有用的候选匹配边只有至多 \(m^2\) 条。于是此时我们仍然可以省略优先队列优化，同样做到 \(O(nm^2)\) 的时间复杂度。
本小节的内容在 这篇文章 中有所体现。
2.3 边权带修改
本节中我们考虑这样一种变式：当我们求解出了本文最开始提到的问题，有了一组分配方案与站点势能后，与某个站点 \(s_x\) 相关的所有匹配权重 \(f(i,x)\) 均被修改，此时能否快速求出一组新的解？
不难发现，这个问题需要我们在残量网络上消去负环。于是类似地构建站点压缩图，不断跑 Bellman-ford 找负环即可。然而这样一来就出现一个问题——我们用 Primal-Dual 精心维护的势能不再有效了。
有没有一种办法让我们能在这个过程中仍然能够使用 Primal-Dual 的优雅性质？这种方法的存在是非常符合直觉的，而实际上它也确实存在。我们只需要对站点 \(s_x\) 的势能做出如下修改：
其中 \(y\) 代表的是除 \(x\) 以外的其它站点，\(\Delta'\) 代表的是修改权重后的新“转移代价”；除 \(x\) 外所有站点的势能保持不变。
进行上述修改之后，Primal-Dual 的性质即得到有效保持。我们可以继续使用 Dijkstra 算法简单地更新残量网络。
这一部分的证明与它的更多算法方面的应用本文不再给出；有兴趣的朋友可以阅读 这篇文章，文中给出了有关这个问题的更多细节。
和 OI 相关性比较强的部分到这里就结束了，下面是关于这个问题的一些扩展应用，以及一些正在进行中的工作。
我认为下文中的应用相当优雅漂亮，也欢迎大家与我交流探讨。
3. 从网络流到图形学
3.1 前置知识
Voronoi Diagram : 给出某一空间中的一个站点集合，以及一个 在空间中连续的 代价函数 \(f(p,s_x)\)（其中 \(p\) 为空间中任意一点，\(s_x\) 为任意站点），Voronoi Diagram 定义了一种空间划分结构，其将空间划分为若干胞体，使得每个胞体内恰好包含一个站点，且对于空间中任意一点 \(p\)，设该点与站点 \(s_x\) 在同一胞体内，则将其与 \(s_x\) 相匹配的代价 \(f(p,s_x)\) 总是小于 \(p\) 到任意其它站点 \(s_y\) 的代价 \(f(p,s_y)\)。
一个最常见的情形是，定义在欧式空间中，代价函数为欧几里得距离的情况。图形学研究中还常见定义于三维模型表面，代价函数为 Geodesic Distance 的情形。
你可以在 Wikipedia 页面 里找到一些图例，以便对这种结构做出更多的了解。有关各种不同情形下 Voronoi Diagram 求解方法的研究在学术界已经很完善了，此处不再过多赘述。
Power Diagram : Power Diagram 可以看作 Voronoi Diagram 的一种变体。在这种变式下，每个站点 \(s_x\) 都被赋予了一个权重 \(w(s_x)\)，代价函数也从 \(f(p,s_x)\) 变成 \(f(p,s_x)-w(s_x)\)。Power Diagram 给出了一种能使得空间中任意点的匹配代价最小的空间划分结构。
你同样可以在 Wikipedia 页面 找到更多关于 Power Diagram 的介绍。
3.2 从离散到连续
回到我们的原始问题：欧氏空间中 有 \(n\) 个散点与 \(m\) 个站点，每个站点 \(s_i\) 有一个容量限制 \(C_i\)（我们不妨设 \(\sum C_i=n\)）；将散点 \(p_i\) 分配给站点 \(s_i\) 的代价为 \(D^2(i,j)\)，即欧氏距离的平方。现需要求出一组分配方案，以最小化总代价；我们不妨设 \(n\gg m\)。
早在 1992 年，Aurenhammer 就在 这篇文章 中证明了一条相当漂亮的结论：当代价函数满足最小二乘性（例如，题目中的“平方欧氏距离”）时，我们总能对这 \(m\) 个站点计算出一组权重，使得按这组权重求出 Power Diagram 后，任意站点 \(s_x\) 所在胞体中的散点数量恰好等于题目中定义的容量限制 \(C_x\)。
在论文中，Aurenhammer 还给出了一组基于几何性质的做法，以求解出最优匹配方案与 Power Diagram 的权重。然而这套算法实在称不上优美——它高度依赖于几何性质，无法扩展至任意代价函数，且需要用到动态凸包等各种复杂的数据结构。
现在我们考虑这样一个问题：既然在本文的第一章中，我们已经给出了一种能求出匹配方案的算法，接下来能否直接用我们的匹配方案，来计算出一套权重，从而避免掉大数据结构？
答案是当然可以。考虑某个散点 \(p_i\)，若它被分配给站点 \(s_x\)，则在 Power Diagram 中，对任意其它站点 \(s_y\)，下式必然成立：
移项后我们得到：
不难发现，这组式子定义了一个 OI 中常见的 差分约束问题。于是我们可以根据用网络流策略计算出的匹配结果，用一次 Bellman-ford 算法来计算出一套权重。
结束了？并没有，一个令人拍案叫绝的事实是，在离散的图论与连续的几何结构之间，还存在着一层更加美妙的联系：
如何证明？只需略微回忆 Primal-Dual 的思想，再观差上面的差分约束不等式，这一结论就呼之欲出了！
这实在是一条太美妙的结论，它让我们在离散的图论模型与连续的几何模型之间建立起了深刻的联系。
3.3 从静态到动态
Centroidal Power Diagram (CPD) : 一种特殊的 Power Diagram，在此类几何结构中，所有站点均位于其所在胞体的集合重心处。这种结构在三维重建、图像采样等领域都有非常重要的应用。
Capacity-Constrained Voronoi Tessellation (CCVT) : 我们可以将它理解成另一种 Power Diagram。CCVT 的每个站点不固定权重，而是对其所在胞体的体积（容量）进行限制（事实上，虽然此处的“容量”是一个连续的概念，但我们同样可以证明，权重大小与容量限制是可以互相转化的）。
下面我们来思考一个问题：给定一个几何空间、空间中的站点数量、各个站点的容量限制，允许我们自定义站点坐标，我们能否在这个空间中构造出一个 CPD？
在图形学圈子里，已经有不少工作研究了这一问题。主流的方法大多采用 L-BFGS 或类似的基于数值计算的迭代法；如果你对这些算法感兴趣的话，可以阅读 这篇文章。
这些基于迭代的算法效率还算不错，但它们有一些不可避免的弊端：首先，它们对空间结构的数学性质有比较强的要求，大部分算法是定义在流形表面的，或者要求空间具有一定的可微性；其次，它们还对代价函数的数学性质有较强的要求，因此如果遇到各向异性或极其特殊的代价函数的话，别说建立连续的空间结构，就连求解离散的匹配方案，它们都很难做到。
因此，如果我们要求解的几何空间是非流形、有孔洞、有多个独立部分或 Triangle Soup 等特殊结构，或者代价函数是非传统的各向异性场，就需要一种全新的模式来求解。
那么到底如何求解呢？答案很简单，将空间采样成离散的点云，然后直接做我们的网络流匹配算法。
诚然，这种做法会损失一些精确度；而且在一些特殊代价函数场景下，我们仍然难以根据匹配方案算出权重，或者刻画出胞体的明确几何边界。但与在这些场景下失效的数值计算方法相比，它已经是一种相当有力的工具了。
此时就出现了一个问题：Power Diagram 的问题解决了，那么我们如何保证站点全部位于胞体重心呢？
这个问题也很好解决——OI 圈子里的大家应该都知道 K-means 吧？在算法最开始的部分，我们先随机选取站点的位置；每一轮 Power Diagram 求解完成后，我们将站点移动到其所在胞体的重心位置，然后进行下一轮的求解，直到总代价在精度范围内不再发生明显变化为止。事实上，这种类似 K-means 的迭代策略还有一个名字，叫作 Lloyd's Method。
至此，我们的问题完美解决了——吗？不，还没有！
在实践中，Lloyd's Method 的收敛是非常迅速的。其实在常数轮的迭代之后，站点的坐标变化就非常小了，散点的匹配方案也几乎没有变化。因此，我们如果在每一轮迭代中重新跑一遍匹配，就会造成巨大的额外时间开销。在实际的工程中，这样的时间代价就显得过于巨大了。
对于这种情况，一个经典的解决策略是，在迭代中增加 Warm-Start。也就是说，在每轮迭代中，我们先继承上一轮的匹配结果，在该结果的基础上进行一些细微的调整，从而避免从零开始计算的巨大时间开销。
还记得我们在 2.3 中介绍过的“边权带修改”问题吗？没错，只需要使用这种策略，我们就能够做到 Warm Start，从而极大地优化迭代过程中的时间开销！
至此，这个问题才算是真正地结束了。
4. 最后的碎碎念
以上就是我想要与大家分享的全部内容了！感谢你能抽时间阅读这些我的小想法，也欢迎与我交流探讨。当然，如果文中有错误或者需要改进的地方，也请无比留言告诉我！
在研究、整理这些内容的时候，我几乎第一次意识到原来 OI 界的算法也可以在科研工作中派上一点实际的用场。这实在是一次相当美妙的体验。
本文中的部分理论以论文形式于 Pacific Graphics 2026 在投中。还有一些是纯粹的小思考，我希望可以把它们整合到后续的工作当中。
浙公网安备 33010602011771号