Slow model construction in JuMP

jacob-roth · March 2, 2023, 4:45pm

I find that when I am modeling a problem in JuMP, the model construction takes a long time. My first thought was that maybe I am writing inefficient code, and I was wondering:

if there might be any quick fixes to my JuMP code
if interfacing with MathOptInterface directly would be faster
if interfacing with Gurobi (or whatever solver I use) directly would be faster

using JuMP, Gurobi
m = 500
n = 100_000
A = randn(m, n)
b = randn(m)
lb = fill(-1.0, n)
ub = fill(+1.0, n)
lb[1] = -Inf
ub[2] = +Inf
julia> @time begin
  model = Model(Gurobi.Optimizer)
  @variable(model, x[1:n])
  @variable(model, y[1:m])
  @variable(model, z[1:m])
  @variable(model, t)
  @objective(model, Min, sum(x[i]^2 for i in 1:n)) # can replace with no objective and still slow
  @constraint(model, y .== A * x .+ b)
  @constraint(model, z .>= t .+ y)
  @constraint(model, sum(z) .>= t)
  @constraint(model, lb .<= x .<= ub)
end                                                                                                                                          
Set parameter Username
Academic license - for non-commercial use only - expires 2023-12-17
GC: pause 67.38ms. collected 3342.926024MB. incr
GC: pause 396.65ms. collected 970.717568MB. full
GC: pause 429.67ms. collected 2750.147009MB. full
481.022531 seconds (3.94 M allocations: 12.238 GiB, 0.19% gc time, 0.02% compilation time)

odow · March 2, 2023, 8:22pm

There are a couple of minor changes you can do to speed things up (using n=10_000, not 100_000):

julia> GC.gc()

julia> @time begin
           model = Model(Gurobi.Optimizer)
           set_silent(model)
           @variable(model, x[1:n])
           @variable(model, y[1:m])
           @variable(model, z[1:m])
           @variable(model, t)
           @objective(model, Min, sum(x[i]^2 for i in 1:n))
           @constraint(model, y .== A * x .+ b)
           @constraint(model, z .>= t .+ y)
           @constraint(model, sum(z) .>= t)
           @constraint(model, lb .<= x .<= ub)
           optimize!(model)
       end
  7.960642 seconds (1.17 M allocations: 1.597 GiB, 1.20% gc time, 1.44% compilation time)

julia> GC.gc()

julia> @time begin
           model = direct_model(Gurobi.Optimizer())
           set_silent(model)
           set_string_names_on_creation(model, false)
           @variable(model, lb[i] <= x[i=1:n] <= ub[i])
           @variable(model, y[1:m])
           @variable(model, z[1:m])
           @variable(model, t)
           @objective(model, Min, sum(x[i]^2 for i in 1:n))
           @constraint(model, [i=1:m], y[i] == sum(A[i, j] * x[j] for j in 1:n) + b[i])
           @constraint(model, [i=1:m], z[i] >= t + y[i])
           @constraint(model, sum(z) >= t)
           optimize!(model)
       end
  7.030907 seconds (25.37 M allocations: 1.170 GiB, 2.17% gc time, 3.04% compilation time)

For example, @constraint(model, lb .<= x .<= ub) is actually a ScalarAffineFunction-in-Interval constraint, which JuMP needs to bridge. You’re much better off passing the bounds explicitly in @variable.

But you’ll notice that it doesn’t do very much.

The real problem is that your A matrix is dense.

500 * 100_000 is 50 million non-zeros. Gurobi is built with the assumption that the constraint matrix is sparse, typically with 1% or fewer non-zeros.

Something seems wrong with the formulation of your model though. It doesn’t make sense. x = 0, and t = very large negative number is the optimal solution?

What is the problem you are actually trying to solve?

jacob-roth · March 2, 2023, 10:45pm

Thanks again for your comments:)

You’re right this is not the real model. The problem I’m considering is a linear-reformulation of a CVaR constraint on y=Ax+b, and I was sloppy about the details here, though it has the main idea. As a first step, I’m testing some different solvers and was comparing them on randomly generated (dense) instances of different sizes.

I’m not at the point of testing on real problems yet with actual A so was using a dense A of a certain size as a rough (bad?) proxy for performance on even larger (but sparse) A matrices.

Do you think it will be faster to work directly in Gurobi’s interface, i.e., GRBnewmodel…?

odow · March 2, 2023, 11:06pm

Using random unstructured dense data as a proxy isn’t very useful. Real problems have sparse structure, and solvers are built to exploit that.

Going via the C API is probably a little faster, but you’ll still have the same issue. You can’t pass a dense matrix to Gurobi, you have to pass the row and column indices, as well as a vector of data. So you end up passing a a few 5e7 length vectors, and the data structures and algorithms inside Gurobi assume sparsity.

For most realistic problems, the cost of building the JuMP model is (much) less than the cost of solving it, so if you’re struggling to build it, Gurobi will likely struggle to solve it.

You also have to factor in the time it takes you to write the C API version, as well as the likelihood of bugs etc . In general, you should probably just use JuMP, and only once all else fails should you consider alternatives.

jacob-roth · March 3, 2023, 1:30am

I tried making/solving an LP from data A,b,c for the problem \max\{c^\top x : Ax \leq b, x \geq 0\} with the GRBnewmodel, GRBaddvars, and GRBaddconstr functions in the dense case (just for fun). To verify my construction, I compare against solving the same problem in JuMP. My code below shows that I am making an error in constructing the problem using the C interface. So I am asking Gurobi to see if they have an example of using the GRBnewmodel, GRBaddvars, and GRBaddconstr functions (or if they can diagnose my below code). I’m posting here as well just in case you’re interested

For m=500 and n=10_000, it seems ~10x faster to use Gurobi directly.

#=
maximize    c' * x
subject to  A * x <= b
                x >= 0
=#

using Gurobi
using JuMP

function solve_lp_direct(A::Matrix{Tf}, b::Vector{Tf}, c::Vector{Tf}) where Tf<:Cdouble
  setup_time = @elapsed begin

    # initialize model
    start_time = time()
    env_p = Ref{Ptr{Cvoid}}()
    error = GRBloadenv(env_p, "lp.log")
    env = env_p[]
    model_p = Ref{Ptr{Cvoid}}()
    error = GRBnewmodel(
      env,
      model_p,
      "lp",
      0,
      C_NULL,
      C_NULL,
      C_NULL,
      C_NULL,
      C_NULL
    )
    model = model_p[]
    m, n = size(A)

    # maximize
    error = GRBsetintattr(model, GRB_INT_ATTR_MODELSENSE, GRB_MAXIMIZE)

    # output
    optimstatus = Ref{Cint}()
    objval = Ref{Cdouble}()
    sol = ones(Tf, n)

    # variables and objective
    error = GRBaddvars(
      model,         # model
      n,             # : numvars
      0,             # : numnz
      C_NULL,        # : *vbeg
      C_NULL,        # : *vind
      C_NULL,        # : *vval
      c,             # : *obj
      zeros(Cdouble, n), # : *lb
      fill(+Inf, n), # : *ub
      C_NULL,        # : *vtype
      C_NULL         # : **varnames
    )

    # constraint
    for i in 1:size(A,1)
      Ai = A[i,:]
      cind = findall(Ai .!= 0.0)
      numnz = Cint(length(cind))
      cval = Ai[cind]
      cind .-= 1
      cind = Cint.(cind)
      sense = GRB_LESS_EQUAL
      rhs = b[i]
      error = GRBaddconstr(
        model,  # : *model
        numnz,  # : numnz
        cind,   # : *cind
        cval,   # : *cval
        sense,  # : sense
        rhs,    # : rhs
        C_NULL, # : *constrname
      )
      @assert(error == 0)
    end
  end # begin
  println("model setup = $(setup_time)")
  solve_time = @elapsed error = GRBoptimize(model)
  println("model solve = $(solve_time)")
  # error = GRBwrite(model, "lp.lp");
  error = GRBgetintattr(model, GRB_INT_ATTR_STATUS, optimstatus);
  error = GRBgetdblattr(model, GRB_DBL_ATTR_OBJVAL, objval);
  error = GRBgetdblattrarray(model, GRB_DBL_ATTR_X, 0, n, sol);
  return sol
end

function solve_lp_jump(A::Matrix{Tf}, b::Vector{Tf}, c::Vector{Tf}) where Tf<:Cdouble
  setup_time = @elapsed begin
    m, n = size(A)
    M = Model(Gurobi.Optimizer)
    @variable(M, x[1:n] .>= 0)
    @objective(M, Max, c'*x)
    @constraint(M, A*x .<= b)
  end
  println("model setup = $setup_time")
  solve_time = optimize!(M)
  println("model solve = $(solve_time)")
  return value.(x)
end

m = 50
n = 100
A = randn(Cdouble, m, n)
b = randn(Cdouble, m)
c = randn(Cdouble, n)
@time sol1 = solve_lp_direct(copy(A), copy(b), copy(c));
@time sol2 = solve_lp_jump(copy(A), copy(b), copy(c));
sum(abs.(sol1-sol2))

odow · March 3, 2023, 1:33am

Let me take a look.

odow · March 3, 2023, 2:10am

It’s probably because you’re not freeing the Gurobi memory at the end. I could trigger the segfault with your original code, but not with my version below. Working with the C API from Julia can be difficult. Gurobi.jl does stuff to ensure that we don’t free memory at the wrong time:

github.com

jump-dev/Gurobi.jl/blob/8c98d0d9988c67005bddfc2cf5b0365357551e5d/src/MOI_wrapper/MOI_wrapper.jl#L120-L131


      
              finalizer(env) do e
                  e.finalize_called = true
                  if e.attached_models == 0
                      # Only finalize the model if there are no models using it.
                      GRBfreeenv(e.ptr_env)
                      e.ptr_env = C_NULL
                  end
              end
              # Even if the loadenv fails, the pointer is still valid.
              _check_ret(env, ret)
              return env
          end

Here’s the code I ran:

using Gurobi
using JuMP

function solve_lp_direct(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    start_time = time()
    env_p = Ref{Ptr{Cvoid}}()
    error = GRBloadenv(env_p, "lp.log")
    env = env_p[]
    error = GRBsetintparam(env, GRB_INT_PAR_OUTPUTFLAG, 0)
    @assert error == 0
    model_p = Ref{Ptr{Cvoid}}()
    error = GRBnewmodel(env, model_p, "lp", 0, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL)
    @assert error == 0
    model = model_p[]
    m, n = size(A)
    error = GRBsetintattr(model, GRB_INT_ATTR_MODELSENSE, GRB_MAXIMIZE)
    @assert error == 0
    error = GRBaddvars(
        model,         # model
        n,             # : numvars
        0,             # : numnz
        C_NULL,        # : *vbeg
        C_NULL,        # : *vind
        C_NULL,        # : *vval
        c,             # : *obj
        zeros(Cdouble, n), # : *lb
        fill(Inf, n),  # : *ub
        C_NULL,        # : *vtype
        C_NULL         # : **varnames
    )
    @assert error == 0
    for i in 1:m
        error = GRBaddconstr(
            model,              # : *model
            n,                  # : numnz
            Cint.(0:(n-1)),     # : *cind
            A[i, :],            # : *cval
            GRB_LESS_EQUAL,     # : sense
            b[i],               # : rhs
            C_NULL,             # : *constrname
        )
        @assert error == 0
    end
    println("model setup = $(time() - start_time)")
    solve_time = @elapsed error = GRBoptimize(model)
    println("model solve = $solve_time")
    ret = GRBfreemodel(model)
    @assert ret == 0
    GRBfreeenv(env)
    return
end

function solve_lp_jump(A::Matrix{Tf}, b::Vector{Tf}, c::Vector{Tf}) where Tf<:Cdouble
    start_time = time()
    m, n = size(A)
    M = Model(Gurobi.Optimizer)
    set_silent(M)
    @variable(M, x[1:n] .>= 0)
    @objective(M, Max, c'*x)
    @constraint(M, A*x .<= b)
    println("model setup = $(time() - start_time)")
    solve_time = @elapsed optimize!(M)
    println("model solve = $(solve_time)")
    return
end

function solve_lp_jump2(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    start_time = time()
    m, n = size(A)
    model = direct_model(Gurobi.Optimizer())
    set_silent(model)
    set_string_names_on_creation(model, false)
    @variable(model, x[1:n] >= 0)
    @objective(model, Max, c' * x)
    @constraint(model, A* x .<= b)
    println("model setup = $(time() - start_time)")
    solve_time = @elapsed optimize!(model)
    println("model solve = $solve_time")
    return
end

m = 500
n = 10_000
A = randn(Cdouble, m, n);
b = randn(Cdouble, m);
c = randn(Cdouble, n);
@time sol1 = solve_lp_direct(A, b, c);
@time sol2 = solve_lp_jump(A, b, c);
@time sol3 = solve_lp_jump2(A, b, c);

and here are the results:

julia> include("/Users/oscar/Desktop/grb.jl")
model setup = 0.19083690643310547
model solve = 4.981149123
  5.296278 seconds (489.36 k allocations: 86.229 MiB, 1.23% gc time, 2.05% compilation time)
model setup = 5.112233877182007
model solve = 10.832208031
 16.768803 seconds (69.07 M allocations: 3.192 GiB, 6.36% gc time)
model setup = 2.414860963821411
model solve = 5.093770094
  7.588198 seconds (43.45 M allocations: 1.522 GiB, 4.90% gc time, 1.67% compilation time)

julia> include("/Users/oscar/Desktop/grb.jl")
model setup = 0.1080930233001709
model solve = 4.50650111
  4.683605 seconds (241.50 k allocations: 71.713 MiB, 1.21% compilation time)
model setup = 2.224179983139038
model solve = 4.767565328
  7.054845 seconds (43.53 M allocations: 1.673 GiB, 4.95% gc time, 0.89% compilation time)
model setup = 2.6410491466522217
model solve = 4.550919043
  7.258149 seconds (43.37 M allocations: 1.518 GiB, 9.25% gc time, 0.91% compilation time)

julia> include("/Users/oscar/Desktop/grb.jl")
model setup = 0.11534619331359863
model solve = 4.978835491
  5.165873 seconds (241.50 k allocations: 71.713 MiB, 1.13% compilation time)
model setup = 2.328436851501465
model solve = 5.220872358
  7.617613 seconds (43.51 M allocations: 1.673 GiB, 5.08% gc time, 0.89% compilation time)
model setup = 2.448287010192871
model solve = 5.000077219
  7.517537 seconds (43.36 M allocations: 1.518 GiB, 5.13% gc time, 0.92% compilation time)

julia> include("/Users/oscar/Desktop/grb.jl")
model setup = 0.11165499687194824
model solve = 5.316515968
  5.502262 seconds (241.50 k allocations: 71.713 MiB, 1.10% compilation time)
model setup = 2.561717987060547
model solve = 5.696499657
  8.333657 seconds (43.54 M allocations: 1.673 GiB, 5.18% gc time, 0.90% compilation time)
model setup = 2.7369630336761475
model solve = 5.279443178
  8.087859 seconds (43.38 M allocations: 1.518 GiB, 6.56% gc time, 0.88% compilation time)

If we ignore the first run because of compilation latency, then yeah, the model setup is much slower via JuMP. But that isn’t surprising because the Gurobi.Optimizer object initializes a whole lot more than our C code. It prepares data structures at the JuMP level, creates a variety of dictionaries to map Cint indices to JuMP variables, normalizes the constraints, etc.

My bigger takeaway is that JuMP took 8 seconds. Gurobi C API took 5. So for the cost of 3 seconds we avoided all of the issues that you’re encountering trying to get the C API to work.

p.s. I posted a link on your Gurobi portal post to this discussion. Note that Gurobi doesn’t provide official support for Julia, so they might not be much help.

The C examples are here: https://www.gurobi.com/documentation/10.0/examples/c_examples.html

jacob-roth · March 3, 2023, 5:39am

Totally agree, I love using JuMP. At the moment, I’m just curious about model construction time, and there the Gurobi API seems quite fast. However, I noticed something odd about the problem constructions: they are giving/solving different problems?

For example, I modified your improved version of my code to return the solution (if one exists) and came across a case where the C API finds an optimal solution to (some) problem but JuMP says the problem is infeasible, so I don’t know what the C API is solving… Which makes me think there’s an issue in how I’m trying to use the C API?

julia> versioninfo()
Julia Version 1.8.3
Commit 0434deb161e (2022-11-14 20:14 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × Intel(R) Xeon(R) W-2145 CPU @ 3.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, skylake-avx512)
  Threads: 1 on 16 virtual cores
Environment:
  LD_LIBRARY_PATH = :/home/myuser/software/gurobi1000/linux64/lib:/home/myuser/software/gurobi1000/linux64/lib

using Gurobi
using JuMP
using LinearAlgebra
using Random

function solve_lp_direct(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    start_time = time()
    env_p = Ref{Ptr{Cvoid}}()
    error = GRBloadenv(env_p, "lp.log")
    env = env_p[]
    error = GRBsetintparam(env, GRB_INT_PAR_OUTPUTFLAG, 1)
    @assert error == 0
    model_p = Ref{Ptr{Cvoid}}()
    error = GRBnewmodel(env, model_p, "lp", 0, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL)
    @assert error == 0
    model = model_p[]
    m, n = size(A)
    error = GRBsetintattr(model, GRB_INT_ATTR_MODELSENSE, GRB_MAXIMIZE)
    @assert error == 0
    error = GRBaddvars(
        model,         # model
        n,             # : numvars
        0,             # : numnz
        C_NULL,        # : *vbeg
        C_NULL,        # : *vind
        C_NULL,        # : *vval
        c,             # : *obj
        zeros(Cdouble, n), # : *lb
        fill(Inf, n),  # : *ub
        C_NULL,        # : *vtype
        C_NULL         # : **varnames
    )
    @assert error == 0
    for i in 1:m
        error = GRBaddconstr(
            model,              # : *model
            n,                  # : numnz
            Cint.(0:(n-1)),     # : *cind
            A[i, :],            # : *cval
            # GRB_LESS_EQUAL,     # : sense
            GRB_EQUAL,          # : sense
            b[i],               # : rhs
            C_NULL,             # : *constrname
        )
        @assert error == 0
    end
    println("model setup = $(time() - start_time)")
    solve_time = @elapsed error = GRBoptimize(model)
    println("model solve = $solve_time")
    sol = ones(T, n)
    ret = GRBgetdblattrarray(model, GRB_DBL_ATTR_X, 0, n, sol);
    @assert ret == 0
    objval = Ref{Cdouble}()
    ret = GRBgetdblattr(model, GRB_DBL_ATTR_OBJVAL, objval);
    @assert ret == 0
    ret = GRBfreemodel(model)
    @assert ret == 0
    GRBfreeenv(env)
    return sol, objval
end

function solve_lp_jump(A::Matrix{Tf}, b::Vector{Tf}, c::Vector{Tf}) where Tf<:Cdouble
    start_time = time()
    m, n = size(A)
    M = Model(Gurobi.Optimizer)
    set_silent(M)
    @variable(M, x[1:n] .>= 0)
    @objective(M, Max, c'*x)
    @constraint(M, A*x .<= b)
    println("model setup = $(time() - start_time)")
    solve_time = @elapsed optimize!(M)
    println("model solve = $(solve_time)")
    return value.(x), JuMP.objective_value(M)
end

function solve_lp_jump2(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    start_time = time()
    m, n = size(A)
    model = direct_model(Gurobi.Optimizer())
    set_silent(model)
    set_string_names_on_creation(model, false)
    @variable(model, x[1:n] >= 0)
    @objective(model, Max, c' * x)
    @constraint(model, A* x .<= b)
    println("model setup = $(time() - start_time)")
    solve_time = @elapsed optimize!(model)
    println("model solve = $solve_time")
    return value.(x), objective_value(model)
end

The problem:

julia> Random.seed!(1234)
TaskLocalRNG()

julia> m = 50
50

julia> n = 100
100

julia> A = randn(Cdouble, m, n);

julia> b = randn(Cdouble, m);

julia> c = randn(Cdouble, n);

julia> @time sol1, objval1 = solve_lp_direct(A, b, c);
Set parameter Username
Academic license - for non-commercial use only - expires 2023-12-17
model setup = 0.02570796012878418
Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (linux64)

CPU model: Intel(R) Xeon(R) W-2145 CPU @ 3.70GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 50 rows, 100 columns and 5000 nonzeros
Model fingerprint: 0xd65f6397
Coefficient statistics:
  Matrix range     [1e-04, 5e+00]
  Objective range  [2e-02, 3e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [5e-02, 2e+00]
Presolve time: 0.01s
Presolved: 50 rows, 100 columns, 5000 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    4.6543195e+31   8.349456e+31   4.654320e+01      0s
      92    5.2160968e+00   0.000000e+00   0.000000e+00      0s

Solved in 92 iterations and 0.01 seconds (0.01 work units)
Optimal objective  5.216096848e+00
model solve = 0.030953253
  0.069868 seconds (7.81 k allocations: 477.694 KiB, 18.47% compilation time)

julia> sum(sol1)
29.62203562836316

julia> @time sol2, objval2 = solve_lp_jump(A, b, c);
Set parameter Username
Academic license - for non-commercial use only - expires 2023-12-17
model setup = 0.0030510425567626953
model solve = 0.003450798
ERROR: Gurobi Error 10005: Unable to retrieve attribute 'ObjVal'
Stacktrace:
  [1] _check_ret
    @ ~/.julia/packages/Gurobi/FliRK/src/MOI_wrapper/MOI_wrapper.jl:368 [inlined]                                                                             [0/1940]
  [2] get(model::Gurobi.Optimizer, attr::MathOptInterface.ObjectiveValue)
    @ Gurobi ~/.julia/packages/Gurobi/FliRK/src/MOI_wrapper/MOI_wrapper.jl:3097
  [3] get(b::MathOptInterface.Bridges.LazyBridgeOptimizer{Gurobi.Optimizer}, attr::MathOptInterface.ObjectiveValue)
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/GIbN0/src/Bridges/bridge_optimizer.jl:1002
  [4] _get_model_attribute(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{Gurobi.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, attr::MathOptInterface.ObjectiveValue)
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/GIbN0/src/Utilities/cachingoptimizer.jl:828
  [5] get
    @ ~/.julia/packages/MathOptInterface/GIbN0/src/Utilities/cachingoptimizer.jl:876 [inlined]
  [6] _moi_get_result(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{Gurobi.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, args::MathOptInterface.ObjectiveValue)
    @ JuMP ~/.julia/packages/JuMP/Xw2ck/src/optimizer_interface.jl:680
  [7] get(model::Model, attr::MathOptInterface.ObjectiveValue)
    @ JuMP ~/.julia/packages/JuMP/Xw2ck/src/optimizer_interface.jl:700
  [8] objective_value(model::Model; result::Int64)
    @ JuMP ~/.julia/packages/JuMP/Xw2ck/src/objective.jl:54
  [9] objective_value
    @ ~/.julia/packages/JuMP/Xw2ck/src/objective.jl:50 [inlined]
 [10] solve_lp_jump(A::Matrix{Float64}, b::Vector{Float64}, c::Vector{Float64})
    @ Main ./REPL[30]:12
 [11] top-level scope
    @ ./timing.jl:262 [inlined]
 [12] top-level scope
    @ ./REPL[55]:0

julia> @time sol3, objval3 = solve_lp_jump2(A, b, c);
Set parameter Username
Academic license - for non-commercial use only - expires 2023-12-17
model setup = 0.023720979690551758
model solve = 0.026847763
ERROR: Gurobi Error 10005: Unable to retrieve attribute 'ObjVal'
Stacktrace:
 [1] _check_ret
   @ ~/.julia/packages/Gurobi/FliRK/src/MOI_wrapper/MOI_wrapper.jl:368 [inlined]
 [2] get(model::Gurobi.Optimizer, attr::MathOptInterface.ObjectiveValue)
   @ Gurobi ~/.julia/packages/Gurobi/FliRK/src/MOI_wrapper/MOI_wrapper.jl:3097
 [3] _moi_get_result(model::Gurobi.Optimizer, args::MathOptInterface.ObjectiveValue)
   @ JuMP ~/.julia/packages/JuMP/Xw2ck/src/optimizer_interface.jl:671
 [4] get(model::Model, attr::MathOptInterface.ObjectiveValue)
   @ JuMP ~/.julia/packages/JuMP/Xw2ck/src/optimizer_interface.jl:700
 [5] objective_value(model::Model; result::Int64)
   @ JuMP ~/.julia/packages/JuMP/Xw2ck/src/objective.jl:54
 [6] objective_value
   @ ~/.julia/packages/JuMP/Xw2ck/src/objective.jl:50 [inlined]
 [7] solve_lp_jump2(A::Matrix{Float64}, b::Vector{Float64}, c::Vector{Float64})
   @ Main ./REPL[31]:13
 [8] top-level scope
   @ ./timing.jl:262 [inlined]
 [9] top-level scope
   @ ./REPL[56]:0

julia>

Thanks! I will try comparing the above with the dense implementation from this example.

odow · March 3, 2023, 5:52am

            # GRB_LESS_EQUAL,     # : sense
            GRB_EQUAL,          # : sense

They aren’t the same problem. JuMP is <= and Gurobi is =.

I’m just curious about model construction time, and there the Gurobi API seems quite fast.

If you just want to measure construction time, set a small time limit of 0.01 seconds and call optimize!. It often isn’t correct to measure the time building with actually calling optimize because solvers may lazily defer construction until the actual solve. (See [2206.03866] JuMP 1.0: Recent improvements to a modeling language for mathematical optimization for a discussion on how to benchmark.)

jacob-roth · March 3, 2023, 5:57am

Oh gosh I missed that as I was playing around with the code, wow just saw, sorry

jacob-roth · March 14, 2023, 8:25pm

Hopefully the last thing. I would like to pass a sparse matrix to Gurobi. Gurobi uses CSR but SparseArrays is CSC. I’m sure that there’s a better way, but I’ve converted my sparse constraint matrix A to CSR with the SparseMatricesCSR package, but GRBXaddconstrs is complaining (gurobi error 1018 = “duplicated indices” error). Is there anything obvious that I’m doing wrong?

using SparseMatricesCSR, Random, Gurobi
Random.seed!(1010101)
m = 1_000
n = 10_000
A = rand(Cdouble, m, n);
b = rand(Cdouble, m);
c = randn(Cdouble, n);
function solve_biglp_direct(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    Acsr = SparseMatricesCSR.SparseMatrixCSR(transpose(sparse(transpose(A))));
    start_time = time()
    env_p = Ref{Ptr{Cvoid}}()
    error = GRBloadenv(env_p, "lp.log")
    env = env_p[]
    error = GRBsetintparam(env, GRB_INT_PAR_OUTPUTFLAG, 1)
    @assert error == 0
    model_p = Ref{Ptr{Cvoid}}()
    error = GRBnewmodel(env, model_p, "lp", 0, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL)
    @assert error == 0
    model = model_p[]
    m, n = size(A)
    error = GRBsetintattr(model, GRB_INT_ATTR_MODELSENSE, GRB_MAXIMIZE)
    @assert error == 0
    error = GRBXaddvars(
        model,             # model
        Base.Csize_t(n),   # : numvars
        0,                 # : numnz
        C_NULL,            # : *vbeg
        C_NULL,            # : *vind
        C_NULL,            # : *vval
        c,                 # : *obj
        zeros(Cdouble, n), # : *lb
        fill(Inf, n),      # : *ub
        C_NULL,            # : *vtype
        C_NULL             # : **varnames
    )
    @assert error == 0
    # for i in 1:m
    #     error = GRBaddconstr(
    #         model,              # : *model
    #         n,                  # : numnz
    #         Cint.(0:(n-1)),     # : *cind
    #         A[i, :],            # : *cval
    #         GRB_LESS_EQUAL,     # : sense
    #         b[i],               # : rhs
    #         C_NULL,             # : *constrname
    #     )
    #     @assert error == 0
    # end
    consense = fill(GRB_LESS_EQUAL, Acsr.m)
    conname = fill(C_NULL, Acsr.m)
    conname = ["$i" for i in eachindex(b)]
    error = GRBXaddconstrs(
        model,                      # : *model
        Acsr.m,                     # : numconstrs
        Base.Csize_t(nnz(Acsr)),    # : numnz
        Base.Csize_t.(Acsr.rowptr.-1), # : *cbeg
        Acsr.colval.-1,                # : *cind
        Acsr.nzval,                 # : *cval
        consense,                   # : sense
        b,                          # : *rhs
        conname,                    # : **constrname
    )
    @assert error == 0
    GRBupdatemodel(model)
    GRBsetdblparam(env, "TimeLimit", 0.0)
    println("model setup = $(time() - start_time)")
    solve_time = @elapsed error = GRBoptimize(model)
    println("model setup = $(time() - start_time)")
    println("model solve = $solve_time")
    sol = ones(T, n)
    ret = GRBgetdblattrarray(model, GRB_DBL_ATTR_X, 0, n, sol);
    @assert ret == 0
    objval = Ref{Cdouble}()
    ret = GRBgetdblattr(model, GRB_DBL_ATTR_OBJVAL, objval);
    @assert ret == 0
    ret = GRBfreemodel(model)
    @assert ret == 0
    GRBfreeenv(env)
    return sol, objval
end

odow · March 14, 2023, 9:21pm

When using the C API, you need to be very, very careful about the types. There are no error checks or warnings. Doing the wrong thing will just give an error, the wrong answer, or a segfault.

In this case, read the documentation:

https://www.gurobi.com/documentation/9.5/refman/c_xaddconstrs.html

The signature is

But your Acsr.colval.-1 is a Vector{Int}, when it needs to be Vector{Cint}. Gurobi will happily accept your argument, but it will reinterpret your Int as two packed Cints.

The fix is probably just Cint.(Acsr.colval.-1), although I haven’t tested.

You should also follow this advice:

Your senses are probably wrong as well. I think it needs to be consense = fill(Cchar(GRB_LESS_EQUAL), Acsr.m).

jacob-roth · March 14, 2023, 9:59pm

Thanks!! It worked with Cint and Cchar Just to document, the working code is below:

using Gurobi
using JuMP
using LinearAlgebra
using Random
using SparseMatricesCSR

function solve_sparselp_direct(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    Acsr = SparseMatricesCSR.SparseMatrixCSR(transpose(sparse(transpose(A))));
    start_time = time()
    env_p = Ref{Ptr{Cvoid}}()
    error = GRBloadenv(env_p, "lp.log")
    env = env_p[]
    error = GRBsetintparam(env, GRB_INT_PAR_OUTPUTFLAG, 1)
    @assert error == 0
    model_p = Ref{Ptr{Cvoid}}()
    error = GRBnewmodel(env, model_p, "lp", 0, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL)
    @assert error == 0
    model = model_p[]
    m, n = size(A)
    error = GRBsetintattr(model, GRB_INT_ATTR_MODELSENSE, GRB_MAXIMIZE)
    @assert error == 0
    error = GRBXaddvars(
        model,             # model
        Base.Csize_t(n),   # : numvars
        0,                 # : numnz
        C_NULL,            # : *vbeg
        C_NULL,            # : *vind
        C_NULL,            # : *vval
        c,                 # : *obj
        zeros(Cdouble, n), # : *lb
        fill(Inf, n),      # : *ub
        C_NULL,            # : *vtype
        C_NULL             # : **varnames
    )
    @assert error == 0
    consense = fill(Cchar(GRB_LESS_EQUAL), Acsr.m)
    conname = fill(C_NULL, Acsr.m)
    error = GRBXaddconstrs(
        model,                         # : *model
        Acsr.m,                        # : numconstrs
        Base.Csize_t(nnz(Acsr)),       # : numnz
        Base.Csize_t.(Acsr.rowptr.-1), # : *cbeg
        Cint.(Acsr.colval.-1),         # : *cind
        Acsr.nzval,                    # : *cval
        consense,                      # : sense
        b,                             # : *rhs
        conname,                       # : **constrname
    )
    @assert error == 0
    GRBupdatemodel(model)
    GRBsetdblparam(env, "TimeLimit", 0.0)
    println("model setup = $(time() - start_time)")
    solve_time = @elapsed error = GRBoptimize(model)
    println("model setup = $(time() - start_time)")
    println("model solve = $solve_time")
    sol = ones(T, n)
    ret = GRBgetdblattrarray(model, GRB_DBL_ATTR_X, 0, n, sol);
    @assert ret == 0
    objval = Ref{Cdouble}()
    ret = GRBgetdblattr(model, GRB_DBL_ATTR_OBJVAL, objval);
    @assert ret == 0
    ret = GRBfreemodel(model)
    @assert ret == 0
    GRBfreeenv(env)
    return sol, objval
end

odow · March 14, 2023, 10:13pm

Perhaps a clarification of what is actually going on will be helpful.

The definition of GRBXaddconstrs is

function GRBXaddconstrs(model, numconstrs, numnz, cbeg, cind, cval, sense, rhs, constrnames)
    ccall((:GRBXaddconstrs, libgurobi), Cint, (Ptr{GRBmodel}, Cint, Csize_t, Ptr{Csize_t}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cchar}, Ptr{Cdouble}, Ptr{Ptr{Cchar}}), model, numconstrs, numnz, cbeg, cind, cval, sense, rhs, constrnames)
end

Julia’s ccall converts arguments that are passed to it to the type given in the corresponding position in the tuple that is the third argument.

It converts arguments using Base.unsafe_convert(T, Base.cconvert(T, arg)) where T is the desired type and arg is the argument.

So take numconstrs. ccall expects a Cint, but you passed it Acsr.m which is a Int64.

julia> y = Int64(1)
1

julia> z = Base.unsafe_convert(Cint, Base.cconvert(Cint, y))
1

julia> typeof(z)
Int32

Bingo. That worked.

But *cind wanted Ptr{Cint}, and you passed it a Vector{Int64}.

julia> x = [1, 2, 3]
3-element Vector{Int64}:
 1
 2
 3

julia> typeof(x)
Vector{Int64} (alias for Array{Int64, 1})

julia> Base.unsafe_convert(Ptr{Cint}, Base.cconvert(Ptr{Cint}, x))
Ptr{Int32} @0x0000000149a046e0

julia> pointer(x)
Ptr{Int64} @0x0000000149a046e0

In this case, Julia converted the pointer of x to the right type, but it didn’t convert the underlying data to the element type Cint.

So as a general rule, if you’re passing scalar, you can be fast and loose with the type and it should convert okay. But if you’re passing a vector, you need to make sure the element type is the bit-exact type expected by the C library.

jacob-roth · March 14, 2023, 10:28pm

Ahh, thanks! Really helpful to see the scalar / vector (pointer) distinction in action. I now see how this explains why it worked with “numconstrs” set to Ascr.m::Int64 even though the C API wants Int32.