简介面向最优传输与距离度量研究者的Matlab/C混合实现资源聚焦非标准化Wasserstein-2UW-2距离的一维与二维高效计算。代码源自Gangbo等2019年论文《Unnormalized Optimal Transport》并封装为C类与Mex接口使Matlab调用便捷高效。压缩包共19个文件包含8个cpp源码一维/二维求解器与Mex入口、7个头文件、2个m示例脚本以及README和LICENSE说明整体仅31KB结构清晰。已有375人学习浏览。资源提供完整求解器类、Mex调用示例及可编译的main入口能够帮助读者深入理解UW-2的算法流程与实现细节并快速迁移到自身实验或工程场景省去从零复现数学推导的成本。包内cpp与h文件构成算法主体m文件提供Matlab调用范例md文档则为编译使用说明便于快速上手尤其适合计算最优传输方向的科研人员与进阶开发者。1. 非标准化 Wasserstein 距离为什么工程里更常见的 EMD 是这个版本做图像块直方图匹配、特征点集对齐或者深度图相似度评估时很多人第一反应是拿标准 EMDEarth Movers Distance来算。但经典 EMD 默认要求两端分布的总质量严格相等而实际场景里这几乎不成立左侧图检测出 3000 个特征点右侧图只有 1200 个两幅图的边缘像素总数不一样同一物体被遮挡后可见面积差一大截。把这种天然失衡的数据硬塞给标准化的最优传输模型要么解不可行要么得先偷偷把总质量归一化成一致——而归一化恰恰丢掉了总质量差异这个关键信息。这份标题指向的 MATLAB 代码库解决的是 unnormalized optimal transport 问题它直接计算非标准化 Wasserstein 距离允许源分布和目标分布质量不相等。你只需要提供两个质量向量和一个代价矩阵就能拿回从源到目标的传输代价和传输流。本文会把这个方案的数学背景、可直接抄的 MATLAB 实现、参数调节和一堆踩坑记录完整拆开让新手能跑通让熟手能看到边界。2. 从连续模型到离散实现非标准化 EMD 的数学拆解2.1 标准 EMD 与 unnormalized 的本质差异总质量不再是常数标准最优传输Standard Optimal Transport问题是这样的给定两个概率分布 P 和 Q总质量都等于 1求一个联合分布 Γ传输计划使得边缘分布分别等于 P 和 Q且传输代价最小。用离散语言写出来就是min Σ_i Σ_j C(i,j) * Γ(i,j)s.t. Σ_j Γ(i,j) P(i) Σ_i Γ(i,j) Q(j) Γ(i,j) 0在标准版本里sum(P) sum(Q) 1所以上面两组约束是自洽的。非标准化版本把约束右端的质量向量换成任意非负实数。于是问题来了如果 sum(P) ≠ sum(Q)这两组等式约束不可能同时严格满足——源总溢出质量或者目标总欠收质量数学上无解。工程上处理这个矛盾的做法很统一引入虚拟节点把质量缺口补平。sum(P) 大于 sum(Q) 时在目标端加一个垃圾场节点把多余的源质量以单位成本 dump_cost 倾倒掉sum(P) 小于 sum(Q) 时在源端加一个无中生有节点以单位成本 dump_cost 凭空制造缺失的质量。这也是标题里 unnormalized-optimal-transport 的核心含义不是把数据归一化后再做 EMD而是允许总质量本身作为度量的一部分参与计算。质量差越大或 dump_cost 定得越高距离值就越大这在有遮挡的面积匹配、有缺帧的轨迹匹配里是真实可解释的语义。2.2 离散化的三个输入质量向量、代价矩阵和传输流落到代码层非标准化 EMD 只需要三样东西质量向量 Pn×1源分布每个 bin/特征点上的质量不必归一化也不必是整数。比如图像直方图 256 个 bin每个 bin 是像素计数这就是质量向量。质量向量 Qm×1目标分布的质量长度可以和 P 不一样总和也可以不一样。这里的不一样正是 unnormalized 的价值所在。代价矩阵 Cn×m把单位质量从源 i 搬到目标 j 的成本。最常用的是坐标欧氏距离也可以是自定义的几何代价、能量代价。输出是标量 cost非标准化 Wasserstein 距离和流量矩阵 flown×m。flow(i,j) 表示实际从源 i 流向目标 j 的质量。cost 按题目的目标函数写成 sum(C .* flow)它衡量的是搬动这些质量花了多少总代价。2.3 为什么这个场景用 MATLAB 写最顺从工程角度MATLAB 的 linprog线性规划求解器可以直接解出上面的传输问题。最优传输本质是一个带等号约束的线性规划变量数是 n*m约束数是 nm对几百个 bin 的问题规模来说linprog 的 dual-simplex 算法几十毫秒内能解完。相比自己写网络单纯形法或者 Sinkhorn 迭代MATLAB 方案的开发周期最短约束矩阵能直接用稀疏形式构造后续还能用 MEX 包一层 C 代码加速。另一个理由是矩阵可视化方便算完 flow 直接用 imagesc 看传输流分布排查问题比 Python 环境下快不少。3. 跑通最小代码非标准化 EMD 主函数与逐段拆解3.1 第一步从坐标构造代价矩阵拿到的是点集坐标而不是现成的代价矩阵时第一步用 pdist2 批量算距离矩阵。注意第 2.2 节说的代价不一定是欧氏距离但欧氏距离是最常见起点。% demo_emd_inputs.m % 源点3 个点目标点2 个点 src_coords [0, 0; 1, 2; 3, 1]; dst_coords [2, 2; 4, 0]; % 源质量与目标质量故意不等体现 unnormalized P [2; 5; 3]; % 总质量 10 Q [4; 2]; % 总质量 6短缺 4 % 代价矩阵欧氏距离 C pdist2(src_coords, dst_coords, euclidean);pdist2 的第一个输入是 n×d 矩阵第二个是 m×d 矩阵返回 n×m 距离矩阵。这里 P 总质量 10Q 总质量 6属于供大于求按 2.1 节方案需要在目标端加一个虚拟垃圾场节点。如果坐标里混入 NaN比如特征点提取失败pdist2 会把整行变成 NaN这一步就要先做过滤。3.2 第二步质量平衡——虚拟节点在哪里加、怎么加这是整个 unnormalized 实现最关键的一步也是最容易翻车的地方。做法是供大于求时在目标端补一列虚拟汇点单位搬运成本 dump_cost供不应求时在源端补一行虚拟源点成本同样是 dump_cost。dump_cost 的含义是单位质量被丢弃或凭空产生的成本。% balance_mass.m function [P2, Q2, C2, ext] balance_mass(P, Q, C, dump_cost) sumP sum(P); sumQ sum(Q); if nargin 4 || isempty(dump_cost) dump_cost max(C(:)); % 常见默认取最大代价 end ext struct(src_added, false, dst_added, false, ... dump_cost, dump_cost); if sumP sumQ % 目标质量多出来需要在源端补一个虚拟源点 deficit sumQ - sumP; P2 [P; deficit]; Q2 Q; C2 [C; dump_cost * ones(1, size(C, 2))]; ext.src_added true; elseif sumP sumQ % 源质量多出来需要在目标端补一个虚拟汇点 excess sumP - sumQ; P2 P; Q2 [Q; excess]; C2 [C, dump_cost * ones(size(C, 1), 1)]; ext.dst_added true; else P2 P; Q2 Q; C2 C; end enddummy 行或列的语义以供大于求为例虚拟汇点列的 flow 值代表从每个源点被丢弃的质量对应成本就是单位丢弃代价。加在哪一侧取决于供需差方向代码里的注释已经写明。注意如果 dump_cost 设得远小于真实搬运代价解算器会倾向于把所有多余质量都丢掉这是参数层面的事本文第 4 章专门展开。3.3 第三步用 linprog 求解线性规划并还原结果平衡之后的总质量一致可以用标准的等号约束线性规划求解。变量是展平后的流量矩阵约束是每个源点流出等于 P(i)、每个目标点流入等于 Q(j)。代码里用稀疏矩阵构造避免大 n、m 时内存炸掉。% emd_unnormalized.m function [cost, flow, info] emd_unnormalized(P, Q, C, dump_cost) % P: n×1 源质量向量不必归一化 % Q: m×1 目标质量向量 % C: n×m 代价矩阵C(i,j) 源i到目标j的单位搬运代价 % dump_cost: 质量差额的单位处理成本可选默认 max(C(:)) % % cost: 非标准化 Wasserstein 距离总传输代价 % flow: n×m 传输流量矩阵flow(i,j) 从i运到j的质量 [P2, Q2, C2, ext] balance_mass(P, Q, C, dump_cost); n length(P2); m length(Q2); % 目标函数最小化 sum(C2 .* flow) f C2(:); % 按列展平与 x 的排列一致 % 源约束每个源点流出之和等于 P2(i) Aeq kron(ones(1, m), speye(n)); beq P2(:); % 目标约束每个目标点流入之和等于 Q2(j) Aeq [Aeq; kron(speye(m), ones(1, n))]; beq [beq; Q2(:)]; lb zeros(n * m, 1); % 流量非负 opts optimoptions(linprog, ... Algorithm, dual-simplex, ... Display, off); [x, val, exitflag] linprog(f, [], [], Aeq, beq, lb, [], [], opts); if exitflag 0 cost Inf; flow []; info struct(exitflag, exitflag, ext, ext); return; end flow_full reshape(x, n, m); if ext.dst_added flow flow_full(:, 1:end-1); % 去掉虚拟汇点列 elseif ext.src_added flow flow_full(1:end-1, :); % 去掉虚拟源点行 else flow flow_full; end cost val; info struct(exitflag, exitflag, ext, ext, ... flow_full, flow_full); end这段代码有四个地方值得展开说明。kron 生成约束矩阵kron(ones(1, m), speye(n))生成 n×(n*m) 稀疏矩阵第 i 行在 x((i-1)m1) 到 x(im) 处为 1正好对应源点 i 到所有目标点的流量之和。第二个 kron 生成 m 行第 j 行在 x(j), x(jn), x(j2n)... 处为 1对应所有源点到目标 j 的流量之和。稀疏构造在 n、m 几百时省内存非常明显。beq 的顺序先拼 P2 再拼 Q2与 Aeq 的行顺序严格对应拼错会直接导致解不可行。exitflag 判断linprog 返回负值时表示不可行或无界这时候强行用 x 会得到 NaN。这里做了保护返回 Inf 让上层逻辑去处理。结果还原linprog 的 x 是按列展平的所以 reshape(x, n, m) 之后 flow_full(i,j) 才是源 i 到目标 j 的流量。去掉虚拟行列后剩下的就是真实点之间的传输流。一个值得注意的点info.flow_full 不要丢掉。尤其在调试阶段看虚拟节点上的流量能立刻判断解是否合理——虚拟汇点流量过大说明丢弃成本太低虚拟源点流量过大说明供需缺口本身太大。4. 参数怎么调bin 数、代价函数、dump_cost 与求解器选项4.1 dump_cost非标准化 EMD 最需要斟酌的参数dump_cost 决定质量差额如何进入最终距离值它直接改变解的结构。下面是不同取值下的行为对比dump_cost 取值行为特征适用场景0多余质量免费丢弃距离值只反映真实搬运部分只关心有效匹配的代价忽略遮挡/缺失面积max(C(:))丢弃单位质量的成本等于最远搬运成本惩罚适中默认值均衡远大于 max(C(:))解会尽量少丢弃哪怕绕远路搬运距离值包含大额惩罚总量失衡本身就是重要的误差信号远小于 min(C(:))解几乎不做远距离搬运直接丢弃多余或伪造缺失点集局部匹配不在乎全局总量这里有个反直觉的坑有人以为 dump_cost 越大越正确但实际上它决定了你是否允许 solver 用丢弃来回避高代价搬运。如果应用场景是遥感图像地物面积对比面积差就该进距离度量dump_cost 要设大如果是特征点匹配点云密度不同是正常的dump_cost 设成 0 反而更干净。我一般会先用 max(C(:)) 跑一遍再对比 dump_cost0 的结果两次的 cost 差值就是质量失衡惩罚占总距离的比例这个比例能帮你判断参数是否合理。4.2 bin 数量与代价函数的设计直方图类数据灰度直方图、色度直方图要决定用多少 bin。bin 太少代价矩阵分辨率不够小位移被平均掉bin 太多每个 bin 质量稀疏、代价矩阵近奇异linprog 的求解时间成倍增长。经验上先从 64 bin 起步算一次看 flow 的稀疏度。如果 flow 里大量非零项集中在相邻 bin 附近说明 bin 密度没问题如果 flow 出现明显横穿的非持续传输从 bin 1 直接到 bin 60说明 bin 量纲和代价函数不匹配。代价函数最常见的三个选择欧氏距离默认首选几何意义直观。曼哈顿距离L1适合网格状街道或图像像素四邻域移动抗噪性比 L2 好。自定义代价比如用颜色差异、梯度方向差的加权和。注意代价函数必须满足非负性和良定义不能出现 NaN 或 Inf。一旦代价矩阵里有一个 NaN整个 linprog 会直接返回空解且 MATLAB 不会给出明确报错这是后面避坑章节要讲的黑匣子问题。4.3 大规模场景下的求解器选项与降维技巧n、m 超过 500 时变量数会到 25 万即使稀疏约束linprog 的求解时间也可能长得不可接受。三个实用手段第一用optimoptions(linprog, Algorithm, dual-simplex)而不是默认的interior-point。dual-simplex 对带非负变量和等号约束的运输问题更稳内存占用也更低。老版本 MATLAB 里这个选项叫ds新版本统一为dual-simplex。第二把质量接近、位置接近的点预聚类。用 kmeans 把 1000 个源点聚成 50 个质量中心每个聚类中心的质量是该簇质量总和坐标取质心。这样把 EMD 的输入从几千维降到几十维代价是丢失簇内结构但作为快速原型完全够用。第三在 C 矩阵上做预剪枝去掉代价超过某个阈值 t 的边设成 Inf 或直接移入虚拟节点。这等于限制只允许短距离搬运在很多图像局部匹配场景下远距离搬运本来就不该发生。剪枝后约束矩阵更稀疏求解速度快一个数量级。5. 非标准化 EMD 最容易翻车的 5 个高频坑5.1 名义上是 unnormalized代码里却偷偷归一化了现象直接拿经典 EMD 教程里的代码开头先P P / sum(P)再丢给求解器。算出的距离在两个场景下失真一是两个分布总质量差距很大但形状相似比如 1000 像素 vs 500 像素的同一物体分割结果归一化后距离很小甚至为 0完全掩盖了面积差异二是遮挡场景中缺失的是哪部分空间结构归一化会把缺失痕迹彻底抹掉。原因归一化隐含了只看分布形状不管总量的假设但 unnormalized EMD 的卖点恰恰是总量差异也要进距离。解决不要在某一个公共预处理函数里偷偷归一化。在代码入口处加断言检查 sum(P) 和 sum(Q) 是否被改动过如果要支持归一化模式显式传一个normalize_flag默认关。5.2 不处理质量差额直接构造等号约束导致不可行现象把 2.1 节那个矛盾问题直接写进 linprogAeq 同时约束源流出和汇流入结果exitflag返回 -2不可行solver 输出空。新手经常在这卡很久以为代码写错了。原因sum(P) ≠ sum(Q) 时两组等式约束在数学上互相矛盾。问题不在求解器而在建模时没加虚拟节点。解决统一走 3.2 节的balance_mass确保进入 linprog 前总质量一致。这是一个 100% 会发生、但只需要三行代码就能规避的问题。5.3 代价矩阵里藏着 NaN/Inflinprog 毫无征兆地返回空现象输入坐标里个别点异常pdist2 在某一行为 NaN或者自定义代价里除法产生 Inf。此时linprog不报错返回空数组xexitflag 经常是 0看起来像没解出来。原因MATLAB 的线性规划接口对 NaN 的处理是直接放弃不会提示你哪一行出问题。解决构造 C 矩阵后加一个断言assert(~any(isnan(C(:))) ~any(isinf(C(:))), ... 代价矩阵存在 NaN 或 Inf先检查坐标输入);同时检查 P、Q 是否有负值或 NaN。这个断言成本极低建议直接写进主线代码别只在调试时开。5.4 用 Sinkhorn 做迭代近似时质量悬殊场景下不收敛的假象现象为了处理大规模数据有人把 linprog 换成 Sinkhorn 迭代熵正则化结果在 sum(P)/sum(Q) 相差 10 倍以上时迭代几百次 cost 还在缓慢漂移输出结果离真实解差很远。原因Sinkhorn 对熵正则项系数 λ 非常敏感质量比越大需要越小的 λ 才能逼近真实解而 λ 越小迭代越慢数值也越不稳。这是算法本身的特性不是参数没调好。解决悬殊场景下别硬调 λ回到 4.3 节的路子——先聚类降维再用 linprog 精确求解或者用 Sinkhorn 的结果做 linprog 的初始点需要把 Sinkhorn 的输出转成可行基工程上麻烦一点但对几百维问题没必要。记住一句话Sinkhorn 适合快速粗糙估计linprog 才是精确答案。5.5 把 flow 直接当作一一对应的匹配关系来用现象算完 flow 矩阵把每一行最大的 flow 值对应的列当作匹配点画连线一看多个源点匹配到同一个目标点还有源点被拆成两半。原因最优传输允许质量拆分。图像直方图场景下这是合理的一个 bin 的颜色扩散到多个相邻 bin但特征点匹配场景下用户期望的是 0-1 配对。解决如果应用要求一一匹配需要额外加整数约束变成 ILP或对 flow 做后处理阈值过滤掉小流量再用匈牙利算法对剩余候选做匹配。还有一个更优雅的做法把 flow(i,j) 归一化成条件概率 P(j|i)用它作为软匹配权重去加权坐标这在多模态配准里比硬匹配更稳。6. 验证与提速怎么确认你算对了再让它跑得更快6.1 三个一眼看穿的验证用例写完代码第一件事不是跑真实数据而是跑三个手工可验的最小例子。% test_emd_unnormalized.m % 用例1单点对单点供不应求 P 5; Q 7; C 3; [cost, ~] emd_unnormalized(P, Q, C); % 总质量差2默认dump3 % 期望真实搬运5*3 虚拟源补齐2*3 21 assert(abs(cost - 21) 1e-9, 用例1失败); % 用例2单点对单点供大于求dump_cost0 [cost, ~] emd_unnormalized(P, 3, C, 0); % 期望真实搬运3*3 9多余2免费丢弃 assert(abs(cost - 9) 1e-9, 用例2失败); % 用例3两个点对两个点且总质量相等 P2 [1; 2]; Q2 [2; 1]; C2 [0, 10; 10, 0]; % 同位置代价0跨位置代价10 % 期望最优是把源1的1送到目标1源2的2送到目标2代价0 [cost, ~] emd_unnormalized(P2, Q2, C2); assert(abs(cost) 1e-9, 用例3失败);用例 1 和用例 2 分别验证供不应求和供大于求两条分支用例 3 验证标准 EMD 退化为 0 的边界。三个例子加起来覆盖 balance_mass 的两个扩展方向和 linprog 的基本正确性。跑通过再上真实数据能省大量排查时间。6.2 提速的三个实用习惯第一个习惯是用tic/toc把 cost 的求解时间打点阈值设为 100ms。超过阈值先看是不是约束矩阵太密——把 C 中超过 3 倍中位代价的边剪掉通常立刻见效。第二个习惯是复用代价矩阵。在批量对比场景下比如一组 100 张图两两算距离代价矩阵只和数据本身有关可以提前算好存成.mat别每次调用都重新 pdist2。第三个习惯是给 linprog 提供一个好的初始可行点。在等号约束下用一个伪随机但满足质量守恒的初始流比如按质量比例均匀打散作为x0传入dual-simplex 的热启动一般能省掉第一阶段的基寻找时间。6.3 从距离走向应用非标准化距离的三个落地方向算出的 cost 除了当相似度阈值还有几个我实际用过的方向一是做图像检索的排序特征cost 比余弦距离在局部微小偏移但总质量接近的场景下区分度高得多二是做点云配准的代价项把非标准化 EMD 的距离梯度反传到位姿参数上比直接 ICP 更鲁棒三是做时间序列的异常检测把每个时间窗口的局部直方图和全局直方图算非标准化距离距离突变点就是异常点。每个方向都只需要把输入换成对应领域的质量和坐标定义主函数不用改一行。说句掏心窝的话这个代码我最初也踩过 5.1 那个归一化的坑被一份看起来更标准的预处理代码坑了两天从那以后我养成了所有输入输出加断言、所有手工用例跑一遍的习惯。也劝你保留那个test_emd_unnormalized.m以后换 MATLAB 版本、改求解器选项跑一遍测试就知道有没有退步。希望帮到你。本文还有配套的精品资源点击获取