圆上所有传输质量的部分最优传输的O(N log N)时间算法

# 圆上所有传输质量的 O(N log N) 部分最优传输
来源：https://arxiv.org/html/2608.23910

Soheil Kolouri  
隶属：联合计算学院  
隶属：范德堡大学，田纳西州纳什维尔  
邮箱：[soheil\.kolouri@vanderbilt\.edu](mailto:[email protected])

###### 摘要
部分最优传输在比较两个测度时，允许部分质量不匹配，这使其对异常值、遮挡和杂乱具有鲁棒性。关注的量通常是整个*轮廓*——即每个传输基数下的最优成本——因为合适的传输量很少能预先知道，而在实数线上，PAWL 算法能在 O(N log N) 时间内返回该轮廓。
许多数据是周期性而非线性的：角度、相位、方向、一天中的时间、色相，以及通过投影到大圆上获得的所有方向。在圆上，相同问题获得了一个全局环流，或等效地获得了一个需要优化的切割。朴素的精确方法通过在每个支撑间隔上运行一次线算法来处理此问题，时间复杂度为 O(N² log N)。
我们证明了这个因子 N 是不必要的。线结构以无切割的形式得以保留，并且一个*自由间隔不变性*在每一步都提供一个切割，使得所有先前的局部更新在该切割上仍然有效的线更新。这产生了 PAWC：一个精确的、O(N log N) 时间、O(N) 内存的算法，一次运行即返回所有 K+1 个成本、嵌套的活动集和计划，同时返回一个对于每个基数都同时最优的单一间隔。通过在大圆上切片，可将其推广到 S^{d-1}。
实证表明，当 N=4096 时，整个轮廓的计算成本为 0.56 毫秒，而一般求解器计算单个传输比例需 1.5 秒；在杂乱遮挡的 mpeg-7 形状上，固定描述符仅改变成本时，它保留了 66% 的干净数据检索分数，而平衡圆形 OT 仅保留 16%；在 S² 上，它将球形切片 Wasserstein 的拟合误差减半，无论是合成数据还是真实污染目标。代码可在 https://github.com/mint-vu/Partial_Wasserstein_on_Circles 获取。

## 1 引言
最优传输 (OT) 提供了一个几何感知的框架来比较测度 (Villani, 2009 (https://arxiv.org/html/2608.23910#bib.bib25); Peyré and Cuturi, 2019 (https://arxiv.org/html/2608.23910#bib.bib6))，但精确的离散求解器通常开销很大。部分传输扩展了该框架，只允许分布的一部分进行匹配：概念上，质量可以从源传输到目标，如果留在源处不匹配则被销毁，或如果在目标处不匹配则被创建，对应有传输成本、销毁成本和创建成本。
在一维情况下，Bai 等人 (2023) (https://arxiv.org/html/2608.23910#bib.bib17) 利用支撑的排序，为固定惩罚选择开发了一个精确的原始-对偶求解器，最坏情况复杂度为 O(n max{m,n})，对于规模相当的分布为 O(N²)。最近，Chapel 和 Tavenard (2025) (https://arxiv.org/html/2608.23910#bib.bib3) 表明，整个由传输质量索引的解路径可以更快地计算：他们的 PAWL 算法在 O(N log N) 时间内返回每个可能传输质量的精确部分传输计划，使用嵌套活动集、相邻候选对、平衡链和链边缘概率的堆。
未匹配质量通常是问题的一部分，而非建模上的麻烦。部分对应关系隔离了两个测度共享的结构，而无需将虚假检测、遮挡特征、位置错误的记录或异常值强制进行几何上不合理的匹配 (Nietert et al., 2022 (https://arxiv.org/html/2608.23910#bib.bib26))。然而，合适的匹配质量量很少是预先已知的，并且可能取决于下游的鲁棒性或覆盖要求。因此，自然的对象是整个部分传输轮廓 k ↦ C_k^∘，其中 C_k^∘ 是恰好匹配 k 对的最小成本。该轮廓量化了覆盖与几何保真度之间的权衡：它显示了随着传输质量的增加，哪些对应关系被保留，以及何时允许额外质量开始需要代价高昂的匹配。
一次性计算所有 C_k^∘ 允许在检查此权衡后选择传输质量，而不是在解决传输问题之前固定 (Chapel and Tavenard, 2025 (https://arxiv.org/html/2608.23910#bib.bib3))。
周期性几何无处不在，而非特例。角度、相位、方向、一天中的时间变量和色相天然地表示在圆上，而球面上定向数据的许多切片构造会产生在大圆上支撑的测度。然而，圆的循环几何阻止了直接应用在线上的有效算法。圆保留了循环顺序，但没有规范的起点、终点或零流边界。对于平衡的 W₁，缺失的边界条件在流公式中表现为一个自由加法环流，或等效地表现为一个切割，必须在该切割上优化一个提升的一维问题 (Delon et al., 2010 (https://arxiv.org/html/2608.23910#bib.bib4); Rabin et al., 2011 (https://arxiv.org/html/2608.23910#bib.bib7))。部分传输将这种拓扑自由度与传输质量耦合起来：随着 k 变化，活动原子以及最优环流（从而最优切割）都可能改变。因此，没有单一的圆展开能恢复完整的轮廓 k ↦ C_k^∘。
一个直接的精确化简在每个支撑间隔处切割圆，在每个生成的线上运行 PAWL，并取生成轮廓的下包络。虽然正确，但这在 O(N) 个切割上重复了 O(N log N) 的计算，因此代价为 O(N² log N)。核心算法挑战是无需为每个可能的切割单独付费，即可恢复整个圆形部分传输轮廓。

##### 贡献。  
线结构以无切割的形式在圆上得以保留（第3节 (https://arxiv.org/html/2608.23910#S3)）：圆形最优活动集可以选择为嵌套的，在每一步激活的两个原子在非活动原子中是循环邻居，并且一个*自由间隔不变性*——每个当前单元总是包含一个尚未被任何选定单元使用的原始间隔——提供了一个切割，使得所有先前的局部更新在该切割上同时是有效的线更新。这给出了一个精确的贪心规则，以及一个对于*每个*基数都同时最优的间隔，它是构造的而非搜索的（定理3 (https://arxiv.org/html/2608.23910#Thmtheorem3) 和 4 (https://arxiv.org/html/2608.23910#Thmtheorem4)），以及一个 O(N log N) 时间、O(N) 内存的算法，一次运行即返回所有成本、嵌套的活动集和计划（第4节 (https://arxiv.org/html/2608.23910#S4)）。在大圆上切片将其推广到 S^{d-1} (第5节 (https://arxiv.org/html/2608.23910#S5))，第6节 (https://arxiv.org/html/2608.23910#S6) 通过与 LP 预言机对比验证了精确性，测量了缩放性，并给出了两个以轮廓为目标的应用——mpeg-7 上形状描述符的部分匹配，以及在 S² 上（包括一个真实的地震目录）的鲁棒拟合。

## 2 背景
### 2.1 圆上的部分最优传输
令 𝕊¹_L = ℝ/Lℤ，其测地距离 d_{𝕊¹}(x,y) = min{ |x-y|, L-|x-y| }，并令 μ = w ∑_{i=1}^n δ_{x_i} 和 ν = w ∑_{j=1}^m δ_{y_j} 为均匀加权的经验测度；记 N = n+m，K = min(n,m)。对于传输质量 s ∈ [0, Kw]，定义圆形部分 1-1 Wasserstein 成本为：
PW_∘(s) := min_{π∈ℝ_+^{n×m}} { ∑_{i=1}^n ∑_{j=1}^m d_{𝕊¹}(x_i, y_j) π_{ij} : ∑_{j=1}^m π_{ij} ≤ w, ∑_{i=1}^n π_{ij} ≤ w, ∑_{i,j} π_{ij} = s }。(1)
因此，PW_∘(s) 是传输恰好 s 单位质量的最小成本，其中每个原子最多参与其全部质量 w。将质量守恒放宽为不等式并固定总传输质量是 Caffarelli 和 McCann (2010) (https://arxiv.org/html/2608.23910#bib.bib2) 的部分问题，他们在成本 h(x-y) 下（h 凸且支撑不相交）建立了存在性和唯一性；Figalli (2010) (https://arxiv.org/html/2608.23910#bib.bib5) 对于二次成本放宽了不相交性，并研究了映射 s ↦ PW(s)。我们计算的正是这个映射，而非其任何单一值。

(a) 𝕊¹_L 上的固定质量部分 OT  
x₁ y₁ x₂ y₂ x₃ y₃ x₄ y₄ w  
d_{𝕊¹}(x₁,y₁)  
w  
d_{𝕊¹}(x₂,y₂)  
w  
∂_d ∂_c ∂_d ∂_c  
实心：已传输；空心：未匹配

(b) 精确的质量和成本核算  
J_k(M) = w ∑_{(i,j)∈M} d_{𝕊¹}(x_i, y_j)_{传输成本} + λ_d (n-k)w_{销毁的源质量} + λ_c (m-k)w_{创建的目标质量}

图 1：圆上的部分传输。(a) 实心弧是选定原子之间的测地传输；空心原子未匹配，其质量在 ∂_d 处销毁或在 ∂_c 处创建。(b) 在固定的传输基数 |M|=k 时，这两个量是常数，因此任何固定的惩罚 λ_d, λ_c 都贡献一个常数，问题等同于基数约束匹配 (2 (https://arxiv.org/html/2608.23910#S2.E2))。

###### 假设 1（常规假设）。
所有原子携带相同质量 w；点集 {x_i} ∪ {y_j} 两两不同；底层成本为 d_{𝕊¹}。

在整数质量 s = kw 时，将可行集按 w 缩放得到基数 k 的二部匹配多面体。其极点是整数的，因此 C_k^∘ := PW_∘(kw) = min_{M: |M|=k} w ∑_{(x_i,y_j)∈M} d_{𝕊¹}(x_i, y_j), k=0,...,K，(2)
其中 M 遍历匹配，因此没有原子出现在多个配对中。我们称 k ↦ C_k^∘ 为*轮廓*，并用 A(M) 表示与 M 关联的原子集合，称为其*活动集*。对于 s ∈ [kw, (k+1)w]，记 λ = (s-kw)/w；该值线性插值，PW_∘(s) = (1-λ) C_k^∘ + λ C_{k+1}^∘，(3)
因此离散轮廓 {C_k^∘}_{k=0}^K 决定了每个 s ∈ [0, Kw] 的 PW_∘(s)；见附录 I (https://arxiv.org/html/2608.23910#A9)。

### 2.2 计算部分传输
OT 及其不平衡变体的精确求解器是立方级的；标准补救措施是熵正则化 (Cuturi, 2013 (https://arxiv.org/html/2608.23910#bib.bib22)) 和切片传输所依赖的一维分位数公式 (Rabin et al., 2012 (https://arxiv.org/html/2608.23910#bib.bib23); Peyré and Cuturi, 2019 (https://arxiv.org/html/2608.23910#bib.bib6))。对于部分传输，Phatak 等人 (2023) (https://arxiv.org/html/2608.23910#bib.bib19) 以 O(n³) 计算完整质量轮廓；Bai 等人 (2023) (https://arxiv.org/html/2608.23910#bib.bib17) 给出一个二次的一维求解器；Bonneel 和 Coeurjolly (2019) (https://arxiv.org/html/2608.23910#bib.bib18) 以拟线性时间解决单射部分分配；Séjourné 等人 (2022) (https://arxiv.org/html/2608.23910#bib.bib20) 以拟线性时间解决一维不平衡 OT，但无法指定传输质量；Chapel 等人 (2021) (https://arxiv.org/html/2608.23910#bib.bib21) 获得一个正则化路径。PAWL (Chapel and Tavenard, 2025 (https://arxiv.org/html/2608.23910#bib.bib3)) 是首个以 O(n log n) 在线上返回完整部分轮廓的算法，是我们扩展的基础。
在圆上，先前的工作是平衡的：Delon 等人 (2010) (https://arxiv.org/html/2608.23910#bib.bib4) 和 Rabin 等人 (2011) (https://arxiv.org/html/2608.23910#bib.bib7) 通过优化切割或全局环流求解 W₁；Díaz Martín 等人 (2024) (https://arxiv.org/html/2608.23910#bib.bib10) 线性化圆形 OT；Bonet 等人 (2023) (https://arxiv.org/html/2608.23910#bib.bib8)；Liu 等人 (2025) (https://arxiv.org/html/2608.23910#bib.bib9) 在球形切片中使用圆形传输。没有一个处理圆上指定质量的部分传输。

##### 为何圆不同于线。
对于一个*固定的*活动集，圆形成本是展开线成本在所有 N 个支撑间隔上的最小值（附录 C (https://arxiv.org/html/2608.23910#A3) 中的命题 1 (https://arxiv.org/html/2608.23910#Thmproposition1)），因此在间隔和基数上的包络是一个精确的 O(N² log N) 参考求解器（推论 2 (https://arxiv.org/html/2608.23910#Thmcorollary2)，也在那里）——我们在第 6 节 (https://arxiv.org/html/2608.23910#S6) 中用作基准，并在测试中用作可信的第二意见。附录 F (https://arxiv.org/html/2608.23910#A6) 中的示例 1 (https://arxiv.org/html/2608.23910#Thmexample1) 显示重用一个固定切割可能相差一个无界因子，第 K.5 节 (https://arxiv.org/html/2608.23910#A11.SS5) 测量了单一切割实际上足够使用的频率。

## 3 圆形部分传输的结构

### 3.1 嵌套活动集与循环邻居：每步 O(N) 个候选
###### 定理 1（嵌套最优扩展；在附录 D (https://arxiv.org/html/2608.23910#A4) 中证明）。
设 M_k 是任意基数 k 的最优匹配，k ≥ 0。则存在一个基数 k+1 的最优匹配 M_{k+1}，使得 A(M_k) ⊆ A(M_{k+1})。

###### 推论 1（每步两个候选）。
设 M_k 如定理 1。设 I = X ∪ Y \ A(M_k) 为非活动原子集。则基数 k+1 的最优扩展 M_{k+1} \ M_k 包含恰好一对 (i,j) ∈ I_x × I_y。此外，(i,j) 在循环顺序下是 I_x 和 I_y 的循环邻居。

**循环邻居的定义**：对于支撑在圆上的两个点集 U 和 V，点对 (u,v) ∈ U × V 称为循环邻居，如果在圆上沿着某个方向（顺时针或逆时针），在 u 和 v 之间没有其他点属于 U 或 V。

###### 定理 2（自由间隔不变性；在附录 D (https://arxiv.org/html/2608.23910#A4) 中证明）。
在定理 1 的设定下，假设 M_k 是基数 k 的最优匹配。那么对于 M_k 的每个最优扩展 M_{k+1}，存在一个间隔 g_k（即支撑 X ∪ Y 中相邻原子之间的圆弧），使得：
1. g_k 包含一个原始间隔（即 X ∪ Y 中至少一对相邻原子之间的弧，该弧不被任何选定单元覆盖）。
2. 对于所有 i ≤ k，M_i 是从 M_k 沿着 g_k 切割并展开后，在线版本上的基数 i 的最优匹配。

自由间隔 g_k 可以构造为：取使得循环单元 C_k(M_k) 包含一个原始间隔的任意间隔。

**关键思想**：该不变性确保在每一步，都存在一个“自由”切割，使得当前和所有之前的最优匹配在该切割下对应于线上的最优匹配。这个间隔是“构造的”而非“搜索的”。

## 4 PAWC 算法
基于第 3 节的结构，我们提出了 PAWC（圆上部分 Wasserstein）算法。

**输入**：圆上两个均匀加权点集 X = {x₁,...,x_n}，Y = {y₁,...,y_m}，权重为 w。假设点按循环顺序排列。
**输出**：对于每个 k=0,...,K（其中 K=min(n,m)），成本 C_k^∘、活动集 A_k 和匹配 M_k，以及一个对每个 k 都同时最优的自由间隔 g。

**算法概述**：
1. 初始化 M₀ = ∅，A₀ = ∅，C₀^∘ = 0，g 为任意间隔。
2. 对于 k = 0 到 K-1：
   a. 确定非活动集 I = X ∪ Y \ A_k。
   b. 根据定理 1 和推论 1，最优扩展 (i,j) 必须是 I_x 和 I_y 的一对循环邻居。
   c. 计算所有这样的候选邻居对 (i,j) 的扩展成本。
   d. 选择使总成本 C_{k+1}^∘ = C_k^∘ + w·d_{𝕊¹}(x_i, y_j) 最小的对 (i^*, j^*)。
   e. 更新：M_{k+1} = M_k ∪ {(i^*, j^*)}，A_{k+1} = A_k ∪ {i^*, j^*}。
   f. 更新自由间隔 g：根据定理 2，确保新的单元 C_{k+1}(M_{k+1}) 包含一个原始间隔（例如，如果添加的边跨越了 g，则更新 g 为新单元中的一个原始间隔）。

**复杂度**：由于在每一步，循环邻居对的数量最多为 O(|I|) = O(N)，并且使用堆等数据结构，可以在 O(log N) 时间内找到最小成本扩展，总时间复杂度为 O(N log N)。内存使用为 O(N) 以存储活动集、匹配和间隔信息。

**正确性**：定理 1、推论 1 和定理 2 保证了算法的贪心步骤是正确的，并且自由间隔不变性确保了所有先前的步骤在线上仍然是最优的，因此整个轮廓是精确计算的。

## 5 推广到 S^{d-1}
圆上的部分传输可以直接推广到球面 S^{d-1}，通过切片方法：将球面测度投影到通过原点的大圆上，得到圆上的一维测度，然后应用 PAWC 算法计算该切片上的部分传输轮廓。对所有大圆切片积分（或求和）即可得到球面上的部分传输轮廓。这类似于经典切片 Wasserstein 距离的构造，但现在是在部分传输的框架下。

## 6 实验与应用
### 6.1 精确性验证与缩放性
- 使用一个线性规划 (LP) 预言机作为基准，在小规模问题上验证 PAWC 的精确性。
- 测量 PAWC 与基准算法（如在 N 个切割上运行 PAWL 的朴素方法）的运行时间缩放。结果显示 PAWC 为 O(N log N)，而朴素方法为 O(N² log N)。

### 6.2 应用一：mpeg-7 形状的部分匹配
- 任务：在 mpeg-7 数据集上进行形状检索，考虑形状有遮挡或杂乱部分。
- 方法：使用形状的轮廓点作为测度。比较 PAWC 计算的整个轮廓与平衡圆形 OT (固定质量传输) 的结果。
- 结果：当保持描述符固定，仅改变匹配成本时，基于 PAWC 的部分匹配方法保留了 66% 的干净数据检索分数，而平衡圆形 OT 仅保留 16%。这证明了部分传输在处理不完整数据时的鲁棒性。

### 6.3 应用二：S² 上的鲁棒拟合
- 任务：将球面模型拟合到可能被污染（包含异常值）的观测数据中。
- 方法：使用球面切片 Wasserstein 距离作为拟合度量。比较使用完整质量传输与使用 PAWC 计算的部分传输轮廓（允许不匹配）的结果。
- 结果：在合成数据和真实地震目录上，使用部分传输轮廓将拟合误差减半，表明它在存在离群点时更具鲁棒性。

## 7 结论
本文提出了 PAWC，一个精确的、O(N log N) 时间的算法，用于计算圆上部分最优传输的整个轮廓。其关键创新是“自由间隔不变性”，它避免了为每个可能的切割单独求解，将复杂度从 O(N² log N) 降低到 O(N log N)。算法通过在大圆上切片推广到球面。实验证明了其精确性、高效性以及在形状匹配和鲁棒拟合等应用中的优势。代码已开源。