ADMM算法MATLAB实现：分布式优化与LASSO求解实战指南

简介MATLAB实现的ADMM算法完整代码包聚焦交替方向乘子法的工程落地适合需要求解大规模线性和非线性优化问题的机器学习、图像处理、信号处理领域开发者尤其适用于分布式与并行计算场景。压缩包共383个文件约17.44MB以171个.m脚本和函数为主另含70个.p预编译文件、22个.mat数据文件、多份.dat/.dat-s测试集以及面向Windows、macOS、Linux的mex编译文件覆盖常见运行环境并附有pdf、txt、md等说明文档。已有1575人浏览学习。资源从问题定义、子问题构建到初始化、迭代流程和终止条件均提供可直接运行的MATLAB实现允许替换目标函数与约束用于稀疏信号恢复、图像去噪、网络优化等任务配合丰富的测试集与跨平台编译版本对深入理解算法细节和开展二次开发很有帮助也适合作为课程设计或科研实验的基准工具。1. ADMM 算法的 MATLAB 实现分布式优化入门为何首选它很多人第一次接触 ADMM 是在某个机器学习项目的卡壳节点数据量大、约束多、直接调用 MATLAB 的通用优化器要不内存不够要不迭代上万次还在震荡。我帮某开发者排查过一个多任务稀疏回归问题三个惩罚项叠在一起梯度类算法跑了一整晚还在波动换成 ADMM 之后十几分钟拿到满意解。ADMM交替方向乘子法的核心思路很简单把一个大优化问题拆成两个子问题轮流更新其中一个变量耦合关系交给乘子迭代消化。这个框架天然契合 MATLAB 这种以矩阵为基本运算单元的“开发语言”向量化写子问题求解代码短、可读性强、跑得快。这篇笔记围绕 MATLAB 实现 ADMM 展开给出可直接套用的迭代骨架、参数调整方法和真实踩坑记录适合作论文复现、算法迁移和集群实验前的单机验证。2. 从增广拉格朗日到三个迭代步骤为什么拆分后子问题都有闭式解2.1 交替方向乘子法的数学结构与投影操作ADMM 处理的标准形式是min f(x) g(z) s.t. Ax Bz cf(x) 一般是光滑损失项比如平方误差 0.5||Ax - b||^2g(z) 一般是不可导正则项比如 L1 范数或核范数。把等式约束放进增广拉格朗日函数得到 L_rho(x, z, y) f(x) g(z) y^T(Ax Bz - c) (rho/2)||Ax Bz - c||^2。迭代就是轮流做三次更新固定 z 和 y 对 x 取极小固定新 x 和 y 对 z 取极小最后让乘子沿残差方向上升。这三步看起来很像增广拉格朗日法关键区别是 x 和 z 完全分离。每个子问题只含一个被优化变量另一个被当成常数所以二次项展开后要么是岭回归形式要么是近端算子形式。L1 正则对应软阈值算子核范数对应奇异值阈值算子都有明确的闭式解MATLAB 两三行就能写完。选型上ADMM 特别适合“损失项光滑、正则项不可导”的组合通用内点法需要二阶信息梯度法需要近似可导正则ADMM 直接把已知的近端算子拼起来绕开求导困难。这也是它成为分布式优化默认算法的原因两个子问题可以分开求解只交换少量中间变量。软阈值算子是最常用的近端算子function z_new soft_threshold(v, kappa) % 软阈值收缩算子向量化实现 % v: 输入向量kappa: 收缩量必须非负 z_new sign(v) .* max(abs(v) - kappa, 0); end逻辑说明sign 和 abs 都是逐元素操作max 把阈值内的元素直接置零阈值外的整体朝原点平移 kappa。参数说明kappa 在 ADMM 流程里通常是 lambda/rho即正则强度与惩罚参数的比例它直接决定稀疏度——kappa 越大输出里零元素越多。这段代码没有循环所有操作都是 MATLAB 内置的向量化函数比逐个元素处理的 for 版本快一个数量级是写近端算子的标准习惯。2.2 用 MATLAB 手写 ADMM 迭代不用工具箱也能跑的骨架代码以 LASSO 为例目标是 min 0.5||Ax - b||^2 lambda||x||_1。引入分裂变量 z令 x - z 0标准形式里 A I、B -I、c 0。三个更新步骤都有闭式解x 更新是岭回归子问题z 更新是软阈值y 更新就是乘子上升。function [x, z, history] admm_lasso(A, b, lambda, rho, max_iter, tol) % ADMM 求解 LASSOmin 0.5*||A*x-b||^2 lambda*||x||_1 % 输入 % A m*n 设计矩阵 % b m*1 观测向量 % lambda 正则系数控制稀疏度 % rho 惩罚参数控制收敛速度 % max_iter 最大迭代次数 % tol 停机阈值 % 输出 % x 原始变量解 % z 分裂变量解 % history 结构体记录主残差 r、对偶残差 s 和目标函数值 [~, n] size(A); AtA A*A; Atb A*b; % 对 (AtA rho*I) 做一次 Cholesky 分解后续迭代复用 [L, p] chol(AtA rho*eye(n), lower); if p ~ 0 error((AA rho*I) 非正定请调大 rho 或用伪逆); end x zeros(n, 1); z zeros(n, 1); y zeros(n, 1); history.r zeros(max_iter, 1); history.s zeros(max_iter, 1); history.obj zeros(max_iter, 1); for k 1:max_iter % x 更新岭回归子问题两次三角回代 x L \ (L \ (Atb rho*z - y)); % z 更新软阈值收缩 zold z; z soft_threshold(x y/rho, lambda/rho); % y 更新乘子沿残差方向上升 y y rho * (x - z); % 残差记录 history.r(k) norm(x - z, 2); history.s(k) rho * norm(z - zold, 2); history.obj(k) 0.5 * norm(A*x - b)^2 lambda * norm(x, 1); % 停机判断主残差和对偶残差同时低于阈值 if history.r(k) tol history.s(k) tol history.iter k; break; end end history.iter k; end逻辑说明主体循环严格按三个更新公式执行。x 更新那行看起来是解线性系统实际是利用 Cholesky 分解把 O(n^3) 的求逆降为两次三角回代这是这段代码最值得留意的性能优化点。z 更新里的 x y/rho 是把对偶变量缩放回原始尺度和纯原始形式等价但数值更稳。参数说明max_iter 一般设 500 到 2000 足够tol 不建议直接用绝对 1e-6对于量级很大的 A 可能过严应该改成下一节的相对判据。2.3 收敛判据与停机条件残差指标选哪个、阈值定多少停机判断最常见的坑是只盯主残差。主残差 r^k ||Ax^k Bz^k - c|| 反映原始可行性但还需要对偶可行性s^k rho||A^TB(z^{k1} - z^k)||。两个残差必须同时低于阈值理论收敛才成立。只盯主残差会导致过早停止解还没走到最优面。工程上普遍使用的判据是相对残差epsilon_pri sqrt(n) * eps_abs eps_rel * max(||Ax||, ||Bz||, ||c||) epsilon_dual sqrt(n) * eps_abs eps_rel * ||A^Ty||eps_abs 取 1e-4 或 1e-5eps_rel 取 1e-3 或 1e-2。sqrt(n) 给残差向量加了维度惩罚避免变量量纲影响。我一般把这个判据抽成独立函数因为它在投影、自适应调参、对比验证里都要用function [stopped, eps_pri, eps_dual] check_stopping(r, s, A, B, x, z, y, c, eps_abs, eps_rel) % 基于相对残差判据判断是否停机 % r: 主残差标量, s: 对偶残差标量 n size(x, 1); eps_pri sqrt(n)*eps_abs eps_rel*max([norm(A*x), norm(B*z), norm(c)]); eps_dual sqrt(n)*eps_abs eps_rel*norm(A*y); stopped (r eps_pri) (s eps_dual); end逻辑说明这个函数把残差压缩成标量和阈值对比。阈值里的相对项跟随变量范数变化保证停机判据的尺度自适应——当 x 的量级是 1e5 时固定绝对阈值会导致永远停不下来。参数说明eps_rel 比 eps_abs 更关键设为 1e-3 意味着允许千分之一的相对误差配合回归场景够用追求更精确时再收紧一个数量级。注意Cholesky 分解在迭代前完成一次只要 rho 不变就全程复用这是线性系统求解开销几乎可忽略的原因也是自适应 rho 实现起来要谨慎的原因——rho 一变分解就得重做。提示绝对残差阈值 tol 和主残差量级不匹配时会出现两种极端情况要么迭代上万次达不到停机要么第一轮就误判收敛。优先使用相对判据。3. 参数怎么设惩罚因子 rho 与松弛因子的调参逻辑3.1 rho 对收敛速度的敏感性一个 LASSO 例子rho 是 ADMM 里最容易出问题的超参数它同时扮演两个角色二次惩罚项的权重和乘子更新的步长。rho 偏小乘子更新慢主残差被拖住rho 偏大二次惩罚项太强变量每轮动一点点前期收敛极慢。还有个维度上的坑rho 必须和 A 的尺度对齐A 的元素在 1e-3 和 1e3 两种量级下同一套 rho 的表现完全不同。用实验验证最直观。固定数据生成方式对同一 LASSO 问题分别用 rho 0.01、0.1、1、10 跑一遍% 生成模拟数据 rng(42); A randn(200, 100); x_true sprandn(100, 1, 0.1); b A * x_true 0.01 * randn(200, 1); lambda 0.1; % 扫描 rho 数量级 rhos [0.01, 0.1, 1, 10]; figure; hold on; for i 1:length(rhos) [x, z, history] admm_lasso(A, b, lambda, rhos(i), 800, 1e-6); semilogy(history.r(1:history.iter), LineWidth, 1.5); end legend(rho0.01, rho0.1, rho1, rho10); xlabel(迭代次数); ylabel(主残差);逻辑说明rng(42) 固定随机种子保证重跑结果一致。sprandn 生成稀疏真解让待估计的 x 有可识别结构。semilogy 用对数纵轴方便观察主残差从 1e0 掉到 1e-6 的全过程。参数说明lambda0.1 是根据噪声水平选的如果数据信噪比明显不同lambda 要同步调整否则稀疏解失真。实验里常见现象是 rho0.01 时主残差掉得快但后期对偶残差震荡rho10 时前 50 轮几乎看不出下降一旦进入线性收敛段却又干净利落。调参第一守则先扫数量级确定一个让主残差平稳下降的区间再微调。3.2 自适应惩罚因子与超松弛技巧手动扫 rho 只是起点工程里更省事的是在迭代中自适应调整。判据基于主残差和对偶残差的比值主残差明显大于对偶残差说明原始可行性差增大 rho反过来说明对偶可行性差减小 rho。mu 10; % 调整阈值一般取 2 或 10 tau 2; % 调整倍率常见 2 左右 if r mu * s rho rho * tau; elseif s mu * r rho rho / tau; end逻辑说明mu 越大触发调整越苛刻rho 变化越稀疏tau 控制调整步长tau 太大容易来回震荡取 2 比较稳妥。改 rho 之后前面 Cholesky 分解的 L、p 必须重新计算所以实际实现里要把分解放进循环并加条件判断只有 rho 变化时才重做。超松弛是另一个可选加速项。ADMM 的迭代可以乘一个松弛因子 alpha常见取值 1.5 到 1.8做法是在变量更新前把当前解和上一步解做 alpha 加权混合。alpha 接近 2 时数值上容易发散超过 2 基本必炸。实测里L1 正则且数据质量较差的问题alpha 取 1.5 通常比不松弛快 1.3 倍左右。松弛因子和自适应 rho 可以同时开一个管乘子更新的幅度一个管变量更新的混合比例。建议先固定 alpha1 调好 rho再逐步引入 alpha避免两个机制互相干扰。3.3 主残差与对偶残差的联动参数调试三板斧我把这套调参流程整理成三板斧拿到代码不知道怎么下手的可以直接照做。第一步固定 rho 扫数量级。用 3.1 的脚本观察主残差曲线是否在 200 轮内进入线性下降。下降速度先快后慢rho 偏小从头到尾都慢rho 偏大。第二步画主残差与对偶残差的比值曲线。比值长期大于 10说明原始可行性滞后rho 要加比值长期小于 0.1说明对偶可行性滞后rho 要减。这个比值就是自适应判据的直接依据。第三步开自适应或者定死一个合适值。问题只跑一次就定死 rho问题要反复求解或输入数据变化大上自适应。下面这张表是我在几个模拟项目里积累的经验值适用条件是 A 为标准正态采样问题场景rho 初值mutaualpha普通 LASSOm/n 在 1~511021.5图像去噪约束较强0.11021.8分布式估计通信受限1021.51.0表里的 rho 初值只是入口跑三轮后要按第二步的观察结果修正。当 m 远大于 n 时A^TA 条件数可能很糟rho 初值取 0.1 以下会等来无尽的迭代这时候优先考虑预处理而不是硬调 rho。4. 工程化配置把 ADMM 写进你自己的 MATLAB 算法库4.1 目标函数拆解与算子符封装前两章的 admm_lasso 把问题写死在主体循环里换一个问题就得改主函数这让 ADMM 的优势打了折扣。工程化第一步是把“问题定义”和“迭代骨架”分开。做法是封装两个函数句柄一个描述损失项 f(x) 的最小化操作一个描述正则项 g(z) 的近端算子。主循环完全不动换问题只换句柄function [x, z, history] admm_solver(x_update, prox_g, rho, max_iter, tol) % 通用 ADMM 骨架问题相关部分通过句柄注入 % x_update: 函数句柄接收 (x_old, z, y, rho) 返回新 x % prox_g: 函数句柄接收 (v, kappa) 返回新 z x zeros(size(x_update([], [], [], rho), 1), 1); % 初始化注意维度 z x; y x; for k 1:max_iter xold x; x x_update(xold, z, y, rho); z prox_g(x y/rho, 1/rho); % 注意第二个参数 y y rho * (x - z); history.r(k) norm(x - z); history.s(k) rho * norm(z - xold); end end逻辑说明x_update 接受当前 x、z、y 和 rho返回子问题的解prox_g 接受 v 和 kappa返回近端算子结果。y 更新里的 rho * (x - z) 是残差方向上升主残差 norm(x - z) 是原始可行性度量。参数说明prox_g 的第二个参数传 1/rho 而不是 lambda/rho因为近端算子定义里含正则强度我们把归一化后的收缩量统一放到 kappa 参数里调用方在构造句柄时再把 lambda 算进去。LASSO 场景下x_update 句柄内部就是 Cholesky 分解后的两次三角回代。这种封装让算法库能同时服务多个目标函数代价是多一层函数调用开销单机仿真完全可接受。4.2 停机条件与历史轨迹输出把 check_stopping 集成进主循环让 history 结构体保存每轮的诊断信息。我习惯至少记录四项主残差 r、对偶残差 s、目标函数值 obj、乘子 y 的范数。y 的范数变化能提前暴露发散苗头目标函数反过来验证算法没有跳出可行域。绘图时用 subplot 把三条曲线叠在一屏figure; subplot(3,1,1); semilogy(1:history.iter, history.r(1:history.iter)); title(主残差); subplot(3,1,2); semilogy(1:history.iter, history.s(1:history.iter)); title(对偶残差); subplot(3,1,3); plot(1:history.iter, history.obj(1:history.iter)); title(目标函数值);逻辑说明前两个子图用对数坐标因为残差跨越多个数量级第三个子图用线性坐标方便观察目标函数的波动细节。历史轨迹的另一个用途是在线调参对偶残差一直偏大而主残差已经很小说明对偶可行性拖后腿把 rho 调小比继续死磕更有效反过来同理没有轨迹记录很难定位这种问题。4.3 与 MATLAB 内置最优化函数的对比验证把 ADMM 的解交给内置求解器对一遍是检查实现正确性性价比最高的办法。LASSO 可以等价改写成二次规划min 0.5||Ax - b||^2 lambda * sum(t_i)约束 -t x t。用内置优化器求参考解% 把 LASSO 改写成 QP 并用内置优化器验证 [~, n] size(A); H A*A; f -A*b; % 展开成 [x; t] 变量 H_full blkdiag(H, zeros(n, n)); f_full [f; lambda*ones(n, 1)]; % 线性约束 -t x t A_lb [eye(n), eye(n)]; % x - t 0 A_ub [-eye(n), -eye(n)]; % -x - t 0 A_con [A_lb; A_ub]; b_con zeros(2*n, 1); opts optimoptions(quadprog, Display, none); q_sol quadprog(H_full, f_full, A_con, b_con, [], [], [], [], [], opts); x_qp q_sol(1:n); t_qp q_sol(n1:end);逻辑说明blkdiag 构造分块 Hessian两条线性约束把绝对值约束编码进辅助变量 t。对比时计算 0.5||Ax - b||^2 lambda||x||_1 的目标值ADMM 和 QP 的差值在 1e-4 量级算正常。参数说明内置优化器在几百维以内验证比较稳妥超过一千维建议直接信 ADMM 的残差判断。这版验证脚本跑通后算法库的入口就齐了通用骨架、近端算子句柄、停机函数、验证脚本下一章聊哪些地方容易翻车。5. 避坑与常见问题运行不下去、收敛缓慢的五个真实案例5.1 现象迭代次数过万但主残差不下降现象history.r(k) 在一个较大值附近横盘压力总压不到停机阈值max_iter 设到 10000 也没用。原因rho 与问题尺度严重不匹配。最常见的是 A 的列范数到达 1e4 量级还在用 rho1。二次惩罚项在增广拉格朗日里权重过小约束残差几乎不受拉扯。解决先打印 norm(A, 2) * norm(x_true, 2) 作为参考量把 rho 提到和该量级同阶比如 rho100 或 1000同时把 tol 从绝对阈值改成相对判据。防止再犯每次换数据先用 3.1 节的小脚本跑 20 轮看主残差下降斜率再定 rho 入口。注意rng 种子只影响随机数据生成不影响调参判断每次换数据都要重新扫 rho 数量级。5.2 现象z 更新直接奇异现象运行到某一步弹出“矩阵接近奇异或缩放错误”指向 z 更新内部的线性求解或者 z 出现 NaN 后整个解崩溃。原因常见于软阈值输入里 kappa 出现负值。lambda/rho 在 lambda 为 0 时无意义或者 soft_threshold 里 max(abs(v) - kappa, 0) 的 kappa 被错误地写成 1/rho 而忘记乘 lambda。解决在 soft_threshold 入口加 assert(kappa 0)。lambda0 建议走另一条分支退化成 Least Squares。若 z 更新用的是通用近端求解器而非闭式算子检查子问题的初值是否代入了解空间必要时给它一个微小扰动比如 z 1e-10避免奇异。5.3 现象收敛了但结果不是最优解现象history 里 r 和 s 双双小于阈值算法声称收敛但目标函数值明显高于内置优化器给出的参考值。原因两个子问题没有严格求解到底。x 更新用的是梯度近似而不是闭式解近端算子用了截断迭代次数不够的内层循环都会让外层残差看上去合格但实际 KKT 条件不成立。解决检查 update 句柄。若 x 更新本身是内层优化把内层迭代阈值收紧一个数量级若用闭式解检查矩阵分解时有没有漏掉参数。我给某开发者排查时发现他把 Cholesky 的 lower 参数漏掉分解结果不完整收敛值自然偏差补上后目标函数差值回到 1e-4 以内。5.4 现象数值溢出与 NaN现象迭代进行到几百轮history.r 突然变成 Inf 或 NaN后面全部失效。原因rho 过大导致乘子 y 以步长 rho 倍增一旦主残差短暂偏大y 冲上天随后 x 更新里的 rho*z - y 两个大数相减精度崩塌。另一种来源是输入数据本身含 NaNMATLAB 的 norm 对 NaN 返回 NaN静默传播。解决脚本开头加 if any(any(isnan(A))) 检查输入数据不清洗不进算法。迭代中给 y 加幅值保护比如 if norm(y) 1e12, rho rho / tau, end。这个保护只影响极端情况不干扰正常收敛。看到 NaN 先怀疑乘子爆炸而不是代码错误方向别搞反。5.5 现象目标函数值不单调曲线反复波动现象history.obj 画出来非单调时而下降时而上升初学者以为是代码写错。原因ADMM 本身不保证目标函数单调下降它是逼近 KKT 条件的过程目标函数小波动是正常现象。真正的问题信号是残差不断放大且伴随发散而不是目标函数一朝一夕的波动。解决把目标函数曲线和残差曲线叠加看。目标函数波动但残差持续下降继续迭代目标函数和残差同步上升回到 5.1 检查 rho。我的习惯是不管波动只看 r、s 和最终目标值三条线的关系。想更稳一点可以在内存里保存目标值最小的一次解极端波动下作为最终输出。6. 进阶玩法把 ADMM 用到稀疏低秩分解与分布式估计上6.1 稀疏低秩矩阵分解拆成两个子问题ADMM 的价值不止 LASSO。稀疏低秩分解求 min ||L||_* lambda||S||_1约束 L S M。两个子问题分别是奇异值阈值和软阈值前者在 MATLAB 里是 svd 后对奇异值做软阈值后者就是前面写过的 soft_thresholdfunction L_new svt(M, kappa) % 奇异值阈值算子低秩成分更新的核心 [U, S, V] svd(M, econ); s diag(S); s_new sign(s) .* max(abs(s) - kappa, 0); L_new U * diag(s_new) * V; end逻辑说明svd 的 econ 模式在矩阵行数远大于列数时节省一半计算量。奇异值收缩到非负后低秩成分的秩可能下降这是低秩分解需要的性质。参数说明kappa 在这里是 lambda/rho控制低秩与稀疏的成分划分kappa 越大 L 的秩越低、S 越稀疏。6.2 从单机到分布式多机器并行更新的折中做法ADMM 的分布式形态是每个节点持有一部分数据块本地求 x_i汇总后更新全局 z。真实 MATLAB 场景里大多数人不真上集群折中做法是用 parfor 把子问题的求解拍平在多核上。先确认数据能按行分块本地 x 更新不依赖其他块只有乘子 y 需要汇总。这样单机多核就能测出分布式形态的收敛曲线验证后再迁移到真集群。6.3 验证脚本与性能对比给一个最小验证合成低秩加稀疏矩阵跑完 ADMM 后输出重构误差和两个成分的支撑集大小。用下面的骨架观察重构误差掉到噪声水平需要多少轮rho 1; % 沿用前面调参经验 M randn(100, 100); truth_l randn(100, 3) * randn(3, 100); truth_s sparse(100, 100); truth_s(1:10:end) 1; M truth_l truth_s; % 构造低秩稀疏数据 lambda 0.05; L zeros(100, 100); S zeros(100, 100); Y zeros(100, 100); for k 1:200 L svt(M - S Y/rho, 1/rho); S soft_threshold(M - L Y/rho, lambda/rho); Y Y rho * (L S - M); end逻辑说明这个演示省掉了收敛判据实际使用要补上。运行后观察 L 的秩在迭代中逐步降到真实秩 3S 的支撑集逐步收敛到构造位置重构误差下降到 1e-6 左右说明两个子问题的耦合被乘子轮次消化掉了。进阶到这一步你对 ADMM 的把握已经从“会跑一个 LASSO”升级成“能设计算子、拆耦合、上并行”。每次新问题到手我习惯先用一周目的最小算例跑通再谈调参和加速从那以后我再没在 ADMM 上翻过车。希望帮到你。本文还有配套的精品资源点击获取