如何用SymPy快速求解超立方体约束下多项式的KKT点？

我正在对n变量低次（如3次）多项式执行梯度下降（GD）实验，定义域受超立方体[-1,1]^n约束，希望将GD的终止点与SymPy计算得到的实际KKT点进行对比。例如，有一个含3个变量a、b、c的SymPy表达式：

0.2*a**3 - 0.1*a**2*b - 0.9*a**2*c + 0.5*a**2 - 0.1*a*b**2 - 0.4*a*b*c + 0.5*a*b - 0.6*a*c**2 - 0.5*a*c - 0.6*a - 0.1*b**3 + 0.4*b**2*c - 0.4*b**2 + 0.7*b*c**2 - 0.4*b*c + 0.9*b + 0.09*c**3 - 0.6*c**2 + 0.4*c + 0.22

约束条件为：

-1 <= a <= 1
-1 <= b <= 1
-1 <= c <= 1

请问在上述约束下，求解该表达式所有KKT点的最优、最快方法是什么？（理想情况是能把搜索范围限定在GD终止点附近来节省计算时间，但不确定怎么严谨实现）

核心求解思路 #

针对超立方体约束下的多项式优化问题，KKT点可拆分为内部临界点和边界/顶点临界点两类处理，结合SymPy的符号运算与数值验证能高效完成：

1. 求解内部KKT点（变量均在(-1,1)区间内） #

内部点的KKT条件等价于梯度为0，直接用SymPy做符号求解：

• 步骤：
  1. 定义符号变量并输入目标表达式；
  2. 计算目标函数对每个变量的偏导，得到梯度方程组；
  3. 解方程组，筛选出所有变量值落在(-1,1)内的解。
• 示例代码片段：

import sympy as sp

a, b, c = sp.symbols('a b c')
f = 0.2*a**3 - 0.1*a**2*b - 0.9*a**2*c + 0.5*a**2 - 0.1*a*b**2 - 0.4*a*b*c + 0.5*a*b - 0.6*a*c**2 - 0.5*a*c - 0.6*a - 0.1*b**3 + 0.4*b**2*c - 0.4*b**2 + 0.7*b*c**2 - 0.4*b*c + 0.9*b + 0.09*c**3 - 0.6*c**2 + 0.4*c + 0.22

# 计算梯度
grad_f = [sp.diff(f, var) for var in (a, b, c)]
# 解梯度为0的方程组
inner_solutions = sp.solve(grad_f, (a, b, c), dict=True)
# 筛选满足区间约束的解
valid_inner = []
for sol in inner_solutions:
    if all(-1 < sp.N(sol[var]) < 1 for var in (a, b, c)):
        valid_inner.append({k: sp.N(v) for k, v in sol.items()})

2. 求解边界上的KKT点（至少一个变量取±1） #

超立方体的边界需分三类处理，遍历所有边界组合后重复「内部+验证」逻辑：

• 单变量取边界：将该变量固定为1或-1，转化为2变量优化问题，求剩余变量的梯度零点并验证区间约束；
• 双变量取边界：固定两个变量为±1，转化为单变量函数，求导找临界点并验证；
• 三变量取边界（顶点）：直接代入原表达式，验证梯度是否满足互补松弛条件（比如极小值问题中，变量取1时对应梯度分量需≥0）。
• 单变量边界处理示例：

# 遍历所有单变量边界情况
boundary_cases = [(a, 1), (a, -1), (b, 1), (b, -1), (c, 1), (c, -1)]
valid_boundary = []

for var, val in boundary_cases:
    # 代入边界值，得到降维后的目标函数
    f_bound = f.subs(var, val)
    # 剩余变量列表
    remaining_vars = [v for v in (a, b, c) if v != var]
    # 计算剩余变量的梯度
    grad_bound = [sp.diff(f_bound, v) for v in remaining_vars]
    # 解梯度为0的方程组
    boundary_solutions = sp.solve(grad_bound, remaining_vars, dict=True)
    
    # 筛选剩余变量在(-1,1)内的解，并补回边界变量值
    for sol in boundary_solutions:
        sol[var] = val
        if all(-1 < sp.N(sol[v]) < 1 for v in remaining_vars):
            # 验证互补松弛条件（以极小值为例：若变量取1，对应梯度分量≥0；取-1则≤0）
            grad_val = sp.N(sp.diff(f, var).subs(sol))
            if (val == 1 and grad_val >= -1e-6) or (val == -1 and grad_val <= 1e-6):
                valid_boundary.append({k: sp.N(v) for k, v in sol.items()})

3. 限定GD终止点附近的搜索（严谨实现） #

若要缩小搜索范围，可基于GD终止点(a0, b0, c0)做以下处理：

• 内部点搜索：在解梯度方程组时，加入不等式约束a0-ε ≤ a ≤ a0+ε（ε取0.3~0.5，根据GD收敛精度调整），用SymPy的solve结合不等式筛选，或转用scipy.optimize.root在局部数值求解；
• 边界点搜索：先判断GD终止点靠近哪个边界（比如a接近1），优先处理该边界的情况，再验证其他边界；
• 关键验证：无论局部搜索得到什么解，必须验证是否满足KKT条件——内部点梯度为0，边界点满足互补松弛。

4. 效率优化建议 #

• 符号解转数值验证：SymPy的符号解可能包含复数或不满足区间的解，用sp.N()转数值后快速筛选；
• 混合符号+数值方法：对复杂的边界方程组，用SymPy推导梯度表达式，再用Scipy数值求解器在局部搜索，比纯符号求解更快；
• 预过滤顶点：直接计算所有8个顶点的函数值和梯度，快速排除不满足互补松弛条件的顶点。

内容的提问来源于stack exchange，提问作者ufghd34