How can I significantly improve the speed of this linear programming problem in JuMP?

Specific Domains Optimization (Mathematical)
,
1

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.

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.
https://github.com/JuliaOpt/Gurobi.jl/issues/279#issue-535364203

2

Hello,
Have you tested with Interior Point Method?

Regards

3

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.

4

If you haven’t yet, I recommend using Julia’s built-in statistical profiler to find out where the time is being spent. You can use ProfileView or StatProfilerHTML to visualize the results.

1 Like
5

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.

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.

2 Likes
6

Hi tkluck,

Thank you for your comment. I will try.

7

Hi odow,

That’s impressive. I will try and let you know. Thank you so much.

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

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.

1 Like
9

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?

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

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
10

You need

@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

model[:con_ub]

instead of just con_ub.

If you don’t know what a “scope” is, you should read the manual, 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.) Are you sure it has to be <= and not ==? In an optimal solution, that inequality should hold with a strict equality.

11

Thanks. In case you miss the reply I posted earlier.

Do you know why when I add the following code, the time is doubled?

    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
12

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.

13

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.