依随深度学习方法的深入发展，人们逐渐意识到最优传输理论所起到的奠基作用。深度学习可以抽象概括为学习流形上的概率分布，因此概率分布之间的变换和衡量概率测度之间的距离是深度学习最为核心的主题。而最优传输理论为解决这些问题提供了理论基础和计算工具。更进一步，最优传输理论为所有概率测度构成的空间定义了黎曼度量和协变微分，从而使得在所有测度构成的空间中进行变分优化成为可能。因此，我们可以预见在深度学习热潮渐落之后，最优传输理论经过大浪淘沙，会凝练成思想精华，全面渗透入工程领域。
但是，描述最优传输映射的偏微分方程具有强烈的非线性，因此精确计算最优传输一直具有挑战性。
目前在工程领域，特别是在深度学习领域，人们都是将概率测度离散化，将Kantorovich能量加上熵正则项，得到光滑的近似，逼近误差用一个参数进行控制，例如常见的Sinkhorn算法，Nesterov近似等等。
我们经过大量实验，发现这些算法有很多内在缺陷。
如果我们希望精确计算连续分布之间的最优传输映射，这些算法需要对源区域和目标区域稠密采样，如此计算复杂度过高，无法用普通硬件实现。
同时这些算法给出了近似解，如果我们希望通过减小控制参数来提高逼近精度，因为算法依赖超越运算（指数对数运算），算法稳定性迅速变差，无法收敛。
如果我们细致地调节参数以保证算法的稳定性，这时得到的近似解与真解相距甚远。
我们一直希望能够对最优传输映射进行深入细致地理论研究，因此要求能够得到高精度的数值解，同时可以控制误差范围。
Sinkhorn算法虽然火热，但是理论上无法真正满足需求。
我们发现，目前几乎不存在高精度最优传输算法的软件工具。
在基础和应用数学领域，虽然大家数十年如一日地研究蒙日-安培方程，但是并没有比较普及的计算工具；
在计算机科学领域，绝大多数的算法都是基于Kantorovich的线性规划及其光滑近似（Sinkhorn）。
这些方法空间复杂度很高，逼近误差较大，更为严重的是，这些方法破坏了最优传输问题本身的理论结构，摒弃了内在的几何特性。
目前工程领域对于这些算法的热衷，很大程度上是用方便来取代有效，并且将最优传输视为一个成型的工具，而没有深入理解其内在的理论结构。
因此，在我们的
课程上，我们着重讲解了最优传输问题的几何求解方法，这一方法与蒙日-安培方程的弱解理论相一致，与Minkowski和Alexandrov的凸几何理论相吻合，从几何角度来理解诠释最优传输理论，给人以强烈的几何直觉和深刻细致的洞察。
关于最优传输映射，我们依然有太多的未知，我们希望几何方法提供的计算工具能够帮助大家更加深入地研究这一理论，对于蒙日-安
培方程具有更
加深刻的理解，从而进一步推动工程应用的发展。
给定欧氏空间的概率测度及其支集
和
，满足总测度相等
，密度函数
和
，
。一个映射
是保持测度的，如果对一切
，都有
记为
。给定传输代价函数
，蒙日问题是求解所有保测度映射中传输代价最小者，
如果传输代价是欧氏距离的平方
，
和
紧致，
为凸区域，那么Brenier定理断言存在一个凸函数
，即所谓的Brenier势能函数，其梯度映射给出了最优传输映射，
，并且这种最优传输映射是惟一的。Brenier势能函数满足蒙日-安培方程（Monge-Ampere）
并且满足边界条件
。由此，求解最优传输映射等价于求解蒙日-安培方程。
我们在
中离散采样，得到采样点集
，用Dirac测度之和来逼近目标测度
，由此来考虑半连续的最优传输问题。由Brenier定理，最优传输映射等于Brenier势能函数的梯度映射，这时Brenier势能函数为分片线性函数，
这里，每个采样点
对应一个支撑平面
，高度
为未知量，
的图是这些支撑平面的上包络(upper envelope)，其勒让德对偶
的图是点
的凸包(convex hull)。
图1. Brenier势能函数与勒让德对偶。
如（1）图所示，Brenier势能函数的投影得到
上的一个Power图（power diagram），
Brenier势能函数的Legendre对偶
的投影，得到
的Weighted Delaunay三角剖分
。通过调节每个
对应支撑平面的高度
，我们可以得到这样一个Power图，每个Power胞腔
的体积等于
，
丘成桐先生，笔者和罗锋【6】曾经证明，所求高度向量
是下面能量的最大化子，
在可容许高度空间
这一能量是严格凹的，最大化子在可容许空间的内部。由定义，我们有能量的梯度分量等于
，能量Hessian矩阵的分量等于Power Voronoi边长与相应的Weighted Delaunay边长之比。由此，我们可以用牛顿法来优化凸能量从而得到唯一解。
图2. 凸包算法的结果。
我们看到，计算所需要的主要概念和算法都来自经典的计算几何。最为主要的是Legendre对偶
，这需要计算
的凸包。我们可以采用Edelsbrunner的经典增量凸包算法。这也等价于计算Weighted Delaunay三角剖分，我们可以用Lawson的边翻转（edge flip）算法来实现。如果图（2）所示。
图3. 上包络算法的结果。
得到凸包
之后，我们计算其相应的上包络，即计算
的勒让德对偶，凸包上每一个顶点
对偶于一个平面
，每条边对偶于两个顶点对偶平面的交线，每个面对偶于三个顶点对偶平面的交点。同时我们添加一个无穷远平面
，如此得到平面族的上包络，如图（3）所示。这里，边界支撑平面被无穷远平面相截。
图4. Sutherland-Hodgman算法。
上包络在
上的投影得到Power图，我们需要计算每个Power胞腔
与
的交集，即用多边形
来剪裁
。平面多边形裁剪有多种算法，比较高效的是Bentley-Ottman的扫描线算法（sweepline），比较简单的是Sutherland-Hodgman的凸多边形剪裁算法。如图（4）所示，我们用凸的红色多边形来剪裁蓝色多边形，每一步我们用红色多边形的一条边来将平面一分为二，保留内侧去除外侧。如此重复，每次切掉一角，直至穷尽红色多边形的所有边。
图5. Brenier势能函数，边界裁剪之后的结果。
图（5）显示了裁剪之后的Brenier势能函数，和相应的Power图。然后，我们计算每个Power胞腔的面积，每条Power Diagram边的长度，构造能量的梯度和Hessian矩阵，计算更新的高度向量
，
然后我们进行试探性的更新：选定一个步长
，
。我们需要判断新的高度向量是否落在可容许空间内部。我们先用新的高度向量计算凸包，如果存在某个
不落在凸包表面，则意味着新的高度向量非法。或者，所有的顶点都落在凸包表面上，但是某个Power胞腔落在
之外，则新的高度向量非法。若新的高度向量非可容许，则我们将步长减半
，重新试探。这种方法我们称之为阻尼算法。几年前，笔者与丘成桐先生和汪徐家教授分别探讨过求解蒙日-安培方程的计算问题，他们都一针见血地指出：求解过程的关键是保证中间解的严格凸性，优化过程一直在可容许解空间
之中。所以，这一步是整个算法的关键。
图6. 算法的收敛过程。
图（6）显示了一个计算实例，由于我们采用了牛顿法来优化凸能量，算法收敛很快。一般情形，如人脸曲面，我们用黎曼映照共形映射到平面圆盘，将共形因子作为概率密度函数
，
为均匀分布。如果在人脸上均匀采样
个点，步长调整为
，则
步迭代就可以达到相对误差（每个Power胞腔的目标面积与现有面积之差处于目标面积）
，绝对平方误差
。计算的稳定性和难度取决于区域
的严格凸性，密度函数和采样点的分布。
