傳輸理論

傳輸理論（tansport theory、transportation theory），又稱為運輸理論，是數學、經濟學等學科中研究最優運輸和資源配置的理論。該問題最早由法國數學家加斯帕爾·蒙日於1781年提出。

1920年代，A·N·托爾斯泰是最早運用數學方法研究傳輸問題的學者之一。1930年，他在蘇聯國家交通部編纂的《運輸規劃》第一卷中發表了題為《尋找太空貨物運輸的最小千公里方法》的論文。

第二次世界大戰期間，蘇聯數學家、經濟學家列昂尼德·坎托羅維奇在該領域取得了重要進展。因此，這一問題有時也被稱為蒙日-坎托羅維奇運輸問題（Monge–Kantorovich transportation problem）。該問題的線性規劃形式也被稱為弗蘭克·勞倫·希區柯克-庫普曼斯運輸問題。

## 背景

### 礦山與工廠

假設有$m$個開採鐵礦石的礦山以及$n$個工廠使用這些礦山生產的鐵礦石，且這些礦山和工廠構成歐幾里得平面$\mathbb{R}^2$中兩個不相交的子集$M$和$F$。同時假設存在一個成本函數$c : \mathbb{R}^2 \times \mathbb{R}^2 \to [0, \infty)$，其中$c(x, y)$表示將一批礦石從$x$運送到$y$的成本。為簡化起見，此處忽略運輸所需的時間。我們還假定每個礦山只能供應一家工廠（不能拆分運輸），並且每家工廠需要恰好一批貨物才能運營（工廠不能半負荷或雙倍負荷運轉）。基於上述條件，一個傳輸計劃可以看作是一個雙射$T: M \to F$。換句話說，每個礦井$m \in M$僅供應一個目標工廠$T(m) \in F$，而每個工廠也只由一個礦山供貨。我們希望找到最優傳輸計劃$T$，使得總成本

$$c(T) := \sum_{m \in M} c(m, T(m))$$

在所有$M$到$F$的傳輸計劃中是最小的。該問題是傳輸問題的一個特例，可看成一個任務分配問題。更具體地說，它等價於在二分圖中尋找最小權重匹配。

### 書籍移動：成本函數的重要性

下面這個簡單的例子說明了成本函數在確定最優傳輸計劃中的重要性。假設我們有$n$本寬度相等的書擺放在書架上（可以看成具像化的實數線），形成連續一排書。我們希望將它們重新排列，在保持其連續性的同時將整體向右移動一本書的寬度。針對這個問題，有兩個顯而易見的最優傳輸候選方案：

- 將所有$n$本書全都向右移動一本書的寬度（許多小步）。
- 將最左側的書向右移動$n$書本的寬度，其他書則保持不動（一大步）。

如果成本函數與歐幾里得距離成正比（即 $c(x, y) = \alpha \|x - y\|$，其中$\alpha > 0$），那麼這兩種候選方案都是最優的。而如果我們選擇與歐幾里得距離的平方成比例的嚴格凸成本函數（即 $c(x, y) = \alpha \|x - y\|^2$，其中$\alpha > 0$），則「許多小步」的方案則是唯一的最優解。

需要注意的是，上述成本函數僅考慮書籍本身移動的水平距離，而沒有考慮拿起每本書並將其移動到位的設備所行進的水平距離。如果考慮後者，那麼在兩種傳輸計劃中，第二種方案對於歐幾里得距離始終是最優的，而第一種方案則對於平方歐幾里得距離是最優的（至少有三本書的情況下）。

## 希區柯克問題

以下傳輸問題的表述由弗蘭克·勞倫·希區柯克提出：

假設有$m$個供應源$x_1, \ldots, x_m$為某一商品供貨，且每個供應源$x_i$處有$a(x_i)$個單位的供應量。同時有$n$個需求點$y_1, \ldots, y_n$需要該商品，每個需求點$y_j$處有$b(y_j)$個單位的需求。若$c(x_i,\ y_j)$表示從$x_i$到$y_j$的單位運輸成本，任務是找到一個流量分配方案，在滿足供應需求的同時最小化運輸成本。這一物流問題由德爾伯特·雷·富爾克森提出，並在他與小萊斯特·倫道夫·福特合著的《網路流》（Flows in Networks）（1962年）一書中得到了闡述。

佳林·庫普曼斯也為運輸經濟學與資源分配問題的表述作出了貢獻。

## 問題的抽象表述

### 蒙日和坎托羅維奇形式

由於黎曼幾何和測度論的發展，在現代或更加技術性的文獻中傳輸問題的表述有所不同。不過，上述礦山與工廠的簡單例子還是可以作為考慮抽象形式時一個有用的參照。此時我們可以考慮並非所有礦山和工廠都營業的情況，並允許一個礦山向多家工廠供貨，而一家工廠也可以從多個礦山接受礦石。

假設$X$和$Y$為兩個可分度量空間，使得$X$（或者$Y$ ) 上的每個機率測度都是拉東測度，並假設$c : X \times Y \to [0, \infty)$是一個博雷爾可測函數。給定$X$上的機率測度$\mu$和$Y$上的機率測度$\nu$，最優傳輸問題的蒙日形式是指尋找一個傳輸映射$T : X \to Y$，使得下式中的下確界成立：

$$\inf \left\{ \left. \int_X c(x, T(x)) \, \mathrm{d} \mu (x) \right| T_* (\mu) = \nu \right\},$$

其中$T_*(\mu)$由$T$推進$\mu$的向前推進算子。如果一個映射$T$達到這個下確界，則該映射被稱為「最優傳輸映射」。

最優傳輸問題的蒙日形式有其局限性，因為有時可能不存在滿足$T_*(\mu) = \nu$的映射$T$。例如，當$\mu$是狄拉克測度而$\nu$不是時，就會出現這種情況。

此時可以通過採用最優傳輸問題的坎托羅維奇形式來克服這一局限，即尋找$X \times Y$上的一個機率測度$\gamma$，使得下式中的下確界成立：

$$\inf \left\{ \left. \int_{X \times Y} c(x, y) \, \mathrm{d} \gamma (x, y) \right| \gamma \in \Gamma (\mu, \nu) \right\},$$

其中$\Gamma (\mu, \nu)$表示$X \times Y$上所有機率測度的集合併滿足邊緣分布$\mu$和$\nu$。可以證明，當成本函數$c$是下半連續並且$\Gamma(\mu, \nu)$是緊測度集合（拉東空間$X$和$Y$蘊含該條件）時，該問題總是存在最小值（另見沃瑟斯坦度量）。西古德·安格嫩特、史蒂文·哈克爾（Steven Haker）與艾倫·坦嫩鮑姆提出了蒙日–坎托羅維奇問題解的梯度下降表述。

### 對偶形式

坎托羅維奇問題的最小值等於

$$\sup \left( \int_X \varphi (x) \, \mathrm{d} \mu (x) + \int_Y \psi (y) \, \mathrm{d} \nu (y) \right),$$

其中上確界遍歷所有有界且連續、滿足

$$\varphi (x) + \psi (y) \leq c(x, y)$$

函數對。

### 經濟學解釋

將以上表述中的符號翻轉更利於從經濟學角度解釋這一問題。假設$x \in X$表示工人特徵的向量，$y \in Y$表示企業特徵的向量，$\Phi(x,y) =-c(x,y)$則表示工人$x$與企業$y$配對所創造的經濟產出。令$u(x) = -\varphi(x)$、$v(y) =-\psi(y)$，蒙日–坎托羅維奇問題可以重新表述為：

$$\sup \left\{ \int_{X\times Y}\Phi(x,y) d\gamma(x,y) ,\gamma \in \Gamma(\mu,\nu) \right\}$$

它的對偶形式為：

$$\inf \left\{ \int_X u(x) \,d\mu(x) +\int_Y v(y) \, d\nu (y) :u(x) +v(y) \geq \Phi(x,y) \right\}$$

其中下確界遍歷所有有界且連續的函數$u:X\to \mathbb{R}$和$v:Y\to \mathbb{R}$。如果對偶問題有解，則有：

$$v(y) =\sup_x \left\{ \Phi(x,y) - u(x)\right\}$$

可以將$u(x)$解釋為$x$類型工人的均衡工資，並將$v(y)$解釋為$y$類型企業的均衡利潤。

## 問題求解

### 一維連續情形

對於$1 \leq p < \infty$，假設$\mathcal{P}_p(\mathbb{R})$表示$\mathbb{R}$上所有$p$階矩有限的機率測度的集合。設$\mu, \nu \in \mathcal{P}_p(\mathbb{R})$且$c(x, y) = h(x-y)$，其中$h:\mathbb{R} \to [0,\infty)$是一個凸函數。

- 如果$\mu$沒有原子 (測度論)，即$\mu$的累積分布函數$F_\mu : \mathbb{R}\to[0,1]$是一個連續函數，則$F_{\nu}^{-1} \circ F_{\mu} : \mathbb{R} \to \mathbb{R}$是一個最優傳輸映射。如果$h$是嚴格凸的，則它是唯一的最優映射。
- 可以得到

$$\min_{\gamma \in \Gamma(\mu, \nu)} \int_{\mathbb{R}^2} c(x, y) \, \mathrm{d} \gamma (x, y) = \int_0^1 c \left( F_{\mu}^{-1} (s), F_{\nu}^{-1} (s) \right) \, \mathrm{d} s.$$

拉切夫（Rachev）與呂申多夫（Rüschendorf）於1998年給出了對此的證明。

### 離散情形與線性規劃

在邊緣分布$\mu$和$\nu$是離散的情形下，令$\mu_x$和$\nu_y$分別是分配給$x\in \mathbf{X}$和$y\in \mathbf{Y}$的機率質量，而$\gamma_{xy}$則是$xy$的分配機率。原始坎托羅維奇問題中的目標函數為

$$\sum_{x\in \mathbf{X},y\in \mathbf{Y}} \gamma_{xy}c_{xy}$$

並滿足約束條件

$$\sum_{y\in \mathbf{Y}}\gamma_{xy}=\mu_x,\forall x\in \mathbf{X}$$

$$\sum_{x\in \mathbf{X}} \gamma_{xy}=\nu_y,\forall y\in \mathbf{Y}.$$

為了將這一問題作為線性規劃問題處理，我們需要將矩陣$\gamma_{xy}$向量化，可以通過堆疊其行或列來完成，我們用$\operatorname{vec}$來表示這一操作。在行主序和列主序的情況下，上述約束條件可改寫為

$$\left( 1_{1\times |\mathbf{Y}| }\otimes I_{|\mathbf{X}|}\right) \operatorname{vec}(\gamma)=\mu$$

$$\left( I_{|\mathbf{Y}\|}\otimes 1_{1\times |\mathbf{X}|}\right) \operatorname{vec}(\gamma)=\nu$$

其中$\otimes$是克羅內克積，$1_{n\times m}$是一個大小為$n\times m$、所有元素為1的矩陣，而$I_{n}$是大小為$n$的單位矩陣。設$z=\operatorname{vec}(\gamma)$，該問題的線性規劃形式為

$$\begin{align} & \text{Minimize } &&\operatorname{vec}(c)^\top z \\[4pt] & \text{subject to:} && z \ge 0, \\[4pt] & && \begin{pmatrix} 1_{1\times |\mathbf{Y}|}\otimes I_{|\mathbf{X}|} \\ I_{|\mathbf{Y}|}\otimes 1_{1\times |\mathbf{X}|} \end{pmatrix} z=\binom{\mu}{\nu} \end{align}$$

這一問題可以很容易地通過大規模線性規劃求解器計算。

### 半離散情形

在半離散情況下，令$X=Y=\mathbb{R}^d$，且$\mu$是$\mathbb{R}^d$上的連續分布，$\nu =\sum_{j=1}^J \nu_j\delta_{y_i}$則是分配機率質量$\nu_j$到$y_j \in \mathbb{R}^d$的離散分布。此時，坎托羅維奇問題的原始形式為：

$$\inf \left\{ \int_X \sum_{j=1}^J c(x,y_j) \, d\gamma_j(x) ,\gamma \in \Gamma(\mu,\nu)\right\}$$

其中$\gamma \in \Gamma(\mu ,\nu)$滿足$\int_X d\gamma_j(x)=\nu_j$和$\sum_j d\gamma_j(x) =d\mu(x)$。

而其對偶形式則為

$$\sup \left\{ \int_X\varphi(x)d\mu (x)+\sum_{j=1}^J\psi_j\nu_j:\psi_j+\varphi (x)\leq c(x,y_j) \right\}$$

還可以寫為：

$$\sup_{\psi \in \mathbb{R}^J}\left\{ \int_X\inf_j\left\{ c(x,y_j)-\psi_j\right\} d\mu (x)+\sum_{j=1}^J\psi_j\nu_j\right\}$$

這是一個有限維凸優化問題，可以通過梯度下降等方法求解。

當$c(x,y)=|x-y|^2/2$時，可以證明分配給特定$j$的$x\in \mathbf{X}$集合是一個凸多面體，而得到的配置稱為冪圖。

### 二次正態情形

假設一個特殊情形$\mu =\mathcal{N}(0,\Sigma_X)$、$\nu =\mathcal{N}(0,\Sigma _{Y})$ 且$c(x,y) =|y-Ax|^2/2$，其中$A$是可逆矩陣。此時有

$$\varphi(x) =-x^\top \Sigma_X^{-1/2}\left( \Sigma_X^{1/2}A^\top \Sigma_Y A\Sigma_X^{1/2}\right) ^{1/2}\Sigma_{X}^{-1/2}x/2$$

$$\psi(y) =-y^\top A\Sigma_X^{1/2}\left( \Sigma_X^{1/2}A^\top \Sigma_Y A\Sigma_{X}^{1/2}\right)^{-1/2} \Sigma_X^{1/2}Ay/2$$

$$T(x) = (A^\top)^{-1}\Sigma_X^{-1/2} \left(\Sigma_X^{1/2}A^\top \Sigma_Y A\Sigma_X^{1/2} \right)^{1/2} \Sigma_X^{-1/2}x$$

加利雄（Galichon）於2016年證明了該情況下的解。

### 可分希爾伯特空間

令$X$是一個可分希爾伯特空間，定義$\mathcal{P}_p(X)$為$X$上所有$p$階矩有限的機率測度的集合，$\mathcal{P}_p^r(X)$則表示其中高斯正則的測度集合，即如果$g$是$X$上任何嚴格正的高斯測度且滿足$g(N) = 0$，則$\mu(N) = 0$也成立。

假設$\mu \in \mathcal{P}_p^r (X)$、$\nu \in \mathcal{P}_p(X)$，並且$c (x, y) = | x - y |^p/p$，其中$p\in(1,\infty), p^{-1} + q^{-1} = 1$。則坎托羅維奇問題存在一個唯一解$\kappa$，並且該解對應一個最優傳輸映射：即存在一個博雷爾映射$r\in L^p(X, \mu; X)$，使得

$$\kappa = (\mathrm{id}_X \times r)_{*} (\mu) \in \Gamma (\mu, \nu).$$

此外，如果$\nu$具有有界支撐，那麼對於$\mu$-幾乎所有的$x\in X$，存在局部利普希茨、 $c$-凹和最大坎托羅維奇勢$\varphi$，使得

$$r(x) = x - | \nabla \varphi (x) |^{q - 2} \, \nabla \varphi (x)$$

其中$\nabla \varphi$表示$\varphi$的加托導數。

## 熵正則化

考慮上述離散問題的一個變體：在原始問題的目標函數中添加一個熵正則化項

$$\begin{align} & \text{Minimize } \sum_{x\in \mathbf{X}, y\in \mathbf{Y}}\gamma_{xy}c_{xy}+\varepsilon \gamma_{xy} \ln \gamma_{xy} \\[4pt] & \text{subject to: } \\[4pt] & \gamma\ge0 \\[4pt] & \sum_{y\in \mathbf{Y}}\gamma_{xy} = \mu_x, \forall x\in \mathbf{X} \\[4pt] & \sum_{x\in \mathbf{X}}\gamma_{xy} = \nu_y, \forall y\in \mathbf{Y} \end{align}$$

相應的對偶問題為

$$\max_{\varphi ,\psi} \sum_{x\in \mathbf{X}} \varphi_x \mu_x + \sum_{y\in \mathbf{Y}} \psi_y v_y - \varepsilon \sum_{x\in \mathbf{X},y\in \mathbf{Y}} \exp \left( \frac{\varphi_x + \psi_y - c_{xy}}{\varepsilon }\right)$$

相較於不含正則化項的問題，原先對偶問題中的硬約束（$\varphi_x + \psi_y - c_{xy}\geq 0$）被替換為了軟約束，即懲罰項$\varepsilon \exp \left( (\varphi _x + \psi_y - c_{xy})/\varepsilon \right)$。對偶問題的最優條件可以表示為

5.1 $\mu_x = \sum_{y\in \mathbf{Y}} \exp \left( \frac{\varphi_x + \psi_y - c_{xy}}{\varepsilon} \right) ~\forall x\in \mathbf{X}$

5.2 $\nu_y = \sum_{x\in \mathbf{X}} \exp \left( \frac{\varphi_x + \psi_y - c_{xy}}{\varepsilon }\right) ~\forall y\in \mathbf{Y}$

令$A$為$|\mathbf{X}| \times |\mathbf{Y}|$的矩陣，其中元素$A_{xy}=\exp \left(-c_{xy} / \varepsilon \right)$。此時對偶問題的求解等價於尋找兩個對角正矩陣$D_1$和$D_2$，它們的大小分別為$|\mathbf{X}|$和$|\mathbf{Y}|$，使得$D_1AD_2 1_{|\mathbf{Y}|}=\mu$、$(D_1AD_2)^\top 1_{|\mathbf{X}|}=\nu$。矩陣$D_1$和$D_2$的存在性是辛克宏定理的推廣，可以使用辛克宏-諾普演算法進行求解。該演算法通過迭代求解5.1中的$\varphi_x$和5.2中的$\psi_y$實現。因此，辛克宏-諾普演算法相當於對偶正則問題的坐標下降法。

## 應用

蒙日-坎托羅維奇運輸問題已廣泛運用於許多領域，例如：

- 圖像配准與變形
- 反射面設計
- 從陰影照相與質子射線成像中提取資訊
- 地震層析成像和反射地震學
- 涉及總替代性質的經濟模型