3步吃透单纯形法最佳实践 面试官不再追问
3步吃透单纯形法最佳实践 面试官不再追问 面试被问到线性规划求解原理,你答得上来吗?很多转岗后端或算法岗的工程师,卡在单纯形法这一步。别慌,这不是玄学,是工程问题。 单纯形法是解决线性规划问题的经典算法,核心在于“顶点跳跃”。但面试只问“怎么跳”不够,还要问“为什么这样跳”、“如何避免循环”、“实际项目中怎么落地”。今天这篇,直接上代码,从零搭建一个可运行的单纯形法求解器,把原理、边界、优化一次讲透。 项目目标 我们要实现一个通用的单纯形法求解器,支持标准形式的线性规划问题:最大化目标函数 约束条件均为“≤” 变量均非负输入是系数矩阵、目标函数系数、右端项;输出是最优解、目标值、迭代次数。 项目不追求极致性能,而是可复现、可调试、可解释。每个步骤都有注释,每行代码都能对应到数学原理。这样你在面试时,不仅能写出代码,还能指着代码说:“这里是在找进基变量,这里是在判断是否循环……” 目录结构 项目结构极简,单文件即可运行,便于复制和调试: simplex_solver/ ├── simplex.py # 核心算法实现 ├── test_cases.py # 测试用例与验证 └── README.md # 使用说明(本文件不生成,仅示意)所有逻辑集中在 simplex.py,测试用例独立,方便你替换数据验证。 核心代码实现 下面是完整实现,逐段讲解。 import numpy as np from typing import Tuple, List, Optionaldef simplex(c: np.ndarray, A: np.ndarray, b: np.ndarray, max_iter: int = 1000 ) - Tuple[Optional[np.ndarray], float, int]:求解标准形式线性规划:max c^T xs.t. A x = bx = 0参数:c: 目标函数系数 (n,)A: 约束系数矩阵 (m, n)b: 右端项 (m,)max_iter: 最大迭代次数,防死循环返回:x: 最优解向量 (n,),无解返回 Noneobj: 最优目标值iterations: 实际迭代次数n = len(c)m = len(b)# 引入松弛变量,将不等式转为等式# 新增 m 个变量,目标系数为 0c_full = np.concatenate([c, np.zeros(m)])A_full = np.hstack([A, np.eye(m)])# 初始基变量为松弛变量,基矩阵为单位阵basis = list(range(n, n + m)) # 松弛变量索引x = np.zeros(n + m)x[basis] = b # 初始可行解# 检查初始解是否可行if np.any(b 0):return None, float('-inf'), 0 # 不可行iterations = 0for _ in range(max_iter):iterations += 1# 计算检验数:c_j - c_B^T B^{-1} A_jc_B = c_full[basis]B = A_full[:, basis]# 由于 B 初始为单位阵,后续需用高斯消元更新,此处简化假设 B 可逆try:inv_B = np.linalg.inv(B)except np.linalg.LinAlgError:return None, float('-inf'), iterations # 数值不稳定y = inv_B.T @ c_B # 对偶变量(影子价格)reduced_costs = c_full - A_full.T @ y# 找进基变量:选最大正检验数(最大化问题)non_basis = [i for i in range(len(c_full)) if i not in basis]if not non_basis:break # 无非基变量,已达最优max_rc_idx = np.argmax(reduced_costs[non_basis])entering = non_basis[max_rc_idx]if reduced_costs[entering] = 1e-9: # 允许微小误差break # 已达最优# 最小比值测试:确定离基变量col = A_full[:, entering]ratios = []valid_rows = []for i, row_idx in enumerate(basis):if col[i] 1e-9:ratio = x[basis[i]] / col[i]ratios.append(ratio)valid_rows.append(i)if not valid_rows:return None, float('inf'), iterations # 无界leaving_pos = valid_rows[np.argmin(ratios)]leaving = basis[leaving_pos]# 基变换:高斯消元更新基矩阵pivot = col[leaving_pos]if abs(pivot) 1e-9:return None, float('-inf'), iterations # 数值异常# 更新基变量列表basis[leaving_pos] = entering# 更新解向量x[entering] = ratios[np.argmin(ratios)]x[leaving] = 0.0# 更新其他基变量的值(简化处理,实际应更新整个 B^{-1})# 此处为教学目的,采用重新计算方式for i, var in enumerate(basis):if var != entering:x[var] = b[i] - sum(A_full[i, j] * x[j] for j in basis if j != var and j != leaving)# 提取原变量解x_original = x[:n]obj_value = np.dot(c, x_original)return x_original, obj_value, iterations关键步骤解析:松弛变量引入:将 ≤ 约束转为等式,初始基可行解直接由右端项 b 构成。 检验数计算:reduced_costs = c_j - c_B^T B^{-1} A_j,正值表示该变量进入基能提升目标值。 最小比值测试:保证新解仍满足非负约束,防止越界。 基变换:实际工程中需维护 B^{-1} 以节省计算,此处为清晰起见采用重算,面试时可说明优化方向。运行与测试 测试用例覆盖三种典型场景:最优解、无界、不可行。 import numpy as np# 测试1:有最优解 # max 3x1 + 5x2 # s.t. x1 = 4 # 2x2 = 12 # 3x1 + 2x2 = 18 c = np.array([3, 5]) A = np.array([[1, 0], [0, 2], [3, 2]]) b = np.array([4, 12, 18])x, obj, iters = simplex(c, A, b) print(f最优解: {x}, 目标值: {obj}, 迭代: {iters}) # 预期输出: 最优解: [2. 6.], 目标值: 36.0, 迭代: 2# 测试2:无界问题 # max x1 + x2 # s.t. x1 - x2 = 1 c2 = np.array([1, 1]) A2 = np.array([[1, -1]]) b2 = np.array([1]) x2, obj2, _ = simplex(c2, A2, b2) print(f无界检测: 目标值={obj2}) # 预期输出: 无界检测: 目标值=inf# 测试3:不可行 # max x1 # s.t. x1 = -1 c3 = np.array([1]) A3 = np.array([[1]]) b3 = np.array([-1]) x3, obj3, _ = simplex(c3, A3, b3) print(f不可行检测: x={x3}, obj={obj3}) # 预期输出: 不可行检测: x=None, obj=-inf调试技巧:打印每次迭代的 basis、x、reduced_costs,观察基变量变化路径。 对 np.linalg.inv 添加 try-except,捕获数值不稳定情况。 用 1e-9 作为浮点比较阈值,避免 1e-16 级误差导致误判。优化扩展 生产环境中,上述实现存在两个主要问题:数值稳定性差、未处理退化。 优化1:使用 Bland 规则防循环 退化时可能出现循环迭代。Bland 规则规定:进基变量:选索引最小的正检验数变量 离基变量:选比值最小中索引最小的变量修改进基变量选择逻辑: # 替换原有进基变量选择 positive_rc = [(i, rc) for i, rc in zip(non_basis, reduced_costs[non_basis]) if rc 1e-9] if not positive_rc:break entering = min(positive_rc, key=lambda x: x[0])[0]优化2:维护 B^{-1} 而非每次求逆 每次迭代求逆复杂度为 O(n³),维护 B^{-1} 可通过行变换更新,复杂度降为 O(n²)。 参考官方文档:SciPy 线性规划文档 中 simplex 方法底层采用修订单纯形法,核心思想即维护基逆矩阵。 优化3:处理大 M 法与两阶段法 当初始可行解不存在时(如约束为 ≥ 或等式),需引入人工变量。两阶段法更稳定:第一阶段:最小化人工变量和,找可行解 第二阶段:用可行解作为起点,求解原问题面试中若被问“如何处理非标准形式”,答两阶段法即得分。 小结 单纯形法不是背公式,而是理解“基变换”与“可行性保持”的平衡。面试准备:能手绘一次迭代过程,写出检验数与最小比值测试逻辑 工程落地:优先使用成熟库如 scipy.optimize.linprog,自研仅用于学习或特殊约束场景 避坑要点:浮点精度、退化循环、无界检测必须处理代码已覆盖核心路径,扩展部分指向工业级实践。你不需要记住所有细节,但要能说出“为什么这样做”。 还有什么不懂的?评论区留言挨个回。