最近很多读者和老顾联系，大家都对深度学习的几何观点有兴趣，询问最优传输理论的自学方法。恰巧老顾的合作者最近有一篇论文将三维的最优传输映射用于医学图像，获得了中国计算机图形学大会2018的最佳论文奖【1】。这篇文章的算法可以直接用于理解最优传输理论。

现代医学技术日益倚重医学图像，通过CT、MRI、PET等现代手段，人们可以获取各种医学图像，窥探人体内部器官，观察心脏脉动，药物扩散等等。通常，CT断层得到一张2维图像，多张CT图像合成一幅3维实体图像，再加上时间，我们可以得到4维图像。每个体素的灰度值代表了机体相应位置的密度。

Fig 1. 右侧膝盖和左侧胫骨的局部放大。

传统的光学放大镜可以放大或者缩小二维图片，对于3维实体或者更高维图像目前不存在相应的光学设备。我们可以采用数字放大镜，将感兴趣的空间区域放大或者缩小。如图1所示，我们用CT设备得到3维图像，左帧用体渲染（volume rendering）方法显示左右腿的膝盖和胫骨，右帧放大了左胫骨和右膝盖。图2显示的是同一种数据，只是改变了视角。图3 是另外一组脊椎和胸廓数据，我们将脊椎放大了8倍。这种体视放大镜允许医生从各个角度观察医学图像，也可以得到断层，详细检查，具有一定的实用价值。

Fig 2. 右侧膝盖和左侧胫骨的局部放大。

Fig 3. 脊椎体数据的放缩。

体视放大镜的几何原理如下：首先，3维体医学图像被放置在单位立方体内，每个体素（voxel）上有相应的灰度值或者伪彩色。用户选取单位立方体内的子区域，我们找到一种映射，将子区域膨胀，其他区域收缩，从而达到体视放大镜的效果。更进一步的描述如下：我们记为单位立方体，欧氏勒贝格测度（即体积元）为，即密度函数恒为1，记成。显然，单位立方体的总测度（体积）为1。用户随便指定另外一个概率测度，其密度函数为，总测度也是1，即

。

我们力图找到一个单位立方体到自身的光滑映射，将空间变形，有的地方膨胀，有的地方收缩，从而引发密度的变化。膨胀的地方得到稀释，密度减小；收缩的地方得到浓缩，密度增加。我们希望稀释浓缩之后，密度函数从变成了。

我们下面考虑这样的映射应该满足的偏微分方程。我们将定义域的参数选取为，值域的参数选取为，映射记为。映射将处的体元映射成处的体元。两个体元的比值，由映射雅克比矩阵（Jacobian Matrix）的行列式给出：

,

这里雅克比矩阵为：

。

我们的目的是将概率测度变成概率测度，这意味着：

,

代入以上各式，我们得到著名的雅克比方程（Jacobi Equation）

.

我们考察一个特例，令源概率测度和目标测度都是标准的勒贝格测度，即密度函数恒为1，那么雅克比方程的解代表保持体积不变的自映射。这种自映射依然有无穷多个，构成一个无穷维的李群，我们称之为单位立方体的特殊微分同胚群 （Special Diffeomorphism Group），记为。

一般而言，满足雅克比方程的解有无穷多个，我们从中择优录取，挑选标准是极小化所谓的传输代价（Transportation Cost），

.

我们在所有满足雅克比方程的映射中，寻找传输代价最小者，即所谓的最优传输映射（Optimal Mass Transportation Map）。

Brenier理论表明，如果源概率测度满足一些非常宽泛的条件，例如绝对连续，二阶矩有限，那么最优传输映射存在并且唯一，更为奇妙的是最优传输映射是某个凸函数的梯度映射，这个凸函数被称为是Brenier势能 。所谓函数是凸的，就是说函数的海森矩阵（Hessian matrix） 

.

是正定的（positive definite）。Brenier势能函数的梯度映射为，将点映成。那么，梯度映射应该满足雅克比方程，由此，我们得到著名的蒙日-安培方程（Monge-Ampere equation），

，

最优传输映射。

至此，体视放大镜的实现等价于求解蒙日-安培方程，从而得到最优传输映射。蒙日-安培方程本身是高度非线性方程，其求解本身需要复杂的计算过程。

根据Brenier极分解定理，所有满足雅克比方程的解都可以分解成一个最优传输映射和一个特殊微分同胚的复合，即，。

蒙日-安培方程最为直观的计算方法来自凸几何，并且所有的步骤都可以用经典的计算几何方法。

假设是欧氏空间的子集，我们说它是一个凸集 （convex set），如果对于其内部的任意两点，连接它们的线段也属于。给定欧氏空间的离散点集，，其凸包（convex hull）是包含点集的最小凸集，换句话说就是包含的所有凸集的交集就是凸包，记为。图4显示了三维空间中一些离散点的凸包。

Fig. 4 三维空间中的凸包。

凸包的计算方法非常成熟，最为普遍的是incremental 算法。比如，我们构造离散点集在三维空间中的凸包，如图4所示。首先，我们任选4个点，构造一个四面体。然后，我们每次增加一个点，逐渐更新凸包。更新方法如下：我们从新加的点处观察当前的凸包，去除凸包上所有能够看得到的面（visible facet），留下所有无法看到的面（invisible facet）。将凸包轮廓线上的每条边和新增加点相连，如此得到更新的凸包。如此重复，直至穷尽所有离散点。

给定维欧氏空间中的一张超平面（hyper-plane），

，

超平面的对偶点为

。

超平面将维欧氏空间分成两个半空间（half space），上半空间和下半空间。上半空间定义为高于平面之上的点集：

，

下半空间类似定义。给定欧氏空间中的一种超平面，这族超平面的上包络（upper envelope）定义成所有超平面对应的上半空间的交集，

 。

超平面族的对偶点集为

,

那么上包络和凸包相互勒让德对偶 (legendre dual)。

所谓的勒让德对偶定义如下：给定函数，其对偶定义为

.

如果我们将上包络视为函数的图（graph），则凸包是的图。我们可以通过计算凸包来计算上包络。

上包络向的投影得到的Power Diagram，凸包向的投影得到的 Weighted Delaunay三角剖分。Power Diagram 和 Weighted Delaunay triangulation彼此相差一个庞加莱对偶。Power Diagram的每个n维胞腔对应着Weighted Delaunay的一个顶点，如果Diagram中的两个n维胞腔共同分享一个（n-1）维的胞腔，则它们在Weighted Delaunay中对偶的两个顶点彼此相连。这意味着Weighted Delaunay中的每条边对偶着Power Diagram中的一个（n-1）维胞腔。

Fig. 5 上包络、凸包、Power Digram、Weighted Delaunay 三角剖分。

图5显示了这4个概念，和它们之间的对偶关系。从计算角度而言，只要能够计算凸包，就可以计算这4个概念。

我们下面考虑最优传输的特殊情形，半离散最优传输（Semi Discrete Optimal Transportation）。给定欧氏空间中的凸区域，源概率密度函数恒等于1，目标概率测度为离散测度，，为狄拉克函数，即质量集中在离散点上。根据Brenier定理，存在一个凸函数（Brenier势能），，其梯度映射就是最优传输映射，

奇妙的是，Brenier定理对于离散测度依然成立，这里Brenier势能函数为分片线性凸函数，其梯度为离散点集，由此我们有离散Brenier势能为

,

即Brenier势能函数的图（graph）等价于平面族的上包络。上包络向的投影得到Power Diagram，支撑平面的投影为胞腔。Brenier势能函数的梯度映射将每个和的交集映射到，。和交集的体积为。我们通过调节每个超平面的截距，使得每个胞腔的体积等于相应像点的目标测度，，那么Brenier势能函数的梯度映射给出了最优传输映射。

Fig.6 亚历山大定理。

由此，问题归结为如何调节截距，使得Power Diagram的胞腔体积等于目标测度。根据Brenier定理，Brenier势能函数存在，并且彼此相差一个常数，最优传输映射唯一。这一结论等价于经典的Alexandrov定理：如图6所示，给定凸多面体每个面的法向量和每个面的投影面积，则凸多面体存在，并且彼此相差一个铅直平移。

我们在【2】中给出了一个基于变分的算法来构造上包络，通过优化下列凸能量，我们能够得到截距，从而得到最优传输映射，

.

这个能量的梯度为。海森矩阵也可以直接构造，Power Diagram中的每个（n-1）维胞腔的体积除以对应的Weighted Delaunay三角剖分中对应的边长，给出了海森矩阵的元素。由此，我们可以用牛顿法直接优化这个能量，得到最优传输映射。

这里，我们用离散概率测度来逼近目标概率测度，从而用离散的上包络来逼近连续Briener势能函数。如果目标概率测度的采样距离为，汪徐家教授的理论工作【3】表明离散和连续的Briener势能函数之差的模为，最优传输映射的误差为。

这里的计算框架对于任意维都成立，理论上可以用于深度学习的生成模型，实践中需要克服维数诅咒（dimension curse）。最优传输理论本身博大精深，我们这里只是从几何角度讲述最为直观的内容，侧重于计算方法。Kontarovich理论更加广泛，目前在深度学习中普遍使用，但是其结构不如Brenier理论丰富，目前收敛阶的理论分析匮乏。

这篇体视放大镜是几位教授在石溪访问期间共同完成。当时大家经常在一起切磋探讨，几乎天天在推导数学、编写程序。在落叶纷飞的枫林中讨论微分拓扑中的叶状结构（foliation），在惊涛拍岸的海边辩论最优传输。夏季的长岛成为旅游胜地，人们在海边游泳嬉闹，纵情欢笑。但是来自武汉大学计算机系的苏科华教授却躺在树荫下独自苦读微分几何，令人深为感动。从这些学者身上，老顾看到了中国年轻一代孜孜以求、奋发图强的精神。相信最优传输理论会在海内外更加广泛地传播，在工程医疗方面得到更多应用。

References

1. Kehua Su， Na Lei， Ming Ma, Li Cui, Xianfeng Gu.  "Focus + Context Visualization Based on Optimal Mass Transportation".  Best Paper Award，Chinagraph 2018.

2. Xianfeng Gu, Feng Luo, Jian Sun, and Shing-Tung Yau. "Variational principles for minkowski type problems, discrete optimal transport", and discrete monge-ampere equations. Asian Journal of Mathematics (AJM), 20(2):383-398, 2016.

3. Haodi Chen， Genggeng Huang and Xu-Jia Wang， “Convergence rate estimates for Aleksandrov's solution to the Monge-Ampere Equation", Accepted by SIAM J. Numerical Analysis.

4. Na Lei,Kehua Su,Li Cui,Shing-Tung Yau,David Xianfeng Gu, "A Geometric View of Optimal Transportation and Generative Model", arXiv:1710.05488. https://arxiv.org/abs/1710.05488