非凸约束下基于Landing的扩散模型：原理、算法与应用

1. 引言：当扩散模型遇见非凸约束

在生成式人工智能的浪潮中，扩散模型以其强大的数据生成能力和坚实的数学基础脱颖而出。其核心思想，简单来说，就是教会一个模型如何“去噪”：先通过一个前向过程，将一张清晰的图片（或任何数据）逐步打乱成纯粹的噪声；然后，再训练一个神经网络去学习如何将这个噪声逆向恢复成原始数据。这个过程在数学上被优雅地描述为求解一个随机微分方程（SDE），其逆向过程则对应着分数匹配和朗之万动力学采样。

然而，现实世界中的许多生成任务并非在无边无际的欧几里得空间中自由翱翔。想象一下，你要生成一个分子的三维结构，这个结构必须满足特定的键长、键角等物理化学约束；或者，你要规划一条机器人的运动轨迹，这条轨迹必须避开障碍物（不等式约束）并精确到达目标点（等式约束）。这些约束共同定义了一个复杂的、可能是非凸的“可行域”或“流形”。传统的扩散模型在这个流形上采样时，就像一架没有起落架的飞机，它可能围绕着目标区域盘旋，却很难精准、稳定地“降落”在跑道上。更糟糕的是，一旦采样点稍微偏离流形，后续的生成过程可能会彻底失控，导致生成无效或物理上不可能的结果。

这正是“非凸约束下基于Landing的高效扩散模型”所要解决的核心难题。Landing方法，或者说“着陆”技术，为这架“飞机”安装了一套智能的自动着陆系统。它通过在朗之万动力学的随机微分方程中，巧妙地引入一个额外的“着陆项”，这个项就像一个强大的引力场，一旦采样点开始偏离约束流形，就会产生一个指向流形的强力，将其指数级快速地拉回正轨。这种方法的美妙之处在于，它不仅在理论上保证了过程几乎必然满足约束，而且在离散化的算法实现中，依然能保持高效的收敛性。

本文将深入拆解这一融合了深度学习和约束优化的前沿方法。我们将从最基础的朗之万动力学与分数匹配讲起，逐步引入流形约束的数学描述，然后重点剖析Landing项的设计原理与威力。接着，我们会详细推导过阻尼（OLLA）和欠阻尼（ULLA）两种约束朗之万动力学的具体形式、它们的反向过程，以及如何将其离散化为可执行的算法。最后，我们会探讨两种关键的训练目标：基于证据下界的离散时间训练和更稳健的条件Wasserstein路径匹配。无论你是希望理解其背后精妙理论的研究者，还是寻求在分子生成、机器人学等领域应用约束生成模型的工程师，这篇文章都将为你提供从理论到实践的完整路线图。

2. 理论基础：从朗之万动力学到流形约束

在深入Landing方法之前，我们必须夯实基础，理解扩散模型、朗之万动力学与流形约束这三者是如何联系在一起的。这不仅仅是公式的堆砌，更是理解整个框架设计动机的关键。

2.1 分数匹配与朗之万动力学：生成的核心引擎

扩散模型的前向过程通常被表述为一个随机微分方程。例如，方差爆炸（VE）SDE可以写作：dX_t = σ(t) dW_t其中W_t是标准布朗运动，σ(t)是随时间变化的噪声尺度。这个过程将数据分布q_0逐渐转化为一个简单的先验分布（通常是高斯分布）q_T。

生成的核心在于反向过程。根据Anderson的反向时间随机微分方程理论，上述SDE的反向过程为：dX_t = -σ(t)^2 ∇ log q_t(X_t) dt + σ(t) dW_t这里出现了关键项∇ log q_t(X_t)，即数据在时间t的分布q_t的分数函数。分数函数指明了概率密度梯度最大的方向，即数据最可能聚集的方向。反向过程沿着分数函数的负方向（即指向高概率区域）进行“去噪”。

然而，我们无法直接得到q_t的分数函数。分数匹配的目标就是训练一个神经网络s_θ(x, t)来逼近这个真实的分数∇ log q_t(x)。一旦我们有了一个训练好的分数估计器s_θ，我们就可以用以下朗之万动力学进行采样，从先验分布p_T ≈ q_T开始，生成数据：dX_t = -1/2 * σ(t)^2 s_θ(X_t, t) dt + σ(t) dW_t这可以看作是在能量函数为-log q_t的势场中，加入布朗运动的粒子运动方程，即过阻尼朗之万动力学。

注意：这里有一个系数1/2的差异需要留意。有些文献（如DDPM）的离散化方案会吸收这个系数到噪声调度中。在本文的框架中，我们显式地保留了它，并与物理中的朗之万方程保持一致。在实现时，务必确保你使用的分数模型s_θ的训练目标与你的采样SDE形式相匹配。

2.2 流形约束的数学表述：我们被限制在何处？

现在，我们引入约束。假设我们的数据点x必须位于一个集合Σ中，这个集合由一组等式和不等式定义：Σ = { x ∈ R^d : h(x) = 0, g(x) ≤ 0 }其中h(x) = [h_1(x), ..., h_m(x)]^T是m个等式约束，g(x) = [g_1(x), ..., g_l(x)]^T是l个不等式约束。

为了后续推导方便，我们定义一个统一的约束函数向量J(x)。对于在点x处“活跃”的不等式约束（即g_i(x) ≥ 0，表示点正在边界上或外部，需要被拉回），我们将其与等式约束合并：J(x) = [h_1(x), ..., h_m(x), g_{i1}(x)+ϵ, ..., g_{i|I_x|}(x)+ϵ]^T这里I_x是在x处活跃的不等式约束索引集，ϵ是一个小的正数，用于将不等式边界g_i(x)=0转换为g_i(x)+ϵ=ϵ>0，使得Landing方法能同样处理。当点严格满足不等式约束（g_i(x) < 0）时，该约束被视为非活跃，不纳入J(x)。

这个约束流形Σ的切空间T_xΣ在满足约束规范的点x处，是所有满足∇J(x)^T v = 0的向量v构成的子空间。而投影到该切空间的算子Π(x)是一个关键工具，它可以将任何向量v投影到切空间上，其形式为Π(x) = I - ∇J(x) G(x)^† ∇J(x)^T，其中G(x) = ∇J(x)∇J(x)^T，†表示伪逆。

2.3 约束朗之万动力学的朴素尝试与困境

一个最直接的想法是：在标准的朗之万动力学中，简单地将梯度（或分数）和噪声都投影到切空间上。对于过阻尼情况，这会产生：dX_t = -1/2 * σ(t)^2 Π(X_t) ∇f(X_t) dt + σ(t) Π(X_t) ◦ dW_t其中∇f(X_t)可以代表分数函数s_θ或其他势场梯度，◦表示Stratonovich积分（在流形上更自然）。

这个方案存在一个根本性问题：边界行为。当粒子（采样点）位于流形边界∂Σ时，上述SDE只能保证其瞬时速度在切空间内，但无法阻止布朗运动的随机冲击将其推出流形之外。一旦推出，投影算子Π(X_t)在流形外可能没有良好定义，或者即使有定义，也无法保证过程能主动返回流形。这会导致采样点偏离约束集，产生无效样本。

我们需要一个机制，不仅能保证在流形上的运动是切向的，还能在偏离流形时，产生一个强力的、指向流形的恢复力。这就是Landing方法的用武之地。

3. Landing方法：指数级收敛的约束执行引擎

Landing方法的灵感来源于航天器的着陆过程：它不仅需要水平方向的控制（在切空间内运动），还需要垂直方向的控制（对抗重力，精准着陆）。我们将这个思想数学化。

3.1 Landing项的设计与直观理解

Landing方法的核心是在朗之万动力学中增加一个额外的漂移项：-α σ(t)^2 ∇J(X_t) G(X_t)^† J(X_t)。让我们拆解这个项：

• J(X_t)：约束函数向量。当X_t在流形Σ上时，J(X_t) = 0（对于等式约束和活跃的不等式约束）。当X_t偏离流形时，J(X_t)的范数度量了偏离的程度。
• ∇J(X_t)：约束函数的雅可比矩阵。它的列张成了流形在X_t处的法空间。∇J(X_t) G(X_t)^†这个组合构成了一个从法空间到自身的映射。
• -α σ(t)^2：一个负的系数，与时间相关的噪声强度σ(t)^2成正比。α > 0是一个可调的超参数，控制着陆力的强度。

因此，整个Landing项的方向是-∇J(X_t)的方向，即沿着约束函数梯度下降最快的方向，也就是直接指向流形Σ的方向。其大小与偏离程度J(X_t)成正比。这就像一个弹簧：偏离越远，拉回力越大。

3.2 指数收敛性：为什么Landing如此有效？

Landing项最强大的理论性质是它能带来指数收敛。考虑约束函数J(X_t)随时间的演化。利用Stratonovich微积分链式法则，我们可以推导出：dJ(X_t) = ∇J(X_t) ◦ dX_t = -α σ(t)^2 J(X_t) dt这是一个关于J(X_t)的确定性常微分方程！布朗运动项dW_t的影响在计算dJ时被神奇地消除了，因为噪声项Π(X_t) ◦ dW_t被∇J(X_t)投影后为零（∇J(X_t)^T Π(X_t) = 0）。

解这个ODE，我们得到：J(X_t) = J(X_0) * exp(-α S(t))，其中S(t) = ∫_0^t σ(s)^2 ds。 这意味着，无论初始点X_0离流形多远，约束违反量J(X_t)都会以指数速度衰减到零。衰减速率由α和累积噪声强度S(t)共同决定。对于不等式约束，一旦g_j(X_t) + ϵ衰减到ϵ（即g_j(X_t)衰减到0），该约束被视为满足，并从活跃集I_x中移除，此后过程将保持在g_j(x) ≤ 0的半空间内。

实操心得：α的选择至关重要。过小的α会导致收敛到流形太慢，在离散时间步长下可能始终存在微小偏差。过大的α会使方程刚性增强，需要更小的时间步长Δt来保持数值稳定。在实践中，α通常在1到100之间，需要根据具体问题和σ(t)的尺度进行调整。一个经验法则是，让α σ(t)^2 Δt在0.1到1的量级，以确保每个时间步都能对约束违反进行有效校正。

3.3 整合：完整的约束朗之万动力学（前向过程）

将Landing项与投影的朗之万动力学结合，我们就得到了完整的约束朗之万方程。

1. 过阻尼朗之万动力学与Landing (OLLA - Overdamped Langevin with Landing):dX_t = -1/2 * σ(t)^2 Π(X_t) ∇f(X_t) dt - α σ(t)^2 ∇J(X_t) G(X_t)^† J(X_t) dt + σ(t) Π(X_t) ◦ dW_t这里∇f(X_t)可以替换为分数估计s_θ(X_t, t)。注意，我们通常使用Stratonovich积分以保证流形上的几何特性。如果转换为Itô积分，会产生一个额外的修正项1/2 σ(t)^2 H(X_t) dt，其中H(X_t)与投影算子的导数有关，计算较为复杂。在离散化时，我们有时会采用拉格朗日乘子法来避免直接计算这个修正项。

2. 欠阻尼朗之万动力学与Landing (ULLA - Underdamped Langevin with Landing):欠阻尼方程引入了动量P_t，通常能实现更快的混合速度。

dX_t = σ(t)^2 P_t dt - α σ(t)^2 ∇J(X_t) G(X_t)^† J(X_t) dt dP_t = Π(X_t) [ σ(t)^2 (-∇f(X_t) - γ P_t) dt + σ(t) √(2γ) ◦ dW_t ] - σ(t)^2 ∇J(X_t)^T G(X_t)^† H_1(X_t, P_t) dt

其中γ是摩擦系数，H_1是一个与动量相关的修正项。欠阻尼方程同样需要Landing项来保证X_t满足约束，并且还需要一个额外的“着陆修正项”来保证动量P_t也保持在切空间内（即∇J(X_t)^T P_t = 0）。

4. 反向过程构建：从噪声到数据的约束生成路径

为了训练扩散模型，我们不仅需要前向过程（从数据加噪到先验），更需要对应的反向过程（从先验去噪生成数据）。在流形上，反向过程的构建需要格外小心。

4.1 流形上SDE的反向过程一般理论

对于一个在流形M上由Stratonovich SDE定义的前向过程：dX_t = Π(X_t) b(X_t, t) dt + κ(t) Π(X_t) ◦ dW_t其对应的反向过程\overleftarrow{X}_t（在时间T-t下）由以下SDE给出：d\overleftarrow{X}_t = Π(\overleftarrow{X}_t) [ -b(\overleftarrow{X}_t, T-t) + κ(T-t)^2 ∇ ln q_{T-t}(\overleftarrow{X}_t) ] dt + κ(T-t) Π(\overleftarrow{X}_t) ◦ d\overline{W}_t其中q_t是前向过程在时间t的分布密度，∇ ln q_t就是我们需要用神经网络学习的分数函数。\overline{W}_t是一个独立的反向时间布朗运动。

这个公式是流形上扩散模型的理论基石。它告诉我们，反向过程不仅需要反转前向漂移b的方向，还需要加上一个分数函数项来引导概率流回到数据分布。

4.2 OLLA与ULLA的反向过程

将上述一般理论应用到我们的OLLA和ULLA前向过程，就可以得到具体的反向过程SDE。

对于OLLA，前向漂移b(x,t) = -1/2 σ(t)^2 ∇f(x)，扩散系数κ(t) = σ(t)。代入公式，并加上对应的Landing项以保证反向过程同样满足边界条件，我们得到反向OLLA SDE (Itô形式)：

d\overleftarrow{X}_t = 1/2 σ(T-t)^2 Π(\overleftarrow{X}_t) [ ∇f(\overleftarrow{X}_t) + 2 ∇ ln q_{T-t}(\overleftarrow{X}_t) ] dt + 1/2 σ(T-t)^2 H(\overleftarrow{X}_t) dt - α σ(T-t)^2 ∇J(\overleftarrow{X}_t)^T G(\overleftarrow{X}_t)^† J(\overleftarrow{X}_t) dt + σ(T-t) Π(\overleftarrow{X}_t) d\overline{W}_t

对于ULLA，过程更复杂，因为它涉及位置和动量两个变量。其反向过程SDE如下：

d\overleftarrow{X}_t = -σ(T-t)^2 \overleftarrow{P}_t dt - α σ(T-t)^2 ∇J(\overleftarrow{X}_t) G(\overleftarrow{X}_t)^† J(\overleftarrow{X}_t) dt d\overleftarrow{P}_t = Π(\overleftarrow{X}_t) [ σ(T-t)^2 ( ∇f(\overleftarrow{X}_t) dt + γ \overleftarrow{P}_t dt + 2γ ∇_p ln q_{T-t}(\overleftarrow{X}_t, \overleftarrow{P}_t) ) dt + σ(T-t)√(2γ) d\overline{W}_t ] + σ(T-t)^2 ∇J(\overleftarrow{X}_t)^T G(\overleftarrow{X}_t)^† [ H_1(\overleftarrow{X}_t, \overleftarrow{P}_t) + α H_2(\overleftarrow{X}_t, \overleftarrow{P}_t) ] dt

这里∇_p ln q_t是分数函数关于动量的部分。ULLA反向过程同样包含了Landing项及其动量修正项H_2。

核心要点：这些反向SDE是训练目标的来源。我们的分数模型s_θ(x, t)（对于OLLA）或s_θ(x, p, t)（对于ULLA）的目标就是去逼近这些方程中的∇ ln q_t项。通过离散化这些SDE，我们可以得到用于训练的具体损失函数。

5. 离散化算法：从连续方程到可执行代码

理论SDE必须被离散化才能在计算机上模拟。离散化的选择直接影响算法的稳定性、准确性和效率。

5.1 时间离散化与噪声调度

我们采用均匀时间网格：t_k = k Δt,k = 0, ..., N, 总时间T = NΔt。噪声调度通常采用线性调度：σ_k = σ_min + (k/N) (σ_max - σ_min)。前向过程从k=0（数据）迭代到k=N（噪声），反向过程则从k=N迭代回k=0。

5.2 OLLA的离散化：欧拉-丸山格式

对于OLLA，我们采用经典的Euler-Maruyama (EM) 格式进行离散化。

前向OLLA (Forward-OLLA):X_{k+1} = X_k - (σ_k^2 Δt / 2) Π(X_k) ∇f(X_k) - α σ_k^2 Δt ∇J(X_k)^T G(X_k)^† J(X_k) + (σ_k^2 Δt / 2) H(X_k) + σ_k √(Δt) Π(X_k) ζ_k其中ζ_k ∼ N(0, I)。H(X_k)是Itô修正项，计算复杂。因此，实践中更常用拉格朗日乘子法。

基于拉格朗日乘子的离散化 (Forward-OLLA-P):X_{k+1} = X_k - (σ_k^2 Δt / 2) Π(X_k) ∇f(X_k) + σ_k √(Δt) Π(X_k) ζ_k + ∇J(X_k)^T λ_k然后求解拉格朗日乘子λ_k，使得J(X_{k+1}) = 0。这通常通过牛顿迭代法完成。这种方法隐式地强制执行了约束，避免了显式计算H(X_k)和Landing项，数值上更稳定。

反向OLLA-P的离散化形式类似，只是漂移项符号相反并加入了分数模型s_θ：X'_k = X'_{k+1} + (σ_{k+1}^2 Δt / 2) Π(X'_{k+1}) [ ∇f(X'_{k+1}) + 2 s_θ^{(k+1)}(X'_{k+1}) ] + σ_{k+1} √(Δt) Π(X'_{k+1}) ζ'_{k+1} + ∇J(X'_{k+1})^T λ_{k+1}其中λ_{k+1}通过求解J(X'_k) = 0得到。X'表示反向过程的变量。

5.3 ULLA的离散化：O-BA分裂与折叠技巧

欠阻尼方程是二阶的，直接离散化容易不稳定。我们采用一种称为“O-BA splitting”的一阶算子分裂方法，并结合“折叠”技巧来节省内存。

核心思想：将ULLA的SDE分裂成三个可积部分：

• O步 (Ornstein-Uhlenbeck)：只处理动量的摩擦和噪声项。在冻结位置X_k和σ_k的情况下，这是一个线性SDE，其解是高斯分布。
• B步 (Boost)：处理势场梯度∇f和与约束相关的动量修正项H_1。
• A步 (Advection)：处理由动量驱动的位置更新和Landing项。

折叠技巧：通过代数变换，将位置X_{k+1}的更新公式用前两步的位置X_k和X_{k-1}来表示，从而消去显式的动量变量P_k。我们定义一个近似的伪动量\tilde{P}_k = Π(X_k) (X_k - X_{k-1}) / (σ_{k-1}^2 Δt)。这带来了两个好处：1) 只需存储位置序列，大幅节省内存，尤其在需要保存整个前向轨迹用于训练时；2) 自动保证了\tilde{P}_k位于切空间T_{X_k} Σ的近似。

前向ULLA (Forward-ULLA) 离散化公式：X_{k+1} = X_k + σ_k^2 Δt Π(X_k) [ a_k \tilde{P}_k - σ_k^2 Δt ∇f(X_k) + √(1 - a_k^2) ζ_k ] - α σ_k^2 Δt ∇J(X_k)^T G(X_k)^† J(X_k) - σ_k^4 Δt^2 ∇J(X_k)^T G(X_k)^† [ H_1(X_k, \tilde{P}_k) - α H_2(X_k, \tilde{P}_k) ]其中a_k = exp(-γ σ_k^2 Δt)。同样，拉格朗日乘子法版本 (Forward-ULLA-P) 用求解λ_k使得J(X_{k+1})=0来替代最后两项复杂的显式项。

反向ULLA的离散化思路类似，但时间方向相反，且漂移项中包含分数模型s_θ。其最终形式为：X'_k = X'_{k+1} + σ_{k+1}^2 Δt √(1 - a_{k+1}^2) Π(X'_{k+1}) ζ'_{k+1} - α σ_{k+1}^2 Δt ∇J(X'_{k+1})^T G(X'_{k+1})^† J(X'_{k+1}) - σ_{k+1}^2 Δt Π(X'_{k+1}) [ a_{k+1} \tilde{P}'_{k+1} + σ_{k+1}^2 ( ∇f(X'_{k+1}) + 2γ ( ∇_p ln q_{k+1}(X'_{k+1}, \tilde{P}'_{k+1}) + \tilde{P}'_{k+1} ) ) Δt ] - σ_{k+1}^4 Δt^2 ∇J(X'_{k+1})^T G(X'_{k+1})^† [ H_1(X'_{k+1}, \tilde{P}'_{k+1}) + α H_2(X'_{k+1}, \tilde{P}'_{k+1}) ]其中\tilde{P}'_{k+1} = Π(X'_{k+1}) (X'_{k+2} - X'_{k+1}) / (σ_{k+2}^2 Δt)。拉格朗日乘子版本 (Backward-ULLA-P) 同样用乘子λ_{k+1}替代后两项。

注意事项：ULLA的离散化公式看起来复杂，但结构清晰。√(1 - a_k^2) ζ_k项来自O步的解析积分，a_k \tilde{P}_k是动量衰减项，- σ_k^2 Δt ∇f(X_k)是势场梯度项，最后一大项是约束相关的修正。在实现时，建议先实现拉格朗日乘子版本，因为它更稳定，无需计算复杂的H_1和H_2。

6. 训练目标：DT-ELBO与条件Wasserstein路径匹配

有了前向和反向过程的离散化，我们需要一个目标函数来训练分数模型s_θ。在流形约束下，经典的证据下界和分数匹配需要重新审视。

6.1 流形上的离散时间证据下界

离散时间证据下界旨在最小化数据分布q_0与模型生成分布p_θ^0之间的KL散度。推导基于这样一个事实：如果我们有从数据x_0到噪声x_N的完整前向轨迹{x_0, ..., x_N}，以及参数化的反向转移核p_θ(x_k | d_k)，我们可以构建一个变分下界。

对于OLLA-P，在假设每个离散步都能通过牛顿法成功将点投影回流形（即x_k ∈ Σ）的前提下，可以推导出反向转移密度：p_θ(x_k | x_{k+1}) ∝ exp( - ||Π(x_{k+1}) (x_k - μ^o_{k+1})||^2 / (2σ_{k+1}^2 Δt) )其中μ^o_{k+1} = x_{k+1} + (σ_{k+1}^2 Δt / 2) Π(x_{k+1}) [ ∇f(x_{k+1}) + 2 s_θ^{(k+1)}(x_{k+1}) ]是反向均值。这里s_θ^{(k+1)}是我们要学习的分数网络。

由此，训练损失（忽略常数项）简化为：L_over(θ) = E_{q(x_0:N)} [ Σ_{k=0}^{N-1} ||Π(x_{k+1}) (x_k - μ^o_{k+1})||^2 / (2σ_{k+1}^2 Δt) ]这非常直观：让网络预测的分数s_θ使得根据它计算出的反向均值μ^o，尽可能接近真实前向轨迹中的上一步x_k（投影到切空间上进行比较）。

对于ULLA-P，推导类似但更复杂，因为转移依赖于连续两个状态(x_{k+1}, x_{k+2})。最终损失为：L_under(θ) = E_{q(x_0:N)} E_{p_N} [ Σ_{k=0}^{N-1} ||Π(x_{k+1}) (x_k - μ^u_{k+1})||^2 / (2σ_{k+1}^4 Δt^2 (1-a_{k+1}^2) ) ]其中μ^u_{k+1} = x_{k+1} + σ_{k+1}^2 Δt Π(x_{k+1}) [ a_{k+1} \tilde{p}_{k+1} + σ_{k+1}^2 Δt ( ∇f(x_{k+1}) + s_θ^{(k+1)}(x_{k+1}, \tilde{p}_{k+1}) ) ]，\tilde{p}_{k+1}是近似动量。

6.2 DT-ELBO的局限性与条件Wasserstein路径匹配

DT-ELBO推导依赖于一个强假设：每个离散步x_k都精确位于流形Σ上。这在实际离散算法中很难绝对保证。牛顿法求解拉格朗日乘子可能失败，或者只能达到数值容差内的解。一旦x_k轻微偏离流形，条件密度p_θ(x_k | ...)可能变得病态，导致训练不稳定，甚至无法计算有效的负对数似然。

为了解决这个问题，论文引入了条件Wasserstein路径匹配框架。它放弃了严格的条件密度建模，转而直接最小化前向过程条件分布q(x_k | c_k)与参数化反向转移分布T^θ_{k+1}(· | d_k)之间的Wasserstein-2距离。

核心思想：通过三角不等式，将数据分布与生成分布之间的Wasserstein距离W_2(q_0, p_θ_0)，分解为每一步条件分布之间距离的和。在一定的正则性假设下，可以证明这个距离的上界是可控的。

关键结论：对于OLLA，CWPM导出的训练损失与DT-ELBO的损失在形式上完全一致：L_over_CWPM(θ) = E_{q(x_0:N)} [ Σ_{k=0}^{N-1} ||Π(x_{k+1}) (x_k - μ^o_{k+1})||^2 / (2σ_{k+1}^2 Δt) ]对于ULLA，也有类似结论。这意味着，即使离散轨迹不完全在流形上，我们仍然可以使用与DT-ELBO形式相同的损失函数进行训练，并且该训练在理论上是合理的，旨在最小化生成分布与数据分布之间的Wasserstein距离。CWPM为在非精确约束执行下使用简单的均方误差损失提供了理论保障。

实操心得：在实际训练中，我们几乎总是使用CWPM框架下的损失函数形式。因为它不要求严格的流形约束，对数值误差更鲁棒。同时，损失函数中的权重1/(2σ_{k+1}^2 Δt)（对于OLLA）或1/(2σ_{k+1}^4 Δt^2 (1-a_{k+1}^2))（对于ULLA）非常重要，它们平衡了不同噪声水平下的损失量级，类似于噪声预测网络中的权重调度。忽略这些权重可能导致训练不稳定或效果不佳。

7. 实现细节与常见问题排查

将理论转化为可运行的代码需要关注许多细节。以下是一些关键的实现要点和常见陷阱。

7.1 核心组件实现清单

1. 约束函数与雅可比计算：需要实现J(x),∇J(x)。对于复杂约束，建议使用自动微分框架（如JAX、PyTorch）来计算雅可比矩阵。注意处理不等式约束的活跃集I_x的逻辑。
2. 投影算子Π(x)：计算G(x) = ∇J(x) ∇J(x)^T，然后求其伪逆G(x)^†。对于小规模问题可以直接用SVD，大规模问题可能需要迭代法。注意数值稳定性，当G(x)接近奇异时，可以添加一个小的正则项δI。
3. 拉格朗日乘子求解（牛顿法）：对于OLLA-P/ULLA-P，需要求解非线性方程组J(x_{new}) = 0。这通常通过牛顿迭代完成：λ_{new} = λ - (∇J^T ∇J)^† J(x_pred)，其中x_pred是未投影的预测值。需要设置最大迭代次数和容差。
4. 分数模型s_θ：网络结构需适应数据形态。对于流形数据，输入通常是欧氏空间坐标，但网络需要学习流形上的分数。可以考虑使用等变或不变网络架构，如果数据具有对称性的话。
5. 训练循环：
   • 从数据集中采样x_0。
   • 运行离散的前向过程（OLLA或ULLA），生成轨迹{x_0, x_1, ..., x_N}。保存所有x_k（对于ULLA还需保存x_{k-1}以计算\tilde{p}_k）。
   • 对于每个k，计算目标μ_k（根据OLLA或ULLA的公式）。
   • 计算损失L = Σ_k weight_k * ||Π(x_{k+1}) (x_k - μ_{k+1})||^2。
   • 反向传播，更新s_θ参数。
6. 采样循环：
   • 从先验分布p_N（通常是切空间上的高斯分布）采样x_N。
   • 运行离散的反向过程（使用训练好的s_θ替换公式中的∇ ln q），从k=N迭代到k=0，得到生成样本x_0。

7.2 常见问题与解决方案速查表

7.3 参数选择经验谈

• 时间步长Δt与总步数N：这是精度与效率的权衡。N通常需要几百到几千步。Δt = T/N，T一般取1。可参考σ_max和σ_min的选择，确保噪声调度能平滑覆盖从数据到先验的整个路径。
• Landing系数α：这是最重要的超参数之一。一个可靠的起点是设置α使得α σ_max^2 Δt ≈ 1。这能保证在最嘈杂的阶段也有强约束力。然后根据约束违反情况微调。
• 噪声调度σ(t)：线性调度σ(t) = σ_min + t/T (σ_max - σ_min)简单有效。σ_min应足够小（如0.01），使得最终先验接近可忽略的噪声。σ_max应足够大，使得数据被充分噪声化（如10或50，取决于数据尺度）。
• 摩擦系数γ(仅ULLA)：与Δt相关。一个经验法则是γ ≈ 1/Δt，这样a_k = exp(-γ σ_k^2 Δt)不会太接近0或1，保持动量的有效性。

我个人在实现中的体会是，从OLLA-P开始是一个更稳妥的选择。它概念相对简单，拉格朗日乘子法提供了严格的约束执行，调试起来也更直观。在OLLA-P工作稳定后，再尝试ULLA-P以追求可能的采样效率提升。另外，对约束违反量\|J(x)\|进行全程监控是必不可少的调试手段，它能第一时间告诉你Landing机制是否在正常工作。