I have a simple MWE, where I try to optimize the selection of the columns of a matrix, so that the columns sum to a minimum number.
using JuMP
using HiGHS
# matrix with 20 columns:
m = rand(1000, 20)
# My objective: find the columns of the matrix that sum to the least
# i.e., given a boolean vector with trues meaining what to choose, minimize:
function objective(chosen)
r = m[:, chosen]
q = sum(sum(c) for c in eachcol(r))
return q
end
# random choice of chosen:
chosen = falses(20)
chosen[[1, 5, 6]] .= true
objective(chosen)
# Attempt with Jump, following the Sudoku tutorial
model = Model(HiGHS.Optimizer)
# The variable chosen is a boolean vector with
# the `true` entries meaning a choice that column
@variable(model, choices[1:20], Bin);
# Constrain the amount of columns that can be chosen
Fmin = 2
Fmax = 3
@constraint(model, Fmin ≤ sum(choices) ≤ Fmax)
@objective(model, Min, objective(choices))
# errors
This example errors because the last line that calls @objective yields: ERROR: ArgumentError: unable to check bounds for indices of type VariableRef.
I’ve searched around in the JuMP docs and even though I found examples similar to the one I need (optimizing recipes for example) I couldn’t find something where the decision variables are used to index an “external” array.
What is the solution I am missing here?
p.s.: In my real use case the objective function is extremely complicated and cannot be expressed in terms of mathematical equations. Nevertheless, it only needs to access the columns of a pre-existing matrix in such a way.
I don’t think you can index directly, but often you don’t need to. For example, you can use them as a mask, e.g.
m = rand(1000, 20)
chosen = rand(Bool, 20) # not the actual variables, just to show
column_sums = [sum(c) for c in eachcol(m)]
selected_sums = column_sums .* chosen # set non-chosen sums to 0
sum(selected_sums) # final sum of selected columns
The reason why you can’t use them directly as indices (as I understand it) is that to solve the whole problem, solvers will often solve relaxations of the problem in which the integer constraints are lifted, meaning those variables may be floating point values during some sub-problems. This can be used to develop bounds on the solution which the solver can use in proving optimality. So the objective needs to be well-defined even when the values are not integers.
Additionally, note the final values chosen by the solver may not be exactly integers either, but rather only within the solver’s allowed tolerance of an integer (e.g. 0.9999 instead of 1.0).
Thanks for the reply, but this is disheartening. I have a decision problem to solve here, not really a summation problem.
I cooked up this trivial MWE to highlight the issue but in the real application I need to choose some of the columns of the matrix and process them in a non trivial way that cannot be written as just a pre-computed sum that is then multiplied by chosen.
In my real example the chosen columns are clustered using DBSCAN, and I want to optimize (maximize) the clustering quality. So the process is not “separatable” nor “linear”: I can’t use all columns to produce a large output and then simply “select” the ouputs.
You could mask the original matrix instead of the sums in the same way (m .* permutedims(chosen)). But if you are trying to propagate JuMP VariableRef’s through an entire clustering process I don’t think that is going to work.
At this level of using JuMP (linear or quadratic programming), you are essentially formulating a problem symbolically in order to describe a problem for the solver to solve. Maybe the non-linear interface could help.
Thanks for the reply, but this is disheartening. I have a decision problem to solve here, not really a summation problem.
I hear you. I think this is a common feeling for people using JuMP for the first time who don’t come from a math programming background. It feels like JuMP is a generic optimization library, but it actually requires a very structured model in the language of mathematical programming. You can’t just thrown an arbitrary Julia function at it.
What is your quality function? I get that it’s not a straight summation, but here’s how I’d write your example above:
using JuMP, HiGHS
m = rand(1000, 20)
model = Model(HiGHS.Optimizer)
@variable(model, choices[1:20], Bin)
@constraint(model, 2 <= sum(choices) <= 3)
@objective(model, Min, sum(m, dims = 1)' * choices)
Thanks for the reply guys. I have a much better understanding now.
My function cannot be written in mathematical terms. It is a clustering operation, whose purpose is to find attractors of chaotic dynamical systems. The output of the function is a real number: the clustering quality multiplied by the number of attractors found. There is no way you can get any analytic expressions there.
But I now know that it is simply not a good problem for JuMP. I will code a brute force random samling loop that samples columns, and I will also use BlackBoxOptim.jl as an alternative, as this can indeed accept any arbitrary Julia function (we have used it before already!).
Ah actually, BlackBoxOptim.jl is wholly unsuitable for this application: it searches in a box of real-valued variables for an optimum. I have instead a high dimensional input space, but it each dimension is binary instead of real valued.
If you are looking for a black-box optimization solver, maybe have a look at NOMAD? It allows you to use binary variables by setting “granularity” parameter and can take any function.
Hi @Datseris, just wondering in the end which package worked for you? I am currently having a very similar problem in which the decision variables (integers/binaries) are used as the indices of the equations.
Would it make sense to try a form of local search? That is a common approach for combinatorial optimization when mathematical programming is either infeasible or intractable.
For instance, start with a vector of chosen variables and then at each iteration try removing or adding every possible variable, then choose the option with the best performance.
Yes this is what I do now but I call this simply the “brute force” approach, as there is no optimization, I just go through all possible combinations. This isn’t scalable but for now I am limiting myself to at most 4 max selections anyways so it works.