简介增广拉格朗日算法Augmented Lagrangian Method的MATLAB源码实现面向优化算法学习者以及机器学习、图像处理、信号处理等领域的工程技术人员用于高效求解带约束的非线性优化问题。该算法在拉格朗日乘子法基础上引入惩罚项与增广项通过迭代中逐步增大的惩罚系数迫使约束逐渐满足同时以乘子更新来平衡目标函数与约束要求。压缩包体积仅2KB包含1个m脚本文件核心代码SolvePALM.m依据近端交替线性化最小化PALM思想完整实现了初始化赋值、决策变量交替更新、拉格朗日乘子修正、惩罚系数动态调整以及收敛判定等关键步骤代码逻辑清晰紧凑很适合逐段研读和二次修改。目前该资源已有3092人学习对于希望快速上手PALM算法、理解增广拉格朗日乘子法工程实现细节或需要在MATLAB中搭建自定义约束优化求解流程的读者都是一份高价值的参考范例无论是算法对比实验还是课程设计都能从中获得直接的代码支撑。1. 增广拉格朗日乘子法为什么罚函数法总是差一口气做约束优化的人多半被罚函数法折磨过ρ调大了Hessian病态迭代半天不动ρ调小了约束又不够紧。增广拉格朗日乘子法Augmented Lagrangian Method就是来解决这个矛盾的——在罚函数里补一个拉格朗日乘子项用有限大小的ρ也能精确收敛到满足约束的解。这篇笔记把增广拉格朗日算法的原理、MATLAB实现和踩坑点从头讲透覆盖等式约束和不等式约束代码可以直接复制跑通。适合正在做稀疏优化、图像重建、最优控制的工程师和研究生尤其是想自己掌控迭代过程、不想被黑匣子求解器困住的人。2. 从罚函数到乘子法增广拉格朗日的数学原理与选型理由2.1 罚函数法的两个痛点病态与收敛慢刚接触约束优化时最容易想到的就是罚函数法把约束违反量加一个二次惩罚项到目标函数里。对等式约束c(x)0问题变成最小化φ(x)f(x)ρ/2·‖c(x)‖²。理论上ρ越大解越靠近可行域但数值上灾难很快就来了在可行方向附近惩罚项贡献的曲率大约是ρ倍目标函数自己的曲率反而被淹没。Hessian矩阵的条件数随ρ线性增长用fminunc或牛顿法迭代时每一步步长都被最小特征值限制收敛速度慢得让人怀疑人生。更麻烦的是罚函数法只有在ρ→∞时才精确满足约束实际你只能取一个有限ρ等于主动接受一个近似解。这就是为什么纯罚函数法在工程里很少单独用。我的第一版图像配准代码就是栽在这里ρ调到1e4约束残差还在10⁻³量级Hessian已经是病态了。后来换成增广拉格朗日乘子法用一个有限的ρ把约束残差压到10⁻⁶附近乘子项把解“拉”回边界才算真正解决问题。需要记住的关键点是罚函数法把乘子信息丢掉了而乘子恰好是约束在最优解处的边际价值丢掉它等于扔掉一半信息。2.2 增广拉格朗日乘子法的推导把拉格朗日项加进罚函数增广拉格朗日函数定义很简单L(x, λ, ρ) f(x) λᵀc(x) ρ/2·‖c(x)‖²对比纯罚函数多出的λᵀc(x)就是拉格朗日乘子项。为什么加这一项就能避免ρ→∞可以从KKT条件看。拉格朗日函数L的驻点条件为∇f(x) ∇c(x)ᵀ(λ ρc(x)) 0把它和原问题最优解处的KKT条件∇f(x*) ∇c(x*)ᵀλ* 0对比会发现λ ρc(x)正在扮演最优乘子λ*的角色。因此每次外迭代最自然的更新就是λ ← λ ρ·c(x)这个更新把所有约束残差“累积”到乘子里。即使ρ不取无穷大乘子λ也会逐步收敛到λ*原始残差c(x)随之被压下去。经验上等式约束凸问题只要ρ 0都能收敛到精确解非凸问题则需要ρ足够大才能稳定。不等式约束稍微绕一点。如果约束是g(x) ≤ 0可以引入非负乘子μ乘子更新变成μ ← max(0, μ ρ·g(x))子问题里构造增广拉格朗日函数时可以省略常数项只用L(x, μ, ρ) f(x) ρ/2·Σ max(0, g_i(x) μ_i/ρ)²这里为什么没有单独的μᵀg(x)项因为把max(0, g_i μ_i/ρ)²展开在起作用的分量上恰好包含μ_i g_i ρ/2·g_i²再加一个与x无关的常数项忽略常数不会影响子问题的最优解。这个写法比传统μᵀg ρ/2‖g‖²更直观也更容易用max实现。2.3 什么时候该选增广拉格朗日而不是ADMM或内点法选算法不能只看名字。我一般先做一张对比表再对号入座方法适用场景优点主要代价纯罚函数约束很软、精度要求低实现最简单ρ→∞数值病态增广拉格朗日中等精度、约束数量不多、子问题可用无约束求解器有限ρ精确收敛外迭代需要解多个无约束子问题ADMM目标可分、大规模分布式近端算子分解、内存低要求问题结构可分离调ρ仍然有玄学内点法中小规模高精度精度可达1e-10Hessian稠密、内存和病态问题明显如果问题只带几个线性或光滑非线性约束目标函数也光滑增广拉格朗日是最好的起点。如果你手里的目标函数含L1范数、全变差TV这类不可导项增广拉格朗日通常要配合近端算子这实际就变成了ADMM。反过来如果要求误差达到1e-10量级继续用ALM会非常吃力内点法这种直接求解KKT系统的方法更合适。工程上“差不多就行”的1e-6精度ALM完全够用而且迭代过程每一步都有明确的物理意义——乘子λ就是约束的影子价格能直接看出哪个约束在拖后腿。3. MATLAB实现增广拉格朗日算法核心框架与最小可运行代码3.1 问题定义等式约束与不等式约束的统一处理动手写代码之前先把问题转换成标准形式目标函数f(x)x是列向量等式约束ceq(x) 0返回列向量不等式约束g(x) ≤ 0返回列向量这里最容易被坑的是不等号方向。MATLAB的fmincon默认线性约束写的是A·x ≤ b但增广拉格朗日的乘子更新公式是按g(x) ≤ 0写的。如果你的原始约束是“大于等于”一定要先取负变成“小于等于”。比如需求是x1x2 ≥ 1写成g(x) 1 - x1 - x2 ≤ 0。方向反了乘子永远不更新约束残差也永远不会下降。初始化时等式乘子λ设为与ceq(x0)同维的全零向量不等式乘子μ设为与g(x0)同维的全零向量。ρ通常从1开始不要贪大。我见过很多新手直接把ρ设成1000理由是要“快速满足约束”结果子问题Hessian病态fminunc一步都走不动。3.2 核心迭代代码用MATLAB写乘子更新与子问题求解先给一个最小可运行的等式约束求解函数可以直接保存为alm_eq.mfunction [x, lam, info] alm_eq(f, ceq, x0, rho, tol, maxit) % ALM_EQ 增广拉格朗日乘子法求解等式约束优化 % min f(x) s.t. ceq(x) 0 % f、ceq 为函数句柄x0 为列向量 % rho 初始惩罚参数tol 收敛容差maxit 最大外迭代次数 x x0(:); lam zeros(size(ceq(x))); for k 1:maxit % 构造当前增广拉格朗日函数 L (xx) f(xx) lam * ceq(xx) (rho/2) * norm(ceq(xx), 2)^2; % 用 fminunc 解无约束子问题热启动初值用上一轮的 x opts optimoptions(fminunc, Display, off, ... OptimalityTolerance, 1e-8, MaxIterations, 500); x_old x; x fminunc(L, x, opts); % 更新乘子λ ← λ ρ·ceq(x) cval ceq(x); lam lam rho * cval; % 收敛判断原始残差和变量变化都足够小 if norm(cval, Inf) tol norm(x - x_old, Inf) tol info.iter k; info.rho rho; return; end end info.iter maxit; info.rho rho; end调用示例求解min (x1-1)^2(x2-2)^2约束x1^2x2^22% 求解 min (x1-1)^2 (x2-2)^2 s.t. x1^2 x2^2 2 f (x) (x(1)-1)^2 (x(2)-2)^2; ceq (x) x(1)^2 x(2)^2 - 2; [x, lam, info] alm_eq(f, ceq, [0; 0], 1, 1e-6, 100); fprintf(x [%.6f, %.6f]\n, x(1), x(2)); fprintf(lambda %.6f\n, lam); fprintf(约束残差 %.6f\n, ceq(x));代码的核心逻辑很直白外层循环解子问题、更新乘子、检查残差。注意fminunc每次的起始点用了上一轮的x这叫热启动能省掉大量重新搜索的迭代。如果初始点每次都重置为x0后半个问题的子问题求解会慢好几倍等于把增广拉格朗日的优势浪费掉。参数OptimalityTolerance我设成1e-8这对外迭代稳定很关键。乘子更新公式里cval来自子问题解如果子问题只求解到1e-4精度乘子更新的误差会被ρ放大外迭代容易来回震荡。所以宁可子问题多迭代几步也不要让乘子被噪声污染。3.3 参数初始化ρ、λ、tol怎么设才能不翻车下面这张参数参考表是我多年调参的默认值适合大多数工程问题参数常见范围默认取值说明ρ0.01 ~ 1001太小约束收敛慢太小子问题病态λ全零即可zeros可以先用松弛问题估计乘子加速μ全零即可zeros不等式乘子非负tol1e-8 ~ 1e-41e-6工程问题足够不要一开始就追1e-8maxit100 ~ 500200超过后先检查子问题精度ρ要不要自适应强烈建议要。我一般每隔几轮检查原始残差‖c(x)‖∞如果连续两轮残差下降都不到5%就把ρ乘以2。但ρ上限不要超过1e6否则子问题Hessian会病态到无法收敛。反过来当残差已经达到tol附近时不要再碰ρ保持住即可。还有一个容易被忽略的坑是量纲。如果f(x)量级在1e6约束残差量级在1e-4那么增广拉格朗日函数里两项的尺度差距极大fminunc的有限差分梯度会全部被大尺度项吞掉。这种情况请先做变量归一化或者手动给目标函数、约束乘一个缩放系数让它们至少在同一个数量级否则后面调什么都白费。4. 把算法套到真实问题上有界约束二次规划与fmincon对比4.1 算例求解带不等式约束的凸二次规划现在看一个带不等式约束的完整例子。问题min f(x)x1²2x2²x1x2≥1, x1≥0.2, x2≥0.3先写成g(x)≤0的标准形式g11-x1-x2, g20.2-x1, g30.3-x2MATLAB代码% 目标函数与不等式约束 f (x) x(1)^2 2*x(2)^2; g (x) [1 - x(1) - x(2); 0.2 - x(1); 0.3 - x(2)]; x [0; 0]; mu zeros(3, 1); rho 1; for k 1:200 % 增广拉格朗日子问题注意这里没有单独的 mu*g 项 L (xx) f(xx) (rho/2) * sum(max(0, g(xx) mu/rho).^2); x_old x; x fminunc(L, x, optimoptions(fminunc, Display, off)); % 乘子投影更新 gval g(x); mu max(0, mu rho * gval); prim norm(max(0, gval), Inf); if prim 1e-6 norm(x - x_old, Inf) 1e-6 break; end end fprintf(x1 %.6f, x2 %.6f\n, x(1), x(2)); fprintf(active constraint value: %.6f\n, g(1)); % 应约等于0运行结果会收敛到x≈(0.6667, 0.3333)。这个点让x1x21也就是第一个约束取等号而x1≥0.2和x2≥0.3都不起作用。对应乘子mu(1)0mu(2)、mu(3)应为0。这个结果正好满足互补松弛条件乘子只在起作用的约束上非零。为什么这里L里可以直接省略mu*g(xx)因为max(0, gmu/rho)平方展开后在起作用的分量上已经包含mu*g rho/2*g²再乘以rho/2后剩下的常数项与x无关子问题极值点不受影响。写成这样既简洁又不容易把符号搞错我后来所有不等式ALM代码都采用这个形式。4.2 与fmincon对比自带求解器在什么情况下更好同样的二次规划用fmincon一行就能写% fmincon 求解同样的线性不等式约束问题 [x_fc, fval_fc] fmincon((x) x(1)^2 2*x(2)^2, [0;0], ... [-1 -1; -1 0; 0 -1], [-1; -0.2; -0.3]);结果完全一致而且fmincon内部用了序列二次规划处理这种光滑小问题非常快。所以如果你的问题规模不大、约束只是简单的线性或边界直接用fmincon不要重复造轮子。什么时候应该自己写增广拉格朗日我遇到过三种场景值得写。第一种是目标函数来自外部仿真比如车辆动力学模型fmincon内部的有限差分会在每个迭代点扰动输入仿真器可能因为扰动越界直接报错ALM把子问题拆出来你可以用无导数的patternsearch或自己设计试探序来做子问题求解。第二种是约束数量很少但状态变量很多比如几十万维的图像重建fmincon每一步都要处理稠密K矩阵内存直接爆掉ALM搭配共轭梯度解子问题内存占用可以降到O(n)。第三种是你在研究对偶变量需要观察乘子λ随迭代的变化规律这只有自己写代码才能看清。4.3 如何判断算法真的收敛了看原始残差、对偶残差和KKT残差只看目标函数值不变就宣布收敛是不行的。ALM真正的停机条件应该包含三部分原始残差‖ceq(x)‖∞和‖max(0, g(x))‖∞都小于tol对偶残差乘子更新量‖λ_new - λ_old‖足够小KKT残差‖∇f(x) Σλ_i·∇ceq_i(x) Σμ_i·∇g_i(x)‖∞足够小第三个才是最优性条件。很多人迭代后约束残差很小但乘子还在慢慢漂移就是因为只看了原始残差。验证代码很简单% 在上面的不等式算例收敛后计算 KKT 残差 grad_f (x) [2*x(1); 4*x(2)]; grad_g (x) [-1 -1; -1 0; 0 -1]; % 每一列是一个约束的梯度注意 g 的定义顺序 kkt_res norm(grad_f(x) grad_g * mu, Inf); fprintf(KKT residual %.6e\n, kkt_res);在x(0.6667,0.3333)处∇f(1.3333,1.3333)grad_g*mu应为(-1.3333,-1.3333)残差接近机器精度。如果这个值不收敛说明子问题精度不够或乘子更新方向写错了。5. 避坑与常见问题增广拉格朗日MATLAB实现中的五个坑5.1 现象迭代发散到NaN迭代几次后x突然变成NaN或者fminunc报错“无法继续因为目标函数返回 Inf 或 NaN”。原因通常有两个一是ρ设置过大子问题Hessian病态fminunc在某一步被数值误差击穿二是约束函数在迭代过程中被送出了定义域比如sqrt(x)、log(x)、外部仿真模型遇到无解状态返回了NaN。解决在子问题求解前后打印norm(x)和ceq(x)定位是哪一轮爆掉的。给x加一个投影或边界保证不会落到定义域外。ρ从1开始不要大于1e4如果必须大ρ改用fmincon的trust-region-reflective算法它对病态Hessian容忍度更高。还可以用OutputFcn在迭代到NaN时直接抛出错误保存断点方便检查现场。5.2 现象乘子λ不收敛反而震荡乘子在外迭代中反复横跳约束残差也下不去。这是ALM最经典的问题。原因十有八九是子问题求解精度不够。乘子更新公式λ←λρ·ceq(x)是一阶累积如果子问题解x存在10⁻⁶量级的随机误差乘子每次都会被放大ρ倍。当ρ1000时一次误差就能把乘子推出理性范围。解决把fminunc的OptimalityTolerance从默认1e-6提高到1e-8甚至1e-10。其次优先使用解析梯度别依赖有限差分。最后检查ρ自适应策略是否过激ρ每次翻倍会让乘子更新量也翻倍震荡几乎是必然。我习惯rho rho * 1.5而不是*2。5.3 现象约束收敛了但结果不是最优解迭代结束后约束残差确实接近0但目标函数值明显不是全局最优甚至从一个初值出发得到了完全不同的解。原因是ALM只能保证局部最优。对于非凸目标或非凸约束子问题本身可能有多个局部极小点fminunc只是一个局部优化器它找到哪个极小点取决于热启动的位置。解决对非凸问题在外循环开始前先用多启动搜索一遍把最好解作为初始点。也可以在外循环的前几轮故意减小ρ让子问题带上更多目标函数自身的曲率降低掉进错误山谷的概率。另外一定要看KKT残差如果KKT残差为0但目标值不是预期值说明它确实收敛到另一个局部解——算法没坏是问题有多峰。5.4 现象不等式约束的乘子出现负值或约束一直违反不等式乘子μ更新后出现负数或者即使迭代很多轮约束残差依然很大。原因几乎都是不等式方向写反。假设你的条件是x1x2≥1如果直接把这行写成g(x)x1x2-1那么g本来就大于0乘子更新式max(0, μρg)会把μ不断推到非负但g依然大于0原始残差永远下不去。正确写法是g(x)1-x1-x2。解决写完之后先算一下初始点的g(x0)。如果约束要求g≤0但初始点处g为正说明方向错了。另外乘子更新必须是μmax(0, μρg)千万不要写μmax(0,μ)ρg后者在μ很大时会越过投影导致负乘子回流。5.5 现象子问题用fminunc报错“objective function returned Inf/NaN”这是MATLAB最常出现的错误之一但很多情况下不是算法问题而是max(0, gμ/ρ)这个投影操作在有限差分求梯度时跨过了不可导点。增广拉格朗日函数在g(x)μ/ρ0处一阶不可微。fminunc的有限差分会在该点左右各采一个点如果左边的梯度是0右边梯度是1差商会给出一个有限值通常没事但如果函数本身定义域超过那里比如经过log或sqrt就会产生NaN。解决对于只带不等式约束的问题可以把不可微点做光滑化替换用0.5*(v sqrt(v^2 eps^2))代替max(0, v)。其中eps取1e-8量级。另一种更省事的方法是子问题改用fmincon的sqp它能处理部分非光滑项且不会像fminunc那样在有限差分上翻车。但如果你追求高性能还是建议用近端算子直接解析求解L1这类项而不是硬塞给通用求解器。6. 进阶技巧用已知解测试和自动微分给增广拉格朗日护航6.1 构造一个已知解的问题验证你的ALM实现每次写完一个新的增广拉格朗日求解器我都会先跑一个能解析算出最优解的小问题。比如min (x1-2)² (x2-3)²s.t. x1x24手工计算最优解是x(1.5, 2.5)最优乘子λ1。把它喂给alm_eqf (x) (x(1)-2)^2 (x(2)-3)^2; ceq (x) x(1) x(2) - 4; [x, lam] alm_eq(f, ceq, [0; 0], 1, 1e-8, 50); % 期望 x [1.5; 2.5], lam 1同时算KKT残差norm(2*(x - [2;3]) lam*[1;1], Inf)。如果这一步都对后续换复杂问题才有底气。这个习惯救了我很多次写错方向或漏掉约束时已知解测试会第一时间暴露问题。6.2 用MATLAB符号工具箱算梯度避免有限差分误差子问题精度不够导致乘子震荡最彻底的解法是给fminunc提供解析梯度。用符号工具箱可以一次性生成梯度的函数句柄syms x1 x2 lam rho f_sym (x1 - 2)^2 (x2 - 3)^2; c_sym x1 x2 - 4; L_sym f_sym lam * c_sym (rho/2) * c_sym^2; grad_L matlabFunction([diff(L_sym, x1); diff(L_sym, x2)], ... Vars, {[x1; x2], lam, rho});然后在每轮迭代构造带梯度的目标函数并开启SpecifyObjectiveGradientL_with_grad (xx) deal(L(xx), grad_L(xx, lam, rho)); opts optimoptions(fminunc, SpecifyObjectiveGradient, true); x fminunc(L_with_grad, x, opts);注意max(0, gμ/ρ)这类不可导项不能直接用符号工具箱求梯度需要分段写。好在工程里最常麻烦我们的是等式约束和光滑不等式约束解析梯度足矣。6.3 rho自适应的一种实用策略最后分享一个固定参数表之外的技巧。每次外迭代记录原始残差p_k‖max(0,g(x))‖∞。如果p_k 0.5 * p_{k-1}说明约束残差没在下降把ρ乘1.5如果p_k已经小于tol则完全冻结ρ。同时设一个硬上限ρ超过1e5后强制提醒自己“该回头检查子问题精度了”。这个策略比从头到尾固定ρ稳得多也比简单粗暴每轮乘2快得多。我现在每次写增广拉格朗日求解器都先跑一遍已知解测试再开启rho自适应最后才上真实数据。希望帮到你。本文还有配套的精品资源点击获取