Scipy minimize constraint: One of two values needs to be zero

Question

I want to minimize the following function :

def objective(B, P, O):
    return - (sum([(p * o - 1) * b for p, o, b in zip(P, O, B)]) /\
    math.sqrt(sum([(1-p) * p * b**2 * o**2 for p, o, b in zip(P, O, B)])))

sol = minimize(objective, x0=bets, args=(P,O,), method='SLSQP', bounds=bnds, constraints=constr)

I want to add the following constraint to B = [1,2,3,4,5,6,....] (B always has even length):

For each pair in the list (1,2), (3,4), (5,6)... (b1,b2), exactly one of the two values should become 0 at the end of optimization. So from a logical standpoint: b1 + b2 = b1 xor b1 + b2 = b2

If I write that as a constraint it looks something like this:

def constb(B):
    for i in range(0, len(B), 2):
        b1, b2 = B[i:i+2]
        if b1 == 0.0:
            return b2 + b1 - b2
        elif b2 == 0.0:
            return b1 + b2 - b1
        else: 
            return 42 

The constraints look like this:

constr = [{'type': 'ineq', 'fun': lambda x: sum(x) - 100},
          {'type': 'eq', 'fun': constb}]

But its not working, since my pairs look like this at the end (My bounds are (0,20) for each value):

20.0 20.0
20.0 20.0
20.0 20.0
20.0 20.0
20.0 20.0

It would seem like the algorithm defaults in the constraint to the else statement. I tried initializing B as 0s but then there arises a MathError because one can't divide by 0.

Is there a way to implement this?

Answer 1

No. That's not possible in general.

A: These constraints are like disjunctions or cardinality-constraints, the typical things which render optimization-problems NP-hard.

B: The solvers in scipy.minimize are based on some strong assumptions and provide local-optima in polynomial-time by design. The core-assumption you are ignoring is:

• constraints + objective are twice differentiable
  • simple rule: if you use if-else stuff in your obj/constraints, it's probably not twice diff!

A and B in conjunction should make it clear, that you are kind of lost.

See also this related introduction from Stanford on Convex-Cardinality Problems which outline the problem. (Keep in mind, that this topic is more general and only related to your problem; for your example, the disjunctive view on the problem is more specialized!)

Either you give up exactness / go for approximation and do those l1-norm tricks as done in machine-learning and co (see lasso-optimization). This is non-trivial to tune = hyper-parameter selection (and also annoying to implement within scipy -> remember => twice-diff => you can't use np.linalg.norm(x, 1) but will need reformulation which also makes the problem much more complex for the solvers within scipy)

Or you give up polynomial-time algorithms and go for mixed-integer mathematical-optimization. This is not part of scipy.

Candidate approaches are then:

• Mixed-Integer Programming (linear)
  • e.g. CoinOR Cbc
• various Convex-Integer Programming
  • QP, SOCP, SDP and co.
• general Nonlinear mixed-integer programming
  • e.g. CoinOR Bonmin (local) / Couenne (global)

I'm too lazy to analyze your objective, but it seems the former is off the table as you have non-linear parts (square-root).