ADMM算法的MATLAB实现：从原理到LASSO与分布式优化

简介一套基于MATLAB实现的ADMM交替方向乘子法算法代码包主要面向机器学习、图像处理、信号处理等领域的算法研究者和工程师用于求解含等式或不等式约束的凸优化问题。压缩包共383个文件内容约17.44MB涵盖171个m源文件、70个p加密文件、22个dat数据、22个mat数据以及面向Windows/Mac/Linux的mexw64、mexmaci64、mexa64预编译文件并附带pdf、md、txt等格式的说明资料目录按源码、数据、文档分类存放便于对照学习。目前已有四千余人浏览学习是实践ADMM的成熟参考资料。代码主体ADMM.m完整实现了问题定义、参数设置、子问题分解、交替更新与乘子更新等核心步骤配套辅助函数用于特定子问题求解大量dat/mat数据文件可供稀疏回归、低秩矩阵恢复等经典场景直接测试mex文件则提供预编译加速便于性能比较。通过阅读代码与pdf文档可以深入理解ADMM的收敛条件、惩罚参数rho的影响及加速策略适合希望快速上手算法复现、二次开发乃至分布式并行化的中高级用户。1. ADMM算法为什么分布式优化和约束问题都绕不开它ADMM算法全称是交替方向乘子法Alternating Direction Method of Multipliers这几年在压缩感知、低秩矩阵恢复、分布式机器学习、图像去模糊这些场景里几乎成了标配解法。你只要手上有一个带约束的优化问题目标函数还能拆成两块那ADMM基本就是最省心的那一个工具。它的聪明之处在于不跟复杂的原问题硬碰硬而是把大问题拆成小步迭代每一步只需要处理一个简单子问题MATLAB里十几行循环就能搭出主干。这份资源给的就是一套可直接运行的MATLAB实现能帮你绕过从论文公式到可运行代码之间的那段空档期。适合正在做调参、算法对比或毕设复现的从业者。2. ADMM的数学原理与MATLAB落地从对偶上升、增广拉格朗日到三段迭代2.1 ADMM的核心思想对偶上升、增广拉格朗日与交替极小化要理解ADMM得先看它从哪里来。最原始的拉格朗日乘子法把约束罚进目标函数然后同时对原变量和对偶变量做梯度上升这就是对偶上升法的思路。问题是对偶上升要求目标函数有比较强的光滑性而且收敛条件偏紧。后来有人加了二次罚项得到增广拉格朗日函数这样对偶上升的收敛范围宽了很多。但仍有一个麻烦增广拉格朗日里的二次罚项让子问题不再可分遇到大规模问题就没法并行。ADMM做的事情其实很直接假设你要解决的问题是min f(x) g(z)约束是 Ax Bz c它不直接对x和z做联合优化而是把x和z分开轮流更新每一轮只解一个子问题。这种交替极小化的好处是两个子问题往往都有闭式解或者能用很便宜的数值方法求解。配合对偶变量的上升更新最终收敛到原问题的最优解。MATLAB里实现这套逻辑核心就是下面三段式子x_new argmin_x ( f(x) (rho/2) * norm(A*x B*z - c u, 2)^2 ); z_new argmin_z ( g(z) (rho/2) * norm(A*x_new B*z - c u, 2)^2 ); u_new u A*x_new B*z_new - c;第一个式子更新原变量x第二个更新分裂变量z第三个是对偶变量u的上升步。rho是惩罚参数控制收敛速度和数值稳定性。代码看起来简单但每一步的求解才是真正的工作量。2.2 MATLAB实现的最小骨架单约束问题三段循环在一个标准的单约束问题里ADMM的MATLAB实现可以缩到最小骨架。下面这段代码展示的是最简单形态function [x, z, u, history] admm_simple(A, B, c, rho, max_iter, tol) % A, B: 约束矩阵 % c: 约束右端项 % rho: 惩罚参数ADMM的核心超参 % max_iter: 最大迭代次数 % tol: 相对终止容差 m size(A, 1); n size(A, 2); p size(B, 2); x zeros(n, 1); z zeros(p, 1); u zeros(m, 1); history zeros(max_iter, 2); for k 1:max_iter % 更新 x这里假设 f(x) ||x||_2^2 即平方正则 % 闭式解来自求解线性方程组 x (A*A rho*eye(n)) \ (A*(c - B*z - u) * rho); % 更新 z这里假设 g(z) indicator_{z0} % 投影到非负象限 z max(0, c - A*x - u); % 对偶更新 u u (A*x B*z - c); % 记录主残差和对偶残差 history(k, 1) norm(A*x B*z - c, 2); history(k, 2) norm(rho * A * B * (z - z_old), 2); end end注意上面这段代码里x的闭式解写成了一个线性方程组的求解。实际使用时要先算(AA rhoI)的分解不要每次迭代都重新求逆。另外z的更新用了非负投影这是g(z)为指示函数时的特例换成别的正则项就要改成对应的近端算子。对偶残差的定义在不同论文里写法略有差异我在这里用了带rho和AB的形式。2.3 为什么用ADMM与你见过的梯度下降和近端梯度的区别很多人在接触ADMM之前已经用过梯度下降或者近端梯度法。梯度下降只适合光滑目标函数近端梯度法能处理一部分非光滑项但它们都要求目标函数整体结构比较规矩。ADMM的出发点不同不试图一次处理整个目标函数而是把目标拆成f和g两块各自有各自的解法只要能算近端算子就能塞进ADMM框架。举个例子LASSO问题里f是平滑的二次损失g是L1范数。用近端梯度做L1近端算子本身很简单但梯度步的步长选择会影响收敛而ADMM里g对应的z-update就是一个软阈值操作没有步长参数调的东西只有rho。这一点在实际使用中的意义很大——少一个超参调参的麻烦就少一截。还有一个实际场景上的区别分布式或联邦式的优化环境里数据可能分布在不同的节点上。梯度下降需要汇总梯度ADMM可以通过交替更新的方式让每个节点独立计算自己的子问题最后只同步对偶变量。这一点在集中式的MATLAB实现里体现不明显但理解了这个背景你就知道ADMM的变量分裂结构不是偶然。3. 以LASSO为例的完整MATLAB实现目标函数拆解与主循环写法3.1 LASSO问题的ADMM重写目标函数拆解与变量分裂LASSO问题的标准形式是min 0.5 * ||A*x - b||_2^2 lambda * ||x||_1要让ADMM上场得先做一步变量分裂把L1范数“隔离”到z上。常见做法是引入z并加约束x - z 0于是问题变成min f(x) g(z)约束 x - z 0其中f(x) 0.5 * ||A*x - b||_2^2g(z) lambda * ||z||_1。这样拆完之后x-update变成了一个带二次正则的最小二乘问题它的解就是线性方程组的解而z-update是对L1范数的近端算子即软阈值操作。拆解的关键是变量分裂这一步的选择。不同的问题分裂方式可以完全不同。比如有的场景是拆出两块目标函数有的场景是拆出约束矩阵A和B。LASSO这种拆法是最直观的把难啃的L1正则项放到一边把光滑的二次损失放到另一边。3.2 主迭代循环x-update、z-update与对偶更新下面给出一个完整的MATLAB函数输入设计矩阵A、观测向量b、正则系数lambda和ADMM参数输出解和收敛历史function [x, history] admm_lasso(A, b, lambda, rho, alpha, max_iter, tol) % admm_lasso - 用ADMM求解LASSO问题 % 输入 % A - m*n 的测量矩阵 % b - m*1 的观测向量 % lambda - L1正则系数非负 % rho - ADMM惩罚参数正数 % alpha - 松弛因子通常取1.0或1.51会加速但可能不稳 % max_iter- 最大迭代次数 % tol - 相对残差容差 % 输出 % x - n*1 的最优解 % history - 迭代历史主残差、对偶残差、目标值 [m, n] size(A); x zeros(n, 1); z zeros(n, 1); u zeros(n, 1); % 预分解x-update需要反复求解(AA rho*I)这个系统 % 先算Cholesky因子后面每次迭代只做回代省掉重复分解 [L, p] chol(A*A rho*eye(n), lower); if p ~ 0 error(矩阵(AA rho*I)必须正定。); end Atb A*b; history zeros(max_iter, 3); for k 1:max_iter % x-update: 求解 (AA rho*I) * x Ab rho*(z - u) % 这一项来自对f(x)求导并令导数为零 x L \ (L \ (Atb rho*(z - u))); % z-update: 软阈值 % 对每个分量执行 shrink(z u, lambda/rho) zold z; z max(0, 1 - lambda/(rho*abs(x u))) .* (x u); % u-update: 对偶变量沿残差方向上升 u u alpha * (x - z); % 收敛历史记录 history(k, 1) norm(x - z, 2); history(k, 2) norm(rho * (z - zold), 2); history(k, 3) 0.5*norm(A*x - b, 2)^2 lambda*norm(x, 1); % 终止判定主残差和原始残差都小于阈值才停 if history(k,1) tol history(k,2) tol break; end end history history(1:k, :); end逻辑说明这段代码的循环体完全对应ADMM的三个基本步骤。x-update使用预分解的Cholesky因子做快速回代z-update实现的是软阈值近端算子u-update则是对偶上升。整体上这个写法与标准ADMM论文里的形式一致只是加入了alpha松弛因子这让代码可以在alpha1.5时获得比alpha1更快的收敛但代价是可能出现轻微震荡。参数说明lambda控制稀疏度lambda越大解越稀疏lambda norm(A*b, Inf)时解全为0这是一个有用的上界rho是ADMM的收敛速度调节器太大导致收敛快但数值精度下降太小导致迭代次数激增alpha是松弛因子通常取1.0加速到1.5之前要先确认残差曲线没有明显的震荡。3.3 收敛参数与输入输出的定义边界输入参数这一块有几个边界条件值得说明。首先是lambda的范围lambda必须是大于等于0的实数如果输入负值软阈值操作里的(1 - lambda/(rho*abs(xu)))可能为负导致z的符号被反向修正整个算法的收敛无从谈起。其次是rhorho不能为0否则增广拉格朗日的二次罚项消失ADMM退化回原始拉格朗日法收敛条件变得极其苛刻。另外要留意的是max_iter和tol的组合。新手容易把tol设得太小比如1e-10但数值精度和机器epsilon有限这种设定往往白跑几千轮还停不下来。合理的做法是先跑到固定迭代次数画残差曲线再决定tol的等级。上面的代码里如果矩阵A的规模较大最好把Cholesky分解换成基于QR或SVD的预处理版本否则每次调用函数时的一次分解开销依然不小。常见做法是在函数外部完成分解把因子传入迭代主体。4. 参数设置与收敛性判定惩罚参数rho、松弛因子alpha与终止条件4.1 惩罚参数rho收敛速度与数值稳定性的天平rho是ADMM里最重要的超参数它的作用在增广拉格朗日函数里看得最清楚二次罚项的系数。rho越大每一步迭代对约束违反的惩罚越重迭代会更快地让Ax Bz - c趋近于零但代价是原变量x的更新步长变小可能导致收敛精度下降甚至在最优解附近来回震荡。rho越小约束违反在早期被容忍迭代过程平滑但很慢。实际调rho的方法我一般习惯先在log空间里试几个值比如0.01、0.1、1、10、100各跑200轮对比主残差收敛曲线然后选曲线下降最快且没有明显震荡的那个量级。注意一个问题rho的选择和问题本身的尺度强相关。如果A是0到1之间的数据rho1通常合理如果A的数值在10^3量级rho可能需要调到10^3以上才有反应。若rho的取值与A的尺度严重不匹配会出现一种极端的现象主残差前期快速下降后期长期不动而z迭代还在缓慢漂移。这时候强行加大rho并不会改善反而是先把数据标准化更有效。4.2 松弛因子alpha加速还是添乱alpha是ADMM的松弛因子标准的ADMM里alpha1就是普通的三段迭代。当alpha1时对偶更新的方向更激进可以在某些问题里把收敛速度提升不少特别是对于二次损失加L1正则这类结构的问题。但alpha不是越大越好。实践中alpha超过1.8以后残差曲线往往出现锯齿状波动甚至发散。原因在于对偶变量u更新里的(1 - alpha)项会引入额外的振荡模式。如果是第一次调这个参数建议只试1.0和1.5两个值对比总的迭代次数不要贪心。很多现成实现里alpha是写死为1.0的。如果你手里的资源代码里有alpha这个参数说明作者考虑到加速场景这是一个加分项。但用的时候要清楚alpha只在目标函数光滑性较好的那段收敛前期有明显收益到了末段收敛区域它的作用会变小。4.3 终止条件怎么设双重残差判定与自适应rho终止条件是ADMM实现里最容易糊弄的部分。只看主残差norm(x - z)小于阈值就停是对偶残差可能还很大算法实际上并没有收敛到鞍点。正确的做法是同时检查主残差和对偶残差两个都小于各自阈值才算完成。阈值的设置也有讲究。相对残差比绝对残差更适合跨问题使用把主残差除以max(norm(Ax), norm(Bz), norm(c))得到一个尺度无关的指标。下面给出一个可用的判定片段eps_pri tol * max([norm(A*x,2), norm(B*z,2), norm(c,2)]) tol_abs; eps_dual tol * norm(A * u, 2) tol_abs; stop (r_pri eps_pri) (r_dual eps_dual);这里tol_abs是绝对容差通常取1e-6量级。tol是相对容差工程上取1e-4到1e-6之间比较常见。如果你发现算法迟迟不终止但残差曲线已经走平先检查是不是tol设得太小再看残差是压在哪个水平走平的——这往往能直接告诉你rho应该增大还是减小。还有一种进阶做法是自适应rho即每隔若干轮根据主残差与对偶残差的比值调整rho的值。常见策略是当主残差大于对偶残差若干倍时加大rho反之减小rho。自适应rho能省去手动搜索的麻烦代价是额外几个参数调整间隔、放大倍数、缩小倍数。下面给一个常规实现片段if mod(k, 10) 0 if r_pri 10 * r_dual rho rho * 2; % 注意改变了rho之后之前预分解的Cholesky因子失效需要重新分解 elseif r_dual 10 * r_pri rho rho / 2; end end这里有个容易翻车的细节rho一变x-update里的(AA rho*I)整个矩阵都变了必须重新做Cholesky分解否则用的是旧因子算出来的x方向不对。每10轮重分解一次代价通常可以承受。5. ADMM避坑指南参数失误、残差误判与实现细节5.1 迭代发散或出现NaNrho与步长设置不当现象运行代码没几轮变量就变成了NaN或者Inf程序直接报错。原因最常见的是rho设得过大导致x-update的线性系统数值上接近奇异或者是数据矩阵A本身的条件数太大比如A里的元素跨越了10^6的量级这时候(AA rho*I)的条件数依然很糟糕解出来的x方向有巨大误差。另一个常见原因是软阈值操作里分母上出现了abs(xu) 0的情况在MATLAB里虽然不会直接报错但结果会变成NaN。解决先把rho往小调比如从10000调到1或者0.1。同时检查A的数值范围如果A的元素量级差异太大先对列做标准化再进ADMM。对于软阈值的数值保护可以在abs(xu)下面加上一个小常数eps1e-12max(0, 1 - lambda/(rho*(abs(xu)eps)))。这一步是血泪经验MATLAB里的0/0结果是NaNNaN会顺着迭代一路传染。5.2 只盯着主残差对偶残差滞后导致的误判现象main残差已经到了1e-6看起来一切正常但最终的解和真实最优解相比有明显偏差。原因主残差衡量的是x和z的一致性对偶残差衡量的是对偶变量的收敛程度。如果只以主残差作为终止条件可能在对偶变量还在漂移的时候就停了。对偶残差反映了当前迭代点到最优鞍点的距离两个残差都小才是真正的收敛。解决终止条件必须写成同时满足主残差和对偶残差的阈值。另外在观察收敛曲线时把两条残差曲线画在同一张图上如果发现主残差已经下降但对偶残差还高出一个量级那就说明当前的rho太小对偶更新太慢把rho调大一点就能拉平两条曲线。5.3 大矩阵场景每次重复求逆缓存分解结果现象矩阵A的规模在几千乘几千代码每跑一轮迭代都要花很长时间但看CPU占用率却不高。原因典型的实现问题每轮迭代都执行(AA rhoI) \ something。MATLAB里反斜杠运算符会先判断矩阵结构并做分解反复调用意味着反复做LU或Cholesky分解而这一步的复杂度是O(n^3)完全淹没了回代的O(n^2)开销。解决进入循环之前先做一次性分解。如果矩阵对称正定用Cholesky如果不保证正定用LU。上面的LASSO代码已经展示了这个技巧。如果rho在迭代中会自适应变化分解必须放在rho更新之后重新做。这里有一个取舍自适应rho的更新间隔不能太短否则频繁重构因子反而拖慢整体速度。5.4 可变rho的震荡更新间隔与边界限制现象使用自适应rho之后残差曲线不是平滑下降而是明显震荡甚至发散。原因rho的更新过于频繁或者调整倍数太大。比如每轮迭代都检查残差比并翻倍rho会让增广拉格朗日的权重大幅摆动收敛过程完全被rho的调整节奏带着走。解决把rho的调整间隔加大至少10轮一次。另外调整倍数不要超过2倍。同时给rho加一个上下界比如rho_min 1e-4rho_max 1e4防止极端情况下rho连续翻倍导致数值溢出。还有一个细节改变rho之后对偶变量u也应该按比例缩放不然旧的u在新尺度下失去了参考意义。标准做法是u u * (rho_old / rho_new)。5.5 不可分目标直接套用用变量分裂改写现象把某个带有非光滑耦合项的目标函数直接套进ADMM框架结果迭代完全不收敛或者z-update根本算不出来。原因ADMM要求目标函数能拆成两个独立部分f(x)和g(z)如果目标函数里有交叉项x.*z这种不可分项把它硬塞进f或者g都是不对的因为子问题不再可分离。解决遇到这种情况先做变量分裂把交叉项用额外约束隔离。比如原问题里有x1*x2这种项可以引入变量z并加约束x1 z让x2和z相互作用。具体做法因问题而异总体上是要保证每个子问题都能解或者能近似解。如果实在无法拆开ADMM不是你该用的工具考虑换成近端梯度或者原始-对偶算法可能更合适。6. 验证ADMM正确性的实用技巧从已知解到KKT残差6.1 构造已知解的反向验证拿到代码之后第一步不应该是直接跑真实数据而是构造一个已知的稀疏解。先生成随机的高斯测量矩阵A用稀疏位置随机生成x_true计算b A*x_true然后跑ADMM。因为b里没有噪声lambda设一个很小的值理论上最优解就是x_true本身。跑完之后比较恢复出的x和x_true看最大误差在什么量级。6.2 用KKT条件残差做定量检查比解误差更严格的是KKT条件残差。ADMM的最优点必须满足0属于A(Ax - b) lambda * sign(x)。这个包含集关系不好直接判断但可以计算每个分量上的次梯度大小来近似。对于非零分量残差应该接近0对于零分量次梯度的绝对值不应该超过lambda。这种检查能暴露主残差达标但解并不正确的问题。6.3 与现成求解器交叉对比如果环境里有优化工具箱可以直接用quadprog或者CVX验证ADMM的解。LASSO可以用CVX一行代码求解把两者的目标函数值对比差值在1e-4以内就算合格。方法设计基础验证脚本%% 基准测试已知稀疏解 rng(2024); % 固定随机种子保证可复现 m 100; n 300; A randn(m, n); % 高斯测量矩阵 x_true zeros(n, 1); support randperm(n, 10); % 随机选择10个支撑位置 x_true(support) randn(10, 1); b A * x_true; % 无噪声观测 lambda 0.01; rho 0.5; % 实际应用中需要与A的尺度匹配 [x_est, history] admm_lasso(A, b, lambda, rho, 1.0, 2000, 1e-6); recovery_err norm(x_est - x_true, 2) / norm(x_true, 2); fprintf(恢复误差: %.4e\n, recovery_err);这种验证脚本的价值在于一旦改动rho、终止条件或预处理方式可以马上看出真实影响而不会被真实数据的噪声干扰判断。从那以后我拿到任何一份优化算法代码第一步永远是先造已知解跑通反向验证再换真实数据已经习惯到改不掉了。希望帮到你。本文还有配套的精品资源点击获取