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 = c

f(x) 一般是光滑损失项，比如平方误差 0.5||Ax - b||^2；g(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 求解 LASSO：min 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('(A''A + 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^{k+1} - 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-5，eps_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('rho=0.01', 'rho=0.1', 'rho=1', 'rho=10'); xlabel('迭代次数'); ylabel('主残差');

逻辑说明：rng(42) 固定随机种子，保证重跑结果一致。sprandn 生成稀疏真解，让待估计的 x 有可识别结构。semilogy 用对数纵轴，方便观察主残差从 1e0 掉到 1e-6 的全过程。参数说明：lambda=0.1 是根据噪声水平选的，如果数据信噪比明显不同，lambda 要同步调整，否则稀疏解失真。实验里常见现象是 rho=0.01 时主残差掉得快但后期对偶残差震荡，rho=10 时前 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 可以同时开，一个管乘子更新的幅度，一个管变量更新的混合比例。建议先固定 alpha=1 调好 rho，再逐步引入 alpha，避免两个机制互相干扰。

3.3 主残差与对偶残差的联动：参数调试三板斧

我把这套调参流程整理成三板斧，拿到代码不知道怎么下手的可以直接照做。

第一步，固定 rho 扫数量级。用 3.1 的脚本，观察主残差曲线是否在 200 轮内进入线性下降。下降速度先快后慢，rho 偏小；从头到尾都慢，rho 偏大。

第二步，画主残差与对偶残差的比值曲线。比值长期大于 10，说明原始可行性滞后，rho 要加；比值长期小于 0.1，说明对偶可行性滞后，rho 要减。这个比值就是自适应判据的直接依据。

第三步，开自适应或者定死一个合适值。问题只跑一次就定死 rho；问题要反复求解或输入数据变化大，上自适应。

下面这张表是我在几个模拟项目里积累的经验值，适用条件是 A 为标准正态采样：

表里的 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(n+1: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 量级，还在用 rho=1。二次惩罚项在增广拉格朗日里权重过小，约束残差几乎不受拉扯。

解决：先打印 norm(A, 2) * norm(x_true, 2) 作为参考量，把 rho 提到和该量级同阶，比如 rho=100 或 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)。lambda=0 建议走另一条分支退化成 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 两个大数相减，精度崩塌。另一种来源是输入数据本身含 NaN，MATLAB 的 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_threshold：

function 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 的秩在迭代中逐步降到真实秩 3，S 的支撑集逐步收敛到构造位置，重构误差下降到 1e-6 左右，说明两个子问题的耦合被乘子轮次消化掉了。进阶到这一步，你对 ADMM 的把握已经从“会跑一个 LASSO”升级成“能设计算子、拆耦合、上并行”。每次新问题到手，我习惯先用一周目的最小算例跑通，再谈调参和加速；从那以后我再没在 ADMM 上翻过车。希望帮到你。

本文还有配套的精品资源，点击获取