使用scipy.optimize SLSQP遇'Singular matrix C in LSQ subproblem'错误求助

Hey there, let's break down why you're hitting that error and walk through concrete fixes. First, let's unpack what that error means: SLSQP solves constrained optimization by formulating a linear least-squares (LSQ) subproblem at each iteration. A singular matrix here tells us the constraint system you've defined is linearly dependent—meaning some constraints are redundant or incorrectly defined, which breaks the solver's ability to compute a valid step.

Root Cause: Closure Trap in Constraint Definition #

The biggest issue in your code is a classic Python closure bug in build_constraint_rows. When you define row_constrain inside the loop, all constraint functions end up referencing the last iteration's vector variable, not the one from their respective loop step.

For your 2-row matrix, this means both equality constraints are actually checking the same row sum (the second row), making your constraint system rank-deficient. SLSQP can't handle this, hence the singular matrix error.

Fixes to Get SLSQP Working #

1. Fix the Closure Trap (Quickest Fix) #

Modify your constraint function definition to capture the current vector for each loop iteration using a default parameter. This breaks the closure reference to the loop variable:

from scipy.optimize import minimize
import numpy as np

def build_objective(ck, sign = -1.00):
    ck = np.concatenate(ck)
    def objective(P):
        return sign*(ck.dot(P))
    return objective

def build_constraint_rows(ck):
    ncol = ck.shape[1]
    nrow = ck.shape[0]
    constrain_dict = []
    for i in range(nrow):
        vector = np.zeros((nrow,ncol))
        vector[i, :] = 1
        vector = np.concatenate(vector)
        # Use a default parameter to capture the current vector
        def row_constrain(P, vec=vector):
            return 1 - vec.dot(P)
        constrain_dict.append({'type': 'eq', 'fun': row_constrain})
    return constrain_dict

# Initialize parameters (unchanged from your code)
c = np.array([[0. , 0. , 0. , 0., 0.], [0. , 20094.21019108, 4624.08079143, 6625.51724138, 3834.81081081]])
P_initial = np.ones(c.shape)*0.01
nrow = c.shape[0]
for i in range(nrow):
    index= np.where(c[i,] == np.max(c[i,]))[0]
    if index.shape[0] > 1:
        index = int(np.random.choice(index, size = 1))
    else:
        index = int(index)
    P_initial[i,index] = 0.96
P_initial = np.concatenate(P_initial)
b = (0.01,0.99)
bnds = (b,)*c.size

# Run optimization with fixed constraints
objective_fun = build_objective(c)
cons = build_constraint_rows(c)
res = minimize(objective_fun, P_initial, method='SLSQP', bounds=bnds, constraints=cons)
print(res)

2. Use Linear Constraint Matrix Form (More Efficient) #

SLSQP supports linear constraints directly via matrix inputs (available in scipy >=1.9.0), which avoids closure issues entirely and is more efficient for linear problems:

from scipy.optimize import minimize
import numpy as np

def build_objective(ck, sign = -1.00):
    ck = np.concatenate(ck)
    def objective(P):
        return sign*(ck.dot(P))
    return objective

# Initialize parameters (unchanged)
c = np.array([[0. , 0. , 0. , 0., 0.], [0. , 20094.21019108, 4624.08079143, 6625.51724138, 3834.81081081]])
P_initial = np.ones(c.shape)*0.01
nrow = c.shape[0]
for i in range(nrow):
    index= np.where(c[i,] == np.max(c[i,]))[0]
    if index.shape[0] > 1:
        index = int(np.random.choice(index, size = 1))
    else:
        index = int(index)
    P_initial[i,index] = 0.96
P_initial = np.concatenate(P_initial)
b = (0.01,0.99)
bnds = (b,)*c.size

# Define linear equality constraints directly as matrix A_eq and vector b_eq
A_eq = np.array([[1,1,1,1,1,0,0,0,0,0], [0,0,0,0,0,1,1,1,1,1]])
b_eq = np.array([1, 1])
cons = {'type': 'eq', 'fun': lambda x: A_eq @ x - b_eq}

# Run optimization
objective_fun = build_objective(c)
res = minimize(objective_fun, P_initial, method='SLSQP', bounds=bnds, constraints=cons)
print(res)

3. Stick with linprog (Best for Linear Problems) #

Since your problem is linear objective + linear constraints + linear bounds, linprog is purpose-built for this scenario—it's faster, more stable, and avoids the pitfalls of using a nonlinear solver like SLSQP for linear tasks. Your original linprog code works perfectly, so there's no need to switch to SLSQP unless you plan to add nonlinear constraints later.

Expected Outcome #

After fixing the constraint issue, SLSQP will successfully converge to the same optimal solution as linprog:

• For the first row (all zeros in c), any valid combination of values summing to 1 within [0.01, 0.99] is optimal (since they don't affect the dot product).
• For the second row, the optimal solution will set the largest c value (20094.21) to 0.96 (since 1 - 4*0.01 = 0.96, which fits within [0.01, 0.99]) and the rest to 0.01, maximizing the dot product.

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