Mpmath是否有quadprog或LSEI的等效实现？

Great question! 我正好在做任意精度优化的项目时遇到过完全一样的需求——用mpmath实现带线性等式/不等式约束的最小二乘/二次规划，而不是用常规的float精度工具（比如quadprog、SLSQP）。下面给你梳理一下可行的思路：

现状说明 #

首先得明确：目前没有官方维护的、直接适配mpmath的LSEI/SLSQP/quadprog类库，但因为mpmath支持任意精度浮点数和完整的线性代数工具，我们可以自己实现核心算法，或者修改现有开源代码适配mpmath。

可行方案 #

1. 手动实现适配mpmath的Active Set方法（推荐） #

quadprog的核心就是Active Set算法，逻辑清晰，很容易改成mpmath版本。核心思路是迭代识别"活动约束"（即等式约束+当前起作用的不等式约束），然后求解对应的KKT系统，直到收敛。

核心步骤：

• 初始化一个满足所有约束的可行点（可以先通过mpmath求解等式约束得到初始点，再调整满足不等式）
• 每次迭代中，构造包含活动约束的KKT矩阵，用mpmath.lu_solve求解搜索方向和拉格朗日乘子
• 根据拉格朗日乘子判断是否需要移除活动约束，或者添加新的约束
• 用线搜索确定步长，更新解，直到满足mpmath级别的收敛条件（比如mp.norm(search_direction) < mp.mpf('1e-40')）

简化代码示例（带等式约束）：

import mpmath as mp

# 设置任意精度，比如50位小数
mp.mp.dps = 50

def qp_active_set(H, c, A_eq, b_eq):
    """
    求解带等式约束的二次规划：min 0.5*x^T H x + c^T x
    约束：A_eq x = b_eq
    """
    # 初始化可行点：先求解等式约束
    x = mp.lu_solve(A_eq, b_eq)
    n = len(x)
    m_eq = len(b_eq)
    
    while True:
        # 构造KKT系统：[H, A_eq.T; A_eq, 0] * [dx; lambda] = [-g; 0]
        g = mp.dot(H, x) + c
        KKT = mp.zeros(n + m_eq, n + m_eq)
        KKT[:n, :n] = H
        KKT[:n, n:] = A_eq.T
        KKT[n:, :n] = A_eq
        
        rhs = mp.zeros(n + m_eq, 1)
        rhs[:n] = -g
        
        # 求解KKT系统
        sol = mp.lu_solve(KKT, rhs)
        dx = sol[:n]
        
        # 收敛判断：搜索方向的范数足够小
        if mp.norm(dx) < mp.mpf('1e-45'):
            break
        
        # 更新解（等式约束下步长为1，因为是可行方向）
        x += dx
    
    return x

如果要扩展到不等式约束，只需要添加活动约束的管理逻辑——每次迭代检查非活动约束是否会被违反，以及活动约束的拉格朗日乘子是否为负（如果是，移除该约束）。

2. 基于mpmath重写SLSQP算法 #

SLSQP是序列二次规划，每一步求解一个带约束的二次子问题，同样可以把常规float版本的SLSQP（比如scipy里的实现）改成mpmath版本：

• 把所有numpy的矩阵运算替换成mpmath的对应函数（比如numpy.dot→mpmath.dot，numpy.linalg.inv→mpmath.inv）
• 用mpmath的精度判断替换常规的收敛阈值
• 保留SLSQP的BFGS近似Hessian的逻辑，只是用mpmath的浮点数计算

这种方法适合需要处理更复杂约束（比如非线性，但你的问题是线性约束，所以Active Set更高效）的场景，但实现起来比Active Set繁琐一些。

3. 用mpmath的符号+数值混合求解KKT条件（适合小规模问题） #

对于小规模的二次规划问题，可以手动推导KKT条件，然后用mpmath.nsolve求解方程组。这种方法不用写完整的优化算法，但需要手动处理互补松弛条件。

示例代码：

import mpmath as mp

mp.mp.dps = 30

# 目标函数：0.5*(x₁² + 2x₂²) + x₁ + x₂
# 约束：x₁ + x₂ ≥ 1，x₁ ≥ 0，x₂ ≥ 0
def kkt_system(x):
    x1, x2, λ1, λ2, λ3 = x
    return [
        # 梯度条件：∇f = λ₁∇c₁ - λ₂∇c₂ - λ₃∇c₃
        x1 + 1 - λ1 + λ2,
        2*x2 + 1 - λ1 + λ3,
        # 互补松弛：λ_i * c_i = 0（这里假设x₁+x₂=1是活动约束，λ₁>0）
        x1 + x2 - 1,
        λ2 * x1,
        λ3 * x2
    ]

# 初始点选可行点[1, 0, 1, 0, 0]
initial_guess = [mp.mpf(1), mp.mpf(0), mp.mpf(1), mp.mpf(0), mp.mpf(0)]
solution = mp.nsolve(kkt_system, initial_guess)

print("最优解：x₁ =", solution[0], "x₂ =", solution[1])

这种方法的缺点是：需要手动推导KKT条件，大规模问题不适用，且nsolve的收敛性严重依赖初始点。

额外建议 #

• 如果你的问题是带约束的最小二乘（即H = AᵀA，c = -Aᵀb），可以简化Active Set方法，因为H是半正定的，数值稳定性更好。
• 任意精度计算的速度比常规float慢很多，所以这种方案适合小规模问题，大规模问题可能需要权衡精度和计算效率。
• 可以参考一些开源的任意精度优化项目，但目前没有成熟的现成库，大多需要自己实现核心逻辑。

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