Performance: Julia/JuMP vs. Python/Pyomo

[Justine]

Dear all,

I could use some good help from the experts here. I was about to translate an optimization model from Python/Pyomo to Julia/JuMP, expecting that Julia/JuMP would outperform Python’s Pyomo. However, model generation in JuMP took an unexpectedly long time.

So I came up with a small performance experiment to compare model generation time in Pyomo and JuMP that I want to share with you and get advice on what I can improve.

Setting:
Windows 11 Intel Core i7-1260P CPU @ 2.1 GHz 48GB RAM

Julia 1.8.3
JuMP 1.5.0
Python 3.8.15
Pyomo 6.4.4

I want to generate the following model with only a single constraint and a constant objective function:

$$\sum_{(i,j,k) \in \mathcal{IJK}} \sum_{(j,k,l) \in \mathcal{JKL}} \sum_{(k,l,m) \in \mathcal{KLM}} x_{i,j,k,l,m} \ge 0 ; \forall ; i \in \mathcal{I}$$

(sorry - only allowed to post one media item)

I generated the data used in the experiment in Python as follows:

def create_random_data(n, m=20):

    I = [f'i{x}' for x in range(1, n + 1)]
    J = [f'j{x}' for x in range(1, m + 1)]
    K = [f'k{x}' for x in range(1, m + 1)]
    L = [f'l{x}' for x in range(1, m + 1)]
    M = [f'm{x}' for x in range(1, m + 1)]

    IJK= pd.DataFrame(np.random.binomial(1, 0.05, size=(len(I)*len(J)*len(K))), 
                       index=pd.MultiIndex.from_product([I,J,K], names=['i', 'j', 'k']), 
                       columns=['value']).reset_index()
    JKL= pd.DataFrame(np.random.binomial(1, 0.05, size=(len(J)*len(K)*len(L))), 
                       index=pd.MultiIndex.from_product([J,K,L], names=['j', 'k', 'l']), 
                       columns=['value']).reset_index()
    KLM= pd.DataFrame(np.random.binomial(1, 0.05, size=(len(K)*len(L)*len(M))), 
                       index=pd.MultiIndex.from_product([K,L,M], names=['k', 'l', 'm']), 
                       columns=['value']).reset_index()

    ijk = [tuple(x) for x in IJK.loc[IJK['value'] == 1][['i', 'j', 'k']].to_dict('split')['data']] 
    jkl = [tuple(x) for x in JKL.loc[JKL['value'] == 1][['j', 'k', 'l']].to_dict('split')['data']] 
    klm = [tuple(x) for x in KLM.loc[KLM['value'] == 1][['k', 'l', 'm']].to_dict('split')['data']] 
    
    return I, J, K, L, M, ijk, jkl, klm

The Julia model implementation:

using JuMP
using BenchmarkTools

function model_ei(I, IJK, JKL, KLM)
    model = Model()

    x_list = [
        (i,j,k,l,m) 
        for (i,j,k) in IJK
        for (jj,kk,l) in JKL if jj == j && kk == k
        for (kkk,ll,m) in KLM if kkk == k && ll == l
    ]

    @variable(model, x[x_list] >= 0)

    @constraint(model, [i in I],  sum(
                                        x[(i,j,k,l,m)] 
                                        for (ii,j,k) in IJK if ii == i
                                        for (jj,kk,l) in JKL if jj == j && kk == k 
                                        for (kkk,ll,m) in KLM if kkk == k && ll == l 
                                    ) >= 0
    )
end

N = [n for n in 50:50:2000]

for n in N
    I, J, K, L, M, IJK, JKL, KLM = read_data(n)
    tmi = @belapsed model_ei($I, $IJK, $JKL, $KLM) samples=7 evals=10
end

I, J, K, L, M, IJK, JKL, KLM are all vectors.

The Pyomo model implementation:

import pandas as pd
import numpy as np
import pyomo.environ as pyo

def pyomo_mi(i, ijk, jkl, klm):    
    model = pyo.ConcreteModel()

    model.IJK = pyo.Set(initialize=ijk)
    model.JKL = pyo.Set(initialize=jkl)
    model.KLM = pyo.Set(initialize=klm)

    model.z = pyo.Param(default=1)

    model.x = pyo.Var([(i,j,k,l,m) for (i,j,k) in ijk 
                       for (jj,kk,l) in jkl if (jj == j) and (kk == k) 
                       for (kkk,ll,m) in klm if (kkk == k) and (ll == l)], 
                      domain=pyo.NonNegativeReals)

    model.OBJ = pyo.Objective(expr=model.z)
    
    model.ei = pyo.Constraint(i, rule=ei_rule)

for n in N:
    I, J, K, L, M, IJK, JKL, KLM = read_data(n)
    tmi = %timeit -o -r 7 -n 10 pyomo_mi(I, IJK, JKL, KLM)

I, J, K, L, M, IJK, JKL, KLM are all list of tuples.

The result for an increasing cardinality of IJK with constant cardinality for J, K, L, M, JKL, KLM:

Can anybody explain why Julia/Jump performs so much worse than Python/Pyomo? Or tell me what I am doing wrong?

Best
Justine

[Christopher_Fisher]

I’m not familiar with JuMP or Pyomo, but I will recommend a few general things to consider. Perhaps there are differences default settings for numeric precision (e.g. Int64 vs Int32), automatic multithreading, or the use of MKL. Another thing you might consider is whether you have any type instabilities. @code_warntype can be helpful for this purpose. Aside from that, I will have to defer to someone who has more expertise, but this might give you a starting point.

[odow]

Hi @Justine, welcome!

In general, I wouldn’t expect JuMP to have this performance difference compared to Pyomo.

I couldn’t run your code because you didn’t provide read_data for the JuMP code.

Did you try this formulation for the constraint?

    @constraint(
        model,
        [i in I], 
        sum(x[k] for k in x_list if k[1] == i) >= 0
    )

You should also try timing how long it takes for the solver to start solving the problem, not just how long it takes to build without the solver. JuMP and Pyomo access solvers very differently, so it might be that Pyomo is quicker to build this part, but slower to get the actual problem to the solver.

but I will recommend a few general things to consider

@Christopher_Fisher provides good advice for general Julia programs, but they’re not relevant in this particular case. JuMP and Pyomo always use int64, and neither build the problem with multithreading.

[Justine]

Hi @odow, thanks for your reply.

I did not provide the read_data for the JuMP code as it is very quick and dirty. I wanted Pyomo and JuMP to work on the exact same data. Thus, I dumped the randomly generated data from Python to JSON to make it available to Julia.

However, your suggestion worked perfectly. My first guess was that it violated the model logic - but it didn’t. I adjusted the Pyomo code using the x_list to make it fair. The adjustment didn’t help Pyomo much in terms of performance. However, it helped JuMP dramatically. The new results look as follows:

I will also follow your advice and look at model generation + model solving in the coming days and will give an update for anyone interested.

[Jack_Chou]

I hope you can share the complete code for this test, I am very interested in it, thx.

[c-ferny]

@Justine
Can you provide complete code for this test. This is a very practical study that would be of great help. Thanks.

[Justine]

Hey folks,

Apologies for the delayed response. I’ve managed to find some time to make progress on the project here and there.

As requested, I have created a GitHub repository where you can find the complete code. You can access it here. I would be delighted if anyone is interested in contributing to it.

I want to mention that the example I chose, which I’ll refer to as the IJKLM model, is primarily a showcase example and doesn’t have any real-world applications. In an attempt to explore the model generation behavior, I decided to now focus on a small-scale model within the supply chain domain. My goal was to compare the model’s performance against both an easy and intuitive Pyomo implementation and a more efficient but less intuitive implementation.

Unfortunately, I haven’t been able to improve the Pyomo performance thus far. I also attempted similar optimizations using JuMP (blue line), but it seems that the proposed improvements that worked for the IJKLM model did not yield the desired results for the supply chain model.

Here is a small snippet from the intuitive JuMP implementation (orange line):

function intuitive_jump(I, L, M, IJ, JK, IK, KL, LM, D, solve)
    model = Model()

    x_list = [
        (i, j, k)
        for (i, j) in IJ
        for (jj, k) in JK if jj == j
        for (ii, kk) in IK if ii == i && kk == k
    ]

    y_list = [
        (i, k, l)
        for i in I
        for (k, l) in KL
    ]

    z_list = [(i, l, m) for i in I for (l, m) in LM]

    @variable(model, x[x_list] >= 0)
    @variable(model, y[y_list] >= 0)
    @variable(model, z[z_list] >= 0)

    @constraint(model, production[(i, k) in IK], sum(
        x[(i, j, k)]
        for (ii, j) in IJ if ii == i
        for (jj, kk) in JK if jj == j && kk == k
    ) >= sum(y[(i, k, l)] for (kk, l) in KL if kk == k)
    )

    @constraint(model, transport[i in I, l in L], sum(
        y[(i, k, l)] for (k, ll) in KL if ll == l
    ) >= sum(z[(i, l, m)] for (ll, m) in LM if ll == l))

    @constraint(model, demand[i in I, m in M], sum(
        z[(i, l, m)] for (l, mm) in LM if mm == m
    ) >= D[i, m])

end

I would greatly appreciate any recommendations or insights from the experts in the community. If you have any suggestions on how to enhance the performance of the supply chain model with JuMP, please feel free to share your thoughts.

Looking forward to your input!

Best regards, Justine

[odow]

That plot doesn’t look right! Let me take a look.

[odow]

So problem is actually Julia-related. You data structures have a lot of Vector{Vector{String}}, where the inner Vector{String} has two or three elements. Then, we spend a lot of time doing for (ii, j) in IJ where we convert the vector ["ii", "j"] into the tuple ("ii", "j").

I get a 3-4x speed-up if I change the data structure to use tuples:

function read_tuple_list(filename)
    return map(JSON.parsefile(filename)) do x
        return tuple(x...)
    end
end

function read_fixed_data()
    N = open(JSON.parse, "supply_chain/data/data_N.json")
    L = open(JSON.parse, "supply_chain/data/data_L.json")
    M = open(JSON.parse, "supply_chain/data/data_M.json")
    JK = read_tuple_list("supply_chain/data/data_JK.json")
    KL = read_tuple_list("supply_chain/data/data_KL.json")
    LM = read_tuple_list("supply_chain/data/data_LM.json")
    return N, L, M, JK, KL, LM
end

function read_variable_data(n)
    I = ["i$i" for i in 1:n]
    IJ = read_tuple_list("supply_chain/data/data_IJ_$n.json")
    IK = read_tuple_list("supply_chain/data/data_IK_$n.json")
    d = read_tuple_list("supply_chain/data/data_D_$n.json")
    D = Dict()
    for x in d
        D[(x[1], x[2])] = x[3]
    end
    return I, IJ, IK, D
end

Your “JuMP” case isn’t really comparable with the others because it’s not doing the same operations. This sum sum(x[p] for p in x_list if p[1] == i && p[3] == k) loops over many terms when most of them are false.

Another issue is that JuMP is adding a lot of 0 >= 0 constraints. It looks like in Pyomo you have

 if lhs or rhs:
        return sum(lhs) >= sum(rhs)
    else:
        return pyo.Constraint.Skip

but there isn’t anything similar in JuMP. You’d have to manually write a for-loop to check for empty summation terms.

there are a couple of other smaller wins. You could use direct_model, you could disable passing names to the solver, and you could avoid creating explicit containers of the constraints. These have a ~20% win for me.

function jump(I, L, M, IJ, JK, IK, KL, LM, D, solve)
    model = if solve == "True"
        direct_model(Gurobi.Optimizer())
    else
        Model()
    end
    set_string_names_on_creation(model, false)
    x_list = [
        (i, j, k)
        for (i, j) in IJ
        for (jj, k) in JK if jj == j
        for (ii, kk) in IK if ii == i && kk == k
    ]
    y_list = [
        (i, k, l)
        for i in I
        for (k, l) in KL
    ]
    z_list = [(i, l, m) for i in I for (l, m) in LM]
    @variable(model, x[x_list] >= 0)
    @variable(model, y[y_list] >= 0)
    @variable(model, z[z_list] >= 0)
    for (i, k) in IK
        @constraint(model, 
            sum(
                x[(i, j, k)]
                for (ii, j) in IJ if ii == i
                for (jj, kk) in JK if jj == j && kk == k
            ) >= sum(y[(i, k, l)] for (kk, l) in KL if kk == k)
        )
    end
    for i in I, l in L
        @constraint(model,
            sum(y[(i, k, l)] for (k, ll) in KL if ll == l) >= 
            sum(z[(i, l, m)] for (ll, m) in LM if ll == l))
    end
    for i in I, m in M
        @constraint(model, sum(z[(i, l, m)] for (l, mm) in LM if mm == m) >= D[i, m])
    end
    if solve == "True"
        set_silent(model)
        optimize!(model)
    end
end

But I’ll let you re-run the benchmarks on your machine so that it’s fair.

[odow]

If you do the same with the IJKL example, you also get a massive speedup:

function read_tuple_list(filename)
    return map(JSON.parsefile(filename)) do x
        return tuple(x...)
    end
end

function read_fixed_data()
    N = open(JSON.parse, "IJKLM/data/data_N.json")
    JKL = read_tuple_list("IJKLM/data/data_JKL.json")
    KLM = read_tuple_list("IJKLM/data/data_KLM.json")
    return N, JKL, KLM
end

function read_variable_data(n)
    I = ["i$i" for i in 1:n]
    IJK = read_tuple_list("IJKLM/data/data_IJK_$n.json")
    return I, IJK
end

[Justine]

Shoutout to @odow again for his suggestions. These are the new plots:

IJKLM

Supply Chain

[odow]

That looks better! The takeaway is that Julia does not like lots of small Vector{Any} datastructures.

Compared with Python you need to be a little more careful with your software engineering. Its easy to write slow code and its easy to write fast code, and the difference can be subtle. If you’re writing a paper on this, then thats a fair downside to point out.

[hellemo]

For high performance in JuMP for sparse problems you might have a look at SparseVariables.jl which was designed for these kinds of problems. There is just a small advantage for the largest problem sizes here, but if you find indexing with tuples cumbersome and/or want to do slicing, it can offer greater advantages both in performance and ergonomics. (see our benchmarks that showcase this here: GitHub - hellemo/SparseVariables.jl: Efficient sparse modelling with JuMP)

See PR Add JuMP version using SparseVariables by hellemo · Pull Request #4 · justine18/performance_experiment · GitHub for details for the supply chain example here.

JuMP done 5600 in 1.09s
Sparse JuMP done 5600 in 1.18s
Intuitive JuMP done 5600 in 1.42s
JuMP done 6700 in 1.47s
Sparse JuMP done 6700 in 1.5s
Intuitive JuMP done 6700 in 2.03s
JuMP done 7900 in 2.02s
Sparse JuMP done 7900 in 1.84s
Intuitive JuMP done 7900 in 2.67s

[trulsf]

With SparseVariables it is possible to use slicing with a minor overhead, leading to a very compact constraint formulation

for (i, k) in IK
    @constraint(model, sum(x[i, :, k]) >= sum(y[i, k, :]))
end

for i in I, l in L
    @constraint(model, sum(y[i, :, l]) >= sum(z[i, l, :]))
end

for i in I, m in M
    @constraint(model, sum(z[i, :, m]) >= D[i,m])
end

[odow]

There’s also this question on StackOverflow; performance - Is there a more efficient implementation for Pyomo models compared to the current intuitive one? - Stack Overflow

The Pyomo folks point out that the choice of data structure overall is poor because of the sparsity. (There are a lot of 0 >= 0 rows.)

But I assumed that was mostly just an artifact of the toy problem and fake data?

[Justine]

I followed the advice given by the Pyomo folks and also implemented your suggested changes and they indeed have a great impact. These are the new plots
IJKLM

Supply Chain

@odow Regarding the sparsity of data: Yes this is an artifact of the IJKLM example. I will implement a more realistic data generator - hopefully soon.

[odow]

Is this model generation time? Interesting the Pyomo is faster for the IJKLM example.

You should really time including solve, perhaps with a time limit of 0.0, so that we time how long it takes for the problem to actually make it to the solver.

[Justine]

I did. The supply chain model did not show any changes. JuMP was faster on this one anyway. Here are the results for the IJKLM model:

Model Generation

Model + Solve

So JuMP is slower than Pyomo if we only measure model generation. If we include the solve statement with a time limit of 0, JuMP is faster again.

[miles.lubin]

EDIT: This thread is cited from a blog post by GAMS: Performance in Optimization Models: A Comparative Analysis of GAMS, Pyomo, GurobiPy, and JuMP. It was originally cited in a way that suggested that JuMP developers collaborated on the blog post. The blog post has since been corrected. I’ve locked the thread to preserve the discussion in its original context.

The JuMP developers’ response: JuMP, GAMS, and the IJKLM model | JuMP.