scipy - 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.)

question from:https://stackoverflow.com/questions/65943996/optimization-with-joint-bounds-symbolic-constraint

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)