Memory consumption and time of @constraint

[Thuener]

I'm having some issues with memory consumption on JuMP. I have a problem that has too many constraints and the JuMP structures for macro @constraint is consuming too much memory.

Example:

using JuMP,CPLEX

function memuse()
 pid = getpid()
 return round(Int,parse(Int,readstring(`ps -p $pid -o rss=`))/1024)
end

memuse()
C = 300000
m = Model(solver=CplexSolver())
N = 100
@variable(m, x[1:N] >= 0)
@variable(m, Θ >= 0)

@objective(m, Max, Θ )
solve(m)

coef = rand(C,N)
m1= memuse()
tic()
@constraint(m, Θ .<= coef*x  )
time = toq()
gc()
m2 = memuse()
print("Memory consumption $(m2-m1) Mb and Time $(time)")

Memory consumption 1022 Mb and Time 11.076951813

Changing the @constraint to CPLEX.add_constrs.

m1= memuse()
tic()
rhs = zeros(C)
coef = hcat(-coef,ones(C));
CPLEX.add_constrs!(m.internalModel.inner, coef, '<', rhs)
time = toq()
gc()
m2 = memuse()

Memory consumption 447 Mb and Time 1.975631208

With JuMP this constraint consumes 2.28x more memory and become 5.6x more slower. Maybe JuMP should have some macro to add constraints direct to the model without storing any information about it.

[joehuchette]

In the first example, you're really seeing the cost of the vectorized syntax. Writing out the loop explicitly gives me about an order of magnitude speedup over the vectorized version:

for i in 1:size(coef,1)
    @constraint(m, Θ <= sum(coef[i,j]*x[j] for j in size(coef,2)))
end

I think the syntax for adding constraints directly to the model, bypassing JuMP, would look a lot like your second code example there :)

[mlubin]

Using vectorized syntax in JuMP is known to consume more memory (because of temporaries) and to be slower than scalar syntax.

[Thuener]

With:

@constraint(m,[i = 1:size(coef,1)], Θ <= sum(coef[i,j]*x[j] for j = 1:size(coef,2)))

Memory consumption 957 Mb and Time 18.123854547

I'm on a different computer, the results for CPLEX.add_constrs on this one are:
Memory consumption 377 Mb and Time 3.32205742

So the scalar syntax consumes 2.53x more memory and is 5.45x more slower then CPLEX.add_constrs. I don't see much difference between this values and the vectorized constraint. So the memory consumption of the vectorized constraint is almost the same as scalar syntax. I think the memory consumption is because of the auxiliary objects created by JuMP.

[mlubin]

Make sure you're timing this inside of a function. Performance at the global scope will be much worse.

[Thuener]

Ok. Using:

function addConstraints(m,x,coef)
 @constraint(m,[i = 1:size(coef,1)], Θ <= sum(coef[i,j]*x[j] for j = 1:size(coef,2)))
end

...

coef = rand(C,N)
m1= memuse()
tic()
addConstraints(m,x,coef)
time = toq()
gc()
m2 = memuse()
print("Memory consumption $(m2-m1) Mb and Time $(time)")

I got: Memory consumption 1180 Mb and Time 11.69368708

To be fair, I also tested the result of the vectorized constraint inside a function.

function addConstraints(m,x,coef)
 @constraint(m, Θ .<= coef*x  )
end

Memory consumption 1035 Mb and Time 14.692361745

[mlubin]

using JuMP,CPLEX

function memuse()
 pid = getpid()
 return round(Int,parse(Int,readstring(`ps -p $pid -o rss=`))/1024)
end

const C = 300000
const N = 100

function t1()
    m = Model(solver=CplexSolver())
    @variable(m, x[1:N] >= 0)
    @variable(m, Θ >= 0)

    @objective(m, Max, Θ )

    coef = rand(C,N)
    gc()
    m1= memuse()
    tic()
    @constraint(m, Θ .<= coef*x  )
    time = toq()
    gc()
    m2 = memuse()
    println("Memory consumption $(m2-m1) Mb and Time $(time)")
end

function t2()
    m = Model(solver=CplexSolver())
    @variable(m, x[1:N] >= 0)
    @variable(m, Θ >= 0)

    @objective(m, Max, Θ )

    coef = rand(C,N)
    gc()
    m1= memuse()
    tic()
    for i in 1:size(coef,1)
        @constraint(m, Θ <= sum(coef[i,j]*x[j] for j = 1:size(coef,2)))
    end
    time = toq()
    gc()
    m2 = memuse()
    println("Memory consumption $(m2-m1) Mb and Time $(time)")
end

println("t1")
t1()
t1()
println("t2")
t2()
t2()

gives

t1
Memory consumption 1359 Mb and Time 6.816507673
Memory consumption 1300 Mb and Time 6.256449671
t2
Memory consumption 649 Mb and Time 1.88936574
Memory consumption 693 Mb and Time 1.878137287

[Thuener]

The code above has some weird behavior. I prefer to create a function to add the constraints, because I'm worried about the overall memory use not the temporary created objects.

using JuMP,CPLEX

function memuse()
 pid = getpid()
 return round(Int,parse(Int,readstring(`ps -p $pid -o rss=`))/1024)
end

const C = 300000
const N = 100

function addconstraints(m,x,Θ,coef,p)
    if p == 1
      @constraint(m, Θ .<= coef*x  )
    elseif p == 2
      for i in 1:size(coef,1)
        @constraint(m, Θ <= sum(coef[i,j]*x[j] for j = 1:size(coef,2)))
      end
    elseif p ==3
      @constraint(m,[i = 1:size(coef,1)], Θ <= sum(coef[i,j]*x[j] for j = 1:size(coef,2)))
    else
      rhs = zeros(C)
      coef = hcat(-coef,ones(C));
      CPLEX.add_constrs!(m.internalModel.inner, coef, '<', rhs)
    
    end

end
function t(p)
    m = Model(solver=CplexSolver(CPX_PARAM_SCRIND=0))
    @variable(m, 0 <= x[1:N] <= 1)
    @variable(m, 0 <= Θ <= 1000)

    @objective(m, Max, Θ )
    solve(m)

    coef = rand(C,N)
    gc()
    m1= memuse()
    tic()
    addconstraints(m,x,Θ,coef,p)
    time = toq()
    gc()
    m2 = memuse()
    println("Memory consumption $(m2-m1) Mb and Time $(time)")
end

println("Vectorized")
t(1)
t(1)

println("Scalar 1")
t(2)
t(2)

println("Scalar 2")
t(3)
t(3)

println("Low-level API")
t(4)
t(4)

Vectorized
Memory consumption 1079 Mb and Time 14.493644756
Memory consumption 1023 Mb and Time 13.870267812
Scalar 1
Memory consumption 1107 Mb and Time 7.674387291
Memory consumption 1119 Mb and Time 8.528182488
Scalar 2
Memory consumption 1130 Mb and Time 8.32419865
Memory consumption 1132 Mb and Time 8.495367176
Low-level API
Memory consumption 327 Mb and Time 3.042267338
Memory consumption 355 Mb and Time 2.844282922

[mlubin]

The inconsistent performance of the scalar case suggests issues in Julia with type inference. Try using @code_warntype.

[Thuener]

The most important problem here is the big difference between memory consumption, the difference in performance should be an consequence of that.

[mlubin]

We store linear constraints as vectors of sparse affine expressions. For each coefficient we additionally store the corresponding Variable (an integer column index and a pointer to the model), so 2.5x more memory than storing the dense matrix on its own is expected. JuMP is designed for sparse problems, and I don't see our data structures changing anytime soon.

Additionally, there will be two copies of the constraint matrix kept in memory: JuMP's internal copy and the solver's internal copy. If this causes a memory bottleneck for you, you should consider not using JuMP.

[Thuener]

Ok, I see that JuMP is not created for that but why not have a special function to add constraints without storing any information about them. I'm not the only one with this problem, many of my collegues have similar problems with memory consumption or perfomance, most of them use with benders decompositions.

[mlubin]

If JuMP has an internal model loaded, then you're free to call MathProgBase.addconstr! to add constraints that JuMP won't keep track of. You should be careful though to not add more constraints through JuMP once you start doing this, because the indices of the constraints will be mismatched and dual values on the constraints will not be correct.

[Thuener]

I added a flag stroreconstr in the JuMP model, that enables to add constraints without storing them. I just set this flag to false after creating the model on this code:

using JuMP,CPLEX

function memuse()
 pid = getpid()
 return round(Int,parse(Int,readstring(`ps -p $pid -o rss=`))/1024)
end

const C = 300000
const N = 100

function addconstraints(m,x,Θ,coef,p)
    m.storeconstr = false
    if p == 1
      @constraint(m, Θ .<= coef*x  )
    elseif p == 2
      for i in 1:size(coef,1)
        @constraint(m, Θ <= sum(coef[i,j]*x[j] for j = 1:size(coef,2)))
      end
    elseif p ==3
      @constraint(m,[i = 1:size(coef,1)], Θ <= sum(coef[i,j]*x[j] for j = 1:size(coef,2)))
    else
      rhs = zeros(C)
      coef = hcat(-coef,ones(C));
      CPLEX.add_constrs!(m.internalModel.inner, coef, '<', rhs)
    
    end

end
function t(p)
    m = Model(solver=CplexSolver(CPX_PARAM_SCRIND=0))
    @variable(m, 0 <= x[1:N] <= 1)
    @variable(m, 0 <= Θ <= 1000)

    @objective(m, Max, Θ )
    solve(m)

    coef = rand(C,N)
    gc()
    m1= memuse()
    tic()
    addconstraints(m,x,Θ,coef,p)
    time = toq()
    gc()
    m2 = memuse()
    println("Memory consumption $(m2-m1) Mb and Time $(time)")
end

println("Scalar 1")
t(2)
t(2)
println("Low-level API")
t(4)
t(4)

has the flowing result:
Scalar 1
Memory consumption 537 Mb and Time 4.433234855
Memory consumption 476 Mb and Time 4.396045515
Low-level API
Memory consumption 445 Mb and Time 2.819263039
Memory consumption 419 Mb and Time 2.84959166

Memory consumption change from 2.65x to 1.17x and performance from 3.49x to 1.55x.

[mlubin]

We already have support in principle for keyword arguments within @constraint. I would consider a well-formed PR that adds a keyword argument to addconstraint(m::Model, c::LinearConstraint) for only passing the constraint to the solver and not storing it in JuMP.

[Thuener]

Thanks, that is all I'm asking for.

[odow]

I don't think this should be added to JuMP. It makes modifications to the JuMP problem too brittle. Unless we implement a JuMP index -> MPB index. Ref JuliaOpt/MathProgBase.jl#125 JuliaOpt/MathProgBase.jl#139

In particular, there are problems with things like updating RHS vectors. Should MathProgBase.setconstrUB!(m, ub) set the first length(ub) constraints? Or should JuMP query the number of constraints in the internalmodel, and be responsible for padding out ub to the correct length?

Users who want this functionality should know the consequences and be forced to do something like

varidx = [v.col for v in [ ... variables ... ]]
coef   = [ ... coefficients ... ]
lb     = -Inf
ub     = 1
push!(m.linconstr, JuMP.LinearConstraint(0, lb, ub)) # dummy constraint for JuMP facing side
MathProgBase.addconstr!(internalmodel(m), varidx, coef, lb, ub)

[odow]

Another stopgap could be

function zeroconstraintmatrix!(m)
    warn("You are about to zero the constraint matrix of the JuMP model. Hopefully you know what this means!")
    for con in m.linconstr
        con.terms = 0
    end
end

function t(p, dozero)
    m = Model(solver=CplexSolver(CPX_PARAM_SCRIND=0))
    @variable(m, 0 <= x[1:N] <= 1)
    @variable(m, 0 <= Θ <= 1000)

    @objective(m, Max, Θ )
    solve(m)

    coef = rand(C,N)
    gc()
    m1= memuse()
    tic()
    addconstraints(m,x,Θ,coef,p)
    dozero && zeroconstraintmatrix!(m)
    time = toq()
    gc()
    m2 = memuse()
    println("Memory consumption $(m2-m1) Mb and Time $(time)")
end

t(2, true)
# Memory consumption 532.87109375 Mb and Time 6.88804248

t(2, false)
# Memory consumption 1048.49609375 Mb and Time 6.702602572

t(4, true)
# Memory consumption 311.9140625 Mb and Time 1.514756458

t(4, false)
# Memory consumption 355.05078125 Mb and Time 1.491463572

[Thuener]

I think memory consumption and performance should be a key factor for JuMP. Many of my colleagues, some that I persuaded to use Julia and JuMP, are coming for a C++/C background and are using Julia because of the simplicity and performance. However, many of them are having problem with memory consumption and performance. JuMP is a wonderful package, but to establish as a package to solve all king of mathematical optimization problems it should be easy for beginners(as it is) and possible to be tweaked by experts.

[mlubin]

@odow has a point. If you're trying to do something that JuMP wasn't designed to do (i.e., not store a copy of the constraints), then your code should look ugly on principle, unless you create your own pretty layer on top. You can also do pop!(m.linconstr) after adding an individual constraint and face the consequences of the JuMP model not matching the solver's model.

We've never claimed that JuMP has similar memory characteristics to hand-tuned solver-specific matrix generators. I think having to keep around two copies of the constraint matrix in memory is a fair price to pay for the added generality. I'm also confused about how push!(m.linconstr, c) is causing such a performance issue.

[joaquimg]

I think JuMP created an awesome environment for mathematical programming, pretty, easy to use, efficient for many situations, we all know that the list is huge!
I never had problems with JuMP for solving MIPs, in my opinion that's because most of the time is spent on the solver anyway, so JuMP bottlenecks disappear...

Writing problem directly in matrix form should faster and lighter. My point is that we could have something in between, not as fast as writing matrices (which are terrible to write LP's and maintain code) nor as slow as some few current JuMP bottlenecks for some algorithms that solve hundreds of thousands of LP´s. This is at least the third time some bottleneck appears while someone is implementing SDDP, myself and @blegat had some other problems, with deleting constraints I think...

I am aware that the current focus is to finish JuMP 1.0 and modifying the design is not an option, I agree with that. However, adding some keyword arguments that could improve algorithm performance at cost of breaking some JuMP functionality would be just ok, in my opinion, as long as the user is aware of that. The pros could be way greater than the cons...

Maybe we could think of some way to make JuMP as awesome for solving millions of similar LP´s as it already is for MIP´s!
Sorry if that is the wrong place for this, but I think we could keep this in mind. Mainly if people are willing to prepare PR´s with improvements.

[joaquimg]

By the way, I liked @odow idea of having something like JuMP index -> MPB index. this could be a way to go...

[mlubin]

I appreciate the enthusiasm for wanting JuMP to do everything that you want it to do and quickly, but I think the discussion will be more productive if we focus on the technical issues here. There are two separate things that are being confused at the moment.

• Memory: There will be two copies of the constraint matrix in memory. I don't think we need to provide a pretty way to avoid this. You're free to go and modify m.linconstr as @odow and I later suggested.
• Performance: Nobody has pointed to a specific thing that's causing seemingly a 2+ times performance hit for no good reason. This merits further investigation and might have a simple fix. I won't be leading this investigation but I'm happy to continue discussing.

[Azizimj]

Hello every body,
I have a problem similar to this one which is with storing variables even before solving the problem. I am using JuMP and Gurobi but the number of variables is large and gives ERROR: OutOfMemoryError().
My code is like this:

@variable(model, w[1:1000,1:1000,1:620], Bin)

I am not that professional in JuMP and Modeling. Can you please help?

[blegat]

A Variable use 16 bytes (8 bytes for a reference to the Model and 8 bytes for a variable id). If you create 620M variables, this makes 9.92 GB of memory. If you have 8GB of RAM, it is not surprising that you get an OutOfMemoryError(). Even if you had enough memory to store so many variable, I doubt that Gurobi could solve a problem of that size.

[Thuener]

The new code presented in JuMP-dev 2019(https://www.youtube.com/watch?v=MLunP5cdRBI):
https://gist.github.com/Thuener/5fd30bda29a84afb126cb5b723574eba

########### Vectorized - WithoutDirect ###########
Memory consumption 1111 Mb and Time 20.904308405
Time to solve model 5.079105861
Memory consumption 1118 Mb and Time 20.27646884
Time to solve model 5.080772802
########### Scalar 1 - WithoutDirect ###########
Memory consumption 997 Mb and Time 12.158498048
Time to solve model 5.005659555
Memory consumption 992 Mb and Time 12.178707011
Time to solve model 5.056525515
########### Scalar 2 - WithoutDirect ###########
Memory consumption 854 Mb and Time 13.162695245
Time to solve model 5.002032299
Memory consumption 1065 Mb and Time 12.788777522
Time to solve model 4.613660219
########### Low-level API ###########
Memory consumption 304 Mb and Time 2.538083196
Time to solve model 4.572178837
Memory consumption 367 Mb and Time 2.641781411
Time to solve model 4.766848752
########### Vectorized - WithDirect ###########
Memory consumption 587 Mb and Time 12.18446959
Time to solve model 5.053543398
Memory consumption 482 Mb and Time 12.490777377
Time to solve model 5.095189864
########### Scalar 1 - WithDirect ###########
Memory consumption 405 Mb and Time 5.850087326
Time to solve model 4.78109017
Memory consumption 456 Mb and Time 5.842091906
Time to solve model 5.22244529
########### Scalar 2 - WithDirect ###########
Memory consumption 440 Mb and Time 6.361990639
Time to solve model 5.057280697
Memory consumption 440 Mb and Time 6.51641373
Time to solve model 4.437310264