0

I am looking for an optimization method in scipy which allows me to minimize object function f(x,y) (returns vector) subject to the constraint g(x,y) < 0.1 and additional bounds on x and y.

I have tried to solve my problem with scipy.optimize.least_squares, scipy.optimize.leastsq and scipy.optimize.minimize. The problem is that leastsq and least_squares allow the object function to be non-scalar, but does not give me the possibility of implementing a constraint (only bounds). minimize on the other hand gives me the possibility of implementing both a constraint and bounds, but f(x,y) must return a scalar. Hence, I am looking for a solution that combines both. Does anyone know whether something like this exists?

The function I want to minimize is

def my_cost(p,f_noise):
    x,y = p[0], p[1]
    f = #some function that returns a 3x1 array
    return (f - fnoise)**2

I did this with the least_squares method.

opti.least_squares(my_cost, p0[:], args = (f_noise,),gtol=1e-2, bounds=bounds)

But here I have the problem that I cannot constrain the variables in p. I need to constrain p so that it fulfils

def constraint(p)
    x = p[0]
    return fy(x) - y <= 0.1 #variable y therefore becomes a function of variable x

To implement the constraint, I tested scipy's minimize function

opti.minimize(my_cost, p0[:], args = (f_noise,), bounds = bounds, constraints={'type': 'eq', 'fun': constraint})

But here I can't seem to find a way to allow my_cost and f_noise to be 3x1 arrays.

For any help I am very grateful! Cheers for your time!

SuperKogito
  • 2,998
  • 3
  • 16
  • 37
Simon
  • 3
  • 1

1 Answers1

1

According to the docs, the objective function must return a float when using scipy.optimize.minimize, whereas with scipy.optimize.least_squares, you cannot use constraints. In this case, you have to be aware of your minimization purpose. Minimizing a difference vector (like f-f_noise) is equivalent to minimizing the element-wise differences and consequently their sum. Therefore a practical solution would be to minimize a defined p-norm of your f(x,y) and g(x). I suggest the square L2-norm as it is very similar to what you are trying in your cost function and it is simple and stable (in comparison to other norms).

enter image description here

You can average the norm and get the Mean Squared Error (MSE):

enter image description here

By applying the previous concepts, you get the following code:

import numpy as np 
from scipy.optimize import minimize

# define fy
def fy(x):
    return x**2 * np.array([[.1],[.2],[.3]])  # some function that returns a 3x1 array

# objective func
def f(p, f_noise):
    x, y = p[0], p[1]
    f    = x * y * np.array([[1],[2],[3]])    # some function that returns a 3x1 array
    return np.linalg.norm(f - f_noise, 2)**2

# constraint
def g(p):
    x         = p[0]
    diff_norm = np.linalg.norm(fy(x) - y) 
    return threshold - diff_norm 

# init 
f_noise   = np.array([[.1],[.2],[.3]])
p0        = np.array([1, 0.5])
bounds    = ((0,2),(0,2))
y         = np.array([[.9],[.7],[.2]])
threshold = 0.1  # make sure to choose an adequate threshold

# minimize
result  =  minimize(f, p0,
                    args        = (f_noise,), 
                    bounds      = bounds, 
                    constraints = {'type': 'ineq', 'fun': g})

# print result
print(result)
SuperKogito
  • 2,998
  • 3
  • 16
  • 37
  • Thanks for your very helpful and quick answer, SuperKogito! That's amazing! I will test it and reply later, wether this solved my problem. – Simon Jun 14 '19 at 09:40
  • Thanks a lot, SuperKogito. This seems to work allthough I am still running some tests. I have one question to this solution: The constraint is select via the threshold. Doesn't this reduce the found solution to the best value where threshold = diff_norm instead of best value for diffnorm < threshold? I.e. in your solution I understand the constraint to be g(x,y) = 0.1 instead of g(x,y) < 0.1. Seems to be some kind of Lagrangian Constraint to me. – Simon Jun 14 '19 at 15:04
  • This depends more on the nature of your problem. It is for you to decide which constraint type (equality or inequality) is more efficient in your case. I think an inequality is probably faster because trivially it is easier to find a value lower than the threshold compared to finding one that is equal to the threshold. – SuperKogito Jun 15 '19 at 19:16
  • Yes, exactly. I think it would be faster. I would like to constrain the the minimization to all points where g(x,y) < 0.1. But as I understand it, the code currently looks for a solution where g(x,y) is closest to 0.1. But I don't know how to implement the inequality constraint. – Simon Jun 15 '19 at 20:10
  • In my example, I am using the inequality constraint. I defined it in this line `constraints = {'type': 'ineq', 'fun': g}`. Therefore, the minimization is looking for a solution while making sure that the return of the constraint is positive: `g(p) >= 0 -> threshold - diff_norm >= 0 -> threshold >= diff_norm`. Please read the constraint section [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) and refer to [this](https://stackoverflow.com/questions/42303470/scipy-optimize-inequality-constraint-which-side-of-the-inequality-is-considere) – SuperKogito Jun 15 '19 at 20:29
  • Ah perfect. I was looking for the 'ineq' command but must have overseen it. Thank you ever so much, Kogito! You helped me a lot! – Simon Jun 15 '19 at 20:42
  • You are welcome :) consider accepting the answer if it helped. – SuperKogito Jun 16 '19 at 09:05