SciPy optimize with additional variables in constraint function

Question

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.

There are tools/solvers that allow a constraint like "matrix A should be pos.def.". See e.g. cvxpy.org/examples/basic/sdp.html. — Erwin Kalvelagen
– Erwin Kalvelagen, Commented May 26, 2022 at 21:25

joni · Accepted Answer · 2022-05-21 11:22:54Z

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

Hi Joni, thanks so much for the great answer (accepted and upvoted). I've started applying these ideas to my full model (different mat[i,j] elements) and am already seeing correct results. Just a quick question: is the non-feasibility of the starting point based on the fact that the constraint is a "minimal eigenvalue of a matrix" constraint? Or does it also depend on the actual contents (the elements) of said matrix? I ask because in my full model, I do believe it does not produce complex eigenvalues at all, so even removing np.real, gives a nice answer.
@user55664123 Yes, it's mainly due to the fact that your matrix is dependent on some optimization variables and needs to be positive definite. If your matrix would have some well-known structure, things would be much easier.

Collectives™ on Stack Overflow

SciPy optimize with additional variables in constraint function

1 Answer 1

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related