活性迹线复杂度界、Moreau–Yosida 正则化、未调整朗之万采样这三个词放在一起说明这篇工作不是简单地把朗之万算法换成“加了正则项的版本”而是想回答一个更偏理论、最终却直接影响计算成本的问题当目标分布带有不可微结构甚至本身存在稀疏性时我能不能用“实际被激活的坐标数量”来估计采样需要付出的总代价先说结论式的阅读体验。标题里的 Active-Trace Complexity Bounds最值得关注的是它把采样算法的运行复杂度从“按全维度 d 算”拉回到“按运行过程中真正碰到的维度规模算”。Moreau–Yosida 在这里起着两个作用把非光滑函数磨光同时保留近端算子的几何结构。Unadjusted Langevin Sampling 则说明整个分析框架不依赖 Metropolis 接受步也就是说理论误差可以从离散化过程和非渐近收敛两个方向同时切入。这篇文章适合谁看如果你在做贝叶斯后验采样尤其是带 L1 正则、稀疏约束、图结构约束或者不等式约束的模型你会关心每一步的梯度计算和近端投影能不能做小。如果你在读优化与采样交叉领域的论文你会关心“复杂度界”到底怎么读、怎么和其他方法对比。还有一个容易被忽略的实用点即使你不碰纯理论这种数学结论也会影响你的算法设计。一旦你知道支撑集在长时间运行中不会到处乱跳就能安全地裁剪维度、稀疏存储、减少通信量。下面我按“先理解问题再落地实验再回到理论最后给避坑建议”的顺序拆开讲。1. Moreau–Yosida 和未调整朗之万采样为什么要绑在一起1.1 采样问题的标准形式在统计计算里很多问题最后都变成已知一个未归一化的概率密度[ \pi(x) \propto \exp(-U(x)) ]要从这个分布中采样。通常把 (U(x)) 叫做势能函数。如果是后验采样(U(x)) 里既有负对数似然又有先验项。最常见、最容易处理的情况是 (U(x)) 强凸且光滑这时朗之万采样器的理论已经很完整。简单说朗之万动力学的离散版本长这样[ x_{k1}x_k-\gamma \nabla U(x_k)\sqrt{2\gamma}\xi_k,\qquad \xi_k\sim N(0,I_d) ]其中 (\gamma) 是步长。这个迭代式就是在“沿着负梯度下降”的基础上加入一个高斯噪声项。没有 Metropolis 接受步的叫做 ULA全称是 Unadjusted Langevin Algorithm加了接受步的叫做 MALA。两者区别不只是名字上的而是直接影响理论分析方式如果每次迭代后都做拒绝接受采样器最终的平稳分布会被修正得很好如果直接跑则需要忍受离散化带来的稳态误差但每一步的计算量更小也没有接受率带来的额外方差。非光滑问题就是从这里冒出来的。(U(x)) 里面如果含有 (|x|_1)、指示函数 (I_C(x))、组稀疏正则项或者分段线性项那么 (\nabla U) 在部分点不存在或者不连续。直接套上式梯度无处定义程序要么报错要么只能在某些区域内“碰运气”跑。1.2 莫罗–约西达包的思路Moreau–Yosida 正则化的想法非常直接既然一个函数在某个点不可微那我就用一个稍微打磨过的函数去逼近它同时要求逼近后的函数仍然保持原函数的最低值和最小值点信息。给定闭凸函数 (h(x))定义它的 Moreau–Yosida 包络[ h_\lambda(x)\min_{u}\left{h(u)\frac{1}{2\lambda}|u-x|^2\right} ]其中 (\lambda0) 是正则化参数。(h_\lambda(x)) 是一个全局光滑的凸函数它的梯度非常简洁[ \nabla h_\lambda(x)\frac{1}{\lambda}\left(x-\operatorname{prox}_{\lambda h}(x)\right) ]这里 (\operatorname{prox}_{\lambda h}(x)) 是近端算子定义为[ \operatorname{prox}{\lambda h}(x)\arg\min{u}\left{h(u)\frac{1}{2\lambda}|u-x|^2\right} ]这个公式是整篇文章最核心的工程支点。如果你只需要记住一个结论那就记这个很多情况下计算 Moreau–Yosida 梯度只需要解一个近端问题而近端问题往往有闭式解或者非常快的迭代器。把 (U gh) 拆成光滑部分 (g) 和非光滑部分 (h) 是常见做法。光滑部分保留显式梯度非光滑部分用 (h_\lambda) 替代得到的势能就是[ U_\lambda(x)g(x)h_\lambda(x) ]这个函数的梯度是可算的、全局 Lipschitz 的。用离散朗之万迭代得到类似这样的更新[ x_{k1}x_k-\gamma \nabla g(x_k)-\frac{\gamma}{\lambda}\left(x_k-\operatorname{prox}_{\lambda h}(x_k)\right)\sqrt{2\gamma}\xi_k ]这就是标题里 Unadjusted Langevin Sampling 加莫罗–约西达的直观含义用一把光滑的小尺度熨斗把函数烫平再用普通的梯度朗之万往上走。由于没有 Metropolis 校正最终采样结果不会完美地等于 (\pi)但是通过选择合适的 (\gamma) 和 (\lambda)可以让误差小到实际可接受。2. “Active-Trace”到底在刻画什么2.1 从维度型复杂度到支撑集型复杂度看复杂度界之前要先明确这里谈的复杂度不是“神经网络 FLOPs”那种工程指标。论文里的复杂度界通常指为了把采样分布逼近到某个 W1 距离或者 W2 距离精度 (\varepsilon)算法所需要的迭代次数或者计算单元总量。如果只是普通的高斯型目标所有 (d) 个维度的地位差不多复杂度写成 (O(d^2))、(O(d\log d)) 也合理。可是一旦涉及稀疏目标或者约束问题很多坐标在绝大多数时间根本没有被真正激活。比如稀疏线性回归的后验分布信号可能只有 (s) 个非零中心在采样过程中大部分维度只是围绕 0 做小波动可能连一次有效更新都不需要。这样还拿 (d) 作为复杂度的主要因子显然很浪费。Active-Trace 想要做的事情就是把“被激活的坐标集合”作为运行过程的一部分记录下来。一个比较典型的设定是在每一步 (k)会有一个激活索引集合 (A_k)它表示这一步中计算近端梯度或者更新坐标时需要真正触碰的坐标。那么运行 (T) 步之后总的代价可以写成[ \operatorname{TraceCost}(T)\sum_{k1}^{T}\left|A_k\right| ]如果每个集合 (A_k) 的大小都接近某个稀疏等级 (s)那么总的代价就是 (O(sT))而不是 (O(dT))。在很多统计模型中(s) 远小于 (d)这就是稀疏带来的真正甜头。Active-Trace Complexity Bounds 研究的就是这种累计激活规模它不只看单步迭代贵不贵还要看运行路径上的活跃集合是否稳定、会不会反复横跳。2.2 为什么“活性迹线”比单纯的非零占比更严格有人可能会说我直接数当前向量有多少非零分量不就行了实际上没那么简单。第一非零分量只是结果不一定等于计算代价。一个向量即使有很多维接近零只要代码仍然显式保存稠密向量并且每次做全量矩阵向量乘复杂度依然是 (O(d^2)) 级别。真正的 Active-Trace 还需要描述“算法在实现时不扫描的那些维度”相当于要求算法在结构上能利用稀疏性。第二活动集合可能比当前支撑集大。优化问题里经常出现活动集的概念指的是当前梯度和约束边界同时起作用的维度。采样问题里噪声始终存在即使某一步的数值非常接近零也不能毫无代价地把这个坐标“冻结”太久否则会引入无法恢复的偏差。所以活动集通常是根据阈值、近端结果和方差一起判断的。第三更麻烦的是活动集可能在支撑集边缘来回跳动。采样器不像确定性优化那样收敛到某个支撑集之后就不动了随机噪声总会在阈值边界附近试探。步长大的时候支撑集翻转可能很频繁步长小的时候翻转少了但要达到同一精度需要更多步。活性迹线类的复杂度界就要在“每步少算点”和“总步数多点”之间做权衡。2.3 追踪的数学符号长什么样虽然不同论文的记号有差异但通常可以这样理解。设原始空间维度为 (d)引入一个支撑迹线[ S_k{i: x_{k,i}\neq 0} ]或者带容差的近似支撑[ S_{k,\epsilon}{i: |x_{k,i}|\epsilon} ]然后统计在整个长度 (T) 的运行过程中出现过的不同坐标集合以及累计激活次数。如果算法能证明存在某个常数 (s)使得[ \sum_{k1}^{T}|S_k|\le sT O(1) ]并且 (s) 不随着总维度 (d) 线性增长那就可以说这个算法具有活性迹线型的复杂度改进。这里的难点不是证明“单步成本低”而是证明“在整个随机路径上翻来覆去激活的集合始终被一个小集合管住”。对于很多带 (\ell_1) 先验、稀疏约束或低秩结构的采样问题这一条在实际中的意义比很多全局复杂度公式都重要。运行时间长的 MCMC 程序最怕的就是偶尔激活一个很难算的坐标把平均成本直接拉起来。3. 用一个小例子把迭代过程落地3.1 为什么选 Lasso 型后验作为最小验证环境要用实际代码感受一下标题里的这些抽象概念最方便的实验对象是线性回归加 L1 先验[ yA\theta\varepsilon,\qquad \varepsilon\sim N(0,\sigma^2I) ]取先验为拉普拉斯先验或者等价的 (\ell_1) 正则后验势能写成[ U(\theta)\frac{1}{2\sigma^2}|y-A\theta|^2\rho|\theta|_1 ]这里 (g(\theta)\frac{1}{2\sigma^2}|y-A\theta|^2) 是光滑的二次项(\nabla g(\theta)\frac{1}{\sigma^2}A^T(A\theta-y)) 可以写成矩阵运算。非光滑项 (h(\theta)\rho|\theta|_1) 对应的近端算子是软阈值收缩[ \operatorname{prox}_{\alpha h}(\theta)\operatorname{soft}(\theta,\alpha\rho) ]其中[ \operatorname{soft}(\theta,\alpha\rho)_i \begin{cases} \theta_i-\alpha\rho, \theta_i\alpha\rho\ 0, |\theta_i|\le\alpha\rho\ \theta_i\alpha\rho, \theta_i-\alpha\rho \end{cases} ]这是几乎所有近端算法里的 Hello World。用它来复现莫罗–约西达 ULA 的实验能让我们在最小范围内看到噪声、阈值和活动集三者之间的互动。3.2 一个最小可运行的验证流程不需要把模型做得很大。常用的流程如下。第一步构造稀疏真值 (\theta^*)。取 (d500)真值只有 20 个位置非零设计矩阵 (A) 可以用随机高斯矩阵并做列归一化。第二步固定步长 (\gamma) 和正则系数 (\rho)。如果采样目标是后验(\rho) 本质上对应先验的尺度不能随便调到完全稀疏。步长可以先从 (1/L) 量级开始其中 (L) 是光滑部分的 Lipschitz 常数。对于二次损失(LL_{\max}(A^TA/\sigma^2)) 一般是可控的。第三步采用最简单的近端朗之万更新import numpy as np def soft_threshold(x, threshold): return np.sign(x) * np.maximum(np.abs(x) - threshold, 0.0) def sample_proximal_ula(A, y, theta_init, rho, sigma, gamma, n_iter, record_every10): d theta_init.shape[0] theta theta_init.copy() trace [] AtA A.T A Aty A.T y inv_sigma2 1.0 / sigma ** 2 for k in range(1, n_iter 1): grad_smooth inv_sigma2 * (AtA theta - Aty) # 先做一个半隐式近端更新先显式走光滑梯度再加噪声 theta theta - gamma * grad_smooth np.sqrt(2.0 * gamma) * np.random.randn(d) # 对非光滑部分做近端收缩 theta soft_threshold(theta, gamma * rho) if k % record_every 0: active np.where(np.abs(theta) 1e-6)[0] trace.append((k, active.copy(), np.linalg.norm(theta - theta_true))) return theta, trace注意这个版本是近端朗之万的一种常见离散化并不是所有论文里一模一样的 Moreau–Yosida 变体。更严格地按照前文的莫罗–约西达梯度更新可以把近端放进梯度公式变成theta theta - gamma * ( inv_sigma2 * (AtA theta - Aty) (1.0 / lam) * (theta - soft_threshold(theta, lam * rho)) ) np.sqrt(2.0 * gamma) * np.random.randn(d)这里额外的参数 (\lambda) 控制莫罗–约西达包扎的程度。使用这两种写法的时候请注意理论性质对应的是哪一种离散化标题里如果出现莫罗–约西达通常是用后一种形式推导光滑势能的梯度。如果你用前一种作为实验代码不要直接断言它能复现论文里的全部结论只能说它是在同一类目标上的近似验证。3.3 观察什么结果才算“跑起来了”跑完之后至少要检查三组输出。第一组是轨迹均值与真值的距离。如果每一步都做软阈值收缩采样均值应该最终落在稀疏向量附近。可以从第 1000 步之后开始算均值不要从第一步就算否则初始点的偏差会污染结果。第二组是活动集合长度随时间的变化。绘制 (|S_k|) 与迭代次数的折线图如果在达到一定步数后稳定在 25 到 40 之间并且没有大范围波动说明活跃坐标集合是可控的。这里不要求严格等于 20因为后验的不确定性会让很多弱信号维度在阈值边缘晃动但是如果活动集跳到 200 以上说明 (\rho) 或者 (\gamma) 没有配对好。第三组是“活性迹线复杂度”的粗略指标累计激活次数total_activation 0 for _, active, _ in trace: total_activation len(active) average_activation total_activation / len(trace) print(average_activation)如果把平均激活数算出来只有 30那么算法每一步的实际计算量就可以围绕这 30 个坐标去优化而不是按 500 维全部更新。这个数字就是 Active-Trace 思想最朴素的体现理论界能不能保证它小是工作标题后半段的重点。4. 判断活动集合稳定性与复杂度界的关键信号4.1 先用小步长把小样本跑通我实际测试时最常踩的坑是一开始就把步长调太大。步长大的优点是混合得快看起来很快就接近真值附近。但大步长同时会让软阈值之后的坐标频繁从零跳成非零。原因很简单高斯噪声的方差正比于步长噪声一旦超过阈值坐标就会被激活。活动集一旦不稳定活性迹线复杂度就会从 (O(sT)) 退化为接近 (O(dT))此时不管理论写得多漂亮实际代码都是白算。所以我的做法是先用 500 到 2000 次迭代跑一个小样本记录活动集长度的平均值和最大值。平均值如果超过模型真实稀疏度的 3 到 5 倍不要急着调采样器结构先降低步长或增大收缩阈值把活动集先压下来。4.2 不要把“当前活动集”和“历史活动集”混淆Active-Trace 关心的通常不是某一瞬间的活动集而是运行过程中曾经被激活过的所有坐标。这两个指标差别巨大。比如真实支撑只有 20 个维度但在前 200 步中随机游走误入了 100 个无关维度。虽然停在平衡状态时每一步的平均活动集可能是 25但历史累计活动集可能已经覆盖了 200 多列。理论上如果算法要保证失败重试或者稀疏存储结构覆盖所有潜在激活坐标就必须按照历史累计最大活动集来预分配空间而不是按均值预留。于是边界条件会写得更保守。所以读论文时要特别留意标题与摘要里的 Active-Trace 是“路径上的集合长度”还是“总状态轨迹的字典长度”。前者相当于每步代价后者相当于内存与总通信代价两者之间通常差一个最大激活维度的因子。4.3 固定时间预算下稀疏采样未必更快如果把矩阵乘法写成稠密形式每一步都是 (O(d^2))那么即使理论活动集很小实际运行时间也不会下降。要享受活性迹线复杂度带来的好处得在实现里真正做稀疏化把近端阈值后的零维从矩阵乘法中剔除用列子集生成梯度或者把稠密向量存储改成稀疏哈希结构。若没有这一步复杂度界就只是理论上的可能性不会自动变成墙钟时间上的加速。这一点在复现论文时常被忽略。一个数学定理说复杂度是 (O(sT))通常默认实现里对零分量不做无效计算。如果你的实验环境是自动微分框架一个全尺寸向量参与矩阵乘法后端的稠密内核会把所有维度全部扫一遍理论上的稀疏结论很难直接体现出来。5. 复杂度界背后的典型条件与可扩展性5.1 光滑化参数、步长、精度三者之间的关系在莫罗–约西达正则化中(\lambda) 不是一个可以随意取小或取大的超参数。(\lambda) 控制包络函数与原函数的接近程度但它也直接控制梯度计算的 Lipschitz 常数。当 (\lambda) 很小时(h_\lambda) 接近原来的不可微函数但 (\nabla h_\lambda(x)\lambda^{-1}(x-\operatorname{prox}_{\lambda h}(x))) 的 Lipschitz 常数会变大导致步长必须取小。步长小要走到同样的时间跨度迭代次数就多。当 (\lambda) 很大时函数磨得太狠近端映射几乎把大部分非光滑细节抹掉采样分布会偏离真正的目标分布。所以这里存在一个三向权衡[ \lambda \to 0 \quad\text{提高忠实度但降低光滑性} ][ \gamma \to 0 \quad\text{降低离散化误差但增加步数} ][ T \to \infty \quad\text{提高遍历度但增加累计成本} ]活性迹线复杂度界往往就是把三者之间的权衡关系整合成更紧凑的信息。比单纯问“多少步能收敛到精度 (\varepsilon)”更进一步的是它还考虑每一步的激活成本会不会随维度增长。强凸假设通常是存在的。更粗略地说如果势能在某个方向上太“平”即使活动集很小采样器也需要走很多步才能探明那个方向。因此很多复杂度界里面会出现条件数、强凸参数、光滑 L 以及莫罗–约西达参数 (\lambda)。看条件时应重点对照本文中的设定看原文有没有限定真值支撑集、限制本征维数、限制近端算子结构。5.2 从有限维到高维的直觉推演在高维情况下最理想的结果是非渐近复杂度里含有的维度因子不是总维度 (d)而是某个比 (d) 小的“活动维度”。典型推演如下。如果目标后验的支撑集合只有 (s) 维真正携带信号其余 (d-s) 维的先验把压缩在零附近那么近端算子能做两件事一是把弱信号向零收缩把真正携带能量的坐标分离出来二是把强耦合矩阵 (A^TA) 的列分成“支撑内”和“支撑外”两个部分。支撑外的列与支撑内列的乘积如果很小就可以在更新时用一个小规模子矩阵代替全矩阵。复杂度界就可以表达成类似于[ \mathcal{O}\left(\frac{s}{\epsilon^2}\cdot \text{某个关于条件数的因子}\right) ]而不是[ \mathcal{O}\left(\frac{d}{\epsilon^2}\cdot \cdots\right) ]这种表达特别适合解释“为什么在超高维环境中还要做 MCMC”。很多工程人员默认 MCMC 在 (d1000) 时不可用一部分原因是他们看到的是依赖全维度的复杂度表达式。如果势能高度稀疏支撑集很小那么有效复杂度可以类比到一个更小的流形或坐标子空间跑起来可能没那么可怕。不过这里要提醒活动集本身也是一个随机变量它可能依赖于数据、先验和迭代历史。证明“活动集长时间受控”通常需要更强的几何假设比如限制强凸性、限制特征值条件或者对近端阈值的确定性判断。当数据矩阵列之间高度相关时活动集的翻动可能非常剧烈活性迹线复杂度界就不一定再是鲜艳的 (O(sT))会有很多额外项。6. 实际跑这种采样器时我最常遇到的三类坑6.1 把莫罗–约西达更新里所有项都塞进自动微分莫罗–约西达梯度里包含一个近端算子而 (\operatorname{prox}_{\lambda h}(x)) 本身是一个优化问题的解。如果直接用自动微分去追踪这个算子内部的迭代梯度过程会非常慢还可能出现高阶不可导问题。正确做法是使用近端算子的闭式解。大部分常见正则项如 (\ell_1)、 (\ell_2) 混合范数、约束指示函数、组稀疏正则都有成熟公式。工程代码里只需要把软阈值公式写对不需要像求解一般优化问题那样反复迭代。如果遇到没有闭式解的正则项可以先检查你是否真的需要莫罗–约西达这种形式还是可以直接套近端-朗之万格式先做显式梯度再做近端收缩后者虽然和理论位点未必完全一致但实现稳定。6.2 只看最后均值不看中间活动集传统 MCMC 诊断通常关注 Trace Plot、ESS、Gelman-Rubin 这些指标。对稀疏采样器来说活动集的结构诊断同样重要。一种常用的快速检查是把采样链分成三段分别统计每一段的活动集覆盖范围。三段都稳定在同一个稀疏集合附近的说明后续估计的方差可信三段活动集差异很大说明链还在不同局部结构之间切换后验概率分布可能是多重模态或者步长太大导致阈值效应不稳定。另一种检查是记录每次近端更新后有多少坐标的数值被精确压到零。这个比例如果变化剧烈说明参数区域从“稀疏空间”到“非稀疏空间”的边界被反复穿过此时活动集复杂度会比理论界大很多。6.3 忘记输出阶段的误差仍然存在未调整朗之万采样不像 MALA 那样每步都有接受/拒绝机制来纠正离散化偏差。无论莫罗–约西达正则化做得多好只要步长 (\gamma0)最后得到的采样样本严格意义上不是来自 (\pi)而是来自另一个与之接近的分布。很多实验里最终估计量与真实值之间的误差来源不只是“采样方差”还包括“离散化偏差”。尤其是之后如果要计算尾部概率、分位数这类对分布形状敏感的量偏差可能比方差还大。一个实用补救措施是最后对样本做很小概率的 Metropolis 调整或者直接把最后若干步样本丢弃只用较小的步长重新再跑一小段。6.4 参数扫描时的优先级最后给一个我常用的顺序先固定 (\lambda)扫描 (\gamma)再固定 (\gamma)扫描 (\rho)最后回到 (\lambda) 精细调整。不要一上来就同时调三个参数否则根本分不清活动集变大是因为阈值太低还是步长太大。正确的诊断逻辑是如果均值离真值远先看中心趋势不一定是步长问题可能是链还没有 burning。如果活动集过大先看 (\rho)(\rho) 增大支撑集会变稀疏这是最直接的控制旋钮。如果活动集大小正常但 ESS 很低再看 (\gamma)通常是步长太小导致移动缓慢。如果活动集反复横跳再看 (\lambda)莫罗–约西达参数取太小会让梯度过于陡峭稳定性下降。最后留几个值得继续追问的点这类研究结果给我的感觉不是告诉你“某个算法一定能跑多快”而是把稀疏模型的内在结构纳入算法设计和复杂度分析。Active-Trace Complexity Bounds 听起来非常理论化但它实际回答的是一个工程人员非常关心的问题能不能只算被激活的部分。真正落地时我建议先把一个小模型跑通记录活动集长度、历史累计覆盖范围和由支撑导致的矩阵子集大小。不要绕过这一层直接去复现大定理里的常数。只有当你亲手看到活动集的规模和模型真实稀疏度之间稳定匹配时再回到论文里的边界条件对照你的数据是否满足那时你才算真正读懂了标题里每个关键词的意思。