简介:交替方向乘子法(ADMM)是求解大规模分布式优化与统计学习问题的经典算法,通过将全局问题分解为若干局部子问题迭代求解,在压缩感知、矩阵分解、Lasso回归等场景应用广泛。本份资源定位为新手友好的MATLAB入门示例,面向正在学习凸优化或最优化方法、希望快速掌握ADMM实现细节的本科生与工程师;压缩包内共含2个文件,均为.m格式的MATLAB源码,资源整体仅2KB,代码量小且结构直观,适合逐行跟读与调试。目前已有1565人学习/下载,说明该示例在同类入门资源中具有一定参考热度。从内容来看,两个脚本分别承担主流程与演示入口的职责:主流程文件实现ADMM的迭代更新、残差计算与收敛判断,演示文件则设置具体的二次规划(quadprog)实例并调用主函数运行;读者通过实际运行可以清晰看到每次迭代中原始残差与对偶残差的下降过程,理解变量分裂与乘子更新的作用。同时,简洁的代码便于自行修改目标函数、约束条件或惩罚参数,延伸到图像复原、分布式计算等实际问题,是一份短小精悍的入门实操素材。
1. ADMM 是什么:分布式优化与稀疏问题里的通用拆解法
在做大规模稀疏回归、图像恢复或分布式训练时,目标函数经常长成「光滑损失项 + 非光滑正则项」的加法,比如最小二乘加上 L1 范数。这类问题整体上不好直接求梯度,但每一部分单独拿出来又都能找到快速解法。ADMM(交替方向乘子法,Alternating Direction Method of Multipliers)做的就是这件事:把一个复杂问题拆成两个或更多子问题,轮换求解,再用一个对偶变量把拆分后的解拉回一致。它的核心卖点有三个:每一轮的子问题往往有闭式解或廉价解法、对目标函数的凸性要求比传统对偶法宽松、天然适合按节点拆分做分布式计算。这篇文章会先讲清楚三个更新式是怎么来的,再给一个能直接跑通的 LASSO 最小实现,然后把惩罚参数、停止准则和各类翻车现场逐一拆开,最后落到自适应参数和热启动这类真正提速的技巧上。适合正在做稀疏优化、图像复原,或者想把训练过程拆到多机上的工程师。
2. 三个更新式的来路:增广拉格朗日为什么能稳住对偶法
2.1 先看一个老办法的毛病:对偶上升遇上不光滑函数就振荡
考虑一个典型的可分离优化问题:min f(x) + g(z),约束是 Ax + Bz = c。直接写出拉格朗日函数 L(x, z, y) = f(x) + g(z) + yᵀ(Ax + Bz − c),最朴素的想法是对偶上升:先固定 y,对 x 和 z 联合求极小,再沿残差方向更新 y。这个思路的麻烦在于两点。
第一是联合极小很难做。Ax + Bz − c 把 x 和 z 缠在一起,对 x 和 z 一起求极小,本质上还是在解一个和原问题一样复杂的优化。第二是对偶函数不光滑时,梯度上升会来回振荡。当 f 或 g 不是强凸时,对偶函数可能是分片线性的,对偶变量 y 的迭代会像台阶一样跳,收敛性没有保证。实际跑起来就是:残差忽大忽小,目标函数缓慢爬升,怎么看都不像一个能用的算法。
增广拉格朗日就是冲着这两个问题去的。把约束违背的平方项加进拉格朗日函数,得到一个带惩罚的版本。关键在于这个二次惩罚项。
2.2 增广拉格朗日的 ρ 做了什么,交替方向又解决了什么
增广拉格朗日的定义是:
L_ρ(x, z, y) = f(x) + g(z) + yᵀ(Ax + Bz − c) + (ρ/2)‖Ax + Bz − c‖²
多出来的二次项有三个直接效果。第一,即使 f 和 g 本身不严格凸,加了这个项之后对 x 和 z 的联合极小问题变成了强凸问题,极小值点唯一存在,不再出现 min 取不到的情况。第二,由于二次项是凸的,对偶函数变得更平滑,对偶变量更新时不会剧烈振荡。第三,ρ 的大小直接决定了约束被满足得多快,但它也把 x 和 z 在联合极小中缠得更紧——这正是问题所在:二次项里 Ax + Bz − c 同时含 x 和 z,联合极小依然难。
ADMM 的处理方式非常直接:不联合极小,改成交替做一次 Gauss-Seidel 式的轮换。写成迭代式就是,令 u = y/ρ 为缩放后的对偶变量:
x_{k+1} = argmin_x f(x) + (ρ/2)‖Ax + Bz_k − c + u_k‖²
z_{k+1} = argmin_z g(z) + (ρ/2)‖Ax_{k+1} + Bz − c + u_k‖²
u_{k+1} = u_k + Ax_{k+1} + Bz_{k+1} − c
x 更新这一步只动 x,z 更新这一步只动 z,u 更新就是简单地把残差累加到对偶变量上。把 u 更新两边乘 ρ 还原成 y,就能看出它等价于步长为 ρ 的对偶上升。所以 ADMM 的本质可以概括为:交替最小化两个子问题,中间夹一步固定步长的对偶上升。这样既避开了联合求极小的复杂度,又靠二次项稳住了对偶更新。
收敛性方面,只要 f、g 是闭的真凸函数,且未增广的拉格朗日存在鞍点,ADMM 就能保证目标函数收敛到最优值,残差收敛到零。最坏情形的收敛速率是 O(1/k),对非光滑项也成立。这个门槛比很多一阶方法宽松,是它被大量用在 L1 类问题上的根本原因。常见的适用信号有三个:目标函数可以分裂成多项且每项有自己的快速更新;约束是线性等式;数据规模大到需要分布式或者矩阵大到没法直接求逆。反过来说,如果两个子问题各自都需要内层迭代几十轮才收敛,那每步代价就很高,不适合直接套 ADMM。
2.3 为什么不用一个大梯度方法直接硬解
一个很自然的疑问是:既然有二次惩罚项,直接用梯度下降不好吗?答案在于可扩展性。原问题里 x 和 z 的维度可能非常高,梯度方法每轮都要处理完全耦合的结构,而 ADMM 把问题拆开后,每个子问题的结构都很简单。以 L1 正则为例,z 子问题得到的恰好是一个软阈值算子,不需要任何迭代;x 子问题如果是二次的,则是一次线性求解。这种「每轮成本极低 + 迭代次数可接受」的组合,才是 ADMM 在工程里被选中的真实理由。
3. 用 ADMM 解 LASSO:变量分裂加最小可复现实现
3.1 第一步永远是变量分裂:把目标塞进标准形式
LASSO 的目标函数是 min (1/2)‖Ax − b‖² + λ‖x‖₁。它天生就是两项相加,但直接把整个目标扔给求解器,你得处理 L1 项带来的非光滑性。常见做法是引入一个副本变量 z,强制 x = z,把问题改写成:
min f(x) + g(z),其中 f(x) = (1/2)‖Ax − b‖²,g(z) = λ‖z‖₁,约束 x − z = 0
这一步叫变量分裂。表面上看多了一个变量和一个等式约束,实际收益很大:f 变成纯二次项,更新时只需解一个线性方程组;g 变成纯 L1 项,更新时只需做一次软阈值。原来要专门设计的近端算子,现在被 ADMM 的两个子问题自动消化掉了。
对照上一章的标准形式,这里 A 矩阵等于单位阵,B 等于负单位阵,c 等于零向量。代入更新式会发现符号很干净:x 子问题是岭回归结构,z 子问题是软阈值,u 子问题就是普通累加。LASSO 因此成为测试 ADMM 实现是否正确的最常用试验场。
3.2 最小可复现代码:软阈值算子、预分解与迭代主循环
下面是完整的 Python 实现,依赖只有数值计算库。代码里每一步都留了注释,可以直接保存运行。
import numpy as np def soft_threshold(v, t): # 软阈值算子:L1 范数的近端算子,逐分量操作 return np.sign(v) * np.maximum(np.abs(v) - t, 0.0) def admm_lasso(A, b, lam, rho=1.0, alpha=1.6, max_iter=1000, eps_abs=1e-6, eps_rel=1e-4): m, n = A.shape AtA = A.T @ A Atb = A.T @ b # 系数矩阵 AtA + rho*I 在整个迭代中不变,提前做 Cholesky 分解 L = np.linalg.cholesky(AtA + rho * np.eye(n)) x = np.zeros(n) z = np.zeros(n) u = np.zeros(n) for k in range(max_iter): # x 更新:解 (AtA + rho*I) x = Atb + rho*(z - u) rhs = Atb + rho * (z - u) x_new = np.linalg.solve(L.T, np.linalg.solve(L, rhs)) # over-relaxation:对 x 做外推,alpha=1.0 时退化为原始 ADMM x_relax = alpha * x_new - (1 - alpha) * x x = x_new # z 更新:软阈值,阈值由 lam/rho 决定 z_prev = z z = soft_threshold(x_relax + u, lam / rho) # u 更新:对偶变量累加约束残差 u = u + x_relax - z # 原始残差与对偶残差,两个都要看 r_prim = np.linalg.norm(x_relax - z) s_dual = np.linalg.norm(rho * (z - z_prev)) # 绝对容差 + 相对容差,分别对应原始和对偶可行性 eps_prim = np.sqrt(n) * eps_abs + eps_rel * max( np.linalg.norm(x_relax), np.linalg.norm(z)) eps_dual = np.sqrt(n) * eps_abs + eps_rel * np.linalg.norm(rho * u) if r_prim <= eps_prim and s_dual <= eps_dual: break return x, z, u, k + 1 if __name__ == "__main__": rng = np.random.default_rng(0) m, n = 100, 300 A = rng.standard_normal((m, n)) x_true = np.zeros(n) x_true[:10] = rng.standard_normal(10) b = A @ x_true + 0.01 * rng.standard_normal(m) x_est, z, u, iters = admm_lasso(A, b, lam=0.1, rho=1.0, alpha=1.6) nz = np.sum(np.abs(z) > 1e-4) print(f"iterations={iters}, nonzero={nz}, primal={np.linalg.norm(x_est - z):.2e}")这段代码的要点有三个。第一是预分解:x 子问题每次都是解同一个系数矩阵 AtA + rhoI,所以 Cholesky 分解只做一次,之后每轮只是两次三角回代,单步成本从 O(n³) 降到 O(n²)。注意只要 rho 大于零,AtA + rhoI 一定正定,Cholesky 不会失败。第二是软阈值算子:它就是 L1 近端算子,推导结果是把每个分量向零压缩,超过阈值的部分保留,不超过的直接置零。阈值取 lam/rho,说明 rho 越大,每轮 z 的收缩越轻,这也是下一章要展开的调参关键。第三是 over-relaxation:alpha 取 1.6 时对 x 做外推,在很多问题上能把收敛轮数减少一半左右,alpha=1.0 时就是教科书标准的 ADMM。
运行这段模拟数据,通常几十轮内原始残差和对偶残差就能同时降到约定的阈值,而且 z 里的非零分量数量能正确恢复出真实的稀疏支撑。这个实验是验证实现是否正确的第一步:如果换几组随机种子都不收敛,问题几乎一定出在 rho 和 lam 的相对尺度上,而不是代码本身。
3.3 想直接抄的分布式变体:consensus 形式
ADMM 能火起来,分布式是一个重要推手。常见的做法是 consensus 形式:把总数据切成 N 份,每份一个工作节点,每个节点维护自己的 x_i 和对偶变量 u_i,中心节点只做一次平均。每轮迭代的骨架如下。
# 工作节点 i 的本地迭代 for epoch in range(max_epoch): # 本地步:解一个带二次惩罚的子问题,f_i 是节点自己的损失 x_i = argmin_local(f_i(x) + (rho / 2) * np.sum((x - z + u_i) ** 2)) # 对偶步:本地残差累积 u_i = u_i + x_i - z # 上报 x_i + u_i 给中心节点 # 中心节点收到所有上报值后 z = np.mean([x_i + u_i for i in range(num_workers)])这个形式的好处是节点之间不需要交换原始数据,只需要上报一个向量,通信量和变量维度成正比。本地子问题如果是二次的,仍然用预分解的线性求解;如果 f_i 很复杂,也可以用几步梯度下降做近似更新,这就是所谓的 inexact ADMM。实际操作中,通信频率和本地计算量的平衡需要按集群规模调,但算法本身非常稳。
4. 惩罚参数 rho 与停止准则:收敛速度和精度的三个必调旋钮
4.1 rho 太小会振荡、太大会拖慢:先理解这个矛盾再动手
rho 在算法里同时扮演两个角色:它是增广拉格朗日里约束惩罚的权重,又是对偶变量更新的步长。这两个角色共同决定了迭代的行为。rho 太小的时候,二次惩罚项很弱,子问题接近不严格凸,对偶更新步长也小,典型表现是原始残差来回振荡,像锯齿一样慢慢磨。rho 太大的时候,每轮都会优先把约束满足掉,但目标函数的下降会被牺牲,典型表现是前期残差掉得飞快,后期目标几乎不动,整个迭代被拖进一种「看起来在收敛、实际在空转」的状态。
我一般会先做一个简单的尺度检查:算出 ‖AᵀA‖ 的谱范数,或者至少看看 A 的最大奇异值。rho 比这个量级小几个数量级时,基本可以断定会振荡。相反,如果 rho 比它大几个数量级,大概率是稳定但慢。实际调参顺序是先把数据归一化,让特征方差在 1 附近,然后从 rho=1 出发,在 1e-2 到 1e2 之间按 10 倍步长扫一遍,看残差曲线的形状。不要一开始就追求精细的 rho 值,先确定量级。
4.2 停止准则别只看目标函数:双残差门槛的来路与参数设置
ADMM 的停止准则在工程上容易犯一个错:盯着目标函数的变化量判断是否收敛。这个做法在 ADMM 里不可靠,因为目标函数经常在残差还很远的时候就已经平了。正确的做法是用两个残差分别约束原始可行性和对偶可行性。
原始残差 r = Ax + Bz − c 衡量约束被满足的程度;对偶残差 s = ρ·Aᵀ·B·(z_{k+1} − z_k) 衡量迭代点离对偶可行性有多远。两个都要低于各自的门槛才能停。工程上常用的门槛是:
ε_prim = √p·ε_abs + ε_rel·max(‖Ax‖, ‖Bz‖, ‖c‖) ε_dual = √n·ε_abs + ε_rel·‖Aᵀy‖
其中 p 是约束个数,n 是变量维度。这两条式子里的绝对值项和相对值项分别对应误差预算:绝对项防止数据量级大时阈值过紧,相对项保证和数据本身的尺度挂钩。对我的 LASSO 代码做简化时,ε_dual 用 ‖ρu‖ 代替 ‖Aᵀy‖,因为缩放形式里 y = ρu,数学上等价。
下面是常用的参数起点,按这个表调一般不会跑偏:
| 参数 | 推荐起点 | 作用 | 调节信号 |
|---|---|---|---|
| rho | 1.0(数据归一化后) | 惩罚权重兼对偶步长 | 残差锯齿调大,目标停滞调小 |
| eps_abs | 1e-6 | 绝对容差 | 数据量级大时适当放大 |
| eps_rel | 1e-4 | 相对容差 | 想要更准就降到 1e-5 |
| alpha | 1.6 | 过松弛外推 | 稳定收敛后可试 1.5~1.8 |
| max_iter | 1000 | 兜底上限 | 到上限未收敛说明参数不合理 |
4.3 观察对象是残差曲线,不是目标函数曲线
判断 ADMM 跑得对不对,最快的方法是把 r 和 s 的范数按迭代轮数画出来。如果两条曲线都在平滑下降,说明 rho 的量级基本合适。如果 r 和 s 交替出现一高一低,这是 rho 失衡的典型信号:r 高 s 低说明原始约束满足得慢,rho 偏小;s 高 r 低说明对偶变量更新过猛,rho 偏大。这个诊断方法比盯目标函数直观得多,也是我在调任何 ADMM 实现时第一个看的东西。
5. 避坑:ADMM 不收敛时的排查路径与 5 条血泪经验
5.1 残差锯齿状振荡却不下落:先把 rho 和数据尺度分开看
现象:原始残差和对偶残差的曲线像锯齿一样,每轮上下跳动,整体下降极其缓慢,甚至先降后升。原因:rho 相对数据尺度太小,二次惩罚项的强凸性不足以稳住子问题,对偶更新步长也过小。解决:第一步永远先做特征标准化,让 A 的列范数在 1 附近;第二步把 rho 按 10 倍步长放大重跑,通常放大两到三个数量级后锯齿就会消失。注意不要一边调 rho 一边动数据,那样你永远分不清问题是出在算法还是出在数据预处理。
5.2 目标函数早就平坦但残差仍高:停止准则只盯一半的典型后果
现象:目标函数值在几百轮后基本不动,但把残差打出来发现还比阈值高两个数量级,程序因为 max_iter 到了才退出。原因:只看目标函数变化量做停止判断,而 ADMM 每轮直接控制的是残差而不是目标函数,目标平滑不代表可行点。解决:把原始残差和对偶残差都打印出来,用双门槛停止;如果残差持续不降,回到 4.1 检查 rho 量级。这条经验的核心是:在 ADMM 里,残差才是算法的真实进度条。
5.3 内层迭代求解器拖后腿:子问题精度与外层容差不匹配
现象:x 子问题没法直接求逆,改用迭代求解器(比如共轭梯度)做内层求解,结果外层 ADMM 怎么都不收敛,或者收敛轮数异常多。原因:内层求解决定的步长精度低于外层容差要求,误差在外层迭代里被不断放大。解决:内层迭代的容差至少要设为外层容差的十分之一,或者直接换用预分解缓存。曾见过一个案例,内层容差设成 1e-3,外层要 1e-6,结果 ADMM 永远停在 1e-4 附近,把内层容差收紧后立刻正常。内层解得太精同样浪费,合适的原则是内层比外层高一个数量级即可,不必追求精确解。
5.4 换一批数据同样代码就翻车:量纲漂移让 rho 的含义变了
现象:同一套代码在数据集 A 上收敛良好,换到数据集 B 后同样参数完全发散,看起来像是代码有 bug。原因:数据量纲变了,等价于 rho 被隐式缩放。A 的特征方差是 1,B 的特征方差是 1000,那么对 B 来说有效的 rho 相当于缩了 1000 倍。解决:把特征标准化写进预处理流程,并且把 rho 的选择和 ‖AᵀA‖ 的谱范数挂钩,而不是拍脑袋定一个固定值。这也是为什么很多工具箱在入口处强制做数据缩放,这不是多此一举,是在保护参数的可迁移性。
5.5 自适应 rho 开得越勤振荡越严重:残差平衡判据需要限速
现象:加了自适应 rho 更新后,原本还算稳定的收敛过程反而振荡加剧,甚至比固定 rho 时更差。原因:残差平衡的判据太敏感,每轮都在改 rho,破坏了交替更新建立的平衡。解决:不要每轮都检查,每 5 到 10 轮检查一次残差比;单次调整幅度限制在 2 倍以内;给 rho 设上下界防止被推到极端值。自适应 rho 不是玄学,但它需要限速和限额,否则就是在给算法加噪声。
6. 进阶:自适应 rho、预条件与热启动,把收敛从「能跑」调到「快跑」
固定 rho 的 ADMM 能跑通大多数问题,但想要更快,三个技巧按优先级排序:自适应 rho、预条件、热启动。自适应 rho 的思路是维持原始残差和对偶残差的平衡。当 r 比 s 大很多时,说明约束满足得慢,加大 rho;反过来则减小 rho。一个稳妥的实现是每 10 轮检查一次,比例超过 10 倍才动,单次变化不超过 2 倍,并且给 rho 设一个上下界。
if k > 0 and k % 10 == 0: old_rho = rho if r_prim > 10 * s_dual: rho = min(rho * 2, rho_max) elif s_dual > 10 * r_prim: rho = max(rho / 2, rho_min) # rho 改变后,缩放后的对偶变量要按比例换算,保持 y 不变 u = u * old_rho / rho # 代价:AtA + rho*I 变了,需要重新做一次分解 L = np.linalg.cholesky(AtA + rho * np.eye(n))这段代码的关键是最后两行。缩放后的对偶变量 u 与真正的对偶变量 y 之间满足 y = ρ·u,rho 变了要维持 y 不变,就得把 u 按 old_rho 和 new_rho 的比例换算。代价是预分解矩阵变了,需要重新分解。所以自适应 rho 不是免费的,在单机小问题上收益有限,在大规模分布式场景里反而要谨慎使用,因为重分解和通信都会跟着变。
预条件处理的思路更简单:在算法启动前对 A 的列做对角缩放,让各列范数接近,相当于从源头上消除量纲差异。我一般把这一步和特征标准化合并做,因为它们是同一件事。做完预条件后,rho=1 通常就是一个非常合理的起点,不需要再花时间扫量级。热启动则是针对一组相关问题的场景,比如给一串 λ 值做 LASSO 路径搜索。用上一个 λ 的解作为下一个 λ 的初值,迭代轮数经常能少一个数量级,代价只是多存一份向量,几乎没有副作用。
我自己调 ADMM 的习惯是:先跑一版固定 rho 的实现,把残差曲线打出来确认量级合适;再决定要不要开自适应;最后才考虑热启动。每换一个数据集,第一件事永远是重新检查数据尺度,而不是急着调 alpha 和 ε。因为绝大多数所谓的「不收敛」,拆开来看都是量纲或 rho 量级的问题。这个排查顺序帮我省掉了大量无效调参的时间,也希望帮到你。
本文还有配套的精品资源,点击获取