2

I'm running a constrained optimisation with scipy.optimize.minimize(method='COBYLA').

In order to evaluate the cost function, I need to run a relatively expensive simulation to compute a dataset from the input variables, and the cost function is one (cheap to compute) property of that dataset. However, two of my constraints are also dependent on that expensive data. So far, the only way I have found to constrain the optimisation is to have each of the constraint functions recompute the same dataset that the cost function already has calculated (simplified quasi-code):

def costfun(x):
    data = expensive_fun(x)
    return(cheap_fun1(data))

def constr1(x):
    data = expensive_fun(x)
    return(cheap_fun2(data))

def constr2(x):
    data = expensive_fun(x)
    return(cheap_fun3(data))

constraints = [{'type':'ineq', 'fun':constr1},
               {'type':'ineq', 'fun':constr2}]

# initial guess
x0 = np.ones((6,))

opt_result = minimize(costfun, x0, method='COBYLA', 
                      constraints=constraints)

This is clearly not efficient because expensive_fun(x) is called three times for every x.

I could change this slightly to include a universal "evaluate some cost" function which runs the expensive computation, and then evaluates whatever criterion it has been given. But while that saves me from having to write the "expensive" code several times, it still runs three times for every iteration of the optimizer:

# universal cost function evaluator
def criterion_from_x(x, cfun):
    data = expensive_fun(x)
    return(cfun(data))

def costfun(data):
    return(cheap_fun1(data))

def constr1(data):
    return(cheap_fun2(data))

def constr2(data):
    return(cheap_fun3(data))


constraints = [{'type':'ineq', 'fun':criterion_from_x, 'args':(constr1,)},
               {'type':'ineq', 'fun':criterion_from_x, 'args':(constr2,)}

# initial guess
x0 = np.ones((6,))

opt_result = minimize(criterion_from_x, x0, method='COBYLA',
                      args=(costfun,), constraints=constraints)

I have not managed to find any way to set something up where x is used to generate data at each iteration, and data is then passed to both the objective function as well as the constraint functions.

Does something like this exist? I've noticed the callback argument to minimize(), but that is a function which is called after each step. I'd need some kind of preprocessor which is called on x before each step, whose results are then available to the cost function and constraint evaluation. Maybe there's a way to sneak it in somehow? I'd like to avoid writing my own optimizer.

One, more traditional, way to solve this would be to evaluate the constraints in the cost function (which has all the data it needs for that, have it add a penalty for violated constraints to the main cost function, and run the optimizer without the explicit constraints, but I've tried this before and found that the main cost function can become somewhat chaotic in cases where the constraints are violated, so an optimizer might get stuck in some place which violates the constraints and not find out again.

Another approach would be to produce some kind of global variable in the cost function and write the constraint evaluation to use that global variable, but that could be very dangerous if multithreading/-processing gets involved, or if the name I choose for the global variable collides with a name used anywhere else in the code: ''' def costfun(x): global data data = expensive_fun(x) return(cheap_fun1(data))

def constr1(x): global data return(cheap_fun2(data))

def constr2(x): global data return(cheap_fun3(data)) '''

I know that some people use file I/O for cases where the cost function involves running a large simulation which produces a bunch of output files. After that, the constraint functions can just access those files -- but my problem is not that big.

I'm currently using Python v3.9 and scipy 1.9.1.

Zak
  • 3,063
  • 3
  • 23
  • 30

2 Answers2

1

You could write a decorator class in the same vein to scipy's MemoizeJac that caches the return values of the expensive function each time it is called:

import numpy as np

class MemoizeData:
    def __init__(self, obj_fun, exp_fun, constr_fun):
        self.obj_fun = obj_fun
        self.exp_fun = exp_fun
        self.constr_fun = constr_fun
        self._data = None
        self.x = None

    def _compute_if_needed(self, x, *args):
        if not np.all(x == self.x) or self._data is None:
            self.x = np.asarray(x).copy()
            self._data = self.exp_fun(x)

    def __call__(self, x, *args):
        self._compute_if_needed(x, *args)
        return self.obj_fun(self._data)

    def constraint(self, x, *args):
        self._compute_if_needed(x, *args)
        return self.constr_fun(self._data)

Followingly, the expensive function is only evaluated once for each iteration. Then, after writing all your constraints into one constraint function, you could use it like this:

from scipy.optimize import minimize

def all_constrs(data):
    return np.hstack((cheap_fun2(data), cheap_fun3(data)))

obj = MemoizeData(cheap_fun1, expensive_fun, all_constrs)
constr = {'type': 'ineq', 'fun': obj.constraint}

x0 = np.ones(6)
opt_result = minimize(obj, x0, method="COBYLA", constraints=constr)
Zak
  • 3,063
  • 3
  • 23
  • 30
joni
  • 6,840
  • 2
  • 13
  • 20
0

While Joni was writing their answer, I found another one, which is admittedly more hacky. I prefer theirs, but for the sake of completeness, I wanted to post this one, too.

It's derived from the material from https://mdobook.github.io/ and the accompanying video tutorials from BYU FLow Lab, in particular this video:

The trick is to use non-local variables to keep a cache of the last evaluation of the expensive function:

import numpy as np

last_x = None
last_data = None

def compute_data(x):
    data = expensive_fun(x)
    return(data)

def get_last_data(x):
    nonlocal last_x, last_data
    if not np.array_equal(x, last_x):
        last_data = compute_data(x)
        last_x = x
    return(last_data)

def costfun(x):
    data = get_last_data(x)
    return(cheap_fun1(data)

def constr1(x):
    data = get_last_data(x)
    return(cheap_fun2(data)

def constr2(x):
    data = get_last_data(x)
    return(cheap_fun3(data)

...and then everything can progress as in my original code in the question.

Reasons why I prefer Joni's class-based version:

  • variable scopes are clearer than with nonlocal
  • If some of the functions allow calculation of their Jacobian, or there are other things worth buffering, the added complexity is held in check better than with
  • Having a class instance do all the work also allows you to do other interesting things, like keeping a record of all past evaluations and the path taken by the optimizer, without having to use a separate callback function. Very useful for debugging/tweaking convergence if the optimizer won't converge or takes too long, but also to visualize or otherwise investigate the objective function or similar.
  • The same ability might actually be really cool for things like constructing a response surface model from the results of previous function evaluations. That could be used to establish a starting guess in case the expensive function is some numerical method that benefits from a good starting point.

Both approaches allow the use of "cheap" constraints which don't require the expensive function to be evaluated, by simply providing them as separate functions. Not sure whether that would help much with compute times, though. I suppose that would depend on the algorithm used by the optimizer.

Zak
  • 3,063
  • 3
  • 23
  • 30