# How can I significantly improve the speed of this linear programming problem in JuMP?
**URL:** https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125
**Category:** Optimization (Mathematical)
**Tags:** question, jump, linear-programming
**Created:** [December 10, 2019, 11:14pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125 "2019-12-10T23:14:10Z")
**Posts on this page:** 13
**Page:** 1
### Author:  [@zxjroger](https://discourse.julialang.org/u/zxjroger)
#### Post date: [December 10, 2019, 11:14pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/1 "2019-12-10T23:14:10Z")
I am doing Bayesian estimation in Julia. For each iteration of Bayesian estimation, there is a for-loop to compute the bounds of some variables in my model. Each bound is computed by a linear programming problem. I am trying to see if the LP problem can be solved in milliseconds.
This is the optimization problem for upper bound.
\begin{align} \max\_{W\_{i^\star j^\star},\lambda}\qquad&W\_{i^\star j^\star}\\ \text{s.t.}\qquad &W\_{ij}=\lambda\_i^1+\lambda\_j^2, \quad\forall ij \quad\text{s.t.}\quad H\_{ij}=1\\ &W\_{ij}\leq\lambda\_i^1+\lambda\_j^2, \quad\forall ij \quad\text{s.t.}\quad H\_{ij}=0 \end{align}
The optimization problem for lower bound is similar, which just replaces \max with \min. Here, the objective is a scalar variable W\_{i^\star j^\star}. Other variables are \lambda\_i^1 and \lambda\_j^2 where i\in\{1,...,N\} and j\in\{1,...,N\}. All W\_{ij} except for W\_{i^\star j^\star} are known constants. H is a known and constant matrix such that its element H\_{ij} is either 0 or 1. In summary, I need to compute the bound of W\_{ij} for each ij one by one in the for-loop.
Let me attach the simplified code I have. First, data on W and H are generated as follows.
```julia-auto
using JuMP
using Gurobi
# Generate fake data and compute matrix H
N = 500
Random.seed!(0)
W = randn(N,N) # Data on W
model_1 = Model(with_optimizer(Gurobi.Optimizer, OutputFlag=0))
@variable(model_1, H_temp[1:N, 1:N] >= 0)
@expression(model_1, row_sum[i = 1:N], sum(H_temp[i, j] for j = 1:N))
@expression(model_1, column_sum[j = 1:N], sum(H_temp[i, j] for i = 1:N))
@constraint(model_1, row_constraint[i = 1:N], row_sum[i] == 1)
@constraint(model_1, column_constraint[j = 1:N], column_sum[j] == 1)
@objective(model_1, Max, sum(W[i, j]*H_temp[i, j] for i = 1:N, j = 1:N))
optimize!(model_1)
if termination_status(model_1) == MOI.OPTIMAL || termination_status(model_1) == MOI.LOCALLY_SOLVED
H = value.(H_temp)
optimal_objective = objective_value(model_1)
else
error("The model_1 was not solved correctly.")
end
H_star = round.(Int, H) # Data on H
```
My main question is how I can improve the speed of the following code, which is to solve for the bounds.
```julia-auto
for i = 1:N # One can change this N to 10 to test the speed
for j = 1:N # One can change this N to 1 to test the speed (so that only ten problems are solved)
if H_star[i,j] == 1 # Compute the lower bound
model_lb = JuMP.direct_model(Gurobi.Optimizer(Presolve=1, Method=0)) # Method=0 is primal simplex
@variable(model_lb, W_lb)
@variable(model_lb, lambda_1[1:N])
@variable(model_lb, lambda_2[1:N])
for i_lb = 1:N
for j_lb = 1:N
if i_lb == i && j_lb == j
@constraint(model_lb, W_lb == lambda_1[i_lb] + lambda_2[j_lb])
elseif H_star[i_lb,j_lb] == 1
@constraint(model_lb, W[i_lb,j_lb] == lambda_1[i_lb] + lambda_2[j_lb])
else H_star[i_lb,j_lb] == 0
@constraint(model_lb, W[i_lb,j_lb] <= lambda_1[i_lb] + lambda_2[j_lb])
end
end
end
@objective(model_lb, Min, W_lb)
optimize!(model_lb)
if termination_status(model_lb) == MOI.OPTIMAL || termination_status(model_lb) == MOI.LOCALLY_SOLVED
W_bound = value.(W_lb)
println("W_lower_bound($i,$j)=$W_bound")
else
error("The model_lb was not solved correctly.")
end
else # H[i,j] == 0
# Compute the upper bound
model_ub = JuMP.direct_model(Gurobi.Optimizer(Presolve=1, Method=0)) # Method=0 is primal simplex
@variable(model_ub, W_ub)
@variable(model_ub, lambda_1[1:N])
@variable(model_ub, lambda_2[1:N])
for i_ub = 1:N
for j_ub = 1:N
if i_ub == i && j_ub == j
@constraint(model_ub, W_ub <= lambda_1[i_ub] + lambda_2[j_ub])
elseif H_star[i_ub,j_ub] == 1
@constraint(model_ub, W[i_ub,j_ub] == lambda_1[i_ub] + lambda_2[j_ub])
else H_star[i_ub,j_ub] == 0
@constraint(model_ub, W[i_ub,j_ub] <= lambda_1[i_ub] + lambda_2[j_ub])
end
end
end
@objective(model_ub, Max, W_ub)
optimize!(model_ub)
if termination_status(model_ub) == MOI.OPTIMAL || termination_status(model_ub) == MOI.LOCALLY_SOLVED
W_bound = value.(W_ub)
println("W_upper_bound($i,$j)=$W_bound")
else
error("The model_ub was not solved correctly.")
end
end
end
end
```
Here in the code, `W_ub` and `W_lb` correspond to W\_{i^\star j^\star} in the optimization problem, and `H_star` corresponds to H in the optimization problem. Since N=500, for each iteration of my Bayesian estimation, I will need to solve the optimization problem for N^2=250000 times. This is why I need to make the computation as fast as possible. Right now, I set the `Method` of optimization to be 0 which is primal simplex and `Presolve` level to be 1. On average, it takes about 1.4 seconds to solve one problem, _i.e._, one bound say W\_{ij}. I would like to know if I can reduce this time to milliseconds. Thank you very much.
* * *
This is the back-story if you are interested. I thank Oscar Dowson for his tremendous help.
[https://github.com/JuliaOpt/Gurobi.jl/issues/279#issue-535364203](https://github.com/JuliaOpt/Gurobi.jl/issues/279#issue-535364203)
---
### Author:  [@bienpierre](https://discourse.julialang.org/u/bienpierre)
#### Post date: [December 12, 2019, 7:04am UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/2 "2019-12-12T07:04:36Z")
Hello,
Have you tested with Interior Point Method?
> [@zxjroger](#):
>
> primal simplex
Regards
---
### Author:  [@zxjroger](https://discourse.julialang.org/u/zxjroger)
#### Post date: [December 13, 2019, 3:55pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/3 "2019-12-13T15:55:53Z")
Thank you for your comment. I tried all optimization methods available in Gurobi, and primal simplex seems to be the best one in term of speed.
---
### Author:  [@tkluck](https://discourse.julialang.org/u/tkluck)
#### Post date: [December 13, 2019, 4:07pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/4 "2019-12-13T16:07:42Z")
If you haven’t yet, I recommend using Julia’s [built-in statistical profiler](https://docs.julialang.org/en/v1/manual/profile/#Profiling-1) to find out where the time is being spent. You can use [ProfileView](https://github.com/timholy/ProfileView.jl) or [StatProfilerHTML](https://github.com/tkluck/StatProfilerHTML.jl) to visualize the results.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [December 13, 2019, 4:17pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/5 "2019-12-13T16:17:26Z")
The issue is that you are rebuilding the model from scratch every time. You can do much better by modifying some of the variable bounds.
```nohighlight
function create_lower_bound_model(N::Int)
model_lb = JuMP.direct_model(Gurobi.Optimizer(OutputFlag = 0))
@variable(model_lb, lambda_1[1:N])
@variable(model_lb, lambda_2[1:N])
@variable(model_lb, lambda_sum[1:N, 1:N])
@constraint(model_lb, [i=1:N, j=1:N], lambda_sum[i, j] == lambda_1[i] + lambda_2[j])
return model_lb
end
function solve_lower_bound(
model_lb::Model, i::Int, j::Int, W::Matrix{Float64}, H_star::Matrix{Int}
)
@objective(model_lb, Min, model_lb[:lambda_sum][i, j])
for i_lb = 1:N, j_lb = 1:N
if i_lb == i && j_lb == j
continue
end
x = model_lb[:lambda_sum][i_lb, j_lb]
set_lower_bound(x, W[i_lb, j_lb])
if H_star[i_lb, j_lb] == 1
set_upper_bound(x, W[i_lb, j_lb])
else
@assert H_star[i_lb, j_lb] == 0
if has_upper_bound(x)
delete_upper_bound(x, Inf)
end
end
end
optimize!(model_lb)
if termination_status(model_lb) == MOI.OPTIMAL
println("W_lower_bound($i, $j) = ", objective_value(model_lb))
else
error("The model_lb was not solved correctly.")
end
end
model_lb = create_lower_bound_model(size(W, 1))
for i = 1:N, j = 1:N
if H_star[i, j] == 1
solve_lower_bound(model_lb, i, j, W, H_star)
else
# TODO
end
end
```
You can do something similar with the upper bound.
---
### Author:  [@zxjroger](https://discourse.julialang.org/u/zxjroger)
#### Post date: [December 13, 2019, 4:19pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/6 "2019-12-13T16:19:38Z")
Hi tkluck,
Thank you for your comment. I will try.
---
### Author:  [@zxjroger](https://discourse.julialang.org/u/zxjroger)
#### Post date: [December 13, 2019, 4:26pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/7 "2019-12-13T16:26:49Z")
Hi odow,
That’s impressive. I will try and let you know. Thank you so much.
---
### Author:  [@zxjroger](https://discourse.julialang.org/u/zxjroger)
#### Post date: [December 13, 2019, 8:09pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/8 "2019-12-13T20:09:22Z")
Hi odow,
Your code indeed improves the speed a lot. I really appreciate it. However, there is a glitch in function `solve_lower_bound`. After I add the following code, the speed is basically the same as my original code, is there another way to write this code?
```julia
if has_upper_bound(model_lb[:lambda_sum][i, j])
delete_upper_bound(model_lb[:lambda_sum][i, j])
end
if has_lower_bound(model_lb[:lambda_sum][i, j])
delete_lower_bound(model_lb[:lambda_sum][i, j])
end
```
Without the above code, it solves 10 problem for N=500 in 9 seconds, after I add the above code, the time is nearly doubled. I am not sure why the new code will take so much time.
The reason I need the above code is as follows. For one iteration, it puts bounds on all `lambda_sum[i_lb, j_lb]` except for `lambda_sum[i, j]`. In the next iteration, say that we are going to solve for `W[i',j']`, so that there should be no bounds on `W[i',j']`. However, the bounds put on `W[i',j']` in the previous iteration are still there. Therefore, I add the following codes to clear bounds before the for-loop in function `solve_lower_bound`.
---
### Author:  [@zxjroger](https://discourse.julialang.org/u/zxjroger)
#### Post date: [December 13, 2019, 11:41pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/9 "2019-12-13T23:41:46Z")
Hi odow,
I wrote the code for upper bound, there is a difference since now I need to change the constraint for i^\star j^\star from equality to inequality. When I run the code within two functions by hand, I got correct answer. When I run the function in the for-loop, however, I got an error saying that `The constraint reference you are trying to delete does not belong to the model.` It seems that Julia did not evaluate the code within the function. Do you no why this is the case?
```julia
function create_upper_bound_model(N::Int)
model_ub = JuMP.direct_model(Gurobi.Optimizer(OutputFlag = 0))
@variable(model_ub, lambda_1[1:N])
@variable(model_ub, lambda_2[1:N])
@variable(model_ub, lambda_sum[1:N, 1:N])
# index the constraint by con_ub[i,j]
con_ub = @constraint(model_ub, [i=1:N, j=1:N], lambda_sum[i, j] == lambda_1[i] + lambda_2[j])
return model_ub
end
```
```julia
function solve_upper_bound(
model_ub::Model, i::Int, j::Int, W::Matrix{Float64}, H_star::Matrix{Int}
)
# change the equality to inequality
delete(model_ub, con_ub[i, j])
con_ub[i, j] = @constraint(model_ub, lambda_sum[i, j] <= lambda_1[i] + lambda_2[j])
@objective(model_ub, Max, model_ub[:lambda_sum][i, j])
if has_upper_bound(model_ub[:lambda_sum][i, j])
delete_upper_bound(model_ub[:lambda_sum][i, j])
end
if has_lower_bound(model_ub[:lambda_sum][i, j])
delete_lower_bound(model_ub[:lambda_sum][i, j])
end
for i_ub = 1:N, j_ub = 1:N
if i_ub == i && j_ub == j
continue
end
x = model_ub[:lambda_sum][i_ub, j_ub]
set_lower_bound(x, W[i_ub, j_ub])
if H_star[i_ub, j_ub] == 1
set_upper_bound(x, W[i_ub, j_ub])
else
@assert H_star[i_ub, j_ub] == 0
if has_upper_bound(x)
delete_upper_bound(x)
end
end
end
optimize!(model_ub)
if termination_status(model_ub) == MOI.OPTIMAL
println("W_upper_bound($i, $j) = ", objective_value(model_ub))
else
error("The model_ub was not solved correctly.")
end
# put back the original constraint for next iteration
delete(model_ub, con_ub[i, j])
con_ub[i, j] = @constraint(model_ub, lambda_sum[i, j] == lambda_1[i] + lambda_2[j])
end
```
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [December 14, 2019, 7:43pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/10 "2019-12-14T19:43:27Z")
You need
```nohighlight
@constriant(model_ub, con_ub[i=1:N, j=1:N], ...)
```
in the first function, and then in the `solve_upper_bound` function, you need to reference objects created in the first function with
```nohighlight
model[:con_ub]
```
instead of just `con_ub`.
If you don’t know what a “scope” is, you should read [the manual](https://docs.julialang.org/en/v1/manual/variables-and-scoping/#Local-Scope-1), which explains which Julia variables can be “seen” by which functions.
Also note that adding and deleting constraints like this is going to cause issues. (In particular, [https://github.com/JuliaOpt/JuMP.jl/issues/1956.](https://github.com/JuliaOpt/JuMP.jl/issues/1956.)) Are you sure it has to be `<=` and not `==`? In an optimal solution, that inequality should hold with a strict equality.
---
### Author:  [@zxjroger](https://discourse.julialang.org/u/zxjroger)
#### Post date: [December 14, 2019, 8:07pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/11 "2019-12-14T20:07:09Z")
Thanks. In case you miss the reply I posted earlier.
> [@How can I significantly improve the speed of this linear programming problem in JuMP?](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/8):
>
> Hi odow, Your code indeed improves the speed a lot. I really appreciate it. However, there is a glitch in function solve\_lower\_bound. After I add the following code, the speed is basically the same as my original code, is there another way to write this code? if has\_upper\_bound(model\_lb[:lambda\_sum][i, j]) delete\_upper\_bound(model\_lb[:lambda\_sum][i, j]) end if has\_lower\_bound(model\_lb[:lambda\_sum][i, j]) delete\_lower\_bound(model\_lb[:lambda\_sum][i, j]) end With…
Do you know why when I add the following code, the time is doubled?
```julia
if has_upper_bound(model_lb[:lambda_sum][i, j])
delete_upper_bound(model_lb[:lambda_sum][i, j])
end
if has_lower_bound(model_lb[:lambda_sum][i, j])
delete_lower_bound(model_lb[:lambda_sum][i, j])
end
```
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [December 14, 2019, 8:16pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/12 "2019-12-14T20:16:12Z")
Is the solve time (i.e., the time spend in Gurobi) doubling? Or the build time?
If you remove bounds, you probably make the model much more difficult to solve.
---
### Author:  [@zxjroger](https://discourse.julialang.org/u/zxjroger)
#### Post date: [December 15, 2019, 7:34pm UTC](https://discourse.julialang.org/t/how-can-i-significantly-improve-the-speed-of-this-linear-programming-problem-in-jump/32125/13 "2019-12-15T19:34:11Z")
It seems that it is more difficult to solve. Your code does not require construct the model for each iteration, so I thought it should save some time compared to my original model.