Consider the following optimization problem: The optimized variable A0,B0,A1,B1 are dimension-2 real symmetric matrices. All of their eigenvalues have absolute value 1 . rho is a given dim-4 matrix. The objective function is max Tr[(A0⊗B0 +A0⊗B1+A1⊗B0 -A1⊗B1)*rho]. Here ⊗ is the Kronecker Product. The background is from the Tsirelson’s bound in quantum information.

I tried to use scipy optimize from Python to solve it, the code is:

import numpy as np
from scipy.optimize import minimize

np.set_printoptions(precision = 8, suppress=True)
paulimatrices = np.zeros((4, 2, 2))

paulimatrices[0] = np.array([[1, 0], [0, 1]]) 
paulimatrices[1] = np.array([[0, 1], [1, 0]])  
paulimatrices[2] = np.array([[0, -1], [1, 0]]) #real
paulimatrices[3] = np.array([[1, 0], [0, -1]])

noise = 1.0


def objective_function(x, epr):
    A0 = x[:4].reshape((2, 2))  
    B0 = x[4:8].reshape((2, 2))
    A1 = x[8:12].reshape((2, 2))
    B1 = x[12:16].reshape((2, 2))
    A0 = A0 + A0.T
    B0 = B0 + B0.T
    A1 = A1 + A1.T
    B1 = B1 + B1.T
    D00 = np.kron(A0,B0)
    D01 = np.kron(A0,B1)
    D10 = np.kron(A1,B0)
    D11 = np.kron(A1,B1)
    return -np.trace((D00+D01+D10-D11) @ epr)
    
rho = 0.25*(np.kron(paulimatrices[0],paulimatrices[0])+noise*np.kron(paulimatrices[1],paulimatrices[1])+noise*np.kron(paulimatrices[2],paulimatrices[2])+noise*np.kron(paulimatrices[3],paulimatrices[3]))
A0_initial = np.eye(2)
B0_initial = np.eye(2)
A1_initial = np.eye(2)
B1_initial = np.eye(2)
initial_guess = np.hstack([A0_initial.flatten(), B0_initial.flatten(),A1_initial.flatten(),B1_initial.flatten()])

def constraint_func1(x):
    A0 = x[:4].reshape((2, 2))
    A0 = A0 + A0.T
    return np.trace(np.dot(A0,A0))-2
    # return np.trace(np.linalg.matrix_power(A0, 2))- 2
def constraint_func2(x):
    A0 = x[:4].reshape((2, 2))
    A0 = A0 + A0.T
    return np.trace(np.dot(np.dot(A0,A0),np.dot(A0,A0))) -2
    # return np.trace(np.linalg.matrix_power(A0, 4))- 2
def constraint_func3(x):
    B0 = x[4:8].reshape((2, 2))
    B0 = B0 + B0.T
    return np.trace(np.dot(B0,B0))-2
    # return np.trace(np.linalg.matrix_power(B0, 2))- 2
def constraint_func4(x):
    B0 = x[4:8].reshape((2, 2))
    B0 = B0 + B0.T
    return np.trace(np.dot(np.dot(B0,B0),np.dot(B0,B0))) -2
    # return np.trace(np.linalg.matrix_power(B0, 4))- 2
def constraint_func5(x):
    A1 = x[8:12].reshape((2, 2))
    A1 = A1 + A1.T
    return np.trace(np.dot(A1,A1))-2
    # return  np.trace(np.linalg.matrix_power(A1, 2))- 2
def constraint_func6(x):
    A1 = x[8:12].reshape((2, 2))
    A1 = A1 + A1.T
    return np.trace(np.dot(np.dot(A1,A1),np.dot(A1,A1))) -2
    # return  np.trace(np.linalg.matrix_power(A1, 4))- 2
def constraint_func7(x):
    B1 = x[12:16].reshape((2, 2))
    B1 = B1 + B1.T
    return np.trace(np.dot(B1,B1))-2
    # return  np.trace(np.linalg.matrix_power(B1, 2))- 2
def constraint_func8(x):
    B1 = x[12:16].reshape((2, 2))
    B1 = B1 + B1.T
    return np.trace(np.dot(np.dot(B1, B1),np.dot(B1,B1))) -2
    # return  np.trace(np.linalg.matrix_power(B1, 4))- 2
constraint1 = {'type': 'eq', 'fun': constraint_func1}
constraint2 = {'type': 'eq', 'fun': constraint_func2}
constraint3 = {'type': 'eq', 'fun': constraint_func3}
constraint4 = {'type': 'eq', 'fun': constraint_func4}
constraint5 = {'type': 'eq', 'fun': constraint_func5}
constraint6 = {'type': 'eq', 'fun': constraint_func6}
constraint7 = {'type': 'eq', 'fun': constraint_func7}
constraint8 = {'type': 'eq', 'fun': constraint_func8}


result = minimize(objective_function, initial_guess, args=(rho,), constraints=[constraint1,constraint2,constraint3,constraint4,constraint5,constraint6,constraint7,constraint8])

A0_optimal = result.x[:4].reshape((2, 2))
B0_optimal = result.x[4:8].reshape((2, 2))
A1_optimal = result.x[8:12].reshape((2, 2))
B1_optimal = result.x[12:16].reshape((2, 2))
A0_optimal = A0_optimal + A0_optimal.T
B0_optimal = B0_optimal + B0_optimal.T
A1_optimal = A1_optimal + A1_optimal.T
B1_optimal = B1_optimal + B1_optimal.T
print("A0:")
print(A0_optimal)
print("B0:")
print(B0_optimal)
print("A1:")
print(A1_optimal)
print("B1:")
print(B1_optimal)
print(-result.fun)

print(np.trace(np.linalg.matrix_power(A0_optimal, 4))- 2)

Here to make the constraints differentiable, I change the constraint that the absolute value of the eigenvalues of A0 is 1 to Tr(A0^2)=2 and Tr(A0^4)=2 (the same for A1,B0,B1). From Cauchy-Schwarz inequality they are equivalent.

However , the result it gives doesn’t satisfy the constraints. Tr(A_0^4) is about 3.97.

My questions: 1. Why is the scipy program unable to give a feasible solution especially the initial point is feasible? Are the degree-4 constraints hard to calculate (but they are smooth functions)?

2. Are there any better ways to use the constraint that all eigenvalues have absolute value 1?