Notions of optimal transport theory and how to implement them on a computer
最优输送理论的概念及如何在计算机上实施它们

Monge’s initial motivation to study this problem was very practical: supposing you have an army of workmen, how can you transform a terrain with an initial landscape into a given desired target landscape, while minimizing the total amount of work ?
蒙日最初研究这个问题的动机非常实际：假设你有一支工人队，如何将初始地形转变为目标地形，同时最大限度地减少总工作量？

Monge’s initial problem statement was as follows:
蒙日最初的问题陈述如下：

$$\begin{array}[]{l}\inf\limits_{T:X\rightarrow X}\int\limits_{X}c(x,T(x))u(x)dx\\[17.07164pt] \mbox{\ \ \ \ subject to: }\\[17.07164pt] \forall B\subset X,\int\limits_{T^{-1}(B)}u(x)dx=\int\limits_{B}v(x)dx\end{array}$$

where $X$ is a subset of ${\mathbb{R}}^{2}$, $u$ and $v$ are two positive functions defined on $X$ and such that $\int_{X}u(x)dx$ = $\int_{Y}v(x)dx$, and $c(\cdot,\cdot)$ is a convex distance (the Euclidean distance in Monge’s initial problem statement).
其中 $X$ 是 ${\mathbb{R}}^{2}$ 的子集， $u$ 和 $v$ 是在 $X$ 上定义的两个正函数，使得 $\int_{X}u(x)dx$ = $\int_{Y}v(x)dx$ ， $c(\cdot,\cdot)$ 是一个凸距离（蒙日初始问题陈述中的欧几里得距离）。

The functions $u$ and $v$ represent the height of the current landscape and the height of the target landscape respectively (symbolized as gray levels in Figure 3). The problem consists in finding (if it exists) a function $T$ from $X$ to $X$ that transforms the current landscape $u$ into the desired one $v$, while minimizing the product of the amount of transported earth $u(x)$ with the distance $c(x,T(x))$ to which it was transported. Clearly, the amount of earth is conserved during transport, thus the total quantity of earth should be the same in the source and target landscapes (the integrals of $u$ and $v$ over $X$ should coincide). This global matter conservation constraint needs to be completed with a local one. The local matter conservation constraint enforces that in the target landscape, the quantity of earth received in any subset $B$ of $X$ corresponds to what was transported here, that is the quantity of earth initially present in the pre-image $T^{-1}(B)$ of $B$ under $T$. Without this constraint, one could locally create matter in some places and annihilate matter in other places in a counterbalancing way. A map $T$ that satisfies the local mass conservation constraint is called a transport map.
函数 $u$ 和 $v$ 分别表示当前地貌的高度和目标地貌的高度(在图 3 中以灰度级表示)。问题在于找到（如果存在）一个从 $X$ 到 $X$ 的函数 $T$ ，将当前地貌 $u$ 转变为期望的地貌 $v$ ，同时最小化通过运输的土壤数量 $u(x)$ 与其运输距离 $c(x,T(x))$ 的乘积。显然，土壤的数量在运输过程中是不变的，因此源地貌和目标地貌中土壤的总量应当相同（在 $X$ 上 $u$ 和 $v$ 的积分应当一致）。这一全局物质守恒约束需要进一步补充一条局部物质守恒约束。局部物质守恒约束要求在目标地貌中，任意子集 $B$ 所接收的土壤数量应与运输到这里的土壤数量对应，也就是最初存在于 $T$ 下的 $B$ 的前像中土壤数量应当一致。如果没有这一约束，人们可以在一些地方本地创造物质，并在其他地方以相等的方式消除物质。 符合局部质量守恒约束的地图称为输运地图。

We now suppose that instead of a “target landscape”, we wish to transport earth (or a resource) towards a set of points (that will be denoted by $Y$ for now on), that represent for instance a set of factories that exploit a resource, see Figure 4. Each factory wishes to receive a certain quantity of resource (depending for instance of the number of potential customers around the factory). Thus, the function $v$ that represents the “target landscape” is replaced with a function on a finite set of points. However, if a function $v$ is zero everywhere except on a finite set of points, then its integral over $X$ is also zero. This is a problem, because for instance one cannot properly express the mass conservation constraint. For this reason, the notion of function is not rich enough for representing this configuration. One can use instead measures (more on this below), and associate with each factory a Dirac mass weighted by the quantity of resource to be transported to the factory.
现在我们假设，我们不是想要将土地（或资源）运输到“目标景观”，而是希望将土地（或资源）运输到一组点（暂时表示为 $Y$ ）。这些点例如代表一组开采资源的工厂，参见图 4。每个工厂都希望获得一定数量的资源（取决于工厂周围潜在客户的数量）。因此，表示“目标景观”的函数 $v$ 被替换为一组有限点上的函数。然而，如果函数 $v$ 在除一组有限点外的所有地方均为零，则其在 $X$ 上的积分也为零。这是一个问题，因为例如不能正确表达质量守恒约束。因此，函数的概念不足以表示这种配置。可以改用度量（下文详述），并将每个工厂与要运往工厂的资源数量加权的 Dirac 质量相关联。

From now on, we will use measures $\mu$ and $\nu$ to represent the “current landscape” and the “target landscape”. These measures are supported by sets $X$ and $Y$, that may be different sets (in the present example, $X$ is a subset of ${\mathbb{R}}^{2}$ and $Y$ is a discrete set of points). Using measures instead of function not only makes it possible to study our “transport to discrete set of factories” problem, but also it can be used to formalize computer objects (meshes) and directly leads to a computational algorithm. This algorithm is very elegant because it is a verbatim computer translation of the mathematical theory (see §7.6). In this particular setting, translating from the mathematical language to the algorithmic setting does not require to make any approximation. This is made possible by the generality of the notion of measure.
从现在开始，我们将使用测度 $\mu$ 和 $\nu$ 来表示“当前景观”和“目标景观”。这些测度由集合 $X$ 和 $Y$ 支持，这些集合可能是不同的（在当前示例中， $X$ 是 ${\mathbb{R}}^{2}$ 的子集， $Y$ 是一组离散的点）。使用测度而不是函数不仅可以研究我们的“运输到离散工厂集合”的问题，还可以用来形式化计算机对象（网格），并直接导致一个计算算法。这个算法非常优雅，因为它是数学理论的逐字计算机转换（见第 7.6 节）。在这种特定情境中，将数学语言转化为算法设置不需要进行任何近似。这是由测度概念的一般性所实现的。

The reader who wishes to learn more on measure theory may refer to the textbook [Tao11]. To keep the length of this article reasonable, we will not give here the formal definition of a measure. In our context, one can think of a measure as a “function” that can be only queried using integrals and that can be “concentrated” on very small sets (points). The following table can be used to intuitively translate from the “language of functions” to the “language of measures” :
希望学习测度论的读者可以参考教材[Tao11]。为了保持本文的长度合理，我们不会在此处给出测度的正式定义。在我们的语境中，可以将测度看作是一个只能通过积分查询并能“聚集”在非常小的集合（点）上的“函数”。以下表格可用于直观地从“函数的语言”翻译成“测度的语言”：

$$\begin{array}[]{|c|c|}\hline\cr\mbox{Function $u$}&\mbox{Measure $\mu$ }\\[5.69054pt] \hline\cr\int_{B}u(x)dx&\mu(B)\mbox{ or }\int_{B}d\mu\\[5.69054pt] \int_{B}f(x)u(x)dx&\int_{B}f(x)d\mu\\[5.69054pt] u(x)&N/A\\ \hline\cr\end{array}$$

(Note: in contrast with functions, measures cannot be evaluated at a point, they can be only integrated over domains).
（注意：与函数相对应，测度不能在一个点被评估，它们只能在域上被积分）。

In its version with measures, Monge’s problem can be stated as follows:
在具有测度的版本中，蒙日问题可以陈述如下：

$$\begin{array}[]{l}\inf\limits_{T:X\rightarrow Y}\int\limits_{X}c(x,T(x))d\mu\ \mbox{ subject to }\ \nu=T\sharp\mu\end{array}$$

(M)

where $X$ and $Y$ are Borel sets (that is, sets that can be measured), $\mu$ and $\nu$ are two measures on $X$ and $Y$ respectively such that $\mu(X)=\nu(Y)$ and $c(\cdot,\cdot)$ is a convex distance. The constraint $\nu=T\sharp\mu$, that reads “$T$ pushes $\mu$ onto $\nu$” corresponds to the local mass conservation constraint. Given a measure $\mu$ on $X$ and a map $T$ from $X$ to $Y$, the measure $T\sharp\mu$ on $Y$, called “the pushforward of $\mu$ by $T$”, is such that $T\sharp\mu(B)=\mu(T^{-1}(B))$ for all Borel set $B\subset Y$. Thus, the local mass conservation constraint means that $\mu(T^{-1}(B))=\nu(B)$ for all Borel set $B$ $\subset$ $Y$.
其中 $X$ 和 $Y$ 是 Borel 集合（即可以被测量的集合）， $\mu$ 和 $\nu$ 分别是 $X$ 和 $Y$ 上的两个测度，使得 $\mu(X)=\nu(Y)$ 和 $c(\cdot,\cdot)$ 是一个凸距离。约束条件 $\nu=T\sharp\mu$ ，“ $T$ 将 $\mu$ 推到 $\nu$ 上”对应于局部质量守恒约束。给定 $X$ 上的一个测度 $\mu$ 和从 $X$ 到 $Y$ 的映射 $T$ ， $Y$ 上的测度 $T\sharp\mu$ ，称为“ $\mu$ 被 $T$ 推送的前向”，满足所有 Borel 集合 $B\subset Y$ 。因此，局部质量守恒约束意味着所有 Borel 集合 $B$ 上 $\subset$ $Y$ 。

The local mass conservation constraint makes the problem very difficult: imagine now that you want to implement a computer program that enforces it: the constraint concerns all the subsets $B$ of $Y$. Could you imagine an algorithm that just tests whether a given map satisfies it ? What about enforcing it ? We will see below a series of transformations of the initial problem into equivalent problems, where the constraint becomes linear. We will finally end up with a simple convex optimization problem, that can be solved numerically using classical methods.
本地质量守恒约束使问题变得非常困难： 现在想象一下，您想要实现一个强制执行它的计算机程序： 约束涉及所有子集 $B$ $Y$ 。 您能想象一个仅测试给定映射是否满足它的算法吗？关于强制执行它呢？我们将在下面看到一系列初始问题转化为等价问题的转换，约束变为线性。 最后，我们将得到一个简单的凸优化问题，可以使用经典方法进行数值解决。

Before then, let us get back to examine the original problem. The local mass conservation constraint is not the only difficulty: the functional optimized by Monge’s problem is non-symmetric, and this causes additional difficulties when studying the existence of solutions for problem (M). The problem is not symmetric because $T$ needs to be a map, therefore the source and target landscape do not play the same role. Thus, it is possible to merge earth (if $T(x_{1})=T(x_{2})$ for two different points $x_{1}$ and $x_{2}$), but it is not possible to split earth (for that, we would need a “map” $T$ that could send the same point $x$ to two different points $y_{1}$ and $y_{2}$). The problem is illustrated in Figure 5: suppose that you want to compute the optimal transport between a segment $L_{1}$ (that symbolizes a “wall of earth”) and two parallel segments $L_{2}$ and $L_{3}$ (that symbolize two “trenches” with a depth that correspond to half the height of the wall of earth). Now we want to transport the wall of earth to the trenches, to make the landscape flat. To do so, it is possible to decompose $L_{1}$ into segments of length $h$, sent alternatively towards $L_{2}$ and $L_{3}$ (Figure 5 on the left). For any length $h$, it is always possible to find a better map $T$, that is a lower value of the functional in (M), by subdividing $L_{1}$ into smaller segments (Figure 5 on the right). The best way to proceed consists in sending from each point of $L_{1}$ half the earth to $L_{2}$ and half the earth to $L_{3}$, which cannot be represented by a map. Thus, the best solution of problem (M) is not a map. In a more general setting, this problem appears each time the source measure $\mu$ has mass concentrated on a manifold of dimension $d-1$ [McC95] (like the segment $L_{1}$ in the present example).
在此之前，让我们回到检查原始问题。当地质量守恒约束并不是唯一的困难时：蒙热问题的优化功能是非对称的，这在研究问题（M）的解存在时导致额外困难。该问题不对称是因为 $T$ 需要是一张地图，因此源地形和目标地形扮演的角色并不相同。因此，可以合并大地（如果 $T(x_{1})=T(x_{2})$ 两个不同点 $x_{1}$ 和 $x_{2}$ ），但不可能分开地球（为此，我们需要一张“地图” $T$ 能将同一点 $x$ 发送到两个不同点 $y_{1}$ 和 $y_{2}$ ）。该问题在图 5 中得到了说明：假设你想要计算一个代表“一堵土墙”的线段 $L_{1}$ 和两个平行线段 $L_{2}$ 和 $L_{3}$ 之间的最佳传输（代表两个深度相当于土墙高度一半的“沟渠”）。现在我们想把土墙运输到沟渠，使地形平坦。为了做到这一点，可以将 $L_{1}$ 分解为长度为 $h$ 的线段，交替地发送到 $L_{2}$ 和 $L_{3}$ （左边的图 5）。 对于任意长度 $h$ ，总是可以通过将 $L_{1}$ 细分为更小的段（右侧第 5 图）而找到一个更好的地图 $T$ ，即(M)中的功能值更低的地图。最佳操作方法是从 $L_{1}$ 的每一点向 $L_{2}$ 发送地球的一半，向 $L_{3}$ 发送地球的一半，而这不能用地图表示。因此，问题(M)的最佳解决方案不是地图。在一个更一般的设定中，每当源测度 $\mu$ 在一个维度 $d-1$ 的流形上有集中质量时[McC95]（例如当前示例中的段 $L_{1}$ ）时，该问题会出现。

In Figure 7, we show a discrete version of the 1D transport between two segments of Figure 6. The measures $\mu$ and $\nu$ are replaced with vectors $U=(u_{i})_{i=1\ldots m}$ and $V=(v_{j})_{i=1\ldots n}$. The transport plan $\gamma$ becomes a set of coefficients $\gamma_{ij}$. Each coefficient $\gamma_{ij}$ indicates the quantity of matter that will be transported from $u_{i}$ to $v_{j}$. The discrete Kantorovich problem can be written as follows:
在图 7 中，我们展示了图 6 中两个区段之间的一维传输的离散版本。测度 $\mu$ 和 $\nu$ 被向量 $U=(u_{i})_{i=1\ldots m}$ 和 $V=(v_{j})_{i=1\ldots n}$ 替换。传输计划 $\gamma$ 成为一组系数 $\gamma_{ij}$ 。每个系数 $\gamma_{ij}$ 表示将从 $u_{i}$ 输送到 $v_{j}$ 的物质数量。离散 Kantorovich 问题可以表示如下：

$$\quad\quad\min\limits_{\gamma}<C,\gamma>\mbox{ subject to }\left\{\begin{array}[]{l}{\bf P_{x}}\gamma=U\\ {\bf P_{y}}\gamma=V\\ \gamma_{i,j}\geq 0\quad\forall i,j\end{array}\right.$$

(2)

where $\gamma$ is the vector of $\mathbb{R}^{m\times n}$ with all coefficients $\gamma_{ij}$ (that is, the matrix $\gamma_{ij}$ “unrolled” into a vector), and $C$ the vector of $\mathbb{R}^{m\times n}$ with the coefficients $c_{ij}$ indicating the transport cost between point $i$ and point $j$ (for instance, the Euclidean cost). The objective function is simply the dot product, denoted by $<C,\gamma>$, of the cost vector $C$ and the vector $\gamma$. The objective function is linear in $\gamma$. The constraints on the marginals (1) impose in this discrete version that the sums of the $\gamma_{ij}$ coefficients over the columns correspond to the $u_{i}$ coefficients (Figure 7-B) and the sums over the rows correspond to the $v_{j}$ coefficients (Figure 7-C). Intuitively, everything that leaves point $i$ should correspond to $u_{i}$, that is the quantity of matter initially present in $i$ in the source landscape, and everything that arrives at a point $j$ should correspond to $v_{j}$, that is the quantity of matter desired at $j$ in the target landscape. As one can easily notice, in this form, both constraints are linear in $\gamma$. They can be written with two matrices ${\bf P_{x}}$ and ${\bf P_{y}}$, of dimensions $m\times mn$ and $n\times mn$ respectively.
其中 $\gamma$ 是具有所有系数 $\gamma_{ij}$ （即矩阵 $\gamma_{ij}$ “展开”为向量）的 $\mathbb{R}^{m\times n}$ 的向量， $C$ 是具有指示点 $i$ 和点 $j$ 之间运输成本 $c_{ij}$ 的 $\mathbb{R}^{m\times n}$ 的向量（例如，欧几里得成本）。目标函数简单地是成本向量 $C$ 和向量 $\gamma$ 的点积，用 $<C,\gamma>$ 表示。目标函数在 $\gamma$ 中是线性的。边际限制（1）在这个离散版本中要求列上系数 $\gamma_{ij}$ 的总和对应于 $u_{i}$ 系数（图 7-B），行上的总和对应于 $v_{j}$ 系数（图 7-C）。直观地看，离开点 $i$ 的一切都应与 $u_{i}$ 相对应，即最初在源景观中存在的物质量，而到达某一点 $j$ 的一切都应与 $v_{j}$ 相对应，即目标景观中 $j$ 处所需的物质量。正如大家可以轻易发现的那样，在这种形式下，这两个约束在 $\gamma$ 中都是线性的。它们可以分别用两个维数分别为 $m\times mn$ 和 $n\times mn$ 的矩阵 ${\bf P_{x}}$ 和 ${\bf P_{y}}$ 来表示。

We introduce, arbitrarily for now, the following function $L$ defined by:
我们介绍，目前随意地，以下由 $L$ 定义的函数。

$$L(\varphi,\psi)=<C,\gamma>-<\varphi,{\bf P_{x}}\gamma-U>-<\psi,{\bf P_{y}}\gamma-V>$$

that takes as arguments two vectors, $\varphi$ in ${\mathbb{R}}^{m}$ and $\psi$ in ${\mathbb{R}}^{n}$. The function $L$ is constructed from the objective function $<C,\gamma>$ from which we subtracted the dot products of $\varphi$ and $\psi$ with the vectors that correspond to the degree of violation of the constraints. One can observe that:
该函数接受两个向量作为参数， $\varphi$ 在 ${\mathbb{R}}^{m}$ 中， $\psi$ 在 ${\mathbb{R}}^{n}$ 中。函数 $L$ 是从目标函数 $<C,\gamma>$ 构建的，我们从中减去了与对应于约束违反程度的向量的点积。可以观察到：

$$\begin{array}[]{lcl}\sup\limits_{\varphi,\psi}[L(\varphi,\psi)]&=&<C,\gamma>\mbox{ if }{\bf P_{x}}\gamma=U\mbox{ and }{\bf P_{y}}\gamma=V\\ &=&+\infty\mbox{ otherwise.}\end{array}$$

Indeed, if for instance a component $i$ of ${\bf P_{x}}\gamma$ is non-zero, one can make $L$ arbitrarily large by suitably choosing the associated coefficient $\varphi_{i}$.
实际上，例如，如果 ${\bf P_{x}}\gamma$ 的一个分量 $i$ 是非零的，则可以通过适当选择相关系数 $\varphi_{i}$ 使 $L$ 随意增大。

Now we consider: 现在我们考虑:

$$\inf\limits_{\gamma\geq 0}\left[\sup\limits_{\varphi,\psi}[L(\varphi,\psi)]\right]=\inf\limits_{\tiny\begin{array}[]{l}\gamma\geq 0\\ {\bf P_{x}}\gamma=U\\ {\bf P_{y}}\gamma=V\end{array}}\left[<C,\gamma>\right].$$

There is equality, because to minimize $\sup[L(\varphi,\psi)]$, $\gamma$ has no other choice than satisfying the constraints (see the previous observation). Thus, we obtain a new expression (left-hand side) of the discrete Kantorovich problem (right-hand side). We now further examine it, and replace $L$ by its expression:
有相等性，因为要使 $\sup[L(\varphi,\psi)]$ 最小化， $\gamma$ 除了满足限制条件之外别无选择（请参见前面的观察）。 因此，我们获得了离散 Kantorovich 问题的新表达式（左侧）， 第 2 项用其表达式替换。

			$\displaystyle\inf\limits_{\gamma\geq 0}\left[\sup\limits_{\varphi,\psi}\left(\begin{array}[]{ll}<C,\gamma>&-<\varphi,{\bf P_{x}}\gamma-U>\\ &-<\psi,{\bf P_{y}}\gamma-V>\end{array}\right)\right]$		(5)
		$\displaystyle=$	$\displaystyle\sup\limits_{\varphi,\psi}\left[\inf\limits_{\gamma\geq 0}\left(\begin{array}[]{ll}<C,\gamma>&-<\varphi,{\bf P_{x}}\gamma-U>\\ &-<\psi,{\bf P_{y}}\gamma-V>\end{array}\right)\right]$		(8)
		$\displaystyle=$	$\displaystyle\sup\limits_{\varphi,\psi}\left[\inf\limits_{\gamma\geq 0}\left(\begin{array}[]{l}<\gamma,C-{\bf P_{x}}^{t}\varphi-{\bf P_{y}}^{t}\psi>+\\[2.84526pt] <\varphi,U>+<\psi,V>\end{array}\right)\right]$		(11)
		$\displaystyle=$	$\displaystyle\sup\limits_{\tiny\begin{array}[]{c}\varphi,\psi\\ {\bf P_{x}}^{t}\varphi+{\bf P_{y}}^{t}\psi\leq C\end{array}}\left[<\varphi,U>+<\psi,V>\right].\quad\quad$		(14)