在贝叶斯推断、机器学习采样和非凸优化的交叉点上有一类算法最近被讨论得越来越多先把目标分布写成一个势函数U(x)再沿着它的负梯度方向加随机扰动做迭代最终让样本点近似服从正比于exp(-U(x))的分布。这条思路从 Langevin 动力学出发经过 Euler–Maruyama 离散化就得到了 ULA即 Unadjusted Langevin Algorithm非校正 Langevin 算法。本文围绕 “Tamed Subgradient ULA beyond Convexity” 这个研究方向展开梳理它的背景、核心原理给出一个可直接运行的 Python 对比实验并总结工程落地时容易踩的坑。适合的读者包括正在研究贝叶斯采样和 Langevin Monte Carlo 的理论方向同学需要在非光滑、非凸目标函数上做采样建模的算法工程师以及想快速理解 Tamed、Subgradient、弱凸这些概念之间关系的入门者。读完本文后你应该能解释 ULA 为什么能采样、直接离散化为什么会发散、次梯度如何进入更新公式、Tamed 为什么能避免梯度爆炸并能在本地跑通一个“双势阱 L1 正则”的完整对比实验。1. 背景与核心概念1.1 从“采样一个非归一化分布”说起很多统计和机器学习任务最终都归结为已知一个分布只给到未归一化的密度π(x) ∝ exp(-U(x))需要从中抽取样本。这里的U通常被称为势函数或者能量函数归一化常数Z ∫ exp(-U(x)) dx在高维情况下几乎不可能精确计算因此我们不能直接做逆变换采样也不能用拒绝采样来硬算。这类问题在贝叶斯后验推断中非常典型后验密度等于似然乘先验再除以证据证据积分通常不可解在生成模型和概率编程中也同样常见。解决思路之一就是构造一条马尔可夫链让它的平稳分布恰好是π然后沿这条链跑足够长的时间把链上的状态当作近似样本。Langevin 算法就是这一族方法里实现最简单、理论最成熟的一种。1.2 什么是 Langevin 算法与 ULALangevin 算法的灵感来自物理中的布朗运动。考虑如下随机微分方程dX_t -∇U(X_t) dt √2 dB_t其中B_t是标准布朗运动。在非常一般的条件下这个扩散过程的平稳分布就是π(x) ∝ exp(-U(x))。可以粗略理解为漂移项-∇U把粒子推向势能低的地方而噪声项√2 dB_t提供了随机探索能力两者平衡之后粒子停留在某个区域的概率密度正比于exp(-U)。为了在计算机上实现我们需要把连续时间离散化。最朴素的做法是用 Euler–Maruyama 格式把上面的 SDE 写成迭代x_{k1} x_k - η ∇U(x_k) √(2η) z_k, z_k ~ N(0, I)其中η是步长。这个迭代格式就是 ULA。所谓 Unadjusted是指它没有像 MALAMetropolis-Adjusted Langevin Algorithm那样在最后加一个 Metropolis-Hastings 接受拒绝步骤而是直接接受每一步的结果。省掉校正步骤让实现更简单代价是离散化会引入偏差步长越大偏差越大步长越小计算成本越高。1.3 次梯度把梯度推广到不可导点∇U存在的前提是U光滑可导。但现实中的目标函数经常不是光滑的例如带 L1 正则的损失函数U(x) f(x) λ||x||_1在x 0处就不可导。这时梯度没有定义我们只能退而求其次使用次梯度subgradient。对于凸函数次梯度的定义是如果对所有y都有f(y) ≥ f(x) g, y - x则g是f在x处的一个次梯度。所有次梯度组成次微分集合∂f(x)。在可导点次微分退化成一个单点集合{∇f(x)}在不可导点它是一个区间或更一般的集合。比如|x|在x 0处次微分是[-1, 1]其中任何一个数都可以作为次梯度参与更新。把 ULA 更新式里的梯度换成任意一个次梯度g_k ∈ ∂U(x_k)就得到 Subgradient ULA。这个替换看起来很简单但理论分析并不是“照抄梯度版本”就能成立的因为次梯度在非光滑点附近的行为更复杂离散化误差也需要重新控制。1.4 Tamed驯化超大漂移项即使有了次梯度ULA 还有一个隐藏的数值问题当∇U或次梯度的增长速度超过线性时例如U(x) (x² - 1)²这种四次多项式势函数其导数4x(x² - 1)在|x|较大时按三次方增长。此时如果迭代点因为噪声偶然跑到远离原点的区域η ∇U(x)会非常大一步更新就可能把点弹到更远的地方形成正反馈最终导致轨迹发散或出现 NaN。这不只是实现细节问题而是 Euler 格式本身在超线性系数下的不稳定性。Tamed驯化的思想最早在随机微分方程数值领域被系统研究代表性工作是 Hutzenthaler、Jentzen 和 Kloeden 提出的 Tamed Euler 格式。核心做法非常朴素在漂移项外面套一个归一化因子让单步位移有界。一种常见形式是x_{k1} x_k - η * g_k / (1 η ||g_k||) √(2η) z_k当η||g_k||很小时η g_k / (1 η||g_k||) ≈ η g_k退化为普通 ULA当η||g_k||很大时漂移位移被压制在常数量级不会出现单步爆炸。相比简单地把梯度截断到某个阈值Tamed 是连续变换数值上更平滑理论分析也更方便。1.5 Beyond Convexity弱凸与非凸势函数早期 ULA 的非渐近收敛分析大多假设U是强凸或凸的。凸假设带来很好的性质目标分布是 log-concave 的梯度指向唯一最小值理论推导可以直接使用强凸不等式。但真实模型很少满足这个条件神经网络后验、双势阱模型、带正则项的鲁棒估计都是典型非凸问题。为了把理论推广到非凸场景研究者引入了弱凸weakly convex假设。函数U称为ρ-弱凸如果U(x) (ρ/2)||x||²是凸函数其中ρ 0。直观理解这个条件允许函数有“负曲率”但负的程度不能超过ρ。弱凸既允许非凸又保留了足够的结构用于分析同时它不要求U可导因此可以自然容纳次梯度。所谓 “beyond convexity”通常就是指把收敛性分析放到弱凸甚至更一般的非凸势函数上。当然全局非凸意味着采样理论通常只能保证收敛到目标分布的某个邻域或者依赖耗散条件把粒子拉回势阱区域这一点在实际使用中要有预期。2. 原理拆解2.1 连续时间 Langevin 动力学的平稳分布为什么dX_t -∇U(X_t)dt √2 dB_t的平稳分布是exp(-U)可以用 Fokker–Planck 方程来理解。设X_t的概率密度为p_t(x)它随时间演化的规律由 Fokker–Planck 方程描述∂p_t/∂t div( p_t ∇U ) Δp_t当p_t不再随时间变化时令右端为零可以验证p_∞(x) ∝ exp(-U(x))是一个平稳解。具体推导只需要求一次散度运算div( exp(-U) ∇U ) Δ( exp(-U) ) exp(-U)( ΔU - ||∇U||² ) exp(-U)( ||∇U||² - ΔU ) 0因此只要我们能够模拟这条 SDE并让轨道跑足够长时间状态分布就会趋于目标分布。连续时间是“精确”的问题只出在离散化上。2.2 ULA 的欧拉离散化Euler–Maruyama 是最直观的离散化在一个小时间步η内认为漂移和噪声都保持不变。这样就得到x_{k1} x_k - η ∇U(x_k) √(2η) z_k这里的噪声方差必须严格等于2η不能随意取其他值因为噪声强度是由 SDE 的扩散系数决定的。如果噪声方差偏大或偏小平稳分布就会偏离exp(-U)。这也是很多自实现采样器常见的错误来源。ULA 是“非校正”的没有 MH 步骤所以每一步都会引入O(η)量级的离散化偏差。理论上的处理方式通常是同时使用两个尺度的步长一个足够小以保证偏差可控另一个足够大以保证快速混合最终在总误差上做平衡。2.3 次梯度 ULA 的更新形式当U不可导时把∇U替换成次梯度g_k ∈ ∂U(x_k)得到x_{k1} x_k - η g_k √(2η) z_k需要注意次梯度的选择不是唯一的。在可导点选哪个都一样在不可导点例如|x|的x 0处次微分是[-λ, λ]选择区间内的不同值会改变单步更新方向但不会改变算法的收敛阶和偏差的本质。实践中通常选0或区间中值因为简单且对称。次梯度 ULA 的直接分析难点在于非光滑点附近的“梯度”是不连续的离散化误差不能再用 Lipschitz 光滑性来控制。一种标准处理是引入 Moreau–Yosida 正则化用光滑近似U_λ代替U再分析近似误差和采样误差的叠加。这类方法也被称为 MYULAMoreau–Yosida ULA。2.4 Tamed 修正的具体作用Tamed 修正的本质是给漂移项加了一个依赖于当前点梯度的缩放。在本文的实现中更新格式写为g_tamed g_k / (1 η ||g_k||) x_{k1} x_k - η g_tamed √(2η) z_k让我们分析两个极端当η ||g_k|| 1即梯度较小或步长很小时g_tamed ≈ g_k算法退化为普通次梯度 ULA保留了原本的收敛性质。当η ||g_k|| 1即梯度巨大时η g_tamed ≈ g_k / ||g_k||单步漂移位移的模大约不超过 1粒子不会因为一步更新而飞出去。与硬截断clip(g, -M, M)相比Tamed 是一个连续函数不会在阈值附近产生跳跃因此对收敛性分析更友好。它相当于让算法在“正常小步长”和“保守有界位移”两种模式之间平滑切换。2.5 收敛性直觉关于 Tamed Subgradient ULA beyond convexity 的理论分析核心工具是建立迭代点分布与目标分布之间的距离上界常用的是 Wasserstein 距离或 KL 散度。典型结论会给出一个形如W_2( μ_k, π ) ≤ 衰减项 步长偏差项的上界。其中衰减项依赖于初始分布与目标的差距以及势函数的耗散条件步长偏差项来自 Euler 离散化和非光滑近似。弱凸假设保证了“负曲率有上界”耗散条件保证了“太远的点会被拉回来”Tamed 保证了“超线性次梯度不会导致轨迹爆炸”。三者合在一起才能得到非凸场景下的非渐近收敛保证。需要强调这类理论通常给出的是“算法与目标之间距离的可控性”而不是“某个局部最优解被找到”。在强非凸、多模态目标上样本可能长时间停留在某个势阱内这是采样的混合问题需要通过多链、退火或更大的噪声来缓解。3. 实验环境与准备工作3.1 运行环境本文示例使用 Python 实现依赖只有numpy和matplotlib不涉及深度学习框架。Python 3.8 numpy 1.20 matplotlib 3.4如果你的环境缺少依赖可以运行pip install numpy matplotlib3.2 示例目标分布为了让实验同时体现“非光滑”和“非凸”两个特点我们选择如下势函数U(x) 0.25 * (x² - 1)² λ|x|其中第一项是一个四次双势阱在x ±1附近有两个势能极小点整体是非凸的而且它的导数x³ - x超线性增长第二项是 L1 正则在x 0处不可导。这个势函数虽然只有一维却能同时触发次梯度、Tamed、非凸分析三个关键点非常适合做概念验证实验。4. 完整实战Python 实现 Tamed Subgradient ULA4.1 构造目标分布与次梯度下面代码定义了势函数和次梯度。np.sign(0)返回 0等价于在[-λ, λ]的次微分集合中选择 0这是最简单且对称的选择。import numpy as np import matplotlib.pyplot as plt LAMBDA 0.5 def potential(x, lamLAMBDA): 势函数 U(x) 0.25*(x^2-1)^2 lam*|x|非凸且非光滑 return 0.25 * (x ** 2 - 1.0) ** 2 lam * np.abs(x) def subgrad(x, lamLAMBDA): 返回一个次梯度光滑部分 x^3 - x加上 lam*sign(x) g x * (x ** 2 - 1.0) g lam * np.sign(x) # np.sign(0)0即次梯度集中选 0 return g4.2 普通 ULA 与 Tamed 版本的更新函数两个更新函数的差异只在漂移项普通 ULA 直接用η * gTamed 版本用η * g / (1 η*|g|)。def ula_step(x, eta, lamLAMBDA): g subgrad(x, lam) return x - eta * g np.sqrt(2.0 * eta) * np.random.randn() def tamed_subgrad_ula_step(x, eta, lamLAMBDA): g subgrad(x, lam) g_tamed g / (1.0 eta * np.abs(g)) return x - eta * g_tamed np.sqrt(2.0 * eta) * np.random.randn()这里有一个需要注意的等价写法也可以把更新写成x - eta * g_tamed展开后是x - eta * g / (1 eta*|g|)与x - (eta * g) / (1 eta*|g|)完全一致都是 Tamed 形式。4.3 轨迹对比检验数值稳定性我们固定一个较远的起点x0 2.0分别用普通 ULA 和 Tamed 版本跑 3000 步观测最大绝对位移。为了防止普通 ULA 发散后污染整个输出先打印统计量再画前若干步轨迹。def compare_trajectory(eta0.3, steps3000, x02.0, seed42): np.random.seed(seed) xs_ula np.empty(steps 1) xs_tamed np.empty(steps 1) xs_ula[0] x0 xs_tamed[0] x0 xa, xt x0, x0 for t in range(steps): xa ula_step(xa, eta) xt tamed_subgrad_ula_step(xt, eta) xs_ula[t 1] xa xs_tamed[t 1] xt return xs_ula, xs_tamed xs_ula, xs_tamed compare_trajectory(eta0.3, steps3000) print(普通 ULA 轨迹的最大绝对值:, np.nanmax(np.abs(xs_ula))) print(Tamed Subgrad ULA 轨迹的最大绝对值:, np.nanmax(np.abs(xs_tamed))) plt.figure(figsize(12, 4)) plt.subplot(1, 2, 1) plt.plot(xs_ula[:300]) plt.title(Ordinary ULA, eta0.3) plt.xlabel(step) plt.ylabel(x) plt.subplot(1, 2, 2) plt.plot(xs_tamed[:300]) plt.title(Tamed Subgrad ULA, eta0.3) plt.xlabel(step) plt.ylabel(x) plt.tight_layout() plt.show()在η 0.3时普通 ULA 很可能在几步之内就飞出[-10, 10]甚至产生无穷大或 NaN而 Tamed 版本会把轨迹限制在势能合理的区域内。这正是 Tamed 修正最直观的收益。4.4 长链采样与密度对比接下来我们验证采样质量用 Tamed Subgradient ULA 跑一条长链丢弃燃烧期样本把经验直方图与理论密度exp(-U)/Z放在一起比较。理论密度用数值积分归一化一维情况下非常可靠。import numpy as np import matplotlib.pyplot as plt def sample_long(step_func, eta, steps120000, burnin20000, x00.0, seed7): np.random.seed(seed) x x0 samples np.empty(steps) for t in range(steps): x step_func(x, eta) samples[t] x return samples[burnin:] eta 0.1 samples sample_long(tamed_subgrad_ula_step, etaeta, steps120000, burnin20000) # 数值计算真实密度 xx np.linspace(-3, 3, 500) true_density np.exp(-potential(xx)) true_density / np.sum(true_density) * (xx[1] - xx[0]) plt.figure(figsize(8, 5)) plt.hist(samples, bins80, densityTrue, alpha0.5, labelTamed Subgrad ULA samples) plt.plot(xx, true_density, r-, linewidth2, labeltrue density exp(-U)/Z) plt.xlabel(x) plt.ylabel(density) plt.title(Sampling comparison, eta%.2f % eta) plt.legend() plt.show()运行后可以看到样本直方图与红色理论密度曲线基本重合说明在步长合适时Tamed Subgradient ULA 能准确逼近目标分布。两个峰的位置大约在±1附近同时由于 L1 正则的存在密度在0附近会有明显的尖角这是非光滑项的典型特征。4.5 结果说明上述实验说明了三件事普通 ULA 在超线性次梯度下可能不稳定Tamed 形式能有效约束单步位移。只要步长合适Tamed Subgradient ULA 仍能保持对目标分布的正确逼近并不会因为 Tamed 修正而显著偏离目标。双势阱示例虽然简单却完整覆盖了非凸势函数、非光滑项、次梯度选择、数值稳定性四个核心问题可以直接作为更大规模实验的测试基线。如果你想要更严谨的结论可以同时跑多组随机种子统计样本的均值、方差和有效样本量并与普通 ULA 做对比表。5. 常见问题与排查思路问题现象常见原因解决思路轨迹出现 NaN 或 inf步长过大或次梯度超线性增长导致单步爆炸减小步长或改用 Tamed / 截断形式直方图与目标分布偏差明显步长太大导致离散化偏差燃烧期太短减小步长丢弃前 10%~50% 样本采样长期困在一个势阱非凸多模态下混合不足增加链长、使用多个随机起点、考虑退火或 tempering结果在不同种子下差异大链长不足样本自相关高增加步数按有效样本量判断收敛在非光滑点处理不当np.sign(0)0只是次梯度集的一种选择选择区间中值或 0 均可不要在同一位置重复切换过大值使用 Tamed 后分布偏移步长仍偏大Tamed 改变了有效漂移尺度适当减小步长再做密度对比下面展开说几个高频问题。第一类问题是发散。如果你看到轨迹数值快速增长首先要检查步长η是否过大。对超线性次梯度即使η 0.1也可能在某些随机种子下发散这是普通 ULA 的已知问题不代表代码写错。改用 Tamed 版本是更稳健的解法。第二类问题是偏差。Tamed 只保证数值稳定不保证离散化偏差消失。η越大Tamed 对漂移的压缩越明显这本身会引入额外偏差。理想做法是用一个小步长做正式采样用稍大步长做烧热或诊断。第三类问题是混合慢。对于双势阱这样的多模态分布单条链可能长时间停留在其中一个势阱。此时不要只增加链长更有效的方法是同时跑多条不同初值的链或者使用并行模拟退火、副本交换等更高级策略。对高维非凸后验采样这是必须正视的问题。6. 工程实践与算法落地建议6.1 步长怎么选步长是 Langevin 类算法最重要的超参数。一个实用的做法是先用不同步长各跑几千步观察轨迹的接受比例、均值和最大位移。由于 ULA 没有 MH 校正无法直接看接受率可以看轨迹的移动范围是否合理移动范围太小说明采样效率低范围剧烈震荡说明步长偏大。理论经验是光滑强凸问题中偏差为O(η)实际中从η 0.1或η 0.01开始调参是常见策略非光滑或高维问题需要进一步缩小。6.2 燃烧期、样本量与自相关马尔可夫链采样必须处理燃烧期。初期样本强烈依赖初始点应当丢弃。经验上至少丢弃前 10%~50% 的样本并用 trace plot 判断“是否已经进入平稳区域”。另外Langevin 链的相邻样本高度相关直接估计均值方差时会低估方差。建议每隔若干步记录一个样本thinning并用有效样本量Effective Sample Size, ESS评估独立样本数。ESS 可以按自相关函数计算ESS ≈ n / (1 2 * sum(rho_k))其中rho_k是 lag-k 自相关系数。6.3 日志与可复现性采样实验的可复现性非常重要。建议在代码中固定随机种子、记录步长、链长、燃烧期和目标函数版本。实际项目中可以设计一个简单的实验配置字典把eta、burnin、steps、seed全部保存到日志中。这样即使换了机器或代码版本也能复现结果并定位问题。6.4 什么时候该用 Tamed不是所有场景都需要 Tamed。如果目标函数的梯度有全局 Lipschitz 界即||∇U(x)||至多线性增长普通 ULA 已经足够稳定Tamed 反而可能增加不必要的偏差。如果势函数是多项式型、指数型或神经网络这类超线性模型特别是你观察到过 NaN 或轨迹飞出的现象Tamed 是成本极低的防御手段。它的实现成本只有一行除法却能显著提升鲁棒性建议在基准测试中同时跑普通版本与 Tamed 版本。6.5 高维与大数据场景的扩展Langevin 类算法的工程价值在高维贝叶斯推断中体现最充分。对于大数据集可以使用随机梯度版本的 SGLDStochastic Gradient Langevin Dynamics用 mini-batch 估计梯度再叠加正确缩放的高斯噪声。对于非光滑正则项可以将 Tamed 思想与近端算子结合每一步先做一次近端映射再处理光滑部分这对应 Proximal Langevin 类算法。高维场景下还要考虑条件数问题必要时加入对角预条件或自适应步长。Tamed 面对的核心矛盾——超线性、非光滑、非凸——在这些扩展中依然存在因此本文的稳定性思路可以迁移过去。7. 总结与学习路线本文从采样问题出发完整介绍了 Tamed Subgradient ULA beyond Convexity 这条技术路线。关键知识点可以概括为ULA 是 Langevin SDE 的欧拉离散化通过“负梯度 噪声”的方式让马尔可夫链的平稳分布逼近exp(-U)当目标不可导时用次梯度替代梯度得到 Subgradient ULA当次梯度超线性增长导致数值不稳定时用 Tamed 形式约束单步漂移当目标非凸时借助弱凸假设和耗散条件完成理论分析。文中的一维双势阱加 L1 正则示例虽然简单却是一个可复现、可调试的最小实验平台。下一步如果你偏理论可以继续学习Langevin 算法的非渐近收敛分析、Wasserstein 距离的工具、Moreau–Yosida 正则化与近端采样、退火与副本交换在多模态采样中的作用。如果你偏工程可以尝试把示例扩展到高维二次型目标、逻辑回归后验采样、稀疏贝叶斯 L1 正则模型用 ESS 和置信区间评估采样质量。动手建议把文中的eta改成0.01、0.1、0.3各跑一遍分别记录普通 ULA 和 Tamed 版本的最大位移、NaN 次数和经验直方图与理论密度的差异。这张小表格做完你对 Tamed 的价值会产生比任何公式都直观的理解。如果本文对你有帮助可以收藏备用后续遇到非光滑非凸采样问题时再回来对照排查。