# Pi (Π) operation in the objective function JuMP

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)

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)

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.

1 Like

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

Thank you! Does mine do it wrongly?

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.

1 Like

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

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.

1 Like

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.

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.

1 Like

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

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

2 Likes

(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.

2 Likes

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.

1 Like

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

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.

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.

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

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!

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

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)
1 Like

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.

Is the logarithm taken here?

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:

2 Likes

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

1 Like

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

1 Like