最优传输逆问题 [inverse-optimal-transport]

2026-10-31
笔记
最优传输
从带噪传输观测到贝叶斯推断、MAP 与 MCMC 的逆最优传输笔记

本篇是基于论文Inverse optimal transport的阅读笔记，结合了论文内容、我自身的理解和 GPT 的解释。

最优传输是在边缘概率分布的约束下，寻找使总运输成本最小的传输计划。就像其他逆问题一样：有了边缘约束和成本矩阵，我们可以通过最优传输算法得到最优传输矩阵；那么，有了最优传输矩阵，我们能否反过来得到对应的成本矩阵呢？这是我们自然会想知道的问题。

另外，为什么这个问题有意义呢？从我的理解来看，成本矩阵包含源与目的地之间的传输代价信息，可能反映物理、经济等机制。因此，反推成本有助于理解产生传输行为的机制。下面是问题的正式表述。

一、最优传输的正问题

给定两个边缘分布 $p,q$ 和成本矩阵 $C$ ，最优传输（Optimal Transport，OT）的目标是求出最优传输计划 $T^*$ 。

主要符号如下：

符号	含义
$p\in\Delta_n,\ q\in\Delta_m$	两个分别为 $n$ 、 $m$ 维的离散概率分布（边缘分布）
$C\in\mathbb R^{n\times m}$	成本矩阵， $C_{ij}$ 表示从位置 $i$ 到位置 $j$ 的运输成本
$T\in\mathbb R_+^{n\times m}$	传输计划矩阵， $T_{ij}$ 表示从 $i$ 向 $j$ 传输的质量
$\mathbf 1_n,\mathbf 1_m$	全 1 的 $n$ 、 $m$ 维列向量
$T^*$	最优传输计划

这里 $\mathbb R_+$ 包括零， $p,q$ 均为列向量； $T_{ij}$ 也可理解为源位置为 $i$ 、目的地为 $j$ 的联合概率。

在离散情况下，概率单纯形可以表示为

\Delta_n=\left\{x\in\mathbb R_+^n:\sum_{i=1}^n x_i=1\right\}.

满足边缘约束的传输计划集合可记为

\Pi(p,q)=\left\{T\in\mathbb R_+^{n\times m}:T\mathbf 1_m=p,\quad T^\top\mathbf 1_n=q\right\}.

那标准离散最优传输问题可以写成

\begin{aligned} T^* &\in\operatorname*{arg\,min}_{T\in\Pi(p,q)}\langle C,T\rangle\\ &=\operatorname*{arg\,min}_{T\in\Pi(p,q)}\sum_{i,j}C_{ij}T_{ij}. \end{aligned}

最优传输问题（康托维奇形式），其实就是用各种方法求出这个计划。计划有了，再计算内积 $\langle C,T^*\rangle$ （逐元素相乘再求和），就能得到最优传输值。

从线性规划的角度看，这个问题可能有唯一解，也可能有多个最优解。 $pq^\top$ 总是可行计划；可行集闭且有界，线性目标因此能取到最小值。如果有两个不同的最优解，由于约束和目标都是线性的，它们的任意凸组合也都是最优解，因此会有无穷多个最优解。

那么，现在可以将求解最优传输计划的过程抽象为一个映射：

T^*=f(p,q,C).

需要注意：最优传输计划不一定唯一。 为了将求解过程写成单值函数 $f$ ，可以固定一个确定性的求解算法以及处理多解的规则，由它选出一个确定的 $T^*$ 。

二、最优传输的逆问题

正问题是从 $(p,q,C)$ 推出 $T^*$ ；逆问题则是：

已知观测到的传输矩阵 $T^{\mathrm{obs}}$ ，反过来推断可能产生它的 $p,q,C$ 。

如果观察到的是没有噪声的真实传输计划，那么边缘分布应满足

T^{\mathrm{obs}}\mathbf 1_m=p, \qquad (T^{\mathrm{obs}})^\top\mathbf 1_n=q.

一种直接思路是：

将 $T^{\mathrm{obs}}$ 的行和、列和当作真实的边缘分布 $p,q$ ；
再通过某种优化方法尝试反推成本矩阵 $C$ 。

最直接的方法是把成本参数限制在有限网格上，遍历候选值，并比较对应的最优传输计划。不过，连续参数不能直接穷举，而且找到相符的成本，也不代表它唯一。

例如，取 $p=q=(1/2,1/2)^\top$ ，成本矩阵

C=\begin{pmatrix}0&1\\1&0\end{pmatrix} \quad\text{与}\quad 2C

都产生同一个唯一最优计划 $T^*=\operatorname{diag}(1/2,1/2)$ ：成本为零要求两个非对角元素都为零，边缘约束再确定两个对角元素；正比例放大成本不改变最小点。这已经说明，仅凭传输计划通常无法唯一恢复成本。

但现实中的 $T^{\mathrm{obs}}$ 往往包含观测噪声，因此由它直接计算出的边缘分布也可能受到噪声影响，无法准确反映真实情况。

例如，在后面的独立高斯噪声假设下，行和、列和的误差分别满足

\sum_{j=1}^m\eta_{ij}\sim\mathcal N(0,m\sigma^2), \qquad \sum_{i=1}^n\eta_{ij}\sim\mathcal N(0,n\sigma^2).

也就是说，直接求和得到的边缘虽然无偏，却仍有随机误差；噪声还可能使它们不满足非负性或总和为 1。

那么，自然地，就要考虑带噪声的观测模型，也就是观测值等于真实值加上噪声

T^{\mathrm{obs}}=f(p,q,C)+\eta.

其中 $\eta$ 是观测噪声矩阵，每个元素都服从所设定的分布。带噪观测 $T^{\mathrm{obs}}$ 本身不必是满足约束的传输计划。多解时，这个模型还把 $f$ 的选择规则纳入了假设；若真实系统选了另一个最优计划，其差异也会被归入观测残差。这个设定在信号处理中很常见：观测信号等于真实信号加上噪声。有电子工程背景的学生，多少都听说过高斯白噪声。

基于这个模型，我们就可以开始推断了，但首先得有个抓手。到底怎么来算呢？没有噪声时， $p,q$ 可以由 $T$ 的行和、列和直接得到，只需尝试反推 $C$ ；现在 $p,q,C$ 都未知，问题就更复杂了。

现在，把 $p,q,C$ 合并记为参数 $\mu=(p,q,C)$ （将 $C$ 展平后，可整体视作一个高维向量），那模型就变成了

T^{\mathrm{obs}}=f(\mu)+\eta.

对于矩阵中的每一个元素，模型为

T_{ij}^{\mathrm{obs}}=f_{ij}(\mu)+\eta_{ij}.

$f$ 是已知的求解算法， $\eta$ 是随机噪声， $T^{\mathrm{obs}}$ 是已经观测到的值，于是就会想到贝叶斯的一套工具：用先验分布表达对参数的偏好，再结合观测推断参数。这里真实参数可以是固定但未知的，贝叶斯分布表达的是我们对它的不确定性。那这个分布怎么建模呢？下一小节来阐述这个问题。

三、高斯噪声假设与似然函数

假设观测噪声独立于参数，各个矩阵元素的噪声相互独立，且服从相同方差的高斯分布，其中 $\sigma>0$ 已知且固定。这是信号处理中常见的建模假设；许多微小扰动叠加时，可在适当条件下由中心极限定理得到高斯近似，但独立、同方差仍需另作假设：

\eta_{ij}\overset{\mathrm{i.i.d.}}{\sim}\mathcal N(0,\sigma^2).

于是，在给定参数 $\mu$ 时，

T_{ij}^{\mathrm{obs}}\mid\mu \sim \mathcal N\bigl(f_{ij}(\mu),\sigma^2\bigr).

这个其实很简单：确定参数 $\mu$ 后， $f(\mu)$ 的值也就确定了，而观测矩阵的每个元素都在对应的 $f_{ij}(\mu)$ 附近服从上述高斯分布。实际推断时，观测值固定，似然则看作参数 $\mu$ 的函数。

根据高斯分布的定义，单个矩阵元素的条件概率密度为

p(T_{ij}^{\mathrm{obs}}\mid\mu) =\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left[ -\frac{\bigl(T_{ij}^{\mathrm{obs}}-f_{ij}(\mu)\bigr)^2}{2\sigma^2} \right].

如果只关注与 $\mu$ 有关的部分，则可以写为

p(T_{ij}^{\mathrm{obs}}\mid\mu) \propto \exp\left[ -\frac{\bigl(T_{ij}^{\mathrm{obs}}-f_{ij}(\mu)\bigr)^2}{2\sigma^2} \right].

由于假设各元素的噪声独立，整个矩阵的条件密度等于各元素条件密度的乘积：

\begin{aligned} p(T^{\mathrm{obs}}\mid\mu) &=\prod_{i=1}^n\prod_{j=1}^m p(T_{ij}^{\mathrm{obs}}\mid\mu)\\[4pt] &\propto\prod_{i=1}^n\prod_{j=1}^m \exp\left[ -\frac{\bigl(T_{ij}^{\mathrm{obs}}-f_{ij}(\mu)\bigr)^2}{2\sigma^2} \right]\\[4pt] &=\exp\left[ -\frac{1}{2\sigma^2}\sum_{i,j} \bigl(T_{ij}^{\mathrm{obs}}-f_{ij}(\mu)\bigr)^2 \right]. \end{aligned}

把观测值固定后，这个概率密度作为参数的函数，就是似然函数，衡量给定参数对当前观测的支持程度。

定义

\Phi(\mu;T^{\mathrm{obs}}) =\frac{1}{2\sigma^2} \left\|T^{\mathrm{obs}}-f(\mu)\right\|_F^2

其中 $\|\cdot\|_F$ 表示矩阵的 Frobenius 范数，即 $\|A\|_F^2=\sum_{i,j}A_{ij}^2$ 。

于是似然函数可简写为

\mathcal L(\mu) =p(T^{\mathrm{obs}}\mid\mu) \propto\exp\bigl[-\Phi(\mu;T^{\mathrm{obs}})\bigr].

预测传输矩阵 $f(\mu)$ 越接近观测矩阵 $T^{\mathrm{obs}}$ ，数据失配 $\Phi$ 越小，对应参数的似然就越大。

四、贝叶斯推断：从似然到后验分布

根据贝叶斯公式，

\pi(\mu\mid T^{\mathrm{obs}}) =\frac{ p(T^{\mathrm{obs}}\mid\mu)\,\pi(\mu) }{p(T^{\mathrm{obs}})}.

其中：

$\pi(\mu\mid T^{\mathrm{obs}})$ ：后验分布；
$p(T^{\mathrm{obs}}\mid\mu)$ ：似然函数；
$\pi(\mu)$ ：先验分布；
$p(T^{\mathrm{obs}})$ ：使后验密度积分为 1 的归一化常数（边际似然）。

这个公式的意思就是：结合先验分布和似然函数，推断在当前观测下，哪些参数更合理。

那来看看我们现在有什么：后验分布的表达式有了，但设置先验时还需要处理参数约束，尤其是 $p,q$ 的非负性与归一化条件。当然，可以直接为它们设置狄利克雷先验；这里稍微处理一下，改用隐变量参数化，便于统一处理约束。

引入隐变量 $u,v,w$ （分别为 $n$ 、 $m$ 维列向量和 $n\times m$ 矩阵），用映射表示原始参数：

p=M_n(u),\qquad q=M_m(v),\qquad C=M_{n\times m}(w).

这里的 $M$ 是参数化或归一化映射，目的是让相应参数满足所设定的约束。例如，对非负且不全为零的向量 $x\in\mathbb R_+^r$ ，可取

M_r(x)=\frac{x}{\sum_{a=1}^r x_a},

使边缘分布自动满足约束； $M_{n\times m}$ 则按具体的成本约束选择。原论文将非负且非零的成本除以全部元素之和；这里保留一般写法，实际计算前还须明确成本映射及参数域。零向量处可另行约定映射值，在下面的连续均匀先验下该情形的概率为零。

定义

\theta=(u,v,w),

以及复合映射

G(u,v,w) =f\bigl(M_n(u),M_m(v),M_{n\times m}(w)\bigr).

于是观测模型简化为

T^{\mathrm{obs}}=G(\theta)+\eta.

作为一种简单的先验选择，可以在有界参数区域 $\mathcal U$ 上对隐变量设置均匀先验。需要注意，隐变量均匀一般不意味着映射后的 $p,q,C$ 均匀。这里直接指定 $\theta$ 的先验，再由映射诱导 $\mu$ 的先验，并非保留某个既定 $\pi(\mu)$ 的密度换元，因此下文不需另加雅可比因子。

例如，当 $\mathcal U=[0,1]^d$ 时，

\theta\sim\mathrm{Uniform}(\mathcal U), \qquad \pi(\theta)=\mathbf 1_{\mathcal U}(\theta),

其中 $\mathbf 1_{\mathcal U}$ 是参数域的指示函数：

\mathbf 1_{\mathcal U}(\theta)= \begin{cases} 1,&\theta\in\mathcal U,\\ 0,&\theta\notin\mathcal U. \end{cases}

这里 $d$ 是隐变量 $\theta=(u,v,w)$ 的总维数；若 $u,v,w$ 的形状分别为 $n$ 、 $m$ 、 $n\times m$ ，则 $d=n+m+nm$ 。

之前推导的似然函数也不白费：只是把 $f(\mu)$ 换成 $G(\theta)$ ，形式仍然一样，而先验更容易处理。相应地，将失配函数写成

\Phi(\theta;T^{\mathrm{obs}}) =\frac{\|T^{\mathrm{obs}}-G(\theta)\|_F^2}{2\sigma^2}, \qquad p(T^{\mathrm{obs}}\mid\theta) \propto\exp\bigl[-\Phi(\theta;T^{\mathrm{obs}})\bigr].

这里沿用 $\Phi$ 表示参数化后的失配函数。

由似然和均匀先验可得

\pi(\theta\mid T^{\mathrm{obs}}) \propto \exp\bigl[-\Phi(\theta;T^{\mathrm{obs}})\bigr] \mathbf 1_{\mathcal U}(\theta).

将归一化常数记为 $Z$ ：

\pi(\theta\mid T^{\mathrm{obs}}) =\frac{1}{Z} \exp\bigl[-\Phi(\theta;T^{\mathrm{obs}})\bigr] \mathbf 1_{\mathcal U}(\theta).

代入 $\theta=(u,v,w)$ ，有

\pi(u,v,w\mid T^{\mathrm{obs}}) =\frac{1}{Z} \exp\left[ -\frac{\|T^{\mathrm{obs}}-G(u,v,w)\|_F^2}{2\sigma^2} \right] \mathbf 1_{\mathcal U}(u,v,w).

其归一化常数为

Z=\int_{\mathcal U} \exp\left[ -\frac{\|T^{\mathrm{obs}}-G(u,v,w)\|_F^2}{2\sigma^2} \right] \,du\,dv\,dw.

当参数维数较高，或者每次求值 $G(\theta)$ 的代价较大时，直接计算 $Z$ 通常很困难，因为高维数值积分可能需要大量函数求值。因此，我们希望使用不必显式求出 $Z$ 的方法研究后验分布。

五、最大后验估计（MAP）

一种求解方式是只寻找后验密度最大的参数，即最大后验估计（Maximum a Posteriori，MAP）：

\theta_{\mathrm{MAP}} \in\underset{\theta\in\mathcal U}{\operatorname{argmax}} \;\pi(\theta\mid T^{\mathrm{obs}}).

在均匀先验、固定噪声方差的假设下，最大化后验密度等价于最小化失配函数：

\theta_{\mathrm{MAP}} \in\underset{\theta\in\mathcal U}{\operatorname{argmin}} \;\frac{\|T^{\mathrm{obs}}-G(\theta)\|_F^2}{2\sigma^2}.

但 MAP 估计通常只提供一个最优点（或一个最优解集合），无法完整描述参数的不确定性。这里求的是隐变量的 MAP；将它映射回 $\mu$ 可得到一个参数估计，但不保证是 $\mu$ 的 MAP，因为密度最大点依赖参数化。

若 $G$ 可微，可以对 $\Phi$ 求导，再用投影梯度下降等方法，在 $\mathcal U$ 内寻找局部最小点。不过，标准线性规划 OT 的解映射通常不光滑，不能直接假定普通梯度法适用；也可以考虑无梯度优化，或改用适当的光滑正则化正问题。

逆问题可能具有多解性： 不同的 $p,q,C$ （或者不同的隐变量 $\theta$ ）有可能产生相同或非常接近的传输矩阵。因此，即便已经观测到 $T^{\mathrm{obs}}$ ，仍可能有多组合理参数。例如，在标准 OT 中，给成本矩阵的每行或每列分别加一个常数，只会使目标函数增加一个由边缘分布决定的常数，不改变最优解集合：若 $C'_{ij}=C_{ij}+a_i+b_j$ ，则 $\langle C',T\rangle=\langle C,T\rangle+a^\top p+b^\top q$ 。变换后的成本是否仍可作为候选，还取决于所选参数域和映射。

产生相同预测的参数具有相同似然，数据无法区分它们；后验还取决于先验，贝叶斯与 MCMC 并不会自动消除这种不可辨识性。

由此自然产生进一步的问题：

哪些参数组合在观测数据下更合理？
后验分布的整体形状是什么？
参数的后验均值和方差是多少？
由这些参数所生成的传输矩阵，元素级的不确定性如何？

这就引出了 MCMC。

六、MCMC：马尔可夫链蒙特卡洛方法

首先来说一下为什么要用 MCMC。我们已经知道后验密度的未归一化形式：

\pi(\theta\mid T^{\mathrm{obs}}) \propto \exp\bigl[-\Phi(\theta;T^{\mathrm{obs}})\bigr] \mathbf 1_{\mathcal U}(\theta).

之前提到过归一化常数 $Z$ 不容易求，但很多马尔可夫链蒙特卡洛（Markov Chain Monte Carlo，MCMC）方法可以只利用未归一化后验密度构造采样算法，而无需显式计算 $Z$ 。在满足收敛条件并经过充分迭代后，可以得到近似服从后验分布的一组样本（样本之间通常相关）：

\theta^{(1)},\theta^{(2)},\ldots,\theta^{(N)} \quad\text{来自}\quad \pi(\theta\mid T^{\mathrm{obs}}).

例如，可用采样均值近似参数的后验均值：

\widehat{\mathbb E}[\theta\mid T^{\mathrm{obs}}] =\frac{1}{N}\sum_{k=1}^N\theta^{(k)}.

进一步，可以根据样本估计各分量的方差、分位数及可信区间，从而研究参数的不确定性。

那具体怎么做呢？可以用均匀分布生成候选参数，但不能直接把它们当作后验样本，还需要接受或拒绝这一步。例如，在 $\mathcal U$ 上采用独立均匀提议，这是 Metropolis–Hastings 的一个简单例子（原论文主要采用随机游走及其分块更新）。先选取 $\theta^{(0)}\in\mathcal U$ ，计算 $G(\theta^{(0)})$ 和失配；每次计算候选的失配，也都要运行一次 OT 求解器。

从当前参数 $\theta^{(k)}\in\mathcal U$ 出发，提出 $\theta'\sim\mathrm{Uniform}(\mathcal U)$ 。
以概率

\alpha=\min\left\{1, \exp\left[\Phi(\theta^{(k)};T^{\mathrm{obs}}) -\Phi(\theta';T^{\mathrm{obs}})\right]\right\}

接受候选参数：另抽 $r\sim\mathrm{Uniform}(0,1)$ ，若 $r\leq\alpha$ ，令 $\theta^{(k+1)}=\theta'$ ；否则保留当前参数，令 $\theta^{(k+1)}=\theta^{(k)}$ 。拒绝后的重复状态也要计入样本。

重复迭代，舍弃预热阶段并检查收敛后，用保留的样本计算统计量。可比较不同初值的多条链的轨迹与统计量，但诊断通过并不保证已充分探索后验。

这个接受概率来自后验密度之比，归一化常数 $Z$ 正好约掉了；独立均匀提议的正反提议密度相同，也会约掉。独立均匀提议容易理解，但在高维或后验集中时可能效率较低，也可以改用随机游走等提议方式；若提议不对称（越界后重抽通常会如此），接受率还需包含反向与正向提议密度之比。给定当前状态后，下一状态的条件分布不再依赖更早状态，这就是马尔可夫性；利用样本平均近似期望，就是蒙特卡洛。上面的均值、方差估计均使用保留的后验样本。

对每个后验样本，再计算

T^{(k)}=G\bigl(\theta^{(k)}\bigr), \qquad k=1,\ldots,N.

于是，对于传输矩阵中的元素 $(i,j)$ ，其后验均值可近似为

\overline T_{ij} =\frac{1}{N}\sum_{k=1}^N G_{ij}\bigl(\theta^{(k)}\bigr).

相应地，无噪声计划的后验方差可用样本方差估计（不包含新的观测噪声）：

\widehat{\operatorname{Var}}(T_{ij}^*\mid T^{\mathrm{obs}}) =\frac{1}{N-1} \sum_{k=1}^N \left(G_{ij}(\theta^{(k)})-\overline T_{ij}\right)^2.

相关样本下，这个估计不保证有限样本无偏；在链满足遍历条件时，它随样本增多趋近后验方差，不能把它当成样本均值的蒙特卡洛误差。

同样可以利用各样本计算相应的

p^{(k)}=M_n(u^{(k)}),\qquad q^{(k)}=M_m(v^{(k)}),\qquad C^{(k)}=M_{n\times m}(w^{(k)}),

从而分析边缘分布及成本矩阵的后验不确定性。应先映射每个样本再求均值，一般不能把隐变量均值直接映射为原参数均值。

七、整体推导链条

整个思路可以概括为：

正问题： 给定 $p,q,C$ ，由最优传输求解器得到 $T^*=f(p,q,C)$ 。
逆问题： 观察到带噪声的矩阵 $T^{\mathrm{obs}}=T^*+\eta$ ，希望反推未知参数。
参数化： 先将原始参数记为 $\mu=(p,q,C)$ ；再引入隐变量 $\theta=(u,v,w)$ ，用复合映射 $G(\theta)$ 表示模型。
似然建模： 假设独立高斯噪声，得到

p(T^{\mathrm{obs}}\mid\theta)\propto \exp\left(-\frac{\|T^{\mathrm{obs}}-G(\theta)\|_F^2}{2\sigma^2}\right).

贝叶斯推断： 加入先验，得到后验分布 $\pi(\theta\mid T^{\mathrm{obs}})$ 。
MAP： 寻找后验最大的参数组合，但这无法描述整体不确定性。
MCMC： 从后验分布采样，计算参数及传输矩阵的均值、方差等统计量。

简写为：

(p,q,C) \xrightarrow{\;f\;}T^* \xrightarrow{\;+\eta\;}T^{\mathrm{obs}} \xrightarrow{\;\text{Bayes}\;} \pi(\theta\mid T^{\mathrm{obs}}) \xrightarrow{\;\text{MCMC}\;} \text{后验样本与不确定性估计}

参考论文：Inverse Optimal Transport（Stuart 和 Wolfram）。