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.

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.

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

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