Gurobi is slow with MOI, likely due to high overhead of model updates

[evanfields]

Here's code for a small working example. The problem is a network flow problem over a complete graph: n points are generated in the unit square, corresponding to vertices of the graph. Vectors x and y of length n and sum 1 correspond to vertex-wise supply and demand. Script:

using JuMP, GLPK, Gurobi, Random

# build a network flow JuMP model
# n (dense) constraints, n^2 variables
function build_jumpmodel(x, y, dm, optimizer)
    n = length(x)
	@assert n == length(y)

	m = JuMP.Model(with_optimizer(optimizer))
	@variable(m, f[1:n, 1:n] >= 0)
	@objective(m, Min, sum(dm[i,j] * f[i,j] for i in 1:n, j in 1:n))


	# correct net change at each voxel to transform x to y
	for i in 1:n
		@constraint(m,
            sum(f[j,i] for j in 1:n)
            - sum(f[i,j] for j in 1:n)
            == y[i] - x[i]
        )
	end
    
    return m
end

function random_network_flow(n, optimizer)
    Random.seed!(123)
    # n random vertices in the unit square
    # random weight vectors, each summing to 1
    vert_locs = rand(2,n)
    dm = zeros(n,n)
    for i in 1:n, j in 1:n
        dm[i,j] = (vert_locs[:,i] - vert_locs[:,j]).^2 |> sum |> sqrt
    end
    x = rand(n); x ./= sum(x)
    y = rand(n); y ./= sum(y)
    
    build_time = @elapsed m = build_jumpmodel(x, y, dm, optimizer)
    @show build_time
    
    solve_time = @elapsed JuMP.optimize!(m)
    println()
    @show solve_time
    
    return JuMP.objective_value(m)
end

Example output:

julia> @profile random_network_flow(150, Gurobi.Optimizer)
Academic license - for non-commercial use only
build_time = 0.288749231

Warning: excessive time spent in model updates.
Consider calling update less frequently.

Optimize a model with 150 rows, 22500 columns and 44700 nonzeros
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [8e-03, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [4e-05, 1e-02]

Concurrent LP optimizer: primal simplex, dual simplex, and barrier
Showing barrier log only...

Presolve removed 1 rows and 150 columns
Presolve time: 0.02s
Presolved: 149 rows, 22350 columns, 44402 nonzeros

Ordering time: 0.00s

Barrier statistics:
 AA' NZ     : 1.103e+04
 Factor NZ  : 1.118e+04 (roughly 9 MBytes of memory)
 Factor Ops : 1.114e+06 (less than 1 second per iteration)
 Threads    : 4

Barrier performed 0 iterations in 0.04 seconds
Barrier solve interrupted - model solved by another algorithm


Solved with primal simplex
Solved in 1445 iterations and 0.04 seconds
Optimal objective  3.122653030e-02

solve_time = 2.96464323
0.031226530301749666

# for context, using GLPK
julia> random_network_flow(150, GLPK.Optimizer)
build_time = 0.112245802

solve_time = 0.216538412
0.031226530301749662

Notice that the overall solve time with Gurobi is order of magnitude worse than with GLPK, but of the ~3 seconds needed for JuMP.optimize! with Gurobi, only ~.04 seconds are spent in Gurobi! Gurobi also warns that excessive time is spent in model updates.

Basic profiling seems a little unreliable here, or I am bad at reading it - reported line numbers and functions don't always agree. That said, looks like lots of time is spent in update_model! of grb_model.jl.

Edit: forgot to mention - this is with Gurobi.jl v0.5.8, LQOI v0.6.0, Gurobi the solver 8.1.

[odow]

looks like lots of time is spent in update_model! of grb_model.jl.

Yes, this will be the culprit. Basically it's because the current wrapper calls update after any mutating operation. Previous versions of JuMP used to call these infrequently, for example by setting all variable lower bounds in the same call, whereas now things get called much more frequently.

There is room in LQOI for improvements (e.g. JuliaOpt/LinQuadOptInterface.jl#73), but given GLPK also uses LQOI that's not the biggest issue.

A very similar issue was previously raised, and my Gurobi specific comments remain valid:
JuliaOpt/LinQuadOptInterface.jl#51 (comment)

The main options are

• changing JuMP to cache calls (for example, calling add_variables(model, N) instead of add_variable(model) N times)
• being clever about when we call update (for example, only get type operations, and only if data has changed)
• being dumb about when we call update (for example, only before optimize!) and forcing the user to instead

Thoughts?

[joaquimg]

changing JuMP to cache calls (for example, calling add_variables(model, N) instead of add_variable(model) N times)

is theoretically possible with CachingOptimizer, but in DIRECT mode it won't work. One options is to allowing the user to pre-allocate a bunch of variables with MOI and the add_variable calls just creates the bindings...

being clever about when we call update (for example, only get type operations, and only if data has changed)

This is the most immediate solution but cumbersome. Someone will ha to sit for a few hours figuring out the right places to call updates.

being dumb about when we call update (for example, only before optimize!) and forcing the user to instead

This could be some sort of argument for gurobi/LQOI that triggers this (non-default) unsafe mode

[odow]

Someone will ha to sit for a few hours figuring out the right places to call updates.

Smartly calling the updates is a pain, because you need to do things like call update after adding a variable in order to use it in a constraint. If the sequence of things is

• add variable
• set bound
• add variable
• set bound
  we're kind of forced to call update after each add variable, although we can hold off after the bound updates... We don't need to call it after adding a constraint or modifying an attribute until a get-type call is made.

This is the reason I just made it "call update after any mutating operation."

This could be some sort of argument for gurobi/LQOI that triggers this (non-default) unsafe mode

Yes, I was thinking of some keyword argument Gurobi.Optimizer(auto_update = true).

[mlubin]

The implementation of copy_to should handle this reasonably. When the whole model can be loaded at once, there's no need to call it after every operation (I haven't checked if this is the case).

Also you don't need to call update after adding variables: https://www.gurobi.com/documentation/8.1/refman/updatemode.html#parameter:UpdateMode