Optimization with joint bounds / symbolic constraint

Question

I am trying to perform least squares optimization on a vector [a1,a2,a3], with the following constraint, where k is a constant:

-k < a3/a2 < k

I will post a simplified version of my code in the hope that this makes it clear enough what I am after.

from scipy import optimize

def loss_function(a_candidate):
    return MyObject(a_candidate).mean_square_error()

def feasibility(a_candidate):
    # Returns a boolean
    k = 1.66711
    a2 = a_candidate[1]
    a3 = a_candidate[2]
    return -k < a3/a2 < k

def feasibility_scipy(a_candidate):
    # As I understand it, for SciPy the constraint has to be a number
    boolean = feasibility(a_candidate)
    if boolean:
        return 0
    else:
        return 1

# Equality constraint forcing feasibility_scipy to be 0.
constraint = optimize.NonlinearConstraint(feasibility_scipy, 0, 0)


optimize_results = optimize.minimize(
    loss_function,
    a_init,  # Initial guess
    constraints=constraint)

Due to how the initial guess a_init is generated, it is in a region where feasibility is False. (The reason we need to use numerical methods in the first place is that an earlier closed-form method returned an infeasible solution. It is possible to supply a very poor feasible guess like (0,0,0), but this will be much further from the true solution).

Since constraint's gradient is zero almost everywhere, the optimization routine cannot find its way out of this infeasible (inadmissible) region, and it does not terminate successfully. With SLSQP it stops after just 1 iteration, with the message Singular matrix C in LSQ subproblem. With the trust-constr solver, it reaches the maximum number of function evaluations, and I believe it has not left the infeasible region because constr_violation is 1.0.

As I understand it, it is not possible in SciPy to provide a 'joint bound' (apologies for invented terminology) on a2 and a3, meaning I am forced to use the NonlinearConstraint approach.

What is the proper way to do this? (A few searches have suggested that I may want to try the mystic package with a symbolic constraint. But before I invest the time to learn this new package I wanted to see if StackOverflow has a SciPy-based solution. Or if you know how to do this in mystic some example code would be very helpful.)

Answer 1

Nevermind, I think the solution may be just:

def a_ratio(a_candidate):
    a2 = a_candidate[1]
    a3 = a_candidate[2]
    return a3/a2

feasibility_constraint = optimize.NonlinearConstraint(a_ratio,-k,k)

Answer 2

Instead of using -k < a3/a2 < k I would make this two linear constraints:

     -k*a2 <= a3
     a3 <= k*a2 

Answer 3

I know the question has been answered already, but you wanted to know how to approach it in mystic. I'm the mystic author. Something like this:

>>> import mystic as my
>>> import numpy as np
>>> 
>>> truth = np.array([1,1,1])
>>> 
>>> 
>>> class MyObject(object):
...   def __init__(self, candidate):
...     self.candidate = candidate
...     self.truth = truth
...   def mean_square_error(self):
...     return ((self.truth - self.candidate)**2).sum()
...
>>>
>>> def loss_function(a_candidate):
...    return MyObject(a_candidate).mean_square_error()
...
>>>

and then the solution:

>>> equations = """
... -k < a3/a2
... a3/a2 < k
... """
>>>
>>> eqn = my.symbolic.simplify(equations, variables=['a1','a2','a3'], locals=dict(k=1.66711), all=False, target='a3')
>>> c = my.symbolic.generate_constraint(my.symbolic.generate_solvers(eqn, variables=['a1','a2','a3']), join=my.constraints.or_)
>>> res = my.solvers.fmin(loss_function, x0=[2,3,4], constraints=c, disp=True, xtol=1e-8)
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 121
         Function evaluations: 229
>>> res
array([1., 1., 1.])

What mystic does differently is that it builds an "operator" c that doesn't allow the optimizer to select infeasible solutions. Not that it rejects bad solutions with a boolean, but that it transforms "bad" candidates into "good" ones, like this:

>>> c([2,10,4])
[2, 10, 4]
>>> c([2,1,4])
[2, 1, 1.6671099999999974]