简介MATLAB 实现的 ADMM 算法资源包面向机器学习、图像处理和信号处理等领域中求解大规模线性和非线性优化问题的开发者与研究人员。资源将复杂优化问题拆解为两个子问题通过松弛步骤、拉格朗日乘子更新和对偶变量更新完成交替迭代适合用于稀疏信号恢复、图像去噪、网络优化等场景。压缩包共 383 个文件约 17.44MB主要包含 .m 源码、封装好的 .p 文件、可直接调用的 .mex 编译接口、.mat/.dat 数据文件以及 pdf、txt 等说明文档兼顾算法学习、二次开发与工程调用需求。已有 1575 人学习下载目录中附有 fap、abalone 等多类数据集可配合算法框架快速复现实验、验证收敛效果并理解参数调优细节。1. 为什么你的 ADMM 在 MATLAB 里总发散先从稀疏恢复场景切入如果你的 ADMM 在 MATLAB 里一跑就出 NaN或者迭代几百次约束残差还不降别急着怀疑算法本身先想想更新顺序、符号方向和正则参数是不是哪里写歪了。matlab实现的 ADMM 算法本质上是围绕增广拉格朗日函数做三轮交替更新x 步、z 步和乘子步。它真正好用的一点是不用为每个新问题重写求解器只要把 f 和 g 的 prox 算子接进来大多数带线性约束的凸优化都能套同一个循环体。这篇文章面向已经知道 ADMM 是什么、但想立刻跑通的人从最小框架开始用 LASSO 做验证再把 rho、停止准则、复数 prox 这些坑一次性讲透你完全可以照着替换成自己的 f 和 g。2. ADMM 的两步分解从目标拆解到最小可运行循环2.1 目标拆解把 f、g 和约束 x - z 0 对上号ADMM 处理的是可分离目标写成通用形式是 min f(x) g(z)s.t. Mx Nz c。但实际写代码时绝大多数落地场景会通过变量代换把约束化简成 x - z 0 的一致性形式。这个形式够覆盖 LASSO、岭回归、全变分去噪、低秩矩阵恢复里最常见的那几种映射。如果你手里的约束是 Mx z那么推导里只需要把单位阵替换成 MMz 步的输入从 x u 换成 Mx u其余结构不变。为什么要把约束特意写成 x z因为 ADMM 的价值在于把原问题拆成两个“单变量子问题”。f 里如果有一个二次数据保真项x 步解出来是一个线性系统g 里如果带 l1 或全变分z 步解出来是一个闭式 prox。两个子问题都比原问题好解这是 ADMM 能落地的根本原因。如果 f 和 g 都复杂到没有可用的 prox那 ADMM 只是把难题换了个位置实际跑起来一样痛苦。三个更新式写下来是这样x 步x^{k1} argmin_x f(x) (rho/2) || x - z^k u^k ||^2z 步z^{k1} argmin_z g(z) (rho/2) || x^{k1} - z u^k ||^2u 步u^{k1} u^k x^{k1} - z^{k1}这里 u 是缩放对偶变量不是原始对偶变量。使用缩放变量之后u 的更新就是简单的“把残差加回去”x 和 z 步里的二次项也不用额外展开。很多人实现翻车就是把缩放变量和原始乘子混着用导致符号对不上。后面避坑章节专门讲这个。2.2 最小可运行循环x、z、u 三个更新在 MATLAB 里长什么样以 LASSO 为例目标函数是 0.5 || Ax - b ||^2 lambda || z ||_1约束 x - z 0。这里 A 是测量矩阵b 是观测f 就是二次项g 就是 l1。x 步的闭式解来自正规方程z 步是软阈值。最小骨架如下% 最小可运行 ADMM 循环骨架 % 问题min 0.5||Ax-b||^2 lambda*||z||_1s.t. x z [m, n] size(A); x zeros(n, 1); z zeros(n, 1); u zeros(n, 1); K A*A rho*speye(n); % x 步法方程矩阵 R chol(K); % 对称正定预分解一次 Atb A*b; for k 1:maxIt zOld z; xNew R \ (R \ (Atb rho*(z - u))); % x 步解正规方程 zNew softThresh(xNew u, lambda/rho); % z 步软阈值 u u (xNew - zNew); % u 步乘子更新 r xNew - zNew; % primal 残差 s rho * (zNew - zOld); % dual 残差 x xNew; z zNew; end这个循环里最关键的是 K 只组装一次。AA rhoI 在迭代中不随 k 变化所以用 chol 做一次分解循环里只做两次三角回代。常见误用是每次迭代都重新算 inv(AA rhoI) 或重新做 mldivide这在几百维时感觉不明显上到几千维就会让你等得怀疑人生。x 步里 Atb rho*(z - u) 是对应右端项。z 步的 softThresh 输入是 xNew u阈值是 lambda/rho。注意不是 lambda因为 z 步的二次项系数是 rho/2prox 点本身已经包含了尺度折算。如果你把这两个子问题中的一个符号写反比如把 xNew u 写成 xNew - u迭代会在前几步看起来正常随后残差开始震荡最终发散。后面避坑章节会给出具体的排查方法。3. 用 LASSO 验证 ADMM软阈值更新与一段完整 MATLAB 脚本3.1 把 LASSO 装进 ADMMprox 算子到底做了什么LASSO 的形式是 min 0.5 ||Ax - b||^2 lambda ||x||_1。它天然就是 f g 的结构不需要额外拆分就能映射到 ADMM。f 的 prox 步等价于求解一个带二次惩罚的最小二乘x 步闭式解为 x (AA rho I)^{-1}(Ab rho(z - u))。g 的 prox 步是 l1 范数的近端算子也就是软阈值function y softThresh(z, t) % 软阈值算子z 可以是实数也可以是复数 % t 是正阈值对应 lambda/rho a abs(z); y max(a - t, 0) ./ max(a, eps) .* z; end软阈值做的事情很简单模长小于 t 的分量直接置零大于 t 的分量按比例收缩。分母里加 eps 是为了避免 a 0 时出现 0/0 的 NaN。这个实现对复数输入也成立因为乘的是 z 而不是 sign(z)相位被完整保留了。后面避坑章节会看到如果写成 max(a - t, 0) .* sign(z)复数信号会直接出问题。lambda/rho 作为一个整体出现在阈值位置这个细节值得多说一句。调 lambda 的时候如果 rho 不跟着变实际阈值会同步变化也就是说你调的不只是正则项强度还在变相调整约束惩罚。我一般习惯先定 rho再在固定的 rho 下去扫 lambda这样能少一个变量。3.2 合成数据实验从生成数据到调用求解器用一个小规模欠定问题验证实现是否合理。生成 200 行 500 列的随机测量矩阵、一个只有 15 个非零元素的稀疏真解再加一点噪声% 生成稀疏恢复数据 rng(20240601); m 200; n 500; A randn(m, n); A A * diag(1 ./ sqrt(sum(A.^2, 1))); % 列归一化 xTrue zeros(n, 1); idx randperm(n); idx idx(1:15); xTrue(idx) randn(15, 1) * 3; b A * xTrue 0.01 * randn(m, 1); % 加小幅噪声 lambda 0.1 * max(abs(A*b)); % 按相关性的量级粗取 rho 1;列归一化这一步不是必须的但能避免某些列范数过小导致尺度失衡。lambda 取 0.1 * max(abs(A*b)) 是一个常见粗估A*b 是观测与字典的相关性最大值反映稀疏信号可能的最大投影强度取十分之一量级作为正则强度通常能保住主要非零位置。如果噪声明显更大这个值要往上调如果信号更稀疏更干净可以往下调。求解器本体可以封装成一个独立函数这样换数据时不用动算法代码function [x, hist] admm_l1(A, b, lambda, rho, opts) if nargin 5, opts struct(); end absTol 1e-6; relTol 1e-4; maxIt 5000; if isfield(opts, absTol), absTol opts.absTol; end if isfield(opts, relTol), relTol opts.relTol; end if isfield(opts, maxIt), maxIt opts.maxIt; end [m, n] size(A); x zeros(n, 1); z zeros(n, 1); u zeros(n, 1); K A*A rho*speye(n); R chol(K); % 预分解一次 Atb A*b; hist.primal zeros(maxIt, 1); hist.dual zeros(maxIt, 1); for k 1:maxIt zOld z; x R \ (R \ (Atb rho*(z - u))); z softThresh(x u, lambda/rho); u u (x - z); r x - z; s rho * (z - zOld); hist.primal(k) norm(r); hist.dual(k) norm(s); epsP sqrt(n)*absTol relTol*max(norm(A*x), norm(z)); epsD sqrt(n)*absTol relTol*norm(rho*u); if hist.primal(k) epsP hist.dual(k) epsD hist.primal hist.primal(1:k); hist.dual hist.dual(1:k); break; end end end这个函数里hist 记录了每一轮的两个残差收敛后截断到实际迭代数。停止条件里 epsP 和 epsD 一个控制原始可行性一个控制对偶可行性。只盯着目标函数值收敛是没用的ADMM 必须两个残差都落到阈值以内才算真正收敛。调用和可视化如下[xHat, hist] admm_l1(A, b, lambda, rho, struct(maxIt, 3000)); figure(Position, [100 100 800 300]); subplot(1,2,1); stem(1:n, xTrue, k, MarkerSize, 3); title(真实信号); subplot(1,2,2); stem(1:n, xHat, r, MarkerSize, 3); title(ADMM 恢复); fprintf(相对误差: %.3e\n, norm(xHat - xTrue)/norm(xTrue));一次典型运行里非零位置基本能对上幅度会有些偏差这是噪声和 lambda 共同作用的结果。如果恢复出的非零位置完全对不上先别调 rho回去看数据生成里列归一化和噪声量级如果位置对了但幅度偏小通常是 lambda 偏大如果出现很多小的伪峰通常是 lambda 偏小。这个排查顺序我从不会反着来。4. rho 和停止准则怎么调参数表、残差曲线与判断方法4.1 三个关键参数rho、容差、maxIter 的职责与调法ADMM 落地时真正需要调的参数不多但每个都不能乱来。rho 控制 x 步里约束项的分量直接决定收敛速度和路径absTol 和 relTol 决定你要求多精确maxIter 只是保险丝不是精度旋钮。四者的关系总结如下参数控制什么调大的典型行为调小的典型行为常用起点rhox 步里约束惩罚的权重primal 残差下降慢x 步钝化dual 残差下降慢可能出现震荡1absTol残差的绝对尺度更早停止精度变低更晚停止容易过拟合噪声1e-6relTol残差相对当前解的尺度相同情况下更宽松更严格迭代变多1e-4maxIter最大迭代轮数防止死循环几乎不影响结果只影响等待时间太小时会在收敛前被截断5000rho 是这里面唯一带玄学气质的参数。rho 太小约束项在 x 步里话语权不足x 和 z 各走各的对偶残差会拉得很高rho 太大x 步直接被约束项主导等于每一步都先把 x 拉到 z 身边再让 f 慢慢修正原问题信息被稀释primal 残差会像钝刀子割肉一样慢慢磨。比较好的做法是先在 0.1、1、10 三个量级各跑一次看收敛曲线选一个再细调。4.2 用残差历史判断收敛方向primal 和 dual 两条线要分开看先跑一次并画出两条残差曲线[xHat, hist] admm_l1(A, b, lambda, rho, struct(maxIt, 3000)); figure(Position, [100 100 800 300]); subplot(1,2,1); loglog(hist.primal, LineWidth, 1.5); hold on; loglog(hist.dual, LineWidth, 1.5); legend(primal, dual, Location, northeast); xlabel(迭代次数); ylabel(残差范数); title(收敛曲线); subplot(1,2,2); semilogy(1:numel(hist.primal), ... 0.5*norm(A*xHat - b)^2 lambda*norm(xHat, 1), k, LineWidth, 1.5); xlabel(迭代次数); ylabel(目标函数值); title(目标值曲线);读图规则很简单primal 降得慢、dual 早就达标是 rho 偏大dual 一直在高位、primal 降到平台是 rho 偏小两条线都快速下降然后同时穿过阈值才是比较健康的曲线。如果目标函数值看起来很平稳但两条残差一直没达标说明算法还在打磨约束可行性这时候停下来得到的结果可能不满足约束直接拿去用会出问题。更系统的做法是扫几个 rho 值看迭代数变化rhoList [0.3, 1, 3]; for j 1:3 [~, histTmp] admm_l1(A, b, lambda, rhoList(j), struct(maxIt, 3000)); fprintf(rho %.1f迭代 %d 次\n, rhoList(j), numel(histTmp.primal)); end同一问题在不同 rho 下的迭代次数可以差好几倍。某个 rho 在这个问题上表现好换一个测量矩阵或换一个噪声水平可能就一般了。所以扫 rho 不能只跑一次就定死最少在三个量级上看趋势。我自己的习惯是写一个小循环把 rho 从 0.3 到 3 各跑一遍选出迭代数最少且残差最后都达标的那个再去做下一步。5. ADMM 避坑指南五个让收敛曲线翻车的实现细节5.1 rho 取值过大primal 残差停在半路现象目标函数值在正常下降但 primal 残差迭代到一两百步后停在一个固定量级比如 1e-3再也不动了dual 残差早就达标。原因rho 被设成 100 甚至更大的值x 步里 rho*I 成了主导项相当于每一次都在把 x 拉到 z 附近f 的信息被严重稀释。x 和 z 互相跟得很紧但离真正的最优解还很远于是 primal 残差被卡住。解决把 rho 往 1 附近调同时保持 lambda/rho 不变。如果目标里 f 是二次函数一个简单经验是看 AA 的特征值量级rho 取接近 AA 最大特征值的十分之一到十分之三通常不会太离谱。改成这个量级后primal 残差会重新开始下降。5.2 停止准则只看绝对容差迭代停不下来现象目标函数已经几乎不变残差曲线也很平但程序一直跑满 maxIter没有任何提前 break。原因absTol 设成 1e-12而这个问题本身的残差量级在 1e-6 就够用了。把精度要求设到比问题尺度还高一个量级以上迭代自然停不下来。relTol 如果设成 0还会让这个情况更严重因为相对尺度完全失效。解决用 absTol 1e-6、relTol 1e-4 起步再看残差曲线的最终水平。如果两条线都落在阈值区间内但你想更精确把 relTol 降一个量级试试如果只是数据噪声本身就有 1e-3 的水平那 absTol 设 1e-6 就是给自己找罪受。先看数据量级再定容差不要反过来。5.3 每次循环都重构矩阵耗时成倍上涨现象迭代轮数不多但每轮都在做矩阵分解500 维的问题还能忍2000 维以上明显变慢运行时间从几十秒涨到几分钟。原因K AA rho*I 在整个迭代过程中不变但代码里把矩阵组装和求解都写进了循环体等于每轮都重新做一次复杂度 O(n^3) 的分解。解决把组装和预分解提到循环外循环内只做三角回代。如果 A 是稀疏的用 ichol 替代 chol配合回代一样能跑if issparse(A) K A*A rho*speye(n); R ichol(K, struct(type, ict, droptol, 1e-3)); else K A*A rho*speye(n); R chol(K); end % 循环内统一用 x R \ (R \ q);这个改动几乎不需要额外思考收益却非常直接尤其在矩阵规模上千之后。我那会儿在一个模拟项目里就是吃了这个亏改完预分解后同样的迭代数跑了不到原来的十分之一时间。5.4 复数信号直接把软阈值写错相位被破坏现象输入是复数观测恢复结果模长看起来对但相位明显不对误差大得离谱。原因软阈值实现成了 max(|z| - t, 0) .* sign(z)。对实数信号这没问题但对复数sign(z) 会把相位直接抹掉或改掉等于把复数当实数在收缩。l1 范数在复数域里定义的是模长prox 应该只压缩模长、保留相位。解决用前面写的 softThresh 实现分母乘 z 而不是乘 sign(z)。这个写法对实数和复数都成立一步到位function y softThresh(z, t) a abs(z); y max(a - t, 0) ./ max(a, eps) .* z; end如果你在做傅里叶变换相关的恢复问题这一步最容易翻车因为中间变量一不留神就变成复数。写完之后可以单独跑一个随机复数向量验证一下输入 z 和输出 y 的 angle(z) 应该完全一致只是模长被收缩了。5.5 对偶变量方向写反残差曲线直接上天现象前几步残差还在下降中间突然开始指数增长然后出现 NaN或者 u 的范数一路飙升不回头。原因把 u 的更新方向写反了。比如 x 步里按 x Atb rho*(z - u) 算但 u 更新却写成 u u - (x - z)两个符号一凑相当于在残差上做负反馈变成正反馈越迭代越放大。解决始终以三个更新式为准。选定 x argmin f (rho/2)||x - z u||^2 之后u 必须更新为 u (x - z)。换一个写法之前先检查符号一致性。一个低成本自检是把 x、z、u 都初始化为零跑一步看 u 的第一个元素是否等于 -x 的对应值再加上 z。如果 u 的第一轮就出现意外的大变化多半是方向写反了。6. 更快的迭代习惯预分解回代与过松弛的取舍验证6.1 预分解循环前的一次分解胜过迭代内一百次回代前面避坑里提到预分解能解决慢的问题但它不只是“减少重复计算”还让 x 步变成两次三角回代。对于 500 维的稠密问题chol 分解一次和每轮回代两次的成本差基本在一个数量级以上对于稀疏大矩阵用 ichol 还能配合稀疏存储进一步压内存。我一般会把这一步作为所有 ADMM 实现的默认选择除非问题规模小到 100 维以下才懒得做。6.2 过松弛一个参数让收敛曲线变得更陡如果残差曲线下降得很稳但偏慢可以试试过松弛。过松弛本质是在 z 步的输入里混合新 x 和旧 z让它更激进地预测下一步位置alphaRelax 1.5; % 常见取值 1.5超过 2 容易发散 zOld z; x R \ (R \ (Atb rho*(z - u))); xTilde alphaRelax * x (1 - alphaRelax) * zOld; z softThresh(xTilde u, lambda/rho); u u (xTilde - z);alpha 取 1 就是标准 ADMM取 1.5 通常能明显减少迭代数但超过 2 之后经常会看到残差开始震荡最后发散。过松弛不能解决 rho 选错的问题它是在 rho 和容差都已经合理的前提下锦上添花。验证方法是把加过松弛前后的收敛曲线叠在一起看迭代数变化不要只凭感觉判断。我也习惯在最后做一步 KKT 抽查确认结果真的站在最优解附近gradF A * (A*xHat - b); signX sign(xHat); kktRes gradF - lambda * signX; fprintf(非零位置 KKT 残差: %.3e\n, norm(kktRes(xHat ~ 0)));这个值如果和提前终止的残差在同一个量级说明停止准则没有骗人。我早年在这个方向上交过学费某个模拟项目里把 rho 设到 100残差停在 1e-2我连着两天在怪数据预处理后来画了残差曲线才发现是 rho 压死了 x 步。之后我给自己定了一条规矩改任何参数之前先看 primal/dual 残差的 log-log 图。这个习惯省下的调试时间比多学几个加速技巧值钱得多。希望帮到你。本文还有配套的精品资源点击获取