ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

增广拉格朗日乘子法实战详解:从数学原理到代码调参

增广拉格朗日乘子法实战详解:从数学原理到代码调参 大多数做工程优化的朋友第一次接触增广拉格朗日乘子法Augmented Lagrangian Method简称 ALM往往是从论文里的某个公式开始的然后在自己上手实现时被一堆细节劝退。这篇文章我想把它彻底讲透包括数学直觉、推导过程、可运行的代码、调参经验以及它和罚函数法、拉格朗日法、ADMM 这些方法之间的区别和联系。内容主要面向有一定优化基础、但想在实际问题里动手用 ALM 的读者也适合被各种教材里抽象推导劝退的人——我会尽量把每个公式背后“为什么长这样”说清楚。1. 为什么需要增广拉格朗日乘子法两种经典方法的困局要理解 ALM 为什么有价值先得从它诞生的背景说起。我们在工程中遇到带约束优化问题时最经典的两条路是二次罚函数法和纯拉格朗日乘子法但这两条路都有各自的硬伤。1.1 二次罚函数法参数越大越病态二次罚函数法的思路很直接把约束违反程度以平方形式加到目标函数上。比如对于等式约束优化问题min f(x) s.t. h(x) 0我们构造罚函数P(x; μ) f(x) (μ/2) * ||h(x)||²然后在一串递增的 μ 下求无约束极小。看起来很简单但实际跑起来问题一堆。当 h(x) 不为零时惩罚项会强烈迫使 x 回到可行域内。理论上当 μ 趋向无穷大时极小点会收敛到原问题的可行解。但在数值上μ 一旦过大Hessian 矩阵的条件数会随 μ 线性恶化梯度下降或牛顿法在这种病态问题上几乎跑不动。我在实际求解带强非线性等式约束的流体参数估计问题时μ 加到 1e6 左右时无约束优化子问题就基本无法收敛了。更重要的问题是纯粹罚函数法的收敛精度受限于 μ 而非机器精度。即便 μ 取到 1e8最终约束残差也可能只到 1e-4 量级这对很多工程问题是不够用的。1.2 纯拉格朗日乘子法的困境对非凸问题缺乏正则化第二种思路是引入拉格朗日乘子 λ直接对拉格朗日函数 L(x, λ) f(x) - λ^T h(x) 做极小化同时更新 λ 来满足约束。对凸问题这本质上是求解 KKT 条件的鞍点理论很漂亮。但对非凸问题拉格朗日函数关于 x 可能既非凸也非凹直接做 min-max 迭代时数值稳定性很差。我在做机构运动学约束的参数辨识时尝试过纯拉格朗日方法λ 迭代经常振荡最终只能靠缩小步长勉强压制但收敛速度又慢得让人无法接受。此外如果 f(x) 本身没有强凸性L(x, λ) 关于 x 可能根本没有下界极小化子问题本身就不良定义。这两种方法的缺陷恰好展示了 ALM 的核心动机用二次正则项提供强凸性稳定子问题同时用乘子迭代保证收敛到精确解。也就是说ALM 同时拿到了两者的优点。2. ALM 的核心数学直觉与公式推导增广拉格朗日乘子法看起来只是在拉格朗日函数后面加了一个二次项但这一步既是保数值稳定性的关键也有很深刻的几何意义。这节我们拆开来看。2.1 增广拉格朗日函数的结构对等式约束问题ALM 的增广拉格朗日函数定义为L_ρ(x, λ) f(x) - λ^T h(x) (ρ/2) * ||h(x)||²其中 ρ 0 是罚参数λ 是乘子向量。与纯拉格朗日函数相比多了最后这个二次正则项。关键点在于L_ρ 关于 x 的 Hessian 近似为 ∇²f ρ * (∇h)(∇h)^T只要 ρ 取得足够大即使 f 高度非凸整个 L_ρ 在 h(x)0 的局部邻域内也可能成为关于 x 的强凸函数。这就是它比纯拉格朗日稳定得多的原因。2.2 乘子更新公式是怎么推出来的先说结论。标准的 ALM 迭代是固定 λ_k求解 x_k argmin_x L_ρ(x, λ_k)更新乘子 λ_{k1} λ_k - ρ * h(x_k)第二个式子初看很突兀为什么要用 ρ 乘以约束残差来更新我们从 KKT 条件出发推一遍。原问题的最优解 x* 必然满足 ∇f(x*) - (∇h(x*))λ* 0且 h(x*) 0。而子问题的极小点 x_k 满足∇f(x_k) - (∇h(x_k))λ_k ρ * (∇h(x_k)) h(x_k) 0整理一下∇f(x_k) - (∇h(x_k)) [λ_k - ρ * h(x_k)] 0对比 KKT 条件可以发现如果令λ_{k1} λ_k - ρ * h(x_k)那么 x_k 恰好是“以 λ_{k1} 为乘子”的拉格朗日函数的一阶驻点条件。换句话说乘子更新就是在用当前约束残差去修正乘子估计使得下一步的极小化更接近真正的拉格朗日驻点。这里还有一个很值得注意的地方从上面的推导来看λ 更新和二次罚项协同工作。即使 ρ 并不太大随着 λ 逼近 λ*x_k 也会逼近 x*而 h(x_k) 会趋于 0。这正是 ALM 能获得精确约束满足的原因——它不依赖 ρ 无穷大而是靠 λ 的自适应逼近。2.3 为什么它叫“精确拉格朗日”一个直观理解 ALM 的角度是增广项改写了约束惩罚的方式。罚函数法中约束残差被“直接压小”所以 μ 需要非常大才能让违反约束的成本足够高。ALM 中乘子 λ 不断修正最小化目标的位置相当于一种模型预测校正机制。二次项起到“局部锚定”的作用防止 λ 震荡太大让迭代路径更平滑。我习惯用一个比喻来向同事解释 ALM纯拉格朗日法像是一个只知道风向的船长每步都把方向调到目标点但风浪一大就原地打转罚函数法像是在船上不断增加压舱物越大越稳但船也被压得动弹不得ALM 是两者结合既知道风向又保留航速同时用适度的压舱物保证不翻船。3. 手写一个 ALM 求解器从推导到可运行代码理论聊完直接上代码。下面我用 Python 实现一个针对等式约束非线性优化问题的 ALM 求解器然后拿一个经典测试问题验证效果。这里的关键是让读者能照着自己的问题替换目标函数和约束表达。3.1 一个简单的测试问题考虑一个有个非凸特性的约束问题这是网上常见教程用的例子也适合检验算法稳定性f(x) exp(x1) x1² 100 * x2² h(x) x1 x2 - 1这个问题的解析解可以算出来方便验证。f 的 Hessian 在 x1 小时的曲率比较小约束与梯度方向几乎垂直对 ALM 来说有足够代表性。3.2 用 scipy 作为子问题求解器ALM 的主循环不需要自己实现牛顿法或梯度下降子问题可以直接用 scipy.optimize.minimize 的 L-BFGS-B 或 trust-constr 方法。示例实现如下import numpy as np from scipy.optimize import minimize def augmented_lagrangian_equality(f, grad_f, h, grad_h, x0, lambda0None, rho1.0, max_iter100, tol1e-7): if lambda0 is None: lambda0 np.zeros(h(x0).shape) lam lambda0.copy() x np.array(x0, dtypefloat) history [] for k in range(max_iter): # 定义增广拉格朗日子问题 def lag_obj(x_inner): h_val h(x_inner) return f(x_inner) - lam h_val (rho / 2) * np.sum(h_val ** 2) def lag_grad(x_inner): h_val h(x_inner) return grad_f(x_inner) - grad_h(x_inner).T lam rho * grad_h(x_inner).T h_val res minimize(lag_obj, x, jaclag_grad, methodBFGS, options{maxiter: 500, gtol: 1e-10}) x res.x h_val h(x) constraint_violation np.linalg.norm(h_val, np.inf) # 乘子更新 lam lam - rho * h_val history.append({iter: k, x: x.copy(), lambda: lam.copy(), violation: constraint_violation}) if constraint_violation tol: print(f收敛于外层迭代 {k}约束残差 {constraint_violation:.2e}) break # 可选的 rho 递增策略后面会细说 if k 5 and constraint_violation / max(history[-2][violation], 1e-12) 0.8: rho * 2 return x, lam, history # 定义目标函数与约束 def f(x): return np.exp(x[0]) x[0]**2 100 * x[1]**2 def grad_f(x): return np.array([np.exp(x[0]) 2 * x[0], 200 * x[1]]) def h(x): return np.array([x[0] x[1] - 1]) def grad_h(x): return np.array([[1.0, 1.0]]) # 1 x 2 if __name__ __main__: x_opt, lam_opt, hist augmented_lagrangian_equality( f, grad_f, h, grad_h, x0np.array([0.0, 0.0]), lambda0np.array([0.0]), rho1.0 ) print(最优解:, x_opt) print(约束值:, h(x_opt)) print(乘子:, lam_opt)这段代码的思路很清晰外层循环维护乘子 λ 和内层子问题的初值每轮用当前 λ 和 ρ 构造一个带二次惩罚的目标函数并执行无约束极小化然后根据约束残差更新乘子。实测在我的机器上这个简单问题大概 6~8 轮外层迭代就收敛到约束残差 1e-8 以下。3.3 代码里容易被忽视的细节第一个是子问题初值的传递。每一步要把上一步的 x_k 作为下一步子问题的初始点。这是因为相邻两次子问题的最优解差异通常很小热启动能显著减少内层迭代次数。有些实现图省事直接用固定的 x0结果外层迭代次数没什么变化但内层每次都要重新收敛浪费时间。第二个是 lag_grad 中 h 对 x 的雅可比矩阵形状。上面例子中 grad_h 返回的是 (1,2) 矩阵但如果你定义多个等式约束grad_h 应返回 (m,n) 矩阵其中 m 是约束数。很多人第一次写这种代码时对不上维度导致梯度计算错误但目标函数看起来还正常问题会在乘子更新时突然放大。第三个是 BFGS 的梯度容差设置。如果不给 minimize 传 gtol默认值可能不够小子问题停在比较粗糙的位置导致约束残差振荡。我建议子问题求解器内部容差比外层停准则高至少两个数量级这样乘子更新才能得到足够好的梯度信息。3.4 一个反直觉的发现我在跑这个简单例子时试过不更新乘子、单纯放大 ρ结果收敛精确度始终上不去罚参数大到 1e8 之后还出现了浮点误差主导的现象。而一旦加入乘子更新即使 ρ 固定在 1.0约束残差也能持续下降到 1e-12。这直观印证了前面的结论ALM 收敛的精确性来自乘子估计不是来自罚参数的不断增大。4. 子问题求解与罚参数调优的实战经验算法结构容易理解真正让 ALM 从玩具代码变成可靠求解器的是参数调整策略。这块没太多教科书内容多半是实际调试中踩坑换来的。4.1 罚参数 ρ 的初始值怎么选初始 ρ 太小会导致子问题非凸极小化可能跑到局部解甚至发散太大会让 L_ρ 的 Hessian 条件数变差子问题收敛变慢。我在实践中一般根据约束残差的量级来估计如果 h(x) 的各分量初始量级在 0.1~1 之间ρ 可以从 1 或 10 开始。如果约束本身做了归一化量级接近 1ρ1 是不错的选择。如果约束残差初始非常大比如 100 量级可以先从 ρ0.1 开始避免一开始二次项压倒目标函数导致子问题把约束压得太快而目标函数恶化过多。总之 ρ 的选取和标度scaling关系密切。一个有效做法是先单独算一下 L_ρ 在初始点附近的 Hessian 条件数选择让条件数小于 1e4 的最小 ρ这能保证内层求解器有较好的收敛速度。4.2 递增策略不要盲目指数增长经典理论分析经常假设 ρ_k 趋向无穷大但实际工程中过度增大 ρ 会带来数值困难。我更推荐一种“按需增加”的策略只有当乘子更新后约束残差下降停滞时才增大 ρ。一个可行判据是连续若干轮约束残差下降率低于阈值。比如在外层迭代里记录 violation 的比值若连续 3 次残差比大于 0.9就把 ρ 乘以 2~5。固定 ρ 缓慢下降的情况完全不需要动 ρ。这种策略兼顾了收敛速度和数值稳定。4.3 终止准则怎么定才靠谱很多人只用约束残差 ||h(x)|| 来判定收敛这其实不够。理论上 ALM 的终止准则应该包含两个部分约束可行度||h(x_k)||∞ ≤ ε_feas乘子变化量||λ_{k1} - λ_k|| / (1 ||λ_k||) ≤ ε_mult后者衡量的是对偶变量是否稳定下来。如果只检查约束残差可能出现这种情况某个子问题求解器精度不足约束残差暂时很小但乘子还在大幅漂移。我的经验是两者同时满足才算收敛ε_feas 通常取 1e-6 到 1e-8ε_mult 可以放宽到 1e-5 到 1e-6因为对偶变量收敛通常比原始残差慢一些。4.4 子问题求解器的选择与容差配置我的经验是约 500 维以内、目标函数和约束梯度好算的问题用 L-BFGS 配合有限差分或解析梯度就够了高维或者有非光滑项的问题则常用 L-BFGS-B 处理简单界约束。真正复杂的大规模问题则会切到 Newton-CG 或 truncated Newton 方法因为增广项带来的曲率信息能在 Newton 型方法中得到充分利用。子问题内部求解时建议内层容差设置为外层停止准则的一半量级以下。比如想达到约束残差 1e-8内层梯度范数至少要达到 1e-10否则乘子更新会掺杂噪声严重时造成后期收敛困难。这类问题很隐蔽表面看算法不停循环实际上深层原因在子问题精度不够。5. 不等式约束与扩展ALM 在真实问题中的用法现实中更多问题是不等式约束比如 g(x) ≤ 0。直接套用等式约束版 ALM 是不行的但通过一些变换可以优雅地复用现有实现。5.1 引入等式约束的转换思路最简单的方法是把不等式约束转为等式约束加非负限制即添加松弛变量 s令g(x) s 0, s ≥ 0然后对等式约束 g(x) s 0 写增广拉格朗日函数同时要求 s 非负。子问题变成带界约束的极小化可以用 L-BFGS-B 这类方法处理。这个变换的代价是问题维度增加了但收益非常明显ALM 的所有现有理论都能沿用而且在 s 上进行投影操作也非常简单只需要把更新后的 s 裁剪到非负区间。5.2 更简洁的广义增广拉格朗日Generalized Augmented Lagrangian还有一种不需要显式引入松弛变量的做法直接把不等式约束嵌入到增广项中。定义L_ρ(x, λ) f(x) (1/(2ρ)) * ( ||max(0, λ - ρ*g(x))||² - ||λ||² )这种形式实质上是把不等式约束的乘子投影到了非负域。更新 λ 后λ_{k1} max(0, λ_k - ρ*g(x_k))。它的好处是问题维度不增加且对于某些非线性规划库更容易实现。我实际测试对比过两种方式对于小规模问题引入松弛变量的做法更稳定因为子问题只处理等式约束间的平衡对于大规模问题广义增广拉格朗日的无扰动版本省掉了松弛变量的存储和更新内存更友好。具体选哪一种取决于你对子问题求解器的熟悉程度和问题规模。5.3 一个带不等式约束的例子最小化 Rosenbrock 函数为了展示实际效果我以经典的 Rosenbrock 函数为例添加一个约束min f(x) (1-x1)² 100*(x2-x1²)² s.t. x1² 3*x2 ≤ 2用广义增广拉格朗日实现时核心就是子问题目标函数包含 max(0, ...) 的二次形式其余和等式约束版本几乎一致。实测这个例子在 ρ1 初始值下约 10 轮迭代内收敛到可行域边界约束残差 1e-7 量级。5.4 ALM 的扩展从 ADMM 到非光滑优化讲 ALM 就不能不提 ADMM因为 ADMM 本质上是把 ALM 应用于一个经过变量分裂的等价问题。具体做法是引入辅助变量 z把原问题转化为min f(x) r(z) s.t. Ax Bz c然后交替极小化 x 和 z再用 ALM 乘子更新。ADMM 的优势是当 f 和 r 都可以独立、廉价地极小化时整个算法比直接在原始变量上做 ALM 快得多。所以可以这么理解如果你的目标函数是可分离结构ADMM 是比 ALM 更好的选择如果没有这种结构ALM 本身更直接。另外对于带 L1 范数正则的非光滑问题ALM 通常需要配合近端算子使用也就是在子问题极小化中加近端项。这样扩展出来的算法常被称为近端增广拉格朗日Proximal ALM它在压缩感知、图像去模糊方面应用很广。6. 它和罚函数法、SQP、内点法的对比与选型建议实战中选优化算法不能只看理论收敛速度还要考虑实现难度、内存开销和可维护性。这里给出我在工程决策时常用的对比框架。6.1 与罚函数法的对比罚函数法实现最简单但精度受限于罚参数。如果约束量级不稳定、问题规模的 Hessian 结构特殊罚函数法很容易遇到病态。ALM 只多维护一个乘子向量实现复杂度增加不多但收敛精度和稳定性明显更好。只要不是只需要一个粗糙可行解的场合我都建议优先使用 ALM 而非纯罚函数法。6.2 与 SQP 的对比SQP 每步需要求解一个带约束二次规划子问题对二阶信息的利用更充分在中小规模、约束光滑的问题上收敛速度常常最快。但 SQP 的代价是每个内层子问题本身就是一个约束优化问题代码实现和调试复杂度较高。ALM 的子问题则是无约束或简单有界约束的极小化可以用更成熟的通用无约束求解器工程实现难度低一大截。如果问题规模不大、约束个数不多我会选择 SQP 来获得更快的末端收敛如果问题规模大、约束结构复杂或者子问题本身已经有很多现成的无约束求解器ALM 往往是更实际的选择。6.3 与内点法的对比内点法在处理线性约束和大规模稀疏问题上非常出色许多商业求解器的底层是内点法。但它对数障碍项的引入比较敏感处理不等式约束时需要路径跟踪对初始化要求高。ALM 不需要求解中心路径对初始点的鲁棒性更好。我的判断是如果问题变量之间耦合很小、约束稀疏内点法大概率更强但如果不是专业数值优化团队ALM 的调参难度和调试成本通常更低。6.4 一张选型表方法子问题难度内存占用收敛精度实现复杂度典型适用场景罚函数法无约束低受限于 μ极低快速得到粗糙可行解纯拉格朗日法无约束低理论上精确低凸性很强的特定问题ALM无约束低精确低一般非线性约束问题SQP约束QP中精确高小型高精度问题内点法线性/非线性系统高精确高大规模稀疏问题ADMM可分裂无约束中中等受步长影响中可分离结构的凸问题6.5 什么情况下不要用 ALMALM 不是万金油。如果约束引入了离散变量比如整数约束 x ∈ {0,1}ALM 完全无能为力需要用分支定界或启发式算法。此外如果目标函数的二阶不可导性非常严重且约束不光滑ALM 的子问题可能出现无界风险最好先做问题变换或在子问题中加近端正则项。再者约束数量极大比如几十万条且大部分非活跃时ALM 每轮对所有约束都计算增广项和梯度浪费很大。这时可能需要先做约束筛选或 active-set 策略本质上又回到 SQP 那套思路上去了。6.6 最后的建议从我做过的大大小小约束优化项目来看ALM 是性价比最高的通用算法之一。它的数学背景深厚但实现门槛很低而且容错能力强。只要先理解乘子更新和罚参数之间的平衡逻辑再注意子问题求解精度大部分非线性约束问题都能在四五个工作日内跑通。希望这篇文章能帮你少踩一些我当年的坑也欢迎在实践中遇到具体问题后进一步交流。
返回列表