二次规划积极集法：从几何直觉到工程实现详解

约束优化的坑我早年踩得最深的一类就是二次规划。表面上只是“目标函数是二次的、约束是线性的”但真上手你会发现直接用无约束最优解去投影不靠谱硬套拉格朗日又容易被不等式条件绕晕。后来把积极集法active set method的完整逻辑捋清楚才算把这类问题看明白。这篇文章不绕弯子直接从几何直觉讲到算法推导再带一个能手算到底的完整例子和一份能跑的Python骨架把二次规划的积极集法讲透。适合正在学最优化、在调MPC控制器、或者被带约束最小二乘折腾过的人参考。1. 先定位二次规划和积极集法解决的是哪类问题1.1 标准二次规划长什么样所有讨论都建立在同一个标准型上[ \min_x \quad \frac{1}{2}x^T H x g^T x ] [ s.t. \quad a_i^T x \le b_i, \quad i1,\dots,m ]这里 (H) 是Hessian矩阵工程里通常要求它是半正定的否则问题本身可能是非凸的积极集法会面临更多麻烦。(g) 是线性项系数(a_i^T x \le b_i) 是一组线性不等式约束。目标函数里的 (1/2) 纯粹是为了求导后系数好看后面所有推导都沿用这个约定。如果约束里还有等式也完全可以并进来把等式视为始终处于积极集合中的特殊约束即可。在实盘中二次规划最常见的来源有三个一是模型预测控制MPC每一步的滚动优化二是带约束的最小二乘问题比如仓位调整、几何配准三是序列二次规划SQP方法内层的子问题。这些场景的共同特点是问题规模不算大但要求实时、可重复、能热启动。1.2 为什么叫“积极集法”“积极集”这个名字英文叫 active set也常翻译成“有效集”或“起作用约束集”。它的含义很直白在所有不等式约束中真正起作用的只有那些让等号成立的约束。比如约束 (x_1x_2\le 2)如果当前点刚好落在 (x_1x_22) 这条线上这个约束就是“积极”的如果当前点离这条线还有距离那它对当前点的行为就没有影响是“消极”的。积极集法的核心假设很朴素如果我们提前知道了最优解处哪些约束是积极的那么不等式约束就可以退化成等式约束剩下的问题就是一个带等式约束的二次规划求解难度直线下降。算法要做的事情就是不断猜测、修正这个积极集直到它和最优解处的真实积极集一致。这个思路放到现实中特别像“摸黑爬山的时候先判断哪面墙是挡着你的然后贴着墙走”。你不必同时考虑所有约束每次只需要盯住当前这面“墙”。1.3 它适合什么场景不适合什么场景先说结论积极集法是中小规模、热启动友好型算法。如果你要解几万个变量的大规模QP或者在线求解不知道初始可行点的通用问题内点法往往更稳。但如果你要在一个控制器里每几十毫秒解一次规模几百维的QP而且上一时刻的积极集在下一时刻大概率还成立那积极集法执行效率经常比内点法高一个量级。积极集法的另一个优势是迭代解释性强。每一次迭代要么朝可行域内部走一步要么向边界靠拢输出的是一个明确的可行点序列而不是内点法那种从可行域内部“蹭”过去的路径。这对工程debug非常有帮助。2. 核心机制拆解工作集、下降方向与乘子增删规则2.1 工作集与真实积极集的差距算法维护一个集合叫工作集working set记为 (W)。它本质上是算法当前猜测的“积极集”。工作集里的约束都要求在当前点处满足等式[ a_i^T x b_i, \quad i \in W ]注意工作集里的约束在数值上必须真的落在边界上否则整个推导会崩。所以每一轮迭代都从当前点出发只考虑工作集内的约束求一个能让目标函数下降、同时不违反工作集等式条件的搜索方向 (d)。搜索方向 (d) 必须满足[ a_i^T d 0, \quad \forall i \in W ]这个条件的意思是沿着 (d) 走一小步工作集内所有约束继续保持等式成立不会撞“墙”。如果 (d0)说明当前点在工作集约束下的最优性条件已经满足这时就需要判断工作集是不是真实积极集于是引出了乘子检验。2.2 乘子符号就是增删约束的判决依据当搜索方向 (d0) 时当前点是满足当前工作集约束的稳定点。接下来要回答一个关键问题工作集里的约束哪些该删掉标准做法是计算KKT乘子。对于形如 (a_i^T x \le b_i) 的不等式约束在最优解处应满足[ Hx g \sum_{i \in W} \lambda_i a_i 0 ]并且对所有积极约束有 (\lambda_i \ge 0)。这个符号约定非常重要后面踩坑部分我会重点讲。如果某个 (\lambda_i 0)说明删掉这个约束后目标函数还能继续下降。实际操作中通常选择最负的乘子对应的约束移出工作集。为什么要移最负的而不是随便移一个从优化角度看最负乘子对应的是“当前最该被释放”的约束释放它对目标函数的收益最大。虽然理论上移任何一个负乘子都能保证算法不终止但工程上从最负的开始迭代次数通常更少。2.3 明确下降方向后怎么确定步长当 (d \ne 0) 时我们把当前点沿 (d) 方向移动一段距离 (\alpha)。但这里有个约束你不能走太远否则会穿破某条不在工作集里的不等式约束。对每条不在工作集里的约束 (a_i^T x \le b_i)沿方向 (d) 行走时它允许的最大步长是[ \alpha_{\max} \min_{i \notin W,\ a_i^T d 0} \frac{b_i - a_i^T x}{a_i^T d} ]这个公式的逻辑是如果 (a_i^T d \le 0)说明沿 (d) 走这条约束会越来越松不用管只有当 (a_i^T d 0) 时约束才是“迎面而来”的需要计算它还能容忍走多远取所有这类约束中的最小值。与此同时我还会计算一个“无约束精确步长”。在纯二次目标函数下沿固定方向 (d) 的目标函数是 (\alpha) 的一元二次函数最优步长可以直接写出[ \alpha^* -\frac{(Hxg)^T d}{d^T H d} ]最终取 (\alpha \min(\alpha^, \alpha_{\max}))。如果 (\alpha^ \alpha_{\max})说明还没撞到新约束就有最优步长直接走 (\alpha^*)工作集不变如果 (\alpha \alpha_{\max})说明走到新边界了要把对应的新约束加入工作集。这一步是整个算法最像“贪心”的地方走一步看一步走到哪堵墙就抱住哪堵墙。2.4 为什么说“符号约定”是新手最容易翻车的地方网上很多资料用的约束写法不同有的写成 (a_i^T x \ge b_i)有的用 (Axb) 的等式形式乘子符号判断标准也跟着变。我自己早期就是吃了这个亏照着A教程推导代码却沿用B教程的符号结果KKT乘子怎么算都反着。我建议统一使用 (a_i^T x \le b_i) 作为标准型并固定KKT条件为[ Hx g \sum_{i \in W} \lambda_i a_i 0, \quad \lambda_i \ge 0 ]这样一来判断规则就一句话乘子为负的约束要删乘子全部非负则停止。如果你换成等式形式 (Axb)符号判断往往会变成“乘子非正停止”很容易绕晕。3. 手推一遍只用三张纸就能算完的完整示例3.1 算法标准流程速览在给出手算例子之前先把算法骨架摆出来给定可行初始点 (x_0)初始化工作集 (W) 为当前点处所有积极约束。计算下降方向 (d)要求工作集内约束满足 (a_i^T d0)。若 (d0)计算KKT乘子全部非负则停止否则移出最负乘子对应的约束回到第2步。若 (d \ne 0)计算允许的最大步长 (\alpha_{\max})以及精确步长 (\alpha^)令 (\alpha\min(\alpha^,\alpha_{\max}))。如果 (\alpha\alpha_{\max})把对应新约束加入工作集。更新当前点 (x \leftarrow x \alpha d)回到第2步。这个流程每一步都有明确含义建议把它抄在纸上手算示例时对照着看。3.2 一个能完整走到底的小问题为了展示算法里“删约束、走边界、加约束、停止”的完整链路我选了这样一个例子[ \min_x \quad \frac{1}{2}\left[(x_1-2)^2 (x_2-2)^2\right] ]三条约束[ -x_1 \le 0,\quad -x_2 \le 0,\quad x_1x_2 \le 2 ]目标函数的最优点本应是 ((2,2))但它不满足 (x_1x_2\le 2)。直观上最优解会被第三条约束压回边界最终落在 ((1,1))。下面用积极集法一步步走。为了和标准型对齐把目标函数展开[ f(x) \frac{1}{2}(x_1^2x_2^2) - 2x_1 - 2x_2 4 ]所以 (HI)(g(-2,-2)^T)。约束三条分别是[ a_1(-1,0)^T,\ b_10 ] [ a_2(0,-1)^T,\ b_20 ] [ a_3(1,1)^T,\ b_32 ]初始点取 (x_0(0,0))。这个点处第1、第2条约束都取等号所以初始工作集 (W{1,2})。3.3 完整的迭代推导手记第一轮迭代在 (x_0(0,0)) 处梯度为[ Hx_0g(-2,-2)^T ]工作集约束要求 (d_10) 且 (d_20)所以 (d0)。这一步说明当前点在工作集约束下已经是最优的了需要检查乘子。解[ (-2,-2)^T \lambda_1(-1,0)^T \lambda_2(0,-1)^T 0 ]得到 (\lambda_1-2)(\lambda_2-2)。两个乘子都为负说明当前两个约束都不该一直卡着。按规则移出最负乘子对应的约束这里两个相等随便选一个我移出第2条约束 (x_2\ge0)工作集变为 (W{1})。第二轮迭代工作集只包含 (-x_1\le0)所以方向 (d) 必须满足 (d_10)。在零空间里求下降方向实际就是最小化[ \frac{1}{2}d_2^2 - 2d_2 ]最优解为 (d_22)所以 (d(0,2)^T)。这个方向沿着 (x_2) 轴往上走。计算允许步长。工作集外的约束有第2条 (-x_2\le0) 和第3条 (x_1x_2\le2)。沿 (d(0,2))第2条约束自动满足不用管第3条要求[ (0,2)\cdot(\alpha d) 2\alpha \le 2 ]所以 (\alpha_{\max}1)。精确步长沿方向最小化[ (-2,-2)^T\cdot(0,2)\alpha \frac12 \alpha^2 (0,2)^T(0,2) -4\alpha 2\alpha^2 ]最优步长也是 (\alpha^*1)。因此走到 (x(0,2))并触发新约束 (x_1x_2\le2)把它加入工作集现在 (W{1,3})。第三轮迭代在 (x(0,2)) 处梯度为[ (-2,0)^T ]工作集两条约束分别是 (-x_10) 和 (x_1x_22)联立让 (d_10) 且 (d_1d_20)只能得到 (d0)。算乘子[ (-2,0)^T \lambda_1(-1,0)^T \lambda_3(1,1)^T 0 ]解出来 (\lambda_30)(\lambda_1-2)。第1条约束乘子为负说明 (-x_1\le0) 这个约束不该留在工作集里把它删掉工作集变为 (W{3})。第四轮迭代在 (x(0,2)) 处工作集只有第3条约束。方向 (d) 要满足 (d_1d_20)不妨设 (d(s,-s)^T)。代入目标函数沿方向的表达式[ \frac{1}{2}(s^2s^2) (-2,0)^T\cdot(s,-s) s^2 - 2s ]最优点是 (s1)所以 (d(1,-1)^T)。这个方向很有意思它沿着 (x_1x_22) 这条边界线往下滑。计算步长。工作集外约束 (-x_1\le0) 和 (-x_2\le0) 中沿 (d) 方向 (x_2) 是减小的所以真正限制的是 (x_2\ge0)[ 2 \alpha(-1) \ge 0 \quad\Rightarrow\quad \alpha \le 2 ]也就是 (\alpha_{\max}2)。精确步长是[ (-2,0)^T\cdot(1,-1)\alpha \frac12\alpha^2(1,-1)^T(1,-1) -2\alpha \alpha^2 ]最优步长 (\alpha^*1)小于2所以走 (\alpha1)到达 (x(1,1))工作集不变仍为 (W{3})。第五轮迭代在 (x(1,1)) 处梯度为[ (-1,-1)^T ]工作集约束要求 (d_1d_20)。同样设 (d(s,-s))沿方向目标变化为[ \frac{1}{2}(s^2s^2) (-1,-1)^T\cdot(s,-s) s^2 0 s^2 ]最优点是 (s0)所以 (d0)。计算乘子[ (-1,-1)^T \lambda_3(1,1)^T 0 ]得到 (\lambda_31 \ge 0)。所有乘子非负满足停止条件。最终最优解就是 ((1,1))目标函数值为[ \frac{1}{2}\left[(1-2)^2(1-2)^2\right]1 ]这个例子虽然只有三维空间里的三条约束但完整展示了移出约束、沿边界前进、加入新约束、最终停止的全部环节。对照代码看这个手推流程基本不会再糊涂。4. 代码骨架与工程化落地的关键细节4.1 一份教学演示用Python实现下面这份代码严格对应上面的算法流程。它不追求极致性能但结构清晰适合用来对照推导理解。核心函数假设约束形式统一为 (Ax \le b)。import numpy as np def null_space(A, tol1e-10): 返回 A 的零空间的一组正交基列向量 if A.size 0: return np.eye(A.shape[1]) u, s, vh np.linalg.svd(A) rank (s tol).sum() return vh[rank:].T def active_set_qp(H, g, A, b, x0, tol1e-8, max_iter100): 求解 min 0.5 x^T H x g^T x s.t. A x b 参数: H: n x n 对称正定矩阵 g: n 维向量 A: m x n 约束矩阵每行一个约束的法向量 b: m 维向量 x0: 初始可行点 返回: 最优解 x m A.shape[0] n H.shape[0] # 初始化工作集所有在当前点处活跃的约束 W [i for i in range(m) if np.isclose(A[i] x0, b[i], atoltol)] x x0.copy() for _ in range(max_iter): grad H x g # 计算下降方向 d满足 A[W] d 0 if len(W) 0: if np.linalg.norm(grad) tol: return x d -grad else: Aw A[W] Z null_space(Aw) if Z.shape[1] 0: d np.zeros(n) else: # 在零空间上求解子问题 Hz Z.T H Z gz Z.T grad v np.linalg.solve(Hz, -gz) d Z v # 如果 d 接近零做乘子检验 if np.linalg.norm(d) tol: if len(W) 0: return x Aw A[W] # 用最小二乘解乘子Aw^T lambda -grad lam, _, _, _ np.linalg.lstsq(Aw.T, -grad, rcondNone) if len(lam) 0 or lam.min() -tol: return x # 移出最负乘子对应的约束 j int(np.argmin(lam)) W.pop(j) continue # 计算沿 d 的最大可行步长 alpha_max np.inf idx_add -1 for i in range(m): if i in W: continue ai_d A[i] d if ai_d tol: alpha_i (b[i] - A[i] x) / ai_d if alpha_i alpha_max: alpha_max alpha_i idx_add i # 无约束精确步长 denom d H d if denom tol: alpha alpha_max if np.isfinite(alpha_max) else 1.0 else: alpha -(grad d) / denom # 收敛到新边界时加入新约束 if np.isfinite(alpha_max): if alpha alpha_max - tol: alpha alpha_max if idx_add ! -1: W.append(idx_add) x x alpha * d raise RuntimeError(达到最大迭代次数未收敛)用这份代码解前面的例子直接传入H np.eye(2) g np.array([-2.0, -2.0]) A np.array([[-1.0, 0.0], [0.0, -1.0], [1.0, 1.0]]) b np.array([0.0, 0.0, 2.0]) x0 np.array([0.0, 0.0]) x_opt active_set_qp(H, g, A, b, x0) print(x_opt) # 应该接近 [1.0, 1.0]代码里的关键设计在于null_space函数。它用SVD求工作集约束矩阵的零空间理论上很稳但每次迭代都做SVD规模大时偏慢。工程化实现中通常会做矩阵分解的低秩更新把上一轮的分解结果通过秩一更新塞进新约束这里就不展开了。4.2 初始可行点怎么找积极集法要求初始点可行否则“当前点处活跃约束”这个说法无从谈起。那遇到初始点不可行怎么办一个实用的做法是先用单纯形法或内点法求解一个“可行性问题”[ \min_{x,s} \quad \sum_i s_i ] [ s.t. \quad a_i^T x - s_i \le b_i,\quad s_i \ge 0 ]如果最优解中所有 (s_i0)就得到了一个可行初始点如果最优值大于0说明原问题本身无可行域这是建模问题而不是算法问题。很多现代QP求解器可以直接输出这个初始可行点不需要你手动处理。4.3 数值问题和性能优化方向代码骨架做好以后真正投入使用前要处理几个工程问题。第一工作集约束的消冗。手算例子中的约束很干净随便哪两条都线性无关。但真实模型里经常出现同一几何边界的微小偏移、重复不等式导致工作集矩阵行秩不足。此时SVD虽然能继续算但乘子里会出现极大的正负项相互抵消数值上非常难看。建议在加入工作集前先检查拟加入约束是不是已有约束的线性组合或者用QR分解配合列主元筛选独立约束。第二热启动。在MPC这类场景上一时刻的解和当前时刻的解变化不大上一时刻的积极集很可能就是当前时刻的好初值。直接用上一次的工作集作为当前的工作集并以上一时刻最优解作为 (x_0)常能节约大量迭代。这也是积极集法在现代MPC里依然有生命力的重要原因。第三大规模下的矩阵更新。真正实用的积极集法不会每轮重新SVD而是维护工作集约束矩阵的QR分解或Cholesky因子通过增删一行一列做低秩更新。这方面成熟开源实现可以参照qpOASES思路就是基于此。5. 避坑实录我实现积极集法时遇到的那些问题5.1 工作集来回切换导致不收敛有一次在调一个带等式和不等式混合约束的QP发现迭代在几个固定工作集之间反复横跳目标函数值不再下降但算法也没报错。后来逐步打印每轮的乘子才发现问题出在终止容差上我的tol设到了1e-12而数值误差导致乘子在小正小负之间波动约束被频繁删了又加。解决办法是把终止容差放宽到1e-8同时给移出约束加一个“滞后”阈值比如只有乘子小于-1e-7才允许移出避免微小数值抖动让约束反复进出工作集。这个技巧听起来很土但真的能救回大量调试时间。5.2 约束尺度差异大导致步长判断失准实际模型里约束的系数尺度可以差出好几个数量级。比如一个是 (1000x_1 x_2 \le 5)另一个是 (0.001x_1 x_2 \le 3)。这时候 (a_i^T d) 的计算结果会受到大尺度约束主导小尺度约束很容易被数值误差淹没。我的处理习惯是提前对约束行做归一化让每行的 (|a_i|_21)同时把对应的 (b_i) 也除以同样的系数。这样步长公式里的分子分母都在同一量级判断谁先触边界才靠谱。5.3 乘子符号写反怎么排查乘子符号问题是最隐蔽的。如果你发现算法总是把本应保留的约束删掉或者明明不是最优却提前停止先别怀疑方向计算大概率是符号约定不一致。我自己的排查技巧是构造一个只有单一不等式约束的二维问题手解出最优解和乘子然后用代码打出来对照。比如约束 (x_1\le1)目标中心在 (x_12)那么最优解在 (x_11)乘子应该为正。如果代码算出负的就把KKT里的正负号整体反转。这种单约束问题几秒钟就能定位问题。5.4 为什么我的迭代次数比别人多积极集法的迭代次数高度依赖初始工作集的猜测质量。如果你总是从空工作集开始每次都从无约束最优方向出发往前走碰到边界再加约束遇到负乘子再删约束那在一堆约束同时生效的问题里迭代次数会明显偏多。除了热启动另一个建议是初始工作集不要只加入严格活跃的约束也可以把那些“当前点离边界非常近”的约束放进去距离阈值设为比如 (10^{-6}) 级别的容差。这样能减少第一轮就大步穿破约束的概率让算法更快进入边界搜索模式。最后再分享一个我在实际工程里的体会积极集法这套东西看起来公式多、符号绕但真正理解价值之后你会发现它其实是“用几何直觉驱动算法”的典范。每一轮操作都有明确的几何含义要么贴着墙走要么换一堵墙要么发现自己其实已经在最优点。这种可解释性在工程调试里太重要了比黑盒求解器出了问题只能干瞪眼要强得多。我建议所有做带约束优化的朋友哪怕最终在工程里用现成库也至少手推一遍类似上面这种小例子把乘子增删、步长截断这些细节真正印在脑子里。等你需要调参、定位数值问题、或者给领导解释“为什么这个解没落在边界上”的时候这些基础会帮你省下大量时间。