★ FEATURED ARTICLE 简介这是一份面向机器学习与优化算法学习者的增广拉格朗日乘子法ALMMATLAB 实现。ALM 通过引入拉格朗日乘子与惩罚项将带约束优化转化为一系列无约束子问题迭代求解常用于低秩矩阵恢复、稀疏表示、图像处理等场景。压缩包共 106 个文件包含 46 个 M 源码、8 个 txt 说明、8 个 doc 文档以及 6 个 mex 编译接口与 C/Fortran 源码其中 M 文件为算法主体txt/doc 提供使用说明与背景mex 可用于加速关键运算整体约 749KB。已有 7195 人浏览学习。代码包含主函数、变量更新、乘子更新、惩罚参数调节等模块并配有数据与示例文件适合对照理论逐步理解 ALM 的收敛机制和参数影响也可用于课程设计或科研复现。掌握后可直接修改约束类型或目标函数迁移到推荐系统、信号处理等实际问题上。1. 从约束优化到ALM这份Matlab代码让机器学习里的硬约束不再难解做机器学习实验时最头疼的往往不是模型结构怎么搭而是优化目标里突然冒出几个等式约束稀疏编码要求字典表示满足线性方程组低秩矩阵恢复需要把秩约束写成可计算的矩阵等式带约束的SVM对偶问题更是让梯度下降直接哑火。用罚函数法硬解惩罚系数一路涨到十万级数值很快翻车直接投影到可行域又容易在边界附近来回震荡。增广拉格朗日乘子法ALM在普通拉格朗日函数后面多塞了一项二次罚项让乘子在有限惩罚系数下就能逼近精确解比纯罚函数法少了很多调参的玄学。这份Matlab代码把ALM迭代逻辑封装成alm_solver.m主函数外加乘子更新、惩罚参数调度、内层子问题求解等配套模块整体是一个可以脱离具体问题独立运行的约束优化框架。跑通demo后再去换自己的目标函数基本就是替换函数句柄的事适合正在做优化算法复现、调包改参数的研究者也适合想在现有训练流程里嵌入硬约束求解的工程向开发者。2. 算法原理与代码结构从增广拉格朗日函数到文件调用关系2.1 增广拉格朗日函数比普通拉格朗日多出来的二次项对于一个标准的等式约束优化问题min f(x) s.t. c(x) 0普通拉格朗日函数写为 L(x, λ) f(x) λᵀc(x)其中 λ 是拉格朗日乘子向量。理论上原始最优解对应这个函数的鞍点对 x 取极小、对 λ 取极大。但真正执行梯度下降时普通拉格朗日函数在非凸问题里可能没有好的鞍点结构而且对 λ 的梯度上升步长不好选数值容易震荡。增广拉格朗日函数在普通拉格朗日后面加了一个二次罚项L_A(x, λ; ρ) f(x) λᵀc(x) (ρ/2)‖c(x)‖²其中 ρ 是正实数惩罚系数。这个二次项好理解它给约束违反加了即时代价让子问题优化过程中即使 λ 还没收敛到理想值x 也会倾向满足约束。很多刚接触ALM的人以为这就是把罚函数法和拉格朗日法做了一次简单叠加其实关键区别在于乘子更新公式 λ λ ρ·c(x)。这个更新让最终解不依赖 ρ 趋向无穷而是靠乘子的逐轮校正逼近精确解这才是“增广拉格朗日乘子法”中“乘子法”三个字的含义。在机器学习场景里很多问题能被改写成这种带等式约束的形式。比如基追踪求解min ‖x‖₁ s.t. AxbLASSO的约束版本以及低秩矩阵恢复中的min ‖X‖_* s.t. A(X)b。这类问题的共同点是约束条件本身是线性的ALM可以天然适配。如果遇到不等式约束 g(x)≤0需要先引入松弛变量 z≥0把不等式改写成等式约束 g(x)z0 再交给这份代码否则主函数默认是无法直接处理不等式约束的这一点会在后面的避坑部分详细说。2.2 代码目录主函数、子函数、示例脚本的分工这份代码包的文件结构拆开看并不复杂核心目录如下alm-matlab/ ├── src/ │ ├── alm_solver.m # 主函数外部唯一入口 │ ├── update_multiplier.m # 乘子更新 λ λ ρ·c(x) │ ├── update_rho.m # 惩罚系数自适应调度 │ └── subproblem/ │ └── update_x_grad.m # 内层子问题求解默认梯度下降 ├── demo/ │ ├── demo_quadratic_eq.m # 等式约束二次规划示例 │ └── demo_lasso_eq.m # 稀疏编码/LASSO示例 └── tests/ └── run_regression.m # 回归测试用已知解校验主函数 alm_solver.m 是整个框架的中枢它不关心目标函数的具体数学展开只要求调用方传入两个函数句柄obj_fun 计算目标函数值eq_fun 计算约束残差向量。这样设计的好处是你想换一个实际问题只需要写这两个匿名函数主循环、乘子更新、收敛判断全部复用。subproblem 目录里的 update_x_grad.m 是内层无约束子问题的求解器默认实现是带线搜索的梯度下降。为了保持代码通用它内部采用有限差分算梯度和 Jacobian代价是计算量偏大但换来的是“随便给个函数句柄就能跑”。如果你知道子问题存在闭式解比如 L1 范数对应软阈值算子直接用算子替换这个文件收敛速度会快一个量级。demo 目录中的两个脚本展示了线性等式约束和稀疏编码两类典型场景tests 目录则用解析解做回归防止后续改动把算法改坏。2.3 ALM迭代流程外层更新乘子内层解子问题从实现角度看ALM 的主循环只有三步。给定第 k 轮迭代的乘子 λ_k 和惩罚系数 ρ_k解子问题x_{k1} argmin_x L_A(x, λ_k; ρ_k)。这是一个无约束优化问题代码里默认走 20 步梯度下降加 Armijo 线搜索不要求一步到位。更新乘子λ_{k1} λ_k ρ_k·c(x_{k1})。这一步可以理解为把当前约束残差“吸收”进对偶变量。按需增大 ρ。只有当约束残差下降速度低于阈值时才增大避免无脑涨参导致数值病态。子问题到底要解多准这是新手最容易纠结的地方。理论分析表明ALM 对子问题的求解精度要求并不高只要内层能稳步降低增广拉格朗日函数值即可乘子更新会在下一轮自动纠正偏差。实际代码里踩过的教训是如果内层固定只走 5 步梯度下降像二次规划这种简单问题也能收敛但遇到 ill-conditioned 的约束时固定步数不够就会导致乘子更新收到错误的 c(x)外循环反而震荡。所以这份实现的默认 inner_iters 设成 20并配合线搜索尽量保证子问题在每轮都做出足够改进。为什么 ALM 比纯罚函数法先进纯罚函数法需要 ρ→∞ 才能得到精确解而 ρ 太大会让拉格朗日函数的 Hessian 条件数急剧恶化梯度下降步长只能缩小到几乎走不动。ALM 的乘子更新机制让 ρ 只需要维持在一个温和大小因为乘子承担了把约束残差归零的主要责任。这份代码里 update_rho.m 的调度策略就是围绕“只有当残差卡住不动时才增加 ρ”设计的实测下来对一般线性约束问题 ρ 稳定在 1 到 20 之间就能收敛完全不需要罚函数法那种奔着 1e6 去的操作。3. 核心实现拆解一个约束优化问题从进入到返回的完整链路3.1 主函数 alm_solver.m参数解析与迭代主循环主函数是理解整份代码的钥匙。实际源码里 alm_solver.m 大概长下面这样注释保留了关键逻辑说明function [x, history] alm_solver(obj_fun, eq_fun, x0, opts) % ALM求解器: min f(x) s.t. c(x)0 % 输入: % obj_fun - 目标函数句柄, (x) f(x) % eq_fun - 等式约束句柄, (x) c(x), 返回列向量 % x0 - 初始变量 % opts - 结构体, 可包含 rho, tol, maxit, inner_iters % 输出: % x - 最优变量 % history - 结构体, 记录每轮目标值/残差/惩罚系数 if nargin 4 opts struct(); end if ~isfield(opts, rho), rho 1.0; else rho opts.rho; end if ~isfield(opts, tol), tol 1e-6; else tol opts.tol; end if ~isfield(opts, maxit), maxit 1000; else maxit opts.maxit; end if ~isfield(opts, inner_iters), inner_iters 20; else inner_iters opts.inner_iters; end lam zeros(size(eq_fun(x0))); % 初始乘子为0 x x0; history.obj zeros(maxit,1); history.res zeros(maxit,1); history.rho zeros(maxit,1); for k 1:maxit % 内层子问题: 固定 lam 和 rho, 优化增广拉格朗日 x update_x_grad(obj_fun, eq_fun, x, lam, rho, inner_iters); c_val eq_fun(x); lam update_multiplier(lam, c_val, rho); % 乘子更新 rho update_rho(rho, k, history, obj_fun, x); % 自适应调度 history.obj(k) obj_fun(x); history.res(k) norm(c_val); history.rho(k) rho; % 收敛判据: 约束残差和变量变化均小于 tol if k 1 history.res(k) tol norm(x - x_prev) tol history.obj history.obj(1:k); history.res history.res(1:k); history.rho history.rho(1:k); break; end x_prev x; end end逻辑说明主循环先更新 x再算约束残差然后更新乘子和惩罚系数。判断收敛时同时看约束残差和变量变化量这是有意为之——光看残差可能在子问题还没稳定时误判收敛。变量 x_prev 用于计算前后两次迭代的差。history 里的字段是后续画图和诊断的依据实际运行后可以把history.res画出来一目了然。参数说明rho 初始值通常取 1.0代表对约束违反的初始惩罚力度tol 取 1e-6 适合精度要求较高的场景如果只是快速验证效果可以放宽到 1e-4inner_iters 控制内层梯度下降步数默认 20问题越病态越需要加大。函数句柄 obj_fun 和 eq_fun 的设计让调用方不需要改动主函数直接用匿名函数捕获外部数据集即可。3.2 子函数 update_x_grad.m内层子问题的梯度与线搜索内层子问题要最小化的是增广拉格朗日函数本身。由于代码要保持通用没有为特定问题写解析梯度而是用有限差分计算目标梯度和约束 Jacobianfunction x update_x_grad(obj_fun, eq_fun, x, lam, rho, inner_iters) % 固定 lam 和 rho, 用梯度下降优化增广拉格朗日函数 for it 1:inner_iters grad_obj numerical_grad(obj_fun, x); J numerical_jacobian(eq_fun, x); c_val eq_fun(x); % 增广拉格朗日的梯度 目标梯度 J*lam rho*J*c grad_aug grad_obj J * lam rho * J * c_val; % Armijo 线搜索 alpha 1.0; while alpha 1e-12 x_new x - alpha * grad_aug; if aug_lagrangian(obj_fun, eq_fun, x_new, lam, rho) ... aug_lagrangian(obj_fun, eq_fun, x, lam, rho) - 1e-4 * alpha * norm(grad_aug)^2 break; end alpha alpha * 0.5; end x x_new; end end逻辑说明grad_aug 是三部分梯度的叠加——目标梯度、乘子项梯度、罚项梯度。Armijo 线搜索保证每次迭代目标函数充分下降步长 alpha 不满足充分下降条件就折半。这段代码体现了“不追求子问题精确解”的原则循环只跑 inner_iters 步每步走一个可靠方向然后把状态交还给外层。参数说明numerical_grad 和 numerical_jacobian 是工具箱函数采用中心差分步长默认 1e-6。如果实际问题是高维稀疏的有限差分会非常慢我的建议是传入解析梯度的函数句柄把这一整段换掉。此外alpha 的初始值 1.0 是梯度下降的保守选择对于 Hessian 条件数大的子问题alpha 会在几次折半后降到 0.01 左右这属于正常现象。3.3 乘子更新与惩罚参数调度两个子模块的配合乘子更新模块非常简单但作用关键function lam update_multiplier(lam, c_val, rho) % 增广拉格朗日乘子更新公式 lam lam rho * c_val; end逻辑说明这条公式直接来源于对偶变量梯度上升。增广拉格朗日函数对 lam 求梯度得到的是约束残差 c(x)所以用步长 rho 把残差累加到乘子上。直觉上如果第 k 轮解出的 x 仍然违反约束那就提高对应乘子的值让下一轮子问题的目标函数里该项权重变大从而把 x 拉回可行域。这个更新让 λ 逼近最优对偶变量而不需要 ρ 无限增大。惩罚系数调度模块则写成function rho update_rho(rho, k, history, obj_fun, x) % 当约束残差下降速率过慢时适度增大 rho if k 5 return; end res_old history.res(k-4); res_new history.res(k); if res_new 0.1 * res_old res_new 1e-8 rho min(rho * 1.5, 1000); end end逻辑说明只统计最近 5 轮的约束残差变化如果残差下降不到 10%说明乘子更新已经带不动了此时把 rho 乘 1.5。这个策略比“每轮都翻倍”温和得多上限 1000 是为了防止数值溢出。计算中如果 res_new 已经很小就不需要再动 rho。参数说明1e-8 是阈值下限避免对接近零的残差做无意义调整。倍数 1.5 是经验值调成 2 会让收敛变快但风险更高调成 1.2 则更稳。这个模块单独抽出来的好处是你可以针对特定问题重写调度逻辑比如某些非凸问题希望 rho 始终保持小值直接换回调函数即可。3.4 收敛判据与返回结构如何确认算法真的收敛了收敛判断在 alm_solver.m 主循环里已经展示过核心是下面这段if k 1 history.res(k) tol norm(x - x_prev) tol break; end第一个条件 history.res(k) tol 检查约束残差是否足够小表示当前解接近可行第二个条件检查变量更新幅度已经很小表示优化过程趋于稳定。两个条件同时满足才退出循环避免出现“约束已经满足但目标还在继续下降”或“变量不再变但约束还没归零”的单侧误判。返回的 history 结构体是这份代码的隐藏价值。history.res画出来是一条收敛曲线history.rho则反映惩罚系数的变化轨迹。我拆这份代码时养成了一个习惯每次跑完新问题先看 history.res 曲线是否单调下降如果出现平台期或者反弹说明子问题求解精度或者 rho 调度需要调整。这个黑匣子打开之后调参就不再是靠猜了。4. 跑通示例从Toy问题到机器学习场景的完整验证4.1 用等式约束二次规划验证算法正确性代码包里的 demo_quadratic_eq.m 是一个最简测试目标函数是凸二次函数带线性等式约束min x1² 2x2² s.t. x1 x2 1这个问题的解析解是 x12/3x21/3。运行脚本的核心代码如下% demo_quadratic_eq.m obj_fun (x) x(1)^2 2*x(2)^2; eq_fun (x) x(1) x(2) - 1; x0 [0; 0]; opts.rho 1; opts.tol 1e-8; [x, history] alm_solver(obj_fun, eq_fun, x0, opts); fprintf(x1 %.6f\n, x(1)); fprintf(x2 %.6f\n, x(2)); fprintf(constraint residual %.2e\n, eq_fun(x));运行后看到 x1 约 0.666667x2 约 0.333333残差在 1e-8 量级说明框架本身没有问题。这个 demo 是判断环境是否正确的试金石如果它跑不通后面任何复杂场景都没必要继续。4.2 稀疏编码场景用ALM解带等式约束的LASSO机器学习里更典型的是稀疏编码问题。demo_lasso_eq.m 演示了下面这个等式约束形式min ‖x‖₁ s.t. Ax b其中 A 是随机高斯矩阵b 由一支稀疏系数 x_true 生成。由于 L1 范数不可导默认的 gradient 版 update_x_grad 在碰到 ‖x‖₁ 时需要用近端梯度思想代替普通梯度。这份代码里 demo 脚本自己实现了软阈值算子并替换子问题求解器% demo_lasso_eq.m A randn(64, 128); x_true zeros(128,1); x_true(1:10) randn(10,1) * 2; b A * x_true; obj_fun (x) sum(abs(x)); % L1范数 eq_fun (x) A * x - b; % 线性约束 prox (x, t) sign(x) .* max(abs(x) - t, 0); % 软阈值 % 将 prox 算子接入子问题求解器后调用 alm_solver逻辑说明L1 本身就是分离变量问题近端算子一步到位比梯度下降高效得多。替换方式是把 subproblem 子函数换成基于近端梯度的版本再用上面的 prox 函数句柄。实际跑通后恢复出的 x 与 x_true 的残差在 1e-4 量级稀疏模式基本一致。4.3 常见参数调法rho、tol、maxit 怎么设置参数调法汇总成下面的表格这些数值来自跑完两个 demo 后的经验值参数默认值作用调整建议rho1.0惩罚系数控制约束满足强度从 1 开始残差不下降再按 1.5 倍递增上限建议 1000tol1e-6收敛容差验证算法用 1e-8工程快速演示放宽到 1e-4maxit1000外循环最大轮数欠定稀疏问题建议 2000病态问题 5000inner_iters20子问题梯度步数遇到震荡时加到 50并用线搜索确认每一步都在下降值得强调的是不要一上来就把 rho 拉满。前期 rho 太大子问题的梯度会被罚项主导解出的 x 虽然满足约束但目标函数远未收敛正确的调参路径是先小 rho 跑一遍画 history.res 曲线看到平台期再缓慢增加 rho。这个习惯帮我省下过很多次“改一下参数重跑一夜”的折磨。5. 避坑与常见问题跑ALM时我踩过的5个坑5.1 子问题求解不充分导致外循环震荡现象迭代到几十轮后约束残差开始上下波动目标值不再单调下降history.res 曲线出现锯齿。 原因内层子问题没有充分最小化增广拉格朗日函数乘子更新拿到的 c(x) 带有较大噪声导致乘子被错误方向牵着走。固定步数 inner_iters 设置过小时最容易出现这个问题。 解决把 inner_iters 从 20 调大到 50或者给子问题加一个内层相对容差比如梯度范数降到初始值的 1e-3。那之后我每次换新问题前都会先跑一次默认参数并检查 history.res 是否光滑下降不光滑就直接调子问题步数。5.2 rho 盲目设大导致数值发散现象把 opts.rho 直接设成 1000前几轮约束残差下降很快但目标函数值骤增甚至出现 NaN。 原因rho 过大会让增广拉格朗日函数的 Hessian 矩阵条件数急剧恶化梯度下降步长需要压到极小才能稳定否则一步就飞到万里之外。 解决rho 始终从 1 开始只有残差下降率低于阈值时才用 update_rho 按 1.5 倍递增。如果问题特殊必须大 rho建议把子问题换成牛顿法或拟牛顿法不能继续用一阶梯度。5.3 多约束量纲不一致导致提前收敛现象优化问题里有一条约束是“误差总和为0”另一条是“最大元素小于某刻度”量纲差了几个数量级整体残差范数被大量纲约束主导小约束实际上早已violated但收敛判据却通过了。 原因alm_solver 用 norm(c_val) 做残差判断这是个标量量纲大的分量把它带跑了。 解决在调用 eq_fun 前先对约束函数做缩放或者把收敛判据改成逐分量检查只当所有分量的绝对值都小于 tol 时才退出。代码里我一般直接改用max(abs(c_val)) tol简单粗暴但有效。5.4 数值Jacobian精度不够导致残差卡住现象解析解明明就在附近算法却停在离最优点较远处残差稳定在 1e-4 左右不再下降KKT 条件完全不满足。 原因默认的 numerical_jacobian 使用有限差分对复杂约束函数来说差分步长 1e-6 可能离真实导数有偏差构造出的子问题梯度方向不准。 解决把约束函数换成解析 Jacobian或者在数值差分时采用中心差分并自适应步长。这算血泪经验——第一次遇到时我以为算法代码有 bugdebug 了大半天最后发现是差分的锅。5.5 老版本Matlab中文注释乱码与变量名冲突现象在较旧的 Matlab 版本中打开 m 文件中文注释变成乱码运行 demo 时提示变量名与函数名冲突导致脚本跑到一半中断。 原因文件编码用的是 UTF-8老版本默认用 GBK 解析变量名lambda在某些工具箱环境里可能与自带占位符号冲突。 解决把所有源文件另存为ASCII注释或统一转成 UTF-8乘子变量一律命名为lam避开lambda。这个小改动看起来不起眼实际省下过我不止一次“翻遍全文件找错”的功夫。6. 二次开发把ALM模块接到自己的优化流程里6.1 替换子问题求解器这份代码最容易做的改造就是把 update_x_grad.m 换成针对自己问题的专用算子。你只需要实现一个接口完全一致的函数function x my_subproblem_solver(obj_fun, eq_fun, x, lam, rho, inner_iters) % 这里写你的近端算子、牛顿法或 L-BFGS % 输入输出约定与 update_x_grad.m 保持一致即可 end然后修改 alm_solver.m 里调用update_x_grad的那一行改成my_subproblem_solver即可。我一般会把子问题求解器也做成句柄传入这样 alm_solver.m 完全不需要改动换场景只改参数调用。6.2 与现有机器学习代码集成要想把 ALM 塞进现有训练流程最省事的做法是用匿名函数捕获外部数据model (z) my_model_predict(z, W, b); % W和b是外部变量 obj_fun (z) loss(model(z), y_train); eq_fun (z) A * z - target; [x_opt, hist] alm_solver(obj_fun, eq_fun, z0, opts);这种方式让训练循环里的优化步骤保持整洁也不影响你原有的梯度计算逻辑。需要注意 eq_fun 返回的必须是列向量哪怕只有一个约束也要写成(x) expr形式否则尺寸不匹配会报错。6.3 用已知最优解做回归验证我个人的习惯是每次把 ALM 接到新问题上先构造一个带有解析解的小规模测试。比如上面的二次规划 demo 就是一个万能模板把目标函数换成你的任务把约束换成任务里真实用的约束再找一个小规模数据的解析解或高精度数值解做对照。验证通过后再增大数据规模否则宁可先不跑。从那以后我每次换到新的约束优化场景都强制自己走一遍“构造解析解→跑demo→画history.res→确认KKT残差”的流程这一步直接劝退了我偷懒的冲动。希望这份拆解能帮你少踩一点调参的坑也希望你改完代码后能回头把这几个经验块补得比我更厚——毕竟约束优化这东西真正跑起来才知道坑在哪里。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
