Support for matrix style model construction (MathProgBase/LinearQuadratic style) in modern Julia + MOI

[lpaizs]

I’m reaching out about support or example workflows for constructing optimization models from raw matrix data in modern Julia, similar to how it was done with MathProgBase’s LinearQuadratic interface.
A bit of background on my situation:
• Years ago I implemented a large scale unit commitment MILP in Julia 0.6.4, using MathProgBase.jl and specifically the LinearQuadratic models interface, together with Gurobi as the solver.
• The model is truly large: on the order of ~2 million variables and constraints. I built the model from big matrices and vectors (A, b, c, etc.) in a Matlab style fashion and handed that directly to the solver.
Now I’d like to move the codebase to a current Julia version, but without changing the overall program structure—i.e., keeping a MathProgBase / LinearQuadratic style formulation rather than rewriting everything in JuMP or another high level modeling layer.
Question / request:
Is there any officially supported or recommended way in the current Julia / MOI ecosystem to construct a large LP/MILP directly from matrices and pass it to a solver, similar to the old MathProgBase LinearQuadratic approach?
Thanks for any pointers or suggestions. Even a small example or outline of how to plug matrix data into MOI, or a note on current limitations, would be very valuable.

[blegat]

We've been advising people to just go with @constraint(model, A * x <= b) and check if that's really a performance issue. We also developed this replacement for that specific use case of MathProgBase in this package: https://github.com/jump-dev/MatrixOptInterface.jl. Let us know if you get into any trouble with it.

[lpaizs]

Hi, thanks so much for the quick reply!

Really appreciate the pointer. I’ll give MatrixOptInterface.jl a try and see how it works with my big model.

[blegat]

Wait, I thought this issue was opened on the MathOptInterface package but I don't see any trace that the issue was transferred. Sorry, then my answer is a bit weird if that's the case ^^

[odow]

Sorry I'm a bit late to this issue. Check out this tutorial: https://jump.dev/MathOptInterface.jl/stable/tutorials/mathprogbase/

[lpaizs]

Hi, sorry for the confusion — the issue is indeed on the MatrixOptInterface side, not MathOptInterface.
We tested MatrixOptInterface in detail for large-scale models and found that the bottleneck is the copy_to step, which scales mainly with the number of constraints and becomes very slow at ~10^6 scale.
Happy to repost benchmark details if useful.

[lpaizs]

@odow Thanks for the link! This looks related to our original MathProgBase transition question from December.
Our current issue is about copy_to performance at large scale, which we described above.

[blegat]

So it's the copy_to from a model in here to Gurobi.Optimizer ? What you pin point a method that is being the bottleneck ? Our maybe share the flamegraph ?

[lpaizs]

yes — the bottleneck I see is specifically the MOI.copy_to(opt, lp) where lp is a MatrixOptInterface.LPForm (matrix form) and opt is Gurobi.Optimizer().

I isolated it by timing components separately (matrix build + LPForm creation are negligible; integrality constraints are basically free; copy_to dominates and scales mostly with the number of constraints).

I don’t yet have a flamegraph, but I can produce one with Julia’s profiler around MOI.copy_to and share it. Is an SVG/HTML flamegraph ok, or would a “top hotspots” text summary be enough?

If helpful, I can also paste the minimal benchmark snippet we used (varying #constraints with fixed nnz/row) so you can reproduce locally.

[lpaizs]

Yes, exactly: bottleneck is MOI.copy_to(opt, lp) where lp is MatrixOptInterface.LPForm and opt is Gurobi.Optimizer().

I profiled a “copy_to-only” reproducer (lp built once; then just MOI.copy_to). The hotspots are in MOI’s generic copy path:
MOI.Utilities.copy.default_copy_to -> pass_nonvariable_constraints_fallback / _copy_constraints.

Inside that, a lot of time is spent (1) extracting rows from SparseMatrixCSC via getindex(A, i, J::UnitRange) (tons of searchsortedfirst), and (2) calling Gurobi.add_constraint(ScalarAffineFunction, LessThan) which ends up in GRBaddconstr (appears in the profile), i.e., constraints are added one-by-one.

This matches our benchmarks: runtime scales mostly with #constraints. I attached a flamegraph SVG.

[lpaizs]

I also generated the same profile using StatProfilerHTML (interactive HTML), which makes the hotspots much clearer. I’ve uploaded that version as well in case it’s helpful.

index.html

[blegat]

It seems this is due to getting rows in SparseMatrixCSC that is slow. Getting a row in a transpose of SparseMatrixCSC would be equivalent to getting a column which might be much faster. Does it help to pass copy(A')' instead of A to the LPForm ? From a quick look, it does not seem that SparseArrays has a shortcut for getting a row form a transpose of sparse matrix so I don't think that it will help :/ We should add a special case for transpose matrices here

[lpaizs]

Thanks! We tried A_fast = copy(transpose(A))' (i.e., passing copy(A')'-style) but for us it got worse:
original copy_to ~43s vs A_fast ~73s at (nvar=ncon=100k, nnz/row=5).

So it seems the Adjoint/transpose trick doesn’t help in our case (maybe it triggers slower generic getindex paths).
Next we’ll try a CSR (row-oriented) sparse matrix for A (e.g., via SparseMatricesCSR.jl) since the profile suggests heavy row extraction.
I’ll report back with numbers.

[lpaizs]

We tested the transpose trick and also a CSR matrix, but both got worse on my side.

At (nvar=ncon=100k, nnz/row=5):

• baseline (SparseMatrixCSC): copy_to ~42.85s
• passing A_fast = copy(transpose(A))' : copy_to ~73s
• passing CSR (SparseMatricesCSR.SparseMatrixCSR(A)): copy_to ~78.5s

So it seems the slowdown isn’t fixed by changing the matrix storage. From the profile it still looks like the dominant cost is MOI’s generic constraint-copy path ending up adding constraints one-by-one (add_constraint -> GRBaddconstr) plus lots of small allocations / index mapping.

Given this, my conclusion is that for million-scale constraints we likely need a bulk-loading approach (Gurobi API GRBaddconstrs / GRBaddvars style) rather than copy_to.
Happy to share more profiling info if useful.

[odow]

Given this, my conclusion is that for million-scale constraints we likely need a bulk-loading approach (Gurobi API GRBaddconstrs / GRBaddvars style) rather than copy_to.

People routinely build and solve problems of 10^6 scale using JuMP without issue.

I think your issue is first creating the sparse CSC, and then doing a copy_to.

My strong suggestion would be that you just use JuMP.

This package has not had much development time put into it.

How long does it take to build your input matrices? Can you share a minimal reproducible example? What is your standard form?

[blegat]

Yes, that aligns with my earlier observation that

From a quick look, it does not seem that SparseArrays has a shortcut for getting a row form a transpose of sparse matrix so I don't think that it will help :/ We should add a special case for transpose matrices here

so until we add a special case for transpose matrices here, it won't help. Until then, I think it will actually be faster to just use JuMP and do @constraint(model, A * x <= b) ?
I can take a look at making it fast with the transpose after #35 is merged. Then it will use MOI.Utilities.MatrixOfConstraints. This issue is a good motivation for using MOI.Utilities.MatrixOfConstraints instead of redefining new structs