Python scipy.optimize求解半正定约束优化出错，如何解决？

问题背景 #

给定4维随机矩阵C，需优化2维对称矩阵A₀和B₀，满足以下约束：

• I+A₀、I-A₀、I+B₀、I-B₀均为半正定矩阵（I为2维单位矩阵）
  目标函数为最大化$\text{Tr}((A₀ \otimes B₀)C)$，其中$\otimes$为Kronecker积。

对应的求解代码如下：

import numpy as np
from scipy.optimize import minimize
np.set_printoptions(precision = 8, suppress=True)
def objective_function(x, C):
    A0 = x[:4].reshape((2, 2))  
    B0 = x[4:8].reshape((2, 2))
    A0 = A0 + A0.T
    B0 = B0 + B0.T
    D00 = np.kron(A0,B0)
    return -np.trace(D00 @ C)
C = np.random.rand(4, 4)

A0_initial = np.zeros((2,2))
B0_initial = np.zeros((2,2))
initial_guess = np.hstack([A0_initial.flatten(), B0_initial.flatten()])

def constraint_func1(x):
    A0 = x[:4].reshape((2, 2))
    A0 = A0 + A0.T
    return  min(np.linalg.eigvals(np.eye(2) - A0))
def constraint_func2(x):
    A0 = x[:4].reshape((2, 2))
    A0 = A0 + A0.T
    return  min(np.linalg.eigvals(np.eye(2) + A0))
def constraint_func3(x):
    B0 = x[4:8].reshape((2, 2))
    B0 = B0 + B0.T
    return  min(np.linalg.eigvals(np.eye(2) - B0))
def constraint_func4(x):
    B0 = x[4:8].reshape((2, 2))
    B0 = B0 + B0.T
    return  min(np.linalg.eigvals(np.eye(2) + B0))
constraint1 = {'type': 'ineq', 'fun': constraint_func1}
constraint2 = {'type': 'ineq', 'fun': constraint_func2}
constraint3 = {'type': 'ineq', 'fun': constraint_func3}
constraint4 = {'type': 'ineq', 'fun': constraint_func4}
result = minimize(objective_function, initial_guess, args=(C,), constraints=[constraint1,constraint2,constraint3,constraint4])
A0_optimal = result.x[:4].reshape((2, 2))
B0_optimal = result.x[4:8].reshape((2, 2))
A0_optimal = A0_optimal + A0_optimal.T
B0_optimal = B0_optimal + B0_optimal.T
print("C",C)
print(A0_optimal)
print(B0_optimal)
print(-result.fun)

遇到的问题 #

运行上述代码后，结果始终为A₀=B₀=零矩阵，显然不符合预期，需解答以下问题：

1. 导致该现象的原因是什么？
2. scipy能否处理这类基于最小特征值的半正定约束？
3. 若scipy不行，Python中还有哪些工具可处理矩阵半正定约束？（已尝试cvxpy，但因Kronecker积导致非凸问题不支持）

问题原因分析 #

1. 局部最优陷阱：初始猜测设为全零矩阵，而scipy.optimize.minimize默认的L-BFGS-B等算法是局部优化器。零矩阵满足所有约束，且目标函数在该点的梯度期望为零（C为随机矩阵时），导致优化器无法跳出局部最优。
2. 约束函数非光滑性：用最小特征值作为约束函数，特征值函数在矩阵发生特征值交叉时不可导，干扰依赖梯度的优化器收敛，难以突破局部最优。
3. 变量冗余：代码将A₀和B₀的4个元素全部作为优化变量，之后再强制对称化，引入了冗余变量（对称矩阵仅需3个独立元素），增加优化维度，降低效率。

scipy对这类约束的支持情况 #

scipy的minimize可以处理基于最小特征值的半正定约束，但存在明显局限性：

• 这类约束属于非光滑约束，依赖梯度的优化器处理时易出现收敛问题，无法保证全局最优。
• 对于本题含Kronecker积的非凸目标函数，scipy的局部优化器仅能找到局部最优，无法确保全局最优解。

替代工具与解决方案 #

1. 现有scipy代码的改进方案 #

调整代码以突破局部最优：

• 减少变量冗余：直接用对称矩阵的独立元素作为优化变量，将总变量数从8降至6。
• 更换初始猜测：采用满足约束的非零初始点，如A₀=0.5I、B₀=0.5I。
• 使用全局优化器：改用scipy.optimize.differential_evolution，无需梯度，可在全局范围搜索最优解，适配非凸问题。

改进后的示例代码：

import numpy as np
from scipy.optimize import differential_evolution

np.set_printoptions(precision=8, suppress=True)

def objective_function(x, C):
    # A0为对称矩阵：[a, b, c] -> [[a,b],[b,c]]
    A0 = np.array([[x[0], x[1]], [x[1], x[2]]])
    # B0为对称矩阵：[d, e, f] -> [[d,e],[e,f]]
    B0 = np.array([[x[3], x[4]], [x[4], x[5]]])
    D00 = np.kron(A0, B0)
    return -np.trace(D00 @ C)

def constraint_func1(x):
    A0 = np.array([[x[0], x[1]], [x[1], x[2]]])
    return min(np.linalg.eigvals(np.eye(2) - A0))

def constraint_func2(x):
    A0 = np.array([[x[0], x[1]], [x[1], x[2]]])
    return min(np.linalg.eigvals(np.eye(2) + A0))

def constraint_func3(x):
    B0 = np.array([[x[3], x[4]], [x[4], x[5]]])
    return min(np.linalg.eigvals(np.eye(2) - B0))

def constraint_func4(x):
    B0 = np.array([[x[3], x[4]], [x[4], x[5]]])
    return min(np.linalg.eigvals(np.eye(2) + B0))

C = np.random.rand(4, 4)
# 变量边界：根据约束，对称矩阵元素范围设为[-1,1]
bounds = [(-1, 1)] * 6
# 全局优化器求解
result = differential_evolution(
    objective_function, bounds, args=(C,),
    constraints=[
        {'type': 'ineq', 'fun': constraint_func1},
        {'type': 'ineq', 'fun': constraint_func2},
        {'type': 'ineq', 'fun': constraint_func3},
        {'type': 'ineq', 'fun': constraint_func4}
    ]
)

A0_optimal = np.array([[result.x[0], result.x[1]], [result.x[1], result.x[2]]])
B0_optimal = np.array([[result.x[3], result.x[4]], [result.x[4], result.x[5]]])

print("C:\n", C)
print("最优A0:\n", A0_optimal)
print("最优B0:\n", B0_optimal)
print("最大化的目标值:\n", -result.fun)

2. 专业优化工具 #

若全局优化器仍无法满足需求，可尝试以下工具：

• Pyomo：支持自定义约束与目标函数，可调用IPOPT等外部非线性规划求解器处理非凸问题。
• CasADi：专注于最优控制与非线性优化，支持自动微分，能高效处理含矩阵运算的非凸问题，可搭配IPOPT、SNOPT等求解器。
• Gurobi/Cplex（商业工具）：若有授权，这类商业求解器对非凸二次/非线性优化的支持更完善，可处理半正定约束与Kronecker积相关的目标函数。

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