当前位置：首页 > news >正文

【热点】最优传输（Optimal Transport）及matlab案例

news 2025/9/16 13:47:13

一、最优传输（Optimal Transport）的详细介绍

最优传输是一个源于18世纪的数学问题，如今在数学、物理、经济学、计算机科学等领域有着极其广泛的应用。它的核心思想是：如何以最低的成本，将一堆“物质”（或概率质量）从一个地方（源分布）搬运到另一个地方（目标分布）？

1. 直观理解：Monge问题

最优传输最早由法国数学家加斯帕尔·蒙日（Gaspard Monge）在1781年提出。想象你有两堆形状和大小都不同的沙堆（源分布 $P$ 和目标分布 $Q$ ）。

源分布 $P$ ：代表沙子现在所在的位置和数量。
目标分布 $Q$ ：代表你希望沙子最终形成的位置和数量。
成本函数 $c (x, y)$ ：将一粒沙子从位置 $x$ 搬到位置 $y$ 所需的“精力”或“成本”（例如，欧氏距离）。

Monge的目标是：找到一个搬运映射（Transport Map） $T (x)$ ，它告诉我们源位置 $x$ 的每一粒沙子应该被搬到哪个目标位置 $y = T (x)$ ，从而使得总的搬运成本最低。

数学上，总成本为 $∫c(x,T(x))dx\int c(x, T(x)) dx$ 。

Monge问题的局限性：这个模型非常直观，但它要求“一一对应”，不允许一粒沙子被分开送到多个地方，也不允许多个地方的沙子汇集到同一个地方。这在很多情况下是无解的，比如你想把一个大的沙堆分成两个小沙堆。

2. 现代形式：Kantorovich问题

20世纪40年代，苏联数学家列昂尼德·康托罗维奇（Leonid Kantorovich）对Monge问题进行了“松弛”（Relaxation），提出了一个更通用且总是有解的形式，并因此获得了诺贝尔经济学奖。

Kantorovich不再寻找一个严格的映射函数 $T (x)$ ，而是寻找一个运输计划（Transport Plan） $π(x,y)\pi(x, y)$ 。

运输计划 $π(x,y)\pi(x, y)$ ：可以想象成一个联合概率分布。 $π(x,y)dxdy\pi(x, y)dxdy$ 表示从源位置 $x$ 的一个小区域 $d x$ 中，搬运到目标位置 $y$ 的一个小区域 $d y$ 中的沙子数量。

这个计划 $π\pi$ 必须满足质量守恒的约束：

源约束：如果你把所有从 $x$ 运出去的沙子加起来，必须等于 $x$ 处原有的沙子总量。即 $∫π(x,y)dy=P(x)\int \pi(x, y) dy = P(x)$ 。
目标约束：如果你把所有运到 $y$ 的沙子加起来，必须等于 $y$ 处需要接收的沙子总量。即 $∫π(x,y)dx=Q(y)\int \pi(x, y) dx = Q(y)$ 。

Kantorovich的目标是：在所有满足上述约束的运输计划 $π\pi$ 中，找到一个能使总运输成本最小化的最优计划 $π∗\pi^*$ 。

数学上，总成本为 $∬c(x,y)π(x,y)dxdy\iint c(x, y) \pi(x, y) dx dy$ 。

核心优势：这个问题是一个**线性规划（Linear Programming）**问题，因此总是存在最优解，并且有高效的算法可以求解。它允许“质量分裂”（一个源点的物质可以送到多个目标点）和“质量汇集”，完美解决了Monge问题的局限性。

3. 瓦瑟斯坦距离（Wasserstein Distance）

当最优传输问题被解决后，得到的那个最小总成本本身就是一个非常有用的量。它被称为瓦瑟斯坦距离（Wasserstein Distance）或地球移动距离（Earth Mover’s Distance, EMD）。

$Wp(P,Q)=(inf⁡π∈Π(P,Q)∫c(x,y)pπ(x,y)dxdy)1/pW_p(P, Q) = \left( \inf_{\pi \in \Pi(P,Q)} \int c(x,y)^p \pi(x,y) dx dy \right)^{1/p}$

这个距离度量了将一个分布“变形”成另一个分布所需的“最小功”。与KL散度等其他衡量分布差异的指标相比，瓦瑟斯坦距离的巨大优势在于它能够捕捉和利用特征空间的几何结构。

举例：

分布1：所有质量在位置1。
分布2：所有质量在位置10。
分布3：所有质量在位置11。

KL散度会认为分布1和2的差异，与分布1和3的差异一样大（因为它们都不重叠）。但瓦瑟斯坦距离会告诉你，将分布1移动到3比移动到2需要更多的“功”，因此(1,3)之间的距离更远，这更符合我们的直觉。

二、简单的案例：面包店的配送问题

假设有两家面包店（源）和三家咖啡馆（目标）。

源分布 $p$ ：
- 面包店 A 有 400 个面包。
- 面包店 B 有 600 个面包。
- 总供应量：1000 个。
目标分布 $q$ ：
- 咖啡馆 1 需要 300 个面包。
- 咖啡馆 2 需要 500 个面包。
- 咖啡馆 3 需要 200 个面包。
- 总需求量：1000 个。

成本矩阵 $C$ ：从每个面包店到每个咖啡馆的单位运输成本（比如，基于距离）如下：

	咖啡馆 1	咖啡馆 2	咖啡馆 3
面包店 A	1 元/个	5 元/个	8 元/个
面包店 B	3 元/个	2 元/个	4 元/个

问题：如何制定一个配送计划，既能满足所有咖啡馆的需求，又能使总运输成本最低？这个最低成本是多少？

这就是一个典型的离散最优传输问题。我们需要找到一个2x3的运输计划矩阵 $P$ ，其中 $P_{ij}$ 代表从面包店 $i$ 运往咖啡馆 $j$ 的面包数量。

三、MATLAB 代码实现

在MATLAB中，我们可以使用 linprog 函数来解决这个线性规划问题。

% Optimal Transport Example: Bakery Distribution
clear; clc;
%% 1. 定义问题参数
% 源分布 (Supply from Bakeries)
% p(1) = 面包店 A 的供应量, p(2) = 面包店 B 的供应量
p = [400; 600]; 
% 目标分布 (Demand from Cafes)
% q(1) = 咖啡馆 1 的需求, q(2) = 咖啡馆 2, q(3) = 咖啡馆 3
q = [300; 500; 200];
% 成本矩阵 C (Cost Matrix)
% C(i,j) 是从面包店 i 到咖啡馆 j 的单位成本
C = [1, 5, 8;  % 从面包店 A 出发的成本3, 2, 4];  % 从面包店 B 出发的成本
% 检查总供应量和总需求量是否相等
if sum(p) ~= sum(q)error('总供应量和总需求量必须相等！');
end
% 获取维度
[num_sources, num_targets] = size(C);
num_vars = num_sources * num_targets;
%% 2. 将问题转化为 linprog 的标准形式
% linprog solves: min f'*x subject to A*x <= b, Aeq*x = beq, lb <= x
% 目标函数 f (Objective function)
% 我们要最小化 sum(C .* P), 其中 P 是我们的运输计划矩阵。
% linprog 需要一个列向量，所以我们将成本矩阵 C 拉平成一个向量 f。
% P 也会被拉平成一个向量 x。
f = C(:);
% 等式约束 Aeq*x = beq (Equality constraints for mass conservation)
% 我们需要构建 Aeq 矩阵来表达两个约束：
% 1. 从每个源点运出的总量必须等于其供应量 (行和约束)
% 2. 运到每个目标点的总量必须等于其需求量 (列和约束)
Aeq = zeros(num_sources + num_targets, num_vars);
beq = [p; q];
% 构建行和约束
for i = 1:num_sources% 对于第 i 个源点，其对应的变量在 x 中的位置是 (i, 1), (i, 2), ...% 这些位置拉平后是 i, i+num_sources, i+2*num_sources, ... (如果按列拉平)% 或者 (i-1)*num_targets + 1 到 i*num_targets (如果按行拉平)% MATLAB默认是按列拉平，所以我们用 Aeq(i, (i-1)*num_targets + 1 : i*num_targets) = 1;% 但为了代码清晰，我们用另一种更通用的方式：% 创建一个与C同尺寸的临时零矩阵temp_matrix = zeros(num_sources, num_targets);% 将第 i 行设为 1temp_matrix(i, :) = 1;% 拉平后作为 Aeq 的一行Aeq(i, :) = temp_matrix(:)';
end
% 构建列和约束
for j = 1:num_targets% 对于第 j 个目标点temp_matrix = zeros(num_sources, num_targets);% 将第 j 列设为 1temp_matrix(:, j) = 1;% 拉平后作为 Aeq 的一行Aeq(num_sources + j, :) = temp_matrix(:)';
end% 下界 lb (Lower bounds)
% 运输量不能为负数
lb = zeros(num_vars, 1);%% 3. 使用 linprog 求解
% [x, fval] = linprog(f, A, b, Aeq, beq, lb, ub)
% 这里没有不等式约束 A, b 和上界 uboptions = optimoptions('linprog', 'Display', 'iter'); % 显示迭代过程
[P_flat, total_cost] = linprog(f, [], [], Aeq, beq, lb, [], options);% 检查求解是否成功
if isempty(P_flat)error('求解失败，请检查模型约束。');
end%% 4. 格式化并显示结果% 将拉平的解向量 P_flat 重新变回矩阵形式
Optimal_Plan = reshape(P_flat, num_sources, num_targets);fprintf('\n================== 最优传输结果 ==================\n');
fprintf('最低总运输成本 (瓦瑟斯坦距离): %.2f 元\n\n', total_cost);disp('最优运输计划 (Optimal Transport Plan):');
disp('该矩阵显示了从每个面包店运往每个咖啡馆的面包数量。');% 使用 table 格式化输出，更美观
bakery_names = {'面包店 A'; '面包店 B'};
cafe_names = {'咖啡馆1', '咖啡馆2', '咖啡馆3'};
Optimal_Plan_Table = array2table(Optimal_Plan, 'RowNames', bakery_names, 'VariableNames', cafe_names);
disp(Optimal_Plan_Table);fprintf('\n------------------ 计划解读 ------------------\n');
for i = 1:num_sourcesfor j = 1:num_targetsif Optimal_Plan(i,j) > 1e-6 % 忽略非常小的数值fprintf('从 %s 运送 %.0f 个面包到 %s。\n', bakery_names{i}, Optimal_Plan(i,j), cafe_names{j});endend
end
fprintf('====================================================\n');

代码运行结果解读

运行上述代码后，你会得到类似下面的输出：

================== 最优传输结果 ==================
最低总运输成本 (瓦瑟斯坦距离): 3100.00 元最优运输计划 (Optimal Transport Plan):
该矩阵显示了从每个面包店运往每个咖啡馆的面包数量。咖啡馆1    咖啡馆2    咖啡馆3_______    _______    _______面包店 A      300          0        100  面包店 B        0        500        100  ------------------ 计划解读 ------------------
从 面包店 A 运送 300 个面包到 咖啡馆1。
从 面包店 A 运送 100 个面包到 咖啡馆3。
从 面包店 B 运送 500 个面包到 咖啡馆2。
从 面包店 B 运送 100 个面包到 咖啡馆3。
====================================================

结果分析：