Pi (Π) operation in the objective function JuMP

[Shayan]

I want to solve a quadratic programming problem with JuMP.
I don’t know if Ipopt can solve it or not. I’m using HiGHS (Although I prefer Ipopt if it’s possible). Here is the model:

\eqalign{ & \max z = \prod\limits_{i = 1}^m {b \cdot {x_i}} \cr & s.t. \cr & \sum\limits_{j = 1}^n {{b_j} = 1} \cr & 0 \le {b_j} \le 1,{\rm{ }}j \in \{ 1,2, \ldots n\} \cr}

where x_i is a vector of size n and b is a vector of size n. Here is the code:

using HiGHS
using JuMP

function example()
  n = m = 3 
  
  # A matrix to store the output of each round
  weights = Matrix{Float64}(undef, n, m)

  # Synthetic data
  data = rand(n, 10);

  # Arbitrary indices
  rnd_idxs = rand(1:3, 3)

  # run the procedure for 3 times
  for iter ∈ 1:3

    # size(x) is (3, 3)
    x = data[:, rnd_idxs]

    # Since HiGHS supports quadratic programming, I chose it.
    model = Model(HiGHS.Optimizer)

    @variables(model, begin
      0<=b[j=1:n]<=1
    end)

    @constraint(model, sum(b[j] for j = 1:n) == 1)

    # Part of the objective function. x[:, i] represents the xᵢ
    h = [sum(b.*x[:, i]) for i ∈ 1:m]

    @NLobjective(model, Max, reduce(*, h))

    optimize!(model)

    # Store the output
    weights[:, iter] = value.(b)
  end

  weights
end;

example()

But it throws the following error:

┌ Warning: Function reduce automatically registered with 2 arguments.
│
│ Calling the function with a different number of arguments will result in an
│ error.
│
│ While you can safely ignore this warning, we recommend that you manually
│ register the function as follows:
│ ```Julia
│ model = Model()
│ register(model, :reduce, 2, reduce; autodiff = true)
│ ```
└ @ MathOptInterface.Nonlinear C:\Users\Shayan\.julia\packages\MathOptInterface\4g9vU\src\Nonlinear\operators.jl:370
ERROR: Unexpected object * of type typeof(*) in nonlinear expression.
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:35
 [2] parse_expression(#unused#::MathOptInterface.Nonlinear.Model, #unused#::MathOptInterface.Nonlinear.Expression, x::Function, #unused#::Int64)
   @ MathOptInterface.Nonlinear C:\Users\Shayan\.julia\packages\MathOptInterface\4g9vU\src\Nonlinear\parse.jl:34
 [3] parse_expression(data::MathOptInterface.Nonlinear.Model, expr::MathOptInterface.Nonlinear.Expression, x::Expr, parent_index::Int64)
   @ MathOptInterface.Nonlinear C:\Users\Shayan\.julia\packages\MathOptInterface\4g9vU\src\Nonlinear\parse.jl:52
 [4] parse_expression
   @ C:\Users\Shayan\.julia\packages\MathOptInterface\4g9vU\src\Nonlinear\parse.jl:14 [inlined]
 [5] set_objective
   @ C:\Users\Shayan\.julia\packages\MathOptInterface\4g9vU\src\Nonlinear\model.jl:43 [inlined]
 [6] set_nonlinear_objective(model::Model, sense::MathOptInterface.OptimizationSense, x::Expr)
   @ JuMP C:\Users\Shayan\.julia\packages\JuMP\pQApG\src\nlp.jl:167
 [7] macro expansion
   @ C:\Users\Shayan\.julia\packages\JuMP\pQApG\src\macros.jl:2688 [inlined]
 [8] example()
   @ Main e:\CORN.jl:139
 [9] top-level scope
   @ e:\CORN.jl:150

The problem is with this line:
@NLobjective(model, Max, reduce(*, h))
How can I solve this problem?

Note that For the example above, for each iteration (iter), the objective function can be written as:
\max z = \left( {{b_1} \times {x_{11}} + {b_2} \times {x_{12}} + {b_3} \times {x_{13}}} \right) \times \cdots \times \left( {{b_1} \times {x_{31}} + {b_2} \times {x_{32}} + {b_3} \times {x_{33}}} \right)

[Shayan]

By replacing @NLobjective(model, Max, reduce(*, h)) with @NLobjective(model, Max, *(h...)) and using Ipopt.Optimizer instead of HiGHS.Optimizer the error is gone. But the sum of each column in the weights variable is bigger than 1!! I mean, it seems that the @constraint(model, sum(b[j] for j = 1:n) == 1) constraint is omitted somehow.

3×3 Matrix{Float64}:
  2.01605e-8   2.01605e-8   2.01605e-8
  1.0          1.0          1.0
 -5.80865e-9  -5.80865e-9  -5.80865e-9

Also, there are negative values among them, which indicates the following command is omitted somehow as well:

@variables(model, begin
      0<=b[j=1:n]<=1
end)

[odow]

Here’s how I’d write your model:

using JuMP, Ipopt
function example()
    n = m = 3 
    weights = zeros(n, m)
    data = rand(n, 10);
    rnd_idxs = rand(1:3, 3)
    for iter in 1:3
        x = data[:, rnd_idxs]
        model = Model(Ipopt.Optimizer)
        @variable(model, 0 <= b[1:n] <= 1)
        @constraint(model, sum(b) == 1)
        @expression(model, h, b' * x)
        @NLobjective(model, Max, prod(h[i] for i in 1:m))
        optimize!(model)
        weights[:, iter] .= value.(b)
    end
    return weights
end
example()

But the sum of each column in the weights variable is bigger than 1!!
Also, there are negative values among them, which indicates the following command is omitted somehow as well:

You should know that Ipopt isn’t an exact solver. Tolerances of ~1e-8 are expected. That applies for both variable bounds, which is why you have a value like -5.80865e-9, and it applies to constraints, which is why some columns do not sum exactly to 1.0.

[Shayan]

@odow

You should know that Ipopt isn’t an exact solver.

Is there any possible way to tackle the inexactness? I.e. by setting etol?

@odow

Here’s how I’d write your model:

Thank you! Does mine do it wrongly?

[odow]

Is there any possible way to tackle the inexactness?

You can change the tol: Ipopt: Ipopt Options

set_attribute(model, "tol", 1e-9)

but you can never achieve a tolerance of exactly 0.0. Rather than trying to find the minimum tolerance, you should step back and ask what level of tolerance is acceptable for your problem. If the variables are weights that sum to 1, then 1e-8 is one part in 100 million. That seems pretty reasonable for practical applications.

Does mine do it wrongly?

No. You can just clean up the syntax a little. No need for sum(b[j] for j = 1:n) etc.

[Shayan]

@odow

If the variables are weights that sum to 1

Yes, exactly. Actually, 4 digits after the decimal point is quite sufficient in my use case!

@odow

No. You can just clean up the syntax a little.

Great! I’ll update mine to what you’ve written. I just wanted to check if I had mistakes on the code. Thank you so much.

[odow]

Actually, 4 digits after the decimal point is quite nice in my use case!

Precisely

You could do something like this to get rid of the small values:

julia> x = [2.01605e-8   2.01605e-8   2.01605e-8
         1.0          1.0          1.0
        -5.80865e-9  -5.80865e-9  -5.80865e-9]
3×3 Matrix{Float64}:
  2.01605e-8   2.01605e-8   2.01605e-8
  1.0          1.0          1.0
 -5.80865e-9  -5.80865e-9  -5.80865e-9

julia> y = round.(x; digits = 4)
3×3 Matrix{Float64}:
  0.0   0.0   0.0
  1.0   1.0   1.0
 -0.0  -0.0  -0.0

Thank you so much.

No problem.

[Shayan]

@odow

You could do something like this to get rid of the small values

Yes, I did it exactly. Also, I used abs to get rid of negative values. However, abs is redundant since round can handle it as well. Thank you for the time you’ve put into this! I appreciate it.

[Shayan]

Now, I remember why I used the abs. To get rid of those -0.0s!

[miles.lubin]

It’s also worth trying out maximizing log(objective) which breaks up into a sum.

[Dan]

@Shayan

the model:

\eqalign{ & \max z = \prod\limits_{i = 1}^m {b \cdot {x_i}} \cr & s.t. \cr & \sum\limits_{j = 1}^n {{b_j} = 1} \cr & 0 \le {b_j} \le 1,{\rm{ }}j \in \{ 1,2, \ldots n\} \cr & }

where x_i is a vector of size n and b is a vector of size n

(sorry for messed up, LaTeX quote)

Orthogonal to the answer about JuMP, mathematically, the problem is in a convex domain (the intersection of a plain and a cube), and the objective is a convex function. The maximum of a convex function in a convex domain is on the boundary, and in this case, the boundary is in {\{0,1\}^n}. So the problem could be solved using HiGHS solver with Bool variables, and then rounding would not be an issue.

Even further, the actual solution is an indicator vector of the highest coefficient of x, and could be found in a closed-form way using argmax. But I may have got confused about the problem or forgotten to dot some i or cross some t.

[mtanneau]

The following assumes that all x_{i} are positive. If some x_{i} have negative entries, feel free to ignore this.

As mentioned by mlubin, maximizing a product is equivalent to maximizing its log, which yields a sum of smaller log terms. Namely, you problem becomes

\begin{array}{rl} \displaystyle \max_{b} \quad & \sum_{i=1}^{m} \log(b^{\top} x_{i})\\ s.t. \quad & \sum_{j} b_{j} = 1\\ & 0 \leq b_{j} \leq 1 \forall j \end{array}

which is equivalent to

\begin{array}{rl} \displaystyle \max_{b} \quad & \sum_{i = 1}^{m} t_{i}\\ s.t. \quad & t_{i} \leq \log(z_{i})\\ & z_{i} = b^{\top} x_{i}\\ & \sum_{j} b_{j} = 1\\ & 0 \leq b_{j} \leq 1, \forall j \end{array}

Finally, this can be formulated as a conic optimization problem with exponential cones, which can be solved by Mosek as well as a number of Julia solvers (e.g. Hypatia).
Conic formulations are typically more numerically stable, faster to solve (because they use more specialized algorithms than Ipopt), and will automatically take care of numerical issues like variables becoming slightly negative.

[Shayan]

Hi Dan,
Sorry for the late response. I’ve not visited Discourse until now.

@Dan

So the problem could be solved using HiGHS solver with Bool variables

Unfortunately, I don’t know how to specify the Boolean variables in JuMP. You’re right, in that case, the rounding wouldn’t be needed.

@Dan

the actual solution is an indicator vector of the highest coefficient of x, and could be found in a closed-form way using argmax.

I don’t know . Considering m=3, and n=2, we can rewrite the obj as:
\eqalign{ & \max z = \left( {\left[ {{b_1},{b_2}} \right] \odot \left[ {{x_{11}},{x_{12}}} \right]} \right) \times \left( {\left[ {{b_1},{b_2}} \right] \odot \left[ {{x_{21}},{x_{22}}} \right]} \right) \times \left( {\left[ {{b_1},{b_2}} \right] \odot \left[ {{x_{31}},{x_{32}}} \right]} \right) \cr & \max z = \left( {{b_1}{x_{11}} + {b_2}{x_{12}}} \right) \times \left( {{b_1}{x_{21}} + {b_2}{x_{22}}} \right) \times \left( {{b_1}{x_{31}} + {b_2}{x_{32}}} \right) \cr & \max z = \left( {{b_1}{x_{11}}{b_1}{x_{21}} + {b_1}{x_{11}}{b_2}{x_{22}} + {b_2}{x_{12}}{b_1}{x_{21}} + {b_2}{x_{12}}{b_2}{x_{22}}} \right) \times \left( {{b_1}{x_{31}} + {b_2}{x_{32}}} \right) \cr & \max z = {b_1}{x_{11}}{b_1}{x_{21}}{b_1}{x_{31}} + {b_1}{x_{11}}{b_1}{x_{21}}{b_2}{x_{32}} + {b_1}{x_{11}}{b_2}{x_{22}}{b_1}{x_{31}} + {b_1}{x_{11}}{b_2}{x_{22}}{b_2}{x_{32}} + {b_2}{x_{12}}{b_1}{x_{21}}{b_1}{x_{31}} + {b_2}{x_{12}}{b_1}{x_{21}}{b_2}{x_{32}} + {b_2}{x_{12}}{b_2}{x_{22}}{b_1}{x_{31}} + {b_2}{x_{12}}{b_2}{x_{22}}{b_2}{x_{32}} \cr}

I think this can lead to a solution that both b_1 and b_2 take values between 0 and 1.

[Shayan]

Hi @mtanneau,
Sorry for the late response. I’ve not visited Discourse until now.

@mtanneau

maximizing a product is equivalent to maximizing its log, which yields a sum of smaller log terms

So cool! This is great. I got it. I would try the rewritten model by utilizing Hypatia. I let you know the result of the experiment!

[Shayan]

I rewrote the model as follows:

using Hypatia

model = Model(Hypatia.Optimizer)
@variables(model, begin
    0<=b[i=1:n]<=1
end)
@constraint(model, sum(b) == 1)
@expression(model, zᵢ, log.(b'*x[:, rnd_idxs]))
@NLconstraint(model, t[i=1:length(zᵢ)], t .<= zᵢ)
@NLobjective(model, Max, sum(t))
optimize!(model)

But, it throws:

ERROR: log is not defined for type GenericAffExpr. Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective. If you are using an `@NL` macro and you encountered this error message, it is because you are attempting to use another unsupported function which calls this method internally.
Stacktrace:

Any thoughts? @odow

[odow]

You need to use a conic formulation:

using JuMP
import SCS
import Ipopt

function model_prod(n, x)
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @variable(model, 0 <= b[1:n] <= 1)
    @constraint(model, sum(b) == 1)
    @expression(model, h, b' * x)
    @NLobjective(model, Max, prod(h[i] for i in 1:m))
    optimize!(model)
    return value.(b)
end

function model_sumlog(n, x)
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @variable(model, 0 <= b[1:n] <= 1)
    @constraint(model, sum(b) == 1)
    @expression(model, h, b' * x)
    @NLobjective(model, Max, sum(log(h[i]) for i in 1:m))
    optimize!(model)
    return value.(b)
end

function model_conic(n, x)
    m = size(x, 2)
    model = Model(SCS.Optimizer)
    set_silent(model)
    @variable(model, 0 <= b[1:n] <= 1)
    @variable(model, t[1:m])
    @constraint(model, sum(b) == 1)
    @expression(model, h, b' * x)
    @constraint(model, [i=1:m], [t[i], 1, h[i]] in MOI.ExponentialCone())
    @objective(model, Max, sum(t))
    optimize!(model)
    return value.(b)
end

function run_model(model_f)
    n = m = 3 
    data = rand(n, 10);
    weights = zeros(n, m)
    for iter in 1:3
        rnd_idxs = rand(1:3, m)
        x = data[:, rnd_idxs]
        weights[:, iter] .= model_f(n, x)
    end
    return weights
end

run_model(model_prod)
run_model(model_sumlog)
run_model(model_conic)

[Shayan]

Thank you so much. You’ve omitted import MathOptInterface as MOI accidentally, I guess. The conic model is so impressive, but I’ve to search for its details.

@odow

@constraint(model, [i=1:m], [t[i], 1, h[i]] in MOI.ExponentialCone())

Is the logarithm taken here?

[odow]

You’ve omitted import MathOptInterface as MOI accidentally, I guess

It gets added to the scope on using JuMP. But yes, you can also use import MathOptInterface as MOI.

but I’ve to search for its details

https://docs.mosek.com/modeling-cookbook/expo.html#logarithm

Is the logarithm taken here?

Yes. I’m adding examples to the JuMP documentation:

• [docs] add more exponential cone tricks by odow · Pull Request #3345 · jump-dev/JuMP.jl · GitHub
• [docs] fix math in tips_and_tricks.jl by odow · Pull Request #3347 · jump-dev/JuMP.jl · GitHub

[Shayan]

Thank you so much for your time and help
It’s great that you keep the doc updated. Thank you for all of these

[odow]

For future reference, here are the docs: Tips and Tricks · JuMP