The goal is to calculate an optimization problem using python scipy.optimize.
Suppose C is a given 4-dimension matrix (in the code I use a random matrix to denote that). The optimized variables are A0 and B0 which are 2-dimension symmetric matrices. The constraints are I+A0, I-A0, I+B0, I-B0 are semidefinite positive, where I is 2-dimension identity matrix. The objective function is max Tr((A0 ⊗ B0)C), where ⊗ is the Kronecker Product.

The code is as follows:

import numpy as np
from scipy.optimize import minimize
np.set_printoptions(precision = 8, suppress=True)

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)

The answer is always A0=B0= zero matrix, which is obviously wrong. Why is the case? Is the scipy able to solve the constraint like this: min(np.linalg.eigvals(np.eye(2) – A0)) > 0? If not, what tools are recommended to solve the semidefinite positve conditions for matrix? (I tried cvxpy, but it doesn’t support the kronecker product of two matrix variables because it is not convex)