基于保凸型Monge–Ampère方程表述的L2最优传输算子分裂方法

《Journal of Computational Physics》：Operator-Splitting Methods for L2 Optimal-Transport Using a Convexity-Preserving Formulation of the Monge–Ampère Equation

　　

针对L2最优传输（Optimal Transport, OT）问题，研究人员提出一种基于保凸型Monge–Ampère方程表述的新型算子分裂（Operator-Splitting）方法，并耦合有限元离散化。在二次代价下，Brenier定理表明最优传输映射为凸势函数u的梯度?u，问题归结为求解带运输边界条件（Transport Boundary Condition, TBC）?u(X)=Y的Monge–Ampère型方程g(?u)det(D2u)=f。由于TBC非标准且直接在全域求解困难，研究人员借助Hessian矩阵特征值恒等式将原方程重构为半线性椭圆方程?Δu=?√(|Δu|2?4det D2u+4f/g(?u))，该形式在稳态自动保持u的严格凸性；引入辅助函数w解耦右端非线性，构造含人工时间的初值问题（IVP），采用Lie分裂在时间上拆为“椭圆演化步”与“松弛耦合步”，并将TBC转化为带法向投影算子N的Neumann条件。空间上采用同一分片线性Lagrange有限元空间近似u、w、?w及D2w，设计梯度恢复算子Gh、Hessian恢复算子Hh（含散度定理内点近似、零Neumann边界修正与Tikhonov正则）及离散投影Nh。该方法可直扩至三维（利用AM-GM不等式推广为?Δu=?d√(Δu)d?dddet D2u+ddf/g(?u)），并通过凸包延拓与Kantorovich对偶处理非单连通源/目标域。二维与三维数值实验在正则及非结构网格上均达到势函数二阶、梯度一阶的最优收敛率，可解带奇性、不同拓扑结构域的OT问题。

研究背景方面，L²最优传输自Monge提出后在网格生成、图像处理、反射天线设计与平均场规划中广泛应用。Brenier定理将二次代价下的最优映射化为凸势梯度，从而归为带运输边界条件（TBC）的Monge–Ampère型方程。Caffarelli与Urbas已给出正密度、凸域下的正则性理论，但数值上直接求解带TBC的方程研究较少：流体法增维代价高，宽模板有限差分或最小二乘优化多限二维Dirichlet条件，Newton型方法需处理斜边界且实现复杂。现有方法难以同时兼顾凸性保持、TBC隐式处理、多维推广及非单连通拓扑。因此研究人员开展本研究，旨在用保凸重构+算子分裂+混合有限元给出易实现、最优收敛且能处理复杂拓扑的算法。论文发表于《Journal of Computational Physics》。

主要关键技术方法上，研究人员首先利用Hessian特征值关系将d维Monge–Ampère方程改写为半线性形式?Δu=?^d√((Δu)^d?d^ddet D²u+d^df/g(?u))，二维取d=2、三维取d=3，该重构在连续层保证λ₂(D²u)>0即保凸；引入辅助w并配γ(w?u)松弛项，得到人工时间初值问题，用Lie分裂拆成“解常系数椭圆方程（Neumann边）+显式松弛”两步；为将TBC ?u(?X)=?Y转为Neumann，设计几何法向位移投影算子N，按点在Y外/上/内分别沿?ν、不动、+ν推进边界Γ=?w(?X)直至贴合?Y；空间离散采用同一P₁有限元空间，通过质量矩阵型恢复得G_h（一阶导）、H_h（三步走：散度定理求内点、零Neumann重算边界点、ε=O(h²) Tikhonov光滑）及N_h，时间步τ~h、α_h~h。非单连通情形通过对源域取凸包X^co延拓密度f^co（退化为障碍问题），对“单连通→非单连通”反向传输则借前向势u用对偶v(y)=max_x∈X(x·y?u(x))求逆映射，样本均来自文内解析密度（高斯、均匀、椭圆/球、环域等）而无外部队列。

研究结果各节浓缩如下。

引言（1. Introduction）指出OT应用广但带TBC的Monge–Ampère数值研究薄弱，研究人员列出四步框架：保凸重构、引辅助函数解耦、Lie分裂+投影处理TBC、混合有限元实现，并宣称可处理非单连通与奇性。

最优传输与Monge–Ampère正则性（2）回顾Brenier定理与Caffarelli定理：目标Y凸、密度有界正时u∈C^2,α(X)，TBC可弱化为?u(?X)=?Y；借Lemma 1给出凸函数边界法向关系ν(?u(x))=(D²u(x))^?1n(x)/∥·∥，为投影算子提供理论依据。

二维分裂策略（3）完成核心推导：3.1用λ_1,2=?(Δu±√(...))证重构式保凸；3.2引w得系统(11)及IVP(12)，γ取?Δ在H₀¹最小特征值；3.3构N，按点相对Y的位置沿±ν移动θ=min(α,|s^Y?s|)；3.4证PDE系统(11)与IVP(12)稳态等价（Theorem 3）；3.5 Lie离散得(16)(17)：向后Euler解椭圆步，显式wⁿ⁺¹=e^?γτwⁿ+(1?e^?γτ)uⁿ⁺¹，并对u作uⁿ⁺¹=u^n+1/2?u^n+1/2(x₀)定常平移。

三维推广（4）用AM-GM将保凸根式写为?Δu=??((Δw)³?27det D²w+27f/g(?w))，指出三维重构本身不保凸，需严格凸初值（如u⁰=?x^TMx, M使MX?Y）并在迭代中可选项目P截断负特征值，文中实验未用投影亦收敛。

非单连通域（5）给出源域取凸包X^co与f^co化归障碍型Monge–Ampère；目标非单连通时借前向u通过对偶(24) v(y)=max_x∈X(x·y?u(x)), r(y)=?v(y)得反向映射，规避分母为零。

混合有限元（6）详述V_h=C⁰∩P₁，G_h由∫D_ih(ψ)?=∫?ψ/?x_i·?恢复；H_h分内点散度型(29)、边界零Neumann重算(33)、Tikhonov(34)三步；N_h在边界节点按ρ(Q_k)相对Y位置沿离散法向±θ_h移；最终格式(37)(38)与(39)(40)分别为二维、三维全离散，初始化用u⁰=w⁰=?x^TMx并平移使u⁰(Q₁)=0。

数值画廊（7）含十个例子：高斯方→矩（表1显L^∞率~1）、均匀方→非均匀方、圆盘→椭圆、断开域正向与反向、环域→盘、凸→非凸、r-adaptive网格生成、方→菱形辨TBC与Neumann差异、立方→长方体、球→椭球；图6示初值近真则Hessian最小特征值全程为正，初值远则先负后正，验证凸性在迭代中恢复。

讨论与结论部分原文指出：研究人员严格建立了(11)与(12)稳态等价，连续层重构保凸但全离散每步不保凸、数值上最小特征值经若干步转正；方法在h加密下对光滑解达势二阶/梯度一阶最优率；τ~h、α_h~h经验选取，过小则停滞；P₁+Hessian恢复比C¹元（如Argyris）自由度少、易扩三维；非单连通通过凸包延拓与Kantorovich对偶融入同一框架，数值实验证实可有效处理奇性与不同拓扑。结论为：所提保凸Monge–Ampère重构+Lie分裂+法向投影N+混合P₁有限元，构成一套易实现、最优收敛、可解二维/三维、单连通/非单连通L²最优传输的统一数值路线，为网格生成与OT应用提供稳健工具。