I am mostly new to Python and StackOverflow so please point out if I am doing anything wrong/need to include more detail e.t.c. My question revolves around minimizing a function, en_sum, which depends on an array of variables en. This optimisation must be subjected to a constraint nlc: the minimum eigenvalue of a matrix which depends on variables en as well as another set of variables, x2, must be positive.
What I would like to do is essentially provide some guess values of en and x2, optimise en_sum subject to nlc and have it produce the resulting optimised en and x2 values.
The minimal code I have written is as follows:
import numpy as np
from scipy import optimize as spo
from scipy.optimize import NonlinearConstraint
def en_min(varia):
#en_min is a function that (aims to) return the optimised values of en and x2.
m = 2
n = 2
s = 2
en = varia[0:m]
x2 = varia[m:len(varia)]
# use_mat(a,b) defines the matrix which depends on en and x2 variables.
def use_mat(a, b):
mat = np.empty([m * s, n * s])
for i in range(0, m * s):
for j in range(0, m * s):
mat[i, j] = (b[j] * a[i - s]) ** 2 + a[j - s] ** 2 + (b[j] + b[i]) ** 2
return mat
# Defining the optimisation function.
def en_sum(q):
return sum(q[iii] ** 2 for iii in range(len(q)))
# This function provides the minimal eigenvalue of a given matrix.
def min_eigenvalue(a):
y = np.linalg.eigvals(a)
return np.amin(y)
# My attempt at setting the constraint (which is that the eigenvalues of the matrix
# defined through use_mat(a,b) are non-negative) is nlc
nlc = NonlinearConstraint(min_eigenvalue(use_mat(en, x2)), 0, np.inf)
# The guess values in opt_res2 are just for the en variables, as en_sum must
# only be a function of those (not x2). However, I do not know how to properly give guess
# values to x2 in the constraint, so I am quite sure this is implemented incorrectly.
opt_res2 = spo.minimize(en_sum(en), np.array([1, 4]), constraints=nlc)
return opt_res2
print(en_min([1, 3, 1.5, 0, 0, 4]))
Naturally, this does not work (I get error TypeError: 'numpy.complex128' object is not callable amongst others), and I wondered if anyone had any tips about where I could be going wrong, and whether the problem is even tractable in this case? My main confusion comes from not knowing how to properly treat the constraint and how to supply the guess values for x2 (is this related to args in the constraint dictionary?). Any ideas are greatly appreciated, thank you.
Interesting problem. IMO, it's a best practice to write the objective function and all the constraints as a function of one (vectorial) optimization variable z. Inside your functions, you can then simply set z = (en, x2) and split the variable z back into en and x2. Thus, your objective function can be written as follows:
def en_sum(z):
en, x2 = z[:m], z[m:]
return np.sum(en**2)
Note also that en_sum(en) and min_eigenvalues(use_mat(en, x2)) aren't functions and thus, aren't callable. Furthermore, your constraint that the minimal eigenvalue has to be positive (the matrix returned by use_mat positive definite), is equivalent to all eigenvalues being positive. This part is a bit messy since it's not trivial to find a feasible starting point such that this constraint is fulfilled. However, we can cheat by only considering the real part of complex eigenvalues. As a consequence, the constraint is fulfilled even if the current point is not feasible. Here, we basically hope that it still converges into the feasible region (where all eigenvalues are real and positive) during the first few iterations.
In code:
import numpy as np
from scipy.optimize import minimize, NonlinearConstraint
def en_min(varia):
#en_min is a function that (aims to) return the optimised values of en and x2.
m = 2
n = 2
s = 2
en = varia[0:m]
x2 = varia[m:len(varia)]
# use_mat(a,b) defines the matrix which depends on en and x2 variables.
def use_mat(z):
a, b = z[:m], z[m:]
mat = np.empty([m * s, n * s])
for i in range(0, m * s):
for j in range(0, m * s):
mat[i, j] = (b[j]*a[i-s])** 2 + a[j-s]**2 + (b[j]+b[i])**2
return mat
# Defining the optimisation function.
def en_sum(z):
en, x2 = z[:m], z[m:]
return np.sum(en**2)
# This function provides the minimal eigenvalue of a given matrix.
def min_eigenvalue(a):
return np.real(np.linalg.eigvals(a))
# constraint
nlc = NonlinearConstraint(lambda z: min_eigenvalue(use_mat(z)), 1e-12, np.inf)
# Bounds on the variables
bounds = [(None, None)]*en.size + [(0, None)]*x2.size
# initial guess
x0 = np.hstack((en, 0.1*np.ones(x2.size)))
opt_res2 = minimize(en_sum, x0=x0, bounds=bounds, constraints=nlc)
return opt_res2
print(en_min(np.array([1, 3, 1.5, 0, 0, 4])))
Alternatively, since a matrix A is positive definite if and only if there's a matrix V such that A = V'V, one could solve an auxiliary problem in order to find a feasible initial guess.
It's also worth mentioning that the above code isn't really efficient and will be quite slow for large optimization problems due to the approximation of the Jacobian (= the derivatives of our constraint function): Currently, each evaluation of the Jacobian requires N*N evaluations of min_eigenvalue, where N is the total number of variables. Thus, we calculate N*N times the eigenvalues per iteration. All in all, for larger problems it's definitely worth providing the exact Jacobian or using algorithmic differentiation