简介面向最优传输计算的研究者与工程人员这一开源仓库提供一维与二维非标准化 Wasserstein-2UW-2距离的完整可运行代码算法实现源自 Gangbo 等论文《Unnormalized Optimal Transport》。核心代码以 C 编写并通过 Mex 封装使 Matlab 可直接调用若编译环境缺少 Mex 编译器代码会自动排除 Matlab 部分仅保留 C 接口兼顾不同使用场景。压缩包体积极小仅 31KB共 19 个文件其中以 8 个 cpp 源文件和 7 个 h 头文件为主另含 2 个 m 脚本、1 个 Markdown 说明和 1 个许可证文件文件分工清晰便于定位求解器主程序、辅助函数或示例入口。已有 375 人学习仓库内分别提供一维与二维两套求解器类UnnormalizedOTSolver1D 和 UnnormalizedOTSolver并配有 main.cpp 示例与 mex 入口脚本演示了从输入数据到距离输出的完整流程。读者可直接编译运行或通过 Matlab 调用免去从零实现底层优化算法的繁琐也便于进一步研究未归一化最优运输的理论与数值实现。1. 非标准化Wasserstein距离为什么你的EMD代码一换数据就翻车写 Matlab 算非标准化 Wasserstein 距离也就是常说的 EMDEarth Movers Distance时最常见的翻车点不是算法本身而是大多数人直接套用了质量守恒的标准 EMD 代码。标准 EMD 假设两个分布总质量相等但实际任务里直方图总质量几乎不可能相等——不同曝光时间的图像直方图、不同面积区域的颜色分布、不完整数据源的概率质量都天然带质量差。把这些数据丢进带质量守恒约束的代码要么 linprog 报约束不可行要么对输入先做归一化、把质量信息直接丢掉。这篇文章从模型讲起给出 LP 精确解和 Sinkhorn 近似解两套可直接落地的 Matlab 代码把约束修正、参数设定和五个真实坑一次讲清新手照着跑熟手对照找边界。2. 从标准 EMD 到非标准化距离两种模型与三套算法路径2.1 标准 EMD 的质量守恒约束与隐藏前提标准 EMD 的离散形式是给定源分布 pnP 个 bin和目标分布 qnQ 个 bin成本矩阵 C(i,j) 表示单位质量从 p 的第 i 个 bin 搬到 q 的第 j 个 bin 的代价目标是找一个耦合矩阵 T满足min C, T 约束T * 1 p T * 1 q T 0这里C, T是 Frobenius 内积T * 1 p表示从每个源 bin 搬出的总质量等于 p(i)T * 1 q表示搬到每个目标 bin 的总质量等于 q(j)。把两个约束分别对 i 和 j 求和立刻推出sum(p) sum(q)。这个等式就是质量守恒约束也是整套模型的隐藏前提。实际工程里绝大多数直方图不满足这个前提。常见做法是悄悄把 p、q 各自除以总和再丢进标准 EMD 代码。这个预处理看似无害实际等于把分布质量信息完全抹掉了。同样的分布形状、总质量相差 10 倍的两组数据算出来的距离完全一样这违背了引入 EMD 度量物理搬运代价的初衷。更要命的是有些代码在质量不等时直接报无可行解初学者一看报错就以为需要归一化形成了错误循环。2.2 松开质量守恒钩子质量缩放模型与非平衡 UOT针对质量不等工程上有两条主流路径。路径 A 是质量缩放模型。保留行和约束T * 1 p把列和约束改成T * 1 lambda * q其中lambda sum(p) / sum(q)。含义是把 q 的总质量按比例折算到与 p 相等再求解最优搬运方案。这个模型的优点是直观、参数少物理意义清晰q 的质量按比例缩小或放大搬运成本随 p 的总质量线性变化。我一般优先用这个模型因为它的结果可解释、可复现而且 linprog 和 Sinkhorn 都能直接落地。路径 B 是真正意义上的 Unbalanced Optimal Transport非平衡最优传输。它允许耦合矩阵 T 的行和、列和同时偏离 p 和 q通过在目标函数里加 KL 散度惩罚项来控制偏离程度min C, T rho1 * KL(T*1, p) rho2 * KL(T*1, q)当 rho1、rho2 趋向无穷大时模型退化为标准 EMDrho 越小允许的质量创建或销毁越多。适合那些两个分布应该部分匹配但有大量异常点或背景噪声的任务比如目标检测里的匹配、基因表达数据的比对。代价是多了两个需要调的惩罚系数调参成本高。2.3 Matlab 选型linprog、网络单纯形、Sinkhorn 怎么选算法路径适用规模精度工程成本linprognP*nQ 小于 1e5 左右精确最低几十行代码网络单纯形 / mex 封装nP*nQ 在 1e5 到 1e7精确需要调用成熟实现别自己造轮子Sinkhorn 迭代nP*nQ 大于 1e7 或需要 GPU近似误差受正则化参数控制中等需要调 eps 和迭代次数选型的核心判断依据是决策变量个数也就是 nP 乘 nQ。nP 和 nQ 各 300 时变量数 9 万linprog 还能在几秒内求解各 1000 时变量数到 100 万linprog 的约束矩阵装进去就几百 MB内存和耗时同时失控。这时候要么换网络单纯形——它利用图结构跳过大量冗余计算要么换 Sinkhorn 做近似。Sinkhorn 每轮迭代只做两次矩阵乘法复杂度 O(nP * nQ)完全避开线性规划是最大规模问题的事实标准。3. 用 linprog 从零实现非标准化 EMD完整函数与参数细节3.1 最小可运行实现质量缩放模型下面这个函数保存为emd_unnorm_lp.m是质量缩放模型的最小实现。输入 p、q 不要求质量相等代码内部自动处理。function [dist, T, lambda] emd_unnorm_lp(p, q, C) % 非标准化Wasserstein距离质量缩放模型 % p: nP×1 源分布列和不需要为1 % q: nQ×1 目标分布列和不需要为1 % C: nP×nQ 成本矩阵C(i,j) 单位质量从p(i)搬到q(j)的代价 % 输出 dist: 非标准化距离T: 最优耦合矩阵lambda: 质量缩放比 p max(p, 0); % 截断负值避免出现负质量 q max(q, 0); lp sum(p); lq sum(q); if lp 0 || lq 0 error(p和q的总质量必须为正); end lambda lp / lq; % 把q折算到与p总质量一致 nP length(p); nQ length(q); nVar nP * nQ; % 组装等式约束 Aeq * x beqx T(:) 按列优先展开 nCon nP nQ; Aeq sparse(nCon, nVar); for i 1:nP % 行和约束对第i行所有列位置在向量化后索引为 i:nP:end Aeq(i, i:nP:end) 1; end for j 1:nQ % 列和约束第j列占据向量化后的一段连续区间 Aeq(nP j, (j-1)*nP1 : j*nP) 1; end beq [p(:); lambda * q(:)]; % 用对偶单纯形法对中等规模稀疏LP更稳 opts optimoptions(linprog, Algorithm, dual-simplex, ... Display, off, MaxIterations, 5000); [x, fval, exitflag] linprog(C(:), [], [], Aeq, beq, ... zeros(nVar, 1), [], opts); if exitflag 0 warning(linprog未收敛exitflag%d请检查成本矩阵和约束, exitflag); end T reshape(x, nP, nQ); dist fval; end这段代码的核心逻辑是把双约束线性规划翻译成 linprog 的标准接口。f 向量是 C 按列展平的结果C(:)和 x 的排列方式一一对应不等式约束留空用[]等式约束交给 Aeq、beq决策变量下界是全 0 向量上界留空。dual-simplex算法对稀疏约束矩阵更稳健比默认的 interior-point 在变量数上万时收敛更快。3.2 变量排列与约束组装的细节最容易写错的就是 Aeq 的索引。Matlab 的T(:)按列优先展平T(i,j)在 x 中的位置是(j-1)*nP i。这意味着行和约束sum_j T(i,j) p(i)固定 i所有 j 对应的索引是i, inP, i2*nP, ...也就是i:nP:end。列和约束sum_i T(i,j) lambda*q(j)固定 j所有 i 对应的索引是(j-1)*nP1 : j*nP是一段连续区间。很多实现把这两个搞反导致行和、列和互换结果 T 的行和等于 lambda*q、列和等于 p距离数值完全错误。组装完之后建议立刻自检assert(norm(sum(T,2) - p, inf) 1e-8, 行和不匹配); assert(norm(sum(T,1) - lambda*q, inf) 1e-8, 列和不匹配);这个自检不仅在开发期有效上线后每次调用也值得保留。我见过太多结果看着差不多的翻车实际是最优解根本没落在可行域上。3.3 linprog 参数设置与数值预处理linprog 有四个参数值得手动调参数建议值作用Algorithmdual-simplex稀疏、中等规模的LP收敛稳定内点法内存高MaxIterations5000 到 20000默认值对大LP可能不够但太小会误报未收敛ConstraintTolerance1e-8 到 1e-9质量比很悬殊时收紧避免约束被容忍度吃掉OptimalityTolerance1e-8距离精度要求高时收紧预处理三招能显著提升数值稳定性。第一p、q 里小于 1e-12 的 bin 直接置 0减少有效变量第二C 整体除以max(C(:))缩放到 [0,1]避免成本量纲跨数量级导致对偶单纯形迭代病态算出的 dist 记得乘回max(C(:))恢复真实量纲第三C 里出现负值时所有元素加一个平移量让最小值变成 0平移量会在后面的 Sinkhorn 章节详细讲恢复公式。4. 大规模场景的非标准化 Sinkhorn 实现从迭代到参数4.1 熵正则化的思想与标准 Sinkhorn 迭代线性规划求 EMD 的瓶颈在于约束优化问题本身。Sinkhorn 的思路是给目标函数加一个熵正则项eps * sum(T .* log(T))让问题变成强凸的然后通过坐标交替更新求解。迭代式非常简洁u p ./ (K * v) v q ./ (K * u)其中 K 是核矩阵K exp(-C / eps)eps 是正则化强度。迭代收敛后最优耦合T diag(u) * K * diag(v)。这里的 eps 控制温度eps 越大T 越平滑、迭代越快但距离偏离真实 EMD 越多eps 越小T 越接近精确最优传输但 K 的数值极不稳定。工程上这是一个必须显式处理的精度-速度权衡。4.2 非标准化质量缩放的 Sinkhorn 代码质量缩放模型的 Sinkhorn 实现分两步先按比例把 q 折算到与 p 质量一致再跑标准 Sinkhorn最后把距离乘以 p 的总质量。这里的关键是记住——熵正则化在单位质量分布上求解恢复时要乘回总质量。function dist emd_unnorm_sinkhorn(p, q, C, eps, tol, maxIter) % 非标准化Wasserstein距离Sinkhorn近似 % p: nP×1q: nQ×1C: nP×nQ成本矩阵非负 % eps: 正则化系数C缩放到[0,1]后建议0.01~0.5 % tol: 残差收敛阈值建议1e-9 % maxIter: 最大迭代次数建议1e4~1e5 if nargin 6, maxIter 50000; end if nargin 5, tol 1e-9; end if nargin 4, eps 0.05; end p max(p, 0); q max(q, 0); lp sum(p); lq sum(q); pn p / lp; qn q / lq; % 单位质量化质量信息保留在lp中 K exp(-C / eps); % 核矩阵C需先缩放和去负值 u ones(length(pn), 1); v ones(length(qn), 1); for it 1:maxIter u pn ./ (K * v); v qn ./ (K * u); % 残差检查两个边际都必须收敛 err max(norm(u .* (K * v) - pn, inf), ... norm(v .* (K * u) - qn, inf)); if err tol, break; end end T (u * v) .* K; % 外积生成耦合矩阵 Wn sum(T .* C, all); % 单位质量平均搬运成本 dist lp * Wn; % 乘回源分布总质量 end注意T (u * v) .* Ku 是 nP×1v 是 nQ×1外积自动广播成 nP×nQ 矩阵元素 u(i)*v(j)再逐点乘 K。这个写法比diag(u)*K*diag(v)省一次稀疏矩阵乘法对大规模友好。残差检查里的u .* (K*v)是当前迭代下 T 的行和和 pn 的无穷范数差就是收敛判据。all参数需要 R2018b 之后的 Matlab老版本改写成sum(sum(T .* C))。4.3 走向非平衡 UOTKL 惩罚下的扩展迭代如果任务需要真正的非平衡最优传输即允许 T 的行和、列和同时偏离 p、q可以在熵正则化基础上叠加两个 KL 惩罚。迭代式变成带幂指数的形式% 非平衡Sinkhorn迭代片段KL-UOT % rho1、rho2: KL惩罚强度越大越接近硬约束 alpha1 rho1 / (rho1 eps); alpha2 rho2 / (rho2 eps); u ones(length(p), 1); v ones(length(q), 1); for it 1:maxIter Ku K * u; v (q ./ max(Ku, eps)) .^ alpha2; % max防除零 Kv K * v; u (p ./ max(Kv, eps)) .^ alpha1; end T (u * v) .* K; dist sum(T .* C, all) ... rho1 * sum(T * ones(length(q),1) .* log(...)) ...这套迭代的本质是把边际约束从必须等于替换成尽量接近。rho 越大T1 越贴近 prho 越小质量可以大比例创建或销毁。注意q ./ Ku遇到 q 中为 0 的 bin 会除零我的做法是分母加eps兜底但更干净的是在迭代前把全为 0 的 bin 直接删掉恢复输出时再补回对应行列。调参建议先固定 eps 在 0.05 左右然后从 rho10 开始往下降观察sum(T(:))偏离sum(p)的程度偏离超过 20% 基本说明 rho 太松了。5. 非标准化 EMD 的 5 个翻车点与排查方法5.1 先归一化再算 EMD丢掉质量信息的经典翻车现象同一个分布形状把 q 放大 10 倍再算距离和原来一模一样。原因用了类似p p / sum(p); q q / sum(q)的预处理EMD 变成了纯粹的形状差异度量质量信息被彻底丢弃。解决先算lp sum(p)把归一化后的标准距离结果乘回lp。我常用的写法是把归一化的代码封装成emd_norm再在外面包一层乘回总质量的壳子禁止在业务代码里裸写归一化。自检方法把 q 缩放 5 倍如果距离不变说明归一化路径在起作用。5.2 约束矩阵索引错误行列和互换现象算出的距离数值明显不合理比如非负成本出现了比上界还大的距离或者最优耦合 T 的行和与 p 对不上。原因T(:)是列优先展平行约束的索引i:nP:end和列约束的连续段(j-1)*nP1:j*nP容易写反结果 T 的行和等于 lambda*q 而不是 p。解决约束组装完后立刻做断言自检把第 3.2 节的两行 assert 放进代码最前面。我第一次写这个函数时索引反了距离值差异不大、T 看起来也是对称的但行和、列和的数值完全错位花了半天才定位到索引问题。5.3 epsilon 太小导致 NaN数值稳定性的坎现象Sinkhorn 迭代在几十轮内全部变成 NaN或者 K 矩阵里出现 Inf。原因K exp(-C/eps)当 C 的最大值接近 1、eps 取 0.001 时exp(-1000)直接下溢为 0对应位置的p ./ (K*v)就是 0 除 0。解决先把 C 缩放到 [0,1]再取 eps 不低于 0.01如果 eps 必须小改用 log-domain 实现——把 C 减去最小值再做核矩阵代价恢复时补回来。具体来说C_shift C - min(C(:))算出的单位质量成本 Wn 要加回min(C(:))最终dist lp * (Wn min(C(:)))因为每个单位质量的成本都被整体平移了一个常数。5.4 大规模 LP 内存爆炸怎么判断该换算法现象nPnQ2000 时linprog 直接卡死或报内存不足。原因决策变量是 nP*nQ 400 万Aeq 虽然是稀疏的但 linprog 的对偶单纯形内部需要维护大规模基矩阵内存和耗时会失控。解决变量数超过 10 万就放弃 linprog转网络单纯形或 Sinkhorn。网络单纯形利用图结构内存占用低一到两个数量级但工程上不建议自己写找成熟的 mex 封装实现更稳。如果追求可控的近似解直接切到第 4 章的 Sinkhorn每轮迭代只有两次矩阵乘内存瓶颈只在 C 和 K 两个 nP×nQ 矩阵本体。5.5 脏数据与负成本距离变负或出现 NaN现象不设任何预处理时距离输出 NaN 甚至负数。原因p 或 q 里有负值比如做了差分留下的残差或者 C 本身是相似度矩阵直接当成本用相似度越高 C 越小目标函数开始追求负成本。解决入口处强制p max(p,0)、q max(q,0)C 如果是相似度越大越好先做C -C再按成本流程走必要时继续平移保证非负。另一个隐蔽来源是 p 里出现 NaN这会在约束 beq 里扩散linprog 直接退出。我习惯在每个输入上做any(isnan(p))检查宁可报错也不带脏数据进求解器。6. 验证非标准化 EMD 代码三个可复跑的检查6.1 一维解析解对照最优传输在一维情形下有闭式解单位质量的 W1 距离等于累积分布差的绝对值积分。用这个性质对照 LP 结果是验证代码正确性最有效的手段。dx 0.01; x (0:dx:1); p exp(-((x - 0.25) / 0.1).^2); q exp(-((x - 0.75) / 0.1).^2); C abs(x - x); % 一维测度 dist emd_unnorm_lp(p, q, C); pn p / sum(p); qn q / sum(q); ref sum(abs(cumsum(pn) - cumsum(qn))) * dx; % 闭式解 assert(abs(dist / sum(p) - ref) 1e-3);这里dist / sum(p)是单位质量平均成本必须和闭式解一致。如果断言不通过优先检查 C 的方向和 Aeq 的索引而不是怀疑闭式解——我栽在这个选择上过一次。6.2 手算 delta 分布测试构造两个极端简单的分布手算出精确期望值p [1; 0]; q [0; 2]; C [0 1; 1 0]; dist emd_unnorm_lp(p, q, C); % 期望lambda 0.5列和约束为[0;1] % 唯一可行耦合 T(1,2) 1距离 1*1 1这个例子虽然只有两个 bin却能同时验证质量缩放逻辑和最优耦合方向。测试通过后再跑一个相同位置但是 q 总质量不同的用例因为 p、q 都落在同一个 bin 上距离必须输出 0否则预处理里的质量缩放就写错了。6.3 质量敏感性测试距离必须随总质量变化非标准化距离的核心特征是质量敏感。写一个脚本确认这一点p [1; 0]; q [0; 1]; C [0 1; 1 0]; d1 emd_unnorm_lp(p, q, C); % 期望1 p2 [10; 0]; d2 emd_unnorm_lp(p2, q, C); % 期望10总质量线性放大 assert(abs(d2 - 10*d1) 1e-10);这个测试的回归意义在于以后有人把归一化逻辑混进代码时质量敏感性测试会立刻失败。这也是我习惯把它留在单元测试集里而不仅仅作为一次性脚本的原因。回头讲一下我踩过的教训第一次实现 Sinkhorn 版本时我忘了把结果乘回 lp所有测试都看起来合理直到质量敏感性测试把 10 倍质量算出 1 倍距离才暴露问题。从那以后任何 EMD 相关代码的验收清单里永远有一行——质量翻倍非标准化距离必须翻倍。希望帮到你。本文还有配套的精品资源点击获取