How do I optimize a multi-vairable scalar function (that involves trace and matrices) in Python?

Question

I am currently trying to optimize a function of the form:

where sigma and x are the optimization variables. We shall say that x is a vector (of the form [x_1, x_2]) that holds the coefficients to X_2 = x_1^2 * A + x_2^2 * B where A and B are positive matrices.

I have attempted to use scipy.optimize but it isn't working. I'm certain that it has to do with the optimization variables being a scalar/matrix. I am attaching my code below:

from scipy.optimize import minimize
import numpy as np

def objective(x):
    
    x_array = x[0]
    sigma = x[1]
    
    fun = np.trace((np.kron(np.diag(np.square(x_array)), rho.T)) @ sigma)
    
    return fun

x0 = np.array([np.array([1 ,1]), np.eye(4)])

res = minimize(objective, x0, method = 'Nelder-Mead')

I am getting an error that says

ValueError: setting an array element with a sequence.

How can I go about solving this optimization problem in Python?

Answer 1

In scipy.optimize the x vector should be 1-dimensional vector of length NUMVAR. You pass on something more complicated. You should do something like:

from scipy.optimize import minimize
import numpy as np


n1 = 2  # number of variables in user x array
n2 = 4  # number of variables in sigma
n = n1 + n2  # number of variables in scipy.optimize x array

def objective(x):
    
    #x_array = x[0]
    #sigma = x[1]
    x_array = x[0:n1]
    sigma = np.diag(x[n1:n])
    
    fun = np.trace((np.kron(np.diag(np.square(x_array)), rho.T)) @ sigma)
    
    return fun

# x0 = np.array([np.array([1 ,1]), np.eye(4)])
# this is a very strange beast. How does one come up with this?
# even numpy thinks this is just bonkers:
#  VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences
# (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. 
#  If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
# let's bring some sanity back and just use a vector of length n=n1+n2
x0 = np.ones((n,))


res = minimize(objective, x0, method = 'Nelder-Mead')

This does not run as we don't have a proper running minimal example. Next time, please have a look at: https://stackoverflow.com/help/minimal-reproducible-example.

Note that I assumed we have 2 variables for x and 4 for sigma and that sigma is a diagonal matrix. If sigma is a full matrix it will have 16 variables. Not clear what the user intended. This is a good example where it is always a good idea to state the problem in clear math first. I don't see any A and B matrices in the code, and in the mathematical model, so the description seems to be out of sync. Of course, we are also missing any constraints. Again, the description is not the same as the math which is not the same as the code.