为何SciPy的minimize函数无法运行而linprog可以？

问题背景 #

• 基于含3787个节点、23292条弧的图结构建模，原使用SciPy的linprog（highs方法）求解线性规划，可快速完成计算；
• 因需添加仅能以函数形式表达的新约束，改用minimize通用优化器，但未添加新约束就触发内存溢出错误，无法分配38.3GiB内存存储50多亿元素的数组。

复现代码 #

import uuid
from random import Random

import numpy as np
from scipy.optimize import linprog, minimize

rand = Random()

nodes = [uuid.uuid4().hex for _ in range(3787)]
temp = [tuple(rand.sample(range(1, 3787), 2)) for _ in range(23292)]
arcs = [(nodes[arc[0]], nodes[arc[1]]) for arc in temp]

pair_utility = {arc: -rand.random() for arc in arcs}

x0 = [1 for _ in enumerate(arcs)]

# Objective coefficients
coefficients = np.array([pair_utility[arc] for arc in arcs])
f = lambda x: coefficients.dot(x)

# Inequality constraints (Ax <= b)
A_incoming = np.zeros((len(nodes), len(arcs)))
A_outgoing = np.zeros((len(nodes), len(arcs)))

for i, node in enumerate(nodes):
    for j, (trip1_id, trip2_id) in enumerate(arcs):
        if node == trip1_id:
            A_outgoing[i, j] = 1
        elif node == trip2_id:
            A_incoming[i, j] = 1

A = np.vstack((A_incoming, A_outgoing))
b = np.ones(2 * len(nodes))
cons = ({"type": "ineq", "fun": lambda x: b - A.dot(x)},)

# Bounds for variables
x_bounds = [(0, 1) for _ in arcs]

# Solve the linear programming problem
result1 = linprog(
    coefficients,
    A_ub=A,
    b_ub=b,
    bounds=x_bounds,
    method="highs",
    options={"disp": True},
)
result2 = minimize(
    f,
    x0,
    constraints=cons,
    bounds=x_bounds,
    options={"disp": True},
)

错误信息 #

Unable to allocate 38.3 GiB for an array with shape (5141509334,) and data type float64
    result2 = minimize(
              ^^^^^^^^^
numpy.core._exceptions._ArrayMemoryError: Unable to allocate 38.3 GiB for an array with shape (5141509334,) and data type float64

错误原因 #

1. 默认方法的自动微分冗余：minimize默认使用SLSQP方法，会自动计算约束函数的雅可比矩阵。你的约束函数b - A.dot(x)的雅可比是-A，但A是稠密矩阵（尺寸7574×23292），自动微分会将其展开为一维数组（共5141509334个元素），直接耗尽内存。
2. 未利用矩阵稀疏性：约束矩阵A本质是稀疏的（每个节点仅对应少量弧，大部分元素为0），但用稠密数组存储既浪费内存，又导致自动微分生成超大稠密雅可比。
3. 初始点不合理：全1的初始点离最优解较远，可能额外增加计算量。

修正方案 #

1. 改用稀疏矩阵存储约束 #

用scipy.sparse构建稀疏约束矩阵，大幅减少内存占用：

from scipy.sparse import lil_matrix, csr_matrix

# 构建稀疏约束矩阵
A_incoming = lil_matrix((len(nodes), len(arcs)))
A_outgoing = lil_matrix((len(nodes), len(arcs)))

for j, (trip1_id, trip2_id) in enumerate(arcs):
    # 找到节点对应的索引
    i_out = nodes.index(trip1_id)
    i_in = nodes.index(trip2_id)
    A_outgoing[i_out, j] = 1
    A_incoming[i_in, j] = 1

# 转换为CSR格式，提升dot运算效率
A = csr_matrix(np.vstack((A_incoming.toarray(), A_outgoing.toarray())))
b = np.ones(2 * len(nodes))

2. 手动提供约束的雅可比矩阵 #

避免自动微分生成超大矩阵，手动指定约束的雅可比：

# 约束函数与雅可比
def constraint_fun(x):
    return b - A.dot(x)

def constraint_jac(x):
    # 约束函数的雅可比是 -A
    return -A

cons = ({
    "type": "ineq",
    "fun": constraint_fun,
    "jac": constraint_jac
},)

3. 优化初始点 #

用linprog的最优解作为minimize的初始点，加速收敛：

x0 = result1.x

4. 指定合适的优化方法 #

明确指定method='SLSQP'，并配置合理迭代参数：

result2 = minimize(
    f,
    x0,
    constraints=cons,
    bounds=x_bounds,
    method="SLSQP",
    options={"disp": True, "maxiter": 1000}
)

后续建议 #

• 优先保留LP求解器：如果新约束能转化为线性形式，继续用linprog，它专门针对线性规划优化，效率远高于通用优化器。
• 利用凸优化库：若新约束是非线性凸约束，可使用cvxpy，它能自动处理稀疏性并选择合适的求解器，无需手动管理雅可比矩阵。
• 避免稠密矩阵：所有大规模矩阵都用scipy.sparse存储，减少内存开销。
• 尝试无导数方法：若无法提供雅可比，可尝试COBYLA方法，但收敛速度会慢很多，适合对精度要求不高的场景。

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