昨天一位做算法岗的朋友跟我聊起线性规划LP说他简历上写着“熟悉最优化”但被问到单纯形法Simplex Method的细节时第一反应是求解器不是都封装好了吗为什么还要知道内部怎么走这话我太熟了。把 LP Simplex 真正理解到位不是为了让你手工去解几十行几百列的规划问题而是为了让你在建模、调参、排查数值问题时脑子里能浮现出一幅清晰的几何画面当前解卡在多面体的哪个顶点下一步沿着哪条边滑为什么滑到某一步就停。这篇文章我从实现者的视角把单纯形法的几何直觉、代数操作、手算流程、工程陷阱和现代变体串起来讲一遍适合刚接触运筹优化的学生也适合用过求解器但一直没搞懂内部逻辑的工程师。1. 一个反直觉的事实指数级算法为什么统治求解器七十年1.1 从一次面试对话说起很多教材在讲单纯形法时都会先强调它的最坏情况复杂度是指数级的。确实有研究者构造过一类特殊的“扭曲立方体”让单纯形法沿着多面体的顶点一个个走走完近乎所有顶点才能找到最优解。理论上说如果问题规模是 n单纯形法最坏可能要跑 2 的 n 次方量级的迭代。按这个说法这套算法早该被淘汰了。但现实完全不是这样。今天几乎所有主流的线性规划求解器内部的默认算法体系里仍然保留着完整的单纯形法家族甚至很多时候它还是最可靠的那一个。成千上万变量和约束的大规模问题单纯形法通常在几十次到几百次迭代内就能收敛。为什么一个理论上“差得离谱”的算法在实际中却如此稳定高效这是我觉得每个用 LP 的人都该停下来想清楚的问题。答案是两个层面一是单纯形法在实际中表现的平均复杂度远好于最坏情况二是它在工程里有两个内点法很难替代的优点——热启动能力和对偶信息的可解释性。理解这一点比背下一堆表格变换步骤有用得多。1.2 真正需要搞懂的三个“为什么”如果给“理解 LP Simplex”划一条最低线我认为是三个问题为什么最优解一定可以在顶点找到而不是可行域内部某个地方为什么沿着多面体的边从一个顶点滑到另一个顶点就足以逼近最优解为什么每次只换入一个变量、换出一个变量目标函数就能不断改善而且最终会在某一步停下来。这三个问题分别对应了单纯形法的几何基础、搜索策略和收敛性判断。后面我会按这个逻辑依次展开。先说结论单纯形法本质上就是一个在凸多面体顶点之间跳跃的搜索算法只是它把几何上的“跳跃”翻译成了代数上的“转轴”pivot操作。把翻译规则搞明白单纯形表就不再是一堆数字的机械变换而是一张写着“现在在哪、下一步去哪、为什么能去”的地图。2. 把LP翻译成几何为什么最优解一定在顶点2.1 标准形式里藏着几何要理解单纯形法先要接受线性规划的标准形式min c^T x s.t. Ax b x 0其中 A 是 m 行 n 列的系数矩阵x 是 n 维变量向量b 是 m 维右侧常数向量。如果原始问题是不等式约束我们可以通过加松弛变量、剩余变量等方式转成等式如果是最大化问题把目标函数取负号就变成最小化。所有这些变形都不是为了刁难人而是让几何结构统一起来。从几何上看Ax b描述的是 n 维空间中的一个 m 维超平面如果 A 行满秩。再加上x 0就等价于在这个超平面上再切一刀只保留非负象限里的部分。最终得到的可行域是一个凸多面体专业上叫 polytope。它可能有界也可能无界可能为空也可能包含无穷多个点。目标函数c^T x是一个线性函数它的等值线是一组平行的超平面。我们要找的就是在所有可行点里让目标函数值最小的那个点。因为目标函数是线性的所以最优解不会出现在某个“圆润”的内部——多面体的边界是平的线性函数在平边界上的变化也是单调的。这个直觉接下来会引导我们走向顶点。2.2 顶点与“基”其实是同一件事多面体的顶点在不退化的正常情况下对应着“恰好有 n 个约束同时取等号”的点。在我们的标准形式里Ax b已经固定了 m 个等式约束还剩下n - m个自由度。要让一个点成为顶点就需要再让n - m个变量恰好等于 0也就是让n - m个非负约束也取等号。这样剩下的 m 个变量就由等式唯一确定了。代数上解释更干净从 A 的 n 列里选出 m 个线性无关的列组成一个可逆矩阵 B对应变量记为 x_B叫基变量其余列组成 N对应变量 x_N叫非基变量。我们强制 x_N 0于是B x_B b解得 x_B B^{-1} b。如果这个解满足 x_B 0那它就是一个基础可行解BFS。几何上每个 BFS 就是一个顶点。我之前带团队做一个产能分配项目时刚开始大家直接用现成求解器跑遇到不可行报告就一头雾水后来我把这个“顶点等于基”的对应关系讲清楚队友们才意识到求解器报infeasible本质是“没有任何一组基变量组合能同时满足等式和非负约束”。这个认知对建模调试非常关键。2.3 高维空间里的“邻居”怎么找顶点知道了搜索策略也顺理成章从当前顶点出发沿着多面体的某条边移动到相邻顶点目标函数值要不减最小化时。所谓“边”在代数上的刻画特别简洁——只让一个非基变量从 0 开始增大其余非基变量保持为 0。因为只有一个自由度在变化这条路径自然是直线而且它一定位于多面体的一条棱上。这一步就是单纯形法名字里“单纯形”的几何来源之一把复杂多面体的搜索拆解成一连串沿着低维单纯形结构顶点、边、面的移动。具体执行时我们需要回答两个问题让哪个非基变量变大入基以及变大多少才不越界出基。这两个问题就是转轴操作的全部内容。3. 核心机制拆解一次转轴到底在决策什么3.1 入基检验数说了算假设当前基变量为 x_B非基变量 x_N 0。用 B^{-1} 左乘Ax b可以把等式改写成x_B B^{-1}b - B^{-1}N x_N如果把目标函数也用非基变量表示就得到z c_B^T B^{-1}b (c_N^T - c_B^T B^{-1}N) x_N括号里那一项对第 j 个非基变量来说就是传说中的检验数reduced costr_j c_j - c_B^T B^{-1} A_j它的含义非常直观在其他非基变量保持为 0 的情况下如果让 x_j 增大一个单位目标函数 z 会变化 r_j。最小化问题中如果 r_j 0说明增大 x_j 能让目标函数变小值得尝试如果所有 r_j 0说明无论放开哪个非基变量都无法再改善目标函数当前 BFS 就是最优解。选择哪个变量入基最常见的策略是选检验数最负的那个对最小化也就是“最值得尝试的方向”。不过工程上这里有个我后来踩过的坑每次都严格选最负检验数在大规模稀疏问题上不一定最快。很多求解器会采用部分定价partial pricing一次只扫描一部分列选其中最好的一个。纯粹按最负检验数走属于教材策略实际实现往往在“搜索质量”和“计算开销”之间做平衡。3.2 出基最小比值规则是保护合法性的选定入基变量后它从 0 开始增大基变量也会跟着变化。我们需要保证所有基变量仍然满足 x_B 0否则就走出了可行域。这一步对应着单纯形表里的最小比值检验ratio test。把入基方向记为 d B^{-1}A_j那么基变量随 x_j 的变化可以写成x_B - d * θ 0其中 θ 是 x_j 增大的步长。只有当 d 的某个分量为正时对应的基变量才会随着 x_j 增大而减小所以步长上限是θ_max min { (x_B)_i / d_i | d_i 0 }比值最小的那一行对应的基变量会最先降到 0因此它要被换出基。如果所有 d_i 0意味着 x_j 可以无限增大而不破坏可行性那么问题是无界的——这在建模时通常意味着漏掉了某个关键约束。最小比值规则听起来只是保护可行性的机械步骤但它有一个容易被忽略的作用它保证从当前顶点移动到的是相邻顶点而不是直接穿过可行域内部。这个“相邻性”是单纯形法搜索效率的基础。如果跳过 ratio test 随便换出可能会得到一个不可行解那就得额外做修复得不偿失。3.3 转轴的本质是一次坐标重写确定入基和出基变量后剩下的就是代数操作把入基变量变成新的基变量把出基变量变成 0。这个操作在表格上表现为高斯消元——以入基变量列和出基变量行的交叉点为主元把主元位置化为 1把该列其他位置化为 0。我更喜欢把转轴理解为一次坐标重写整个线性方程组没变我们只是换了一组变量作为“当前坐标系”原来用 s1、s2 等基变量表达的量现在改用新的基变量来表达。这一步计算量不大但它把几何上的“移动到相邻顶点”完整地翻译成了代数上的“换一组基”。很多同学学到这里会陷入机械消元的泥潭忘了观察目标函数行的变化。我建议每次转轴后都做两件事一是确认 RHS 全部非负二是看目标函数值是否严格朝期望方向变化。如果 RHS 出现负数说明转轴出错了如果目标函数值完全不变化那可能是遇到退化这个后面会专门讲。4. 手算一个完整例子跟着单纯形表走一遍4.1 构造模型与初始表空谈不好消化我们用一个简单但完整的例子手算一遍。假设某生产计划问题要最大化两种产品的收益max z 5x1 4x2 s.t. x1 2x2 6 3x1 2x2 12 x1, x2 0引入松弛变量 s1、s2 转成标准形式max z 5x1 4x2 s.t. x1 2x2 s1 6 3x1 2x2 s2 12 x1, x2, s1, s2 0初始基选择 s1、s2因为它们的系数列正好是单位矩阵。初始基础可行解就是 s1 6s2 12x1 x2 0目标值为 0。表格如下最后一行的数字是各变量的检验数这里是 z 相对于每个变量的方向基x1x2s1s2RHSs112106s2320112z54000这里我用的是最大化问题的直接检验数形式z 5x1 4x2检验数越正越值得入基。4.2 第一次转轴从 (0,0) 滑到 (4,0)当前 z 行的检验数里x1 为 5x2 为 4选最大的 5所以 x1 入基。接着做最小比值检验。x1 列的数字是 1 和 3分别对应 s1 行和 s2 行s1 行比值6 / 1 6s2 行比值12 / 3 4最小比值在第 s2 行所以 s2 出基。主元位置是 s2 行、x1 列交叉处值为 3。执行高斯消元把主元化为 1x1 列在 s1 行和 z 行的位置化为 0。得到新表基x1x2s1s2RHSs104/31-1/32x112/301/34z02/30-5/320这个表的几何含义是我们从原点 (0,0) 沿着第二约束的边界滑到了顶点 (4,0)。此时 x1 4x2 0z 20s1 2s2 0。注意 z 行里 x2 的检验数是 2/3仍然是正数说明目标还有改善空间。4.3 第二次转轴从 (4,0) 滑到 (3,1.5)现在唯一正检验数是 x2因此 x2 入基。x2 列的数字是 4/3 和 2/3s1 行比值2 / (4/3) 1.5x1 行比值4 / (2/3) 6最小比值是 1.5所以 s1 出基主元位置是 s1 行、x2 列交叉处值为 4/3。消元后得到基x1x2s1s2RHSx2013/4-1/43/2x110-1/21/23z00-1/2-3/221此时 x1 3x2 1.5z 21s1 s2 0。z 行的两个检验数都是负数没有可以改善的方向了所以这就是最优解。几何上这个点正好是第一约束x1 2x2 6与第二约束3x1 2x2 12的交点两条资源都被完全利用。4.4 从最终表读出更多信息单纯形表的最后一行不只是判断最优性的工具它还携带对偶信息。在这个例子里s1 和 s2 的检验数分别是 -1/2 和 -3/2。对最大化问题检验数的相反数就是对应约束的影子价格s1 代表第一种资源的影子价格为 1/2s2 代表第二种资源的影子价格为 3/2。它的实际意义是如果第一种资源的右侧常数 b1 增加一个单位目标函数大约会增加 1/2第二种资源增加一个单位目标函数大约增加 3/2。在做资源采购或容量规划时这个数字比单纯的最优解值更有决策价值。可惜很多工程团队跑完求解器只看最优解白白浪费了这张表。我自己的习惯是任何 LP 模型跑完之后都会顺手把最终的检验数导出看一眼。很多时候模型里的“瓶颈约束”和“剩余富余约束”一眼就能分辨出来这比反复调模型参数要快得多。5. 工程里的暗坑退化、循环与数值稳定性5.1 退化是什么为什么让目标函数“卡住”单纯形法有一个理想化假设每一步转轴目标函数值都会严格改善。但实际模型中经常出现退化——某个基变量恰好等于 0。这时候转轴前后顶点在几何上可能没有变化目标函数值也不变算法就像原地踏步。退化通常发生在顶点处有“多余”的约束同时生效也就是说一个顶点实际上可以由超过 n 个约束共同确定。代数上表现为某个 x_B 分量为 0。如果连续多步退化的转轴单纯形法可能绕着同一个顶点打转极端情况下会出现循环永远达不到最优。我实际遇到过的最典型情况是重资产调度模型里大量“刚好持平”的约束。比如某个班组的工时刚好排满某些通道的流量刚好达到上限这些数值让基变量变成 0单纯形法的迭代进程就慢下来了。不是模型错了但确实需要算法层面的处理来绕过。5.2 循环与反循环策略针对循环问题理论上有几种经典方案。最朴素的是 Blands 规则入基和出基都选择下标最小的候选变量。它能从理论上保证不循环实际算起来数字也不复杂。但工程实现很少直接用这个规则因为“最小下标”完全忽略了问题本身的数值结构会导致迭代次数偏多。更实用的做法是扰动法把所有约束右侧常数轻微随机扰动一下破坏退化顶点的“病态切面”跑完之后再把扰动去掉从相邻顶点继续。另一种是字典序规则本质上给基向量定义一个字典序作为比较依据保证每次转轴后某些指标严格改善。现代求解器为了性能通常会组合使用默认放宽零阈值判断配合动态扰动只有在检测到长时间无改善时才切换反循环策略。这个细节提醒我们单纯形法不是只看教材上那套标准流程就能直接上生产环境的。真正稳定的求解器必然包含一堆关于“什么时候该打破常规”的启发式判断。5.3 没有初始基怎么办大M法与两阶段法标准形式里如果 A 里没有现成的单位矩阵我们就没有初始基不能直接开始转轴。经典方案有两个大M法和两阶段法。大M法的思路是给每个没有单位列的约束加一个非负人工变量并在目标函数里对人工变量施加重大的惩罚系数 M最大化时减去 M 乘以人工变量之和。只要 M 足够大最优解就会把人工变量压到 0从而还原原问题。听起来很简单但 M 太难选了太小会被原始目标“淹没”导致最优解里混入人工变量太大会在数值计算中引入巨大的条件数让基矩阵的病态程度加剧。我后来基本不在自己的代码里用大M法只在教学示例里演示。两阶段法更稳健第一阶段先忽略原目标函数只求人工变量之和的最小值。如果最小值是 0说明可以找到原问题的可行基于是把人工变量删掉带着第一阶段得到的基进入第二阶段恢复原目标函数继续迭代。工程求解器内部往往还会先做预求解presolve通过约束化简直接消掉一部分变量很多时候根本不需要显式的人工变量。5.4 主元选择里的数值洁癖单纯形法核心操作是高斯消元每一步都会产生浮点误差。如果一个主元非常接近 0消元时会把行里的数值放大好几个量级导致精度灾难。这正是为什么工程实现里不能强行把所有“小于某个阈值”的数字全当真值处理。常见的保护措施有三个维度一是设置主元阈值绝对值太小的候选主元直接跳过换到同一列里其他更大的元素二是做行缩放把约束矩阵各行、各列的数值量级拉齐避免某个约束的单位是“万亿”而另一个是“0.001”三是利用稀疏矩阵技术尽量选那些消元后填充fill-in少的转轴顺序降低计算量和误差累积。我年轻的时候写过一个很糙的单纯形法 Demo用的是全稠密矩阵加完全消元跑几百行约束的小模型都能卡到无语。后来才发现问题不只在于复杂度更在于没有做任何数值保护。改用阈值选主元和混合精度处理后同一批问题瞬间顺畅了。纯数学推导里所有数都是精确的工程实现里最怕的恰恰是这些“差不多”。6. 从教材到工业求解器单纯形法家族的现代分工6.1 Revised Simplex不重算整张表原始单纯形法每次转轴都要更新整张表复杂度是 O(mn)在大规模稀疏问题上浪费严重。**修正单纯形法Revised Simplex**只维护基矩阵的逆 B^{-1}每次迭代只计算入基列B^{-1}A_j和检验数c_B^T B^{-1}A_j转轴时用高效更新的方式修正 B^{-1}而不是重算整表。这样做有两个好处一是大幅降低单次迭代的计算量尤其适合变量多、约束中等但矩阵稀疏的问题二是数值上更稳定因为可以定期从头计算 B^{-1} 来校正累积误差。今天几乎所有工业级求解器实现的单纯形法都是修正版本而不是教材里那张大表。理解这一点对我个人的帮助是当我自己写实验性求解器时从“整表消元”切到“维护 B^{-1}”是一次质的飞跃。如果只是学习和验证思路用整表法没问题一旦要处理超过几百行约束的模型必须走修正单纯形法的路子。6.2 Dual Simplex为什么现代求解器青睐它单纯形法还有一个重要变体叫对偶单纯形法它求解的不是原问题而是对偶问题。单纯形法原版要求从“原可行”的解出发而对偶单纯形法从“对偶可行”的解出发允许原变量暂时不可行在迭代中逐步恢复可行性。为什么工业界如此依赖对偶单纯形法最重要的原因是它和整数规划的割平面法、分支定界法配合极好。在每个分支节点父节点的最优解常常只是轻微违反新加的约束对偶单纯形法可以借着这个“几乎可行”的初始点做热启动快速收敛。这一点原版单纯形法和内点法都很难做到。我记得在一个项目排产模型里我们反复往模型里添加“防同时开工”之类的逻辑约束每次重新求解求解器用的几乎都是对偶单纯形法的热启动路径求解时间从秒级压到了毫秒级。如果当初只用冷启动的算法整棵分支树的节点数会让人崩溃。6.3 单纯形还是内点实际项目怎么挑到了真实项目里最常被问的问题是我到底用单纯形法还是内点法没有绝对答案但我给一个根据问题形态挑选的方向。先看几点经验约束数量不算大、但是变量很多同时需要反复求解修改后的模型优先考虑单纯形法尤其是对偶单纯形法因为它能热启动问题和网络流、多商品流、电力调度这种带有强烈稀疏结构的大规模模型通常内点法更稳迭代次数少且不容易受退化影响如果只求一次静态解且矩阵规模巨大、机器内存充足内点法往往更快如果求解时间很敏感且后续会加约束单纯形法有机会从上一个解继续这是内点法最不擅长的地方。一个实际的工程判断流程是先用预求解把模型大幅化简然后根据稀疏程度和热启动需求选算法。很多商业求解器默认会同时开启多条路径实验谁先跑完就用谁的结果。我在做大规模多周期调度时经常看到单纯形法和内点法在同一个模型上表现完全相反——内点法冷启动飞快单纯形法热启动飞快。关键不是迷信某个算法而是让工具按场景去选。最后说一个我在实际项目里反复体会到的点单纯形法带给你的不只是“一个最优解”它每一次转轴都在告诉你当前解在多面体上的位置。当你需要向业务方解释为什么某个原料库存“刚好卡死”产量时检验数和对偶信息是最有力的工具。这也是为什么即使内点法在纯计算速度上频频占优单纯形法家族依然牢牢占据工业求解器核心位置的原因。理解它不是理解一个过时算法而是理解了整个优化建模世界里最常用的一套“空间语言”。
