0

I have this kind of data :

import random
data=random.sample(range(1, 100), 5)
x= [-1,1,1,-1,1]

def f(x,data):
        prod=[a * b for a, b in zip(data, x)]
        result=abs(sum(prod))
        return result

I Would like to find the best x composed of -1 or 1 to minimize the value of f(x)

Maybe we can use scipy.minimise() but how can we add the -1 or 1 as a constrain on the value inside of x ? Does somebody have an idea ?

Thomas LESIEUR
  • 408
  • 4
  • 14

2 Answers2

2

You want to solve a mixed-integer linear programming problem (MILP), which aren't supported yet by scipy.optimize.

However, you can use a modelling package like PuLP to formulate your MILP and pass it to a MILP solver. Note that your MIP can be formulated as

(P)

min   |f(x)| = |d_0 * x_0 + ... + d_n * x_n|
s.t. x_i ∈ {-1, 1}  ∀ i = 0,...,n

which is the same as

(P')

min   |f(x)| = |d_0 * (2*x_0 - 1) + ... + d_n * (2*x_n - 1)|
s.t. x_i ∈ {0, 1}  ∀ i = 0,...,n

and can be implemented like this

min abs_obj
s.t. f(x) <= abs_obj
     f(x) >= -1.0*abs_obj
     x_i   ∈ {0, 1}  ∀ i = 0,...,n

In code:

import pulp
import random

data = random.sample(range(1, 100), 5)

# pulp model
mdl = pulp.LpProblem("our_model", sense=pulp.LpMinimize)

# the binary variables x
x = pulp.LpVariable.dicts("x", range(5), cat="Binary")
# the variable that stores the absolute value of the objective
abs_obj = pulp.LpVariable("abs_obj")

# set the MIP objective
mdl += abs_obj

# Define the objective: |f(x)| = abs_obj
mdl += pulp.lpSum((2 * x[i] - 1) * data[i] for i in range(5)) <= abs_obj
mdl += pulp.lpSum((2 * x[i] - 1) * data[i] for i in range(5)) >= -1.0*abs_obj

# solve the problem
mdl.solve()

# your solution
signs = [1 if var.varValue > 0 else -1 for var in x.values()]

Alternatively, if you don't want to use another package, you can use scipy.optimize.minimize and implement a simple penalty method. Thereby you solve the problem (P') by solving the penalty problem

min |f(x)| + Ɛ * (x_0 * (1 - x_0) + ... + x_n * (1 - x_n))
with 0 <= x_i <= 1

where Ɛ is a given penalty parameter. Here, the idea is that the right penalty term equals zero for an integer solution.

Note that as the case may be that you need to solve a sequence of penalty problems to achieve convergence to an integer solution. Thus, I'd highly recommend sticking to a MILP solver instead of implementing a penalty method on your own.

joni
  • 6,840
  • 2
  • 13
  • 20
  • Thank you very much ! That help a lot. pulp solution seems great ! Do you have an idea on how i can implement an absolute value on the lpSum ? – Thomas LESIEUR Sep 07 '21 at 08:55
  • 1
    @timebis Sure, I totally missed the absolute value in your question. I edited my answer. – joni Sep 07 '21 at 09:21
  • Thanks a lot. Your help and explanation are amazing :) Have a nice day ! – Thomas LESIEUR Sep 07 '21 at 09:28
  • Do you know how can I add a square of my optimisation problem ? smthg like this : mdl += ( pulp.lpSum((2 * x[i] - 1) * data[i] for i in range(5)) )^2 <= abs_obj I tried many things but nothing works ^^ Is there a function like pulp.lpmult() or smth ? – Thomas LESIEUR Sep 07 '21 at 11:21
  • 1
    This would turn the problem into a *quadratic* problem. Unfortunately, PuLP only supports *linear* problems. However, you could use other modelling libraries like [Pyomo](http://www.pyomo.org/) or [pySCIPOpt](https://github.com/scipopt/PySCIPOpt). Both support quadratic problems but the latter is more beginner-friendly, IMHO. Feel free to ask a separate question in case you need help. – joni Sep 07 '21 at 11:33
  • Thanks, I will have a look – Thomas LESIEUR Sep 07 '21 at 11:40
0

Yes, you can do it using scipy.optimize.minimize:

from scipy.optimize import minimize
minimize(f, [0] * len(data), args=data, bounds=[(-1, 1)] * len(data))

This call minimizes f which you defined in the original post. It passes a zero array as an initial guess for the minimization problem. The argument f requires is 'data' which is specified by the argument 'args'. The constraints you want are specified by the argument 'bounds' as a list of min/max tuples with the length of the input data.

Roman Tz
  • 89
  • 4
  • Thanks for the help ! But i don't want to set bounds : I want the value in x to be -1 or 1 (not between -1 and 1) – Thomas LESIEUR Sep 06 '21 at 20:42
  • Sorry, didn't realize you need specific integer solutions. Generally, SciPy does not support integer programming. See here https://stackoverflow.com/questions/39101137/how-can-i-get-integer-solutions-with-scipy-optimize-linprog If the function f in your example is the real use case, than in order to get the minimized solution you just need to take the sign of every element in 'data' and multiply it by -1, like so: numpy.sign(data) * -1 – Roman Tz Sep 06 '21 at 21:46
  • Please add further details to expand on your answer, such as working code or documentation citations. – Community Sep 06 '21 at 21:54
  • No, it's not the real use case, f is just a small example. – Thomas LESIEUR Sep 06 '21 at 22:04
  • Then can you please elaborate your specific use case? what are the requirements for your input and output? – Roman Tz Sep 06 '21 at 22:31