Do we know why open source MIP solvers are significantly slower than commercial ones?

[SanPen]

I develop power system software, and part of that involves solving large MIPs.

While testing different solvers, I noticed that the open source ones (Highs, CBC, Scip) perform about the same, while the commercial ones (CPLEX, XPRESS, and Gurobi) perform about the same among them, but are about two orders of magnitude faster than the open source ones for the same MIP, sometimes yielding different solutions that also solve the problem.

This seems to be due to an algorithmic advantage rather than to programming techniques.

¿Do we know what commercial solvers are doing to be much faster?

I've noticed that for LP problems of the kind that I solve, all solvers perform comparably.

[jajhall]

Commercial MIP solvers benefit from person-decades of development investment - funded by high license fees. They have parallel tree search, and all sorts of tricks to handle different classes of problem. You may have seen "two orders of magnitude" difference in performance on your instances, but it's about one order of magnitude for Mittlemann's benchmarks between HiGHS and Gurobi - and not much more for SCIP and Cbc. HiGHS is looking to parallelise its tree search in 2024.

Again YMMV, but for LP on Mittlemann's benchmarks the best open-source solver - HiGHS - is 20 times slower than the best commercial solver (COPT). That's down to the idiosyncratic IPM implementation in HiGHS - which is why it's better than other open-source LP solvers - and parallelism. HiGHS is writing a new solver with the aim of reducing the gap.

[timellemeet]

Was a typo. Yes fully understand! Happy to help test some real world cases when we get there. I do column generation and solve the subproblems already on a GPU (using JAX). So it would be the ideal application I think. We are running many larger and smaller "jobs" each day by now in production, all powered by HiGHS and JAX. Happy to tell you more about if manage to make it to the workshop in June.

[jajhall]

Sounds exotic. It would be interesting to have you present your work at the workshop

[kutlay]

@timellemeet This sounds very exciting, is there any demo/presentation anywhere? Are you planning to pursue this as an open source project?

[raulhuatuco]

Was a typo. Yes fully understand! Happy to help test some real world cases when we get there. I do column generation and solve the subproblems already on a GPU (using JAX). So it would be the ideal application I think. We are running many larger and smaller "jobs" each day by now in production, all powered by HiGHS and JAX. Happy to tell you more about if manage to make it to the workshop in June.

sound interesting

[sk-surya]

@timellemeet, very interesting approach. I thought the cuPDLP-C works best for very large LPs, never thought of using it to solve many many small subproblems parallelly. How many subproblems are you solving at the same time? And what's the speedup you get when compared to solving the subproblems using CPU?

[anonyco]

I was browsing around the topic of optimization software when I stumbled upon this

Throwing in my two cents as someone with significant experience in high performance computing, what I saw browsing around this project was:

• No SIMD vectorization anywhere. (On the upside, at least HiGHs doesn’t use non-portable SIMD intrinsics. That’d be even worse than no SIMD at all.)
• Use of SparseArrays for computation. SparseArrays only make sense for storage and never net a profit for computation, which becomes extremely apparent in simd vectorization (which just isn’t possible with sparse arrays.) The only valid usecase I’ve seen for SparseArrays was in a long-running computation where storing intermediate calculations densely quickly gobbled up gigabytes of RAM and crashed; the solution in this case was to pack the intermediates into SparseArrays when they’re not in use and unpack them in parallel with the previous calculation, which netted a 230x speedup over the baseline non-SIMD SparseArray-based code just from those two optimizations. In all my other experiences with SparseArrays, I found the best solution to be redesigning the algorithm to cycle between sort+dedupe and intersect+accumulate, which can be functionally equivalent to SparseArrays for many use-cases.
• Deleteriously long memory dependency chains that completely obstruct the critical path everywhere.
• Ridiculously excessive use of static inline (which significantly pessimises performance more often than not), seemingly no consideration for structuring control flow to minimize non-tailcalls, and not one use anywhere of the most powerful performance optimization in existence: noinline-coerced tailcall chaining.
• No specialization for small integers (which are very common); it was all uint64_t all around
• switch is used inappropriately all over the place as some kind of performance cure all, when it’s actually one of the least efficient ways to do dispatch. Almost every switch I saw browsing through this project would benefit significantly being replaced with an indirect call into an array of function pointers.
• Extremely inefficient handling of various low-level integer routines such as for gcd and hashing
• Inefficient modular arithmetic that uses significantly more steps than necessary AND fails to leverage 128-bit multiplication cpu instructions.
• Liberal use of slow division by the same number and never computing the constant of inverse multiplication.
• Inappropriate usage of stdlib math functions in hot functions as most of these can’t be inlined with simd due to signaling nan and errno side effects, making any function that uses them incredibly slow.
• Excessive function call arguments in sooo many places and seemingly not one instance of recycling hoisted argument packing classes optimization.
• Extremely inefficient control flow in general, e.g. poorly organized excessive error checking complexity, poorly nested loops, unhoisted unpredictable branches, lack of branchless conditionals, etc
• General lack of debugging separation makes countless places significantly slower than necessary in production
• Generally excessively complex branchy check-all-conditions; very few places apply the concise (and much more optimizable due to nonobvious differences in unprovable behavior) bitwise logic equivalents.
• No collectivized memory allocation discipline. If you solve all the above performance issues, HiGHs will still be agonizingly slow even for just this one reason alone.
• Etc and etc

Based on my experiences optimizing other math-heavy software (albeit in different problem domains), I’d ballpark a 100-1000x potential speedup from these micro-optimizations, maybe even as much as a 10000x speedup depending on how heavily the sparse array is (ab)used. I’m not even joking about this: HiGHs looks extremely poorly optimized.

The insistence of @jajhall on algorithmic improvements to HiGHs in spite of the potential 100-10000x speedup of better micro-optimized C++ only serves to further reinforce my point that theres no magic rhyme or reason the performance of commercial solvers are better.

P.s. Also I think I happened to notice a concurrency bug in

HiGHS/highs/parallel/HighsSplitDeque.h

Line 213 in 3c26f33

splitRq = splitRequest.load(std::memory_order_relaxed);

on line 213. I believe it should be memory_order_acquire and that both are missing the second access of the double checked locking idiom (if my guess is right that that’s what you’re doing there.)

[jajhall]

Thanks for your extensive comments. If the best of commercial linear optimization software is 10x faster than HiGHS - much of which is due to reasons that we understand and are working on - and there is a potential 100-10000x speedup of better micro-optimized C++ in HiGHS, this implies that there's a 10-1000x speedup for the commercial solvers. Do you think that the developers of commercial solvers could also benefit from your expertise from "optimizing other math-heavy software... in different problem domains"?

SparseArrays only make sense for storage and never net a profit for computation

How sparse do you think the arrays are in HiGHS? Constraint matrix columns typically have 2-10 nonzeros in a vector that of dimension up to 10^7: are you suggesting that we scatter these into full-length vectors and then use SIMD vectorization? Similarly the solution vectors of linear systems in simplex have O(10) nonzeros in vectors of dimension up to 10^7. Should these be scattered and SIMD vectorization then used?

I'm not suggesting that everything you suggest lacks insight into the computational reality of linear optimization, and I'll share your observations with our developers (@filikat, @Opt-Mucca, @fwesselm), in particular the concurrency bug in line 213 of 3c26f33 in HiGHS/highs/parallel/HighsSplitDeque.h

[anonyco]

To answer your first question about expertise in commercial solvers, let me cite my recent experience with Lua. Lua is 33 year old software, so I’d expect its bytecode interpreter to be pretty efficient (and the only gains to be made in JIT like LuaJIT.) Infact, after just a few hours of tinkering and rewriting Lua’s dispatch loop to chain noinline-coerced tailcalls, I was able to realize a 1.4x boost in speed. That increased to 1.7x when I applied the preserve_none convention and unrolled instruction argument preloading 2 ahead. This further increased to 1.9x when I rewrote Lua’s primitive number ops to do all possibilities in parallel and cmov branchless select at the end. Honestly, I don’t understand how people can write such slow software (and for it to go unnoticed for so long) as it seems so intuitive to me how to make it fast.

Do you think that the developers of commercial solvers could also benefit from your expertise from "optimizing other math-heavy software... in different problem domains"?

I went to trade school (didn’t finish college) and got a computerized manufacturing degree. I cut pieces of metal to the specifications of blueprints for a living and it’s significantly better work than my expertise in software engineering, CI/CD, VLSI, systems architecture, Linux sysadmin, fullstack, DevOps, embedded, and HPC ever landed me, ha ha. So, no, I don’t do any paid tech work. I do, however, volunteer for FOSS.

How sparse do you think the arrays are in HiGHS? Constraint matrix columns typically have 2-10 nonzeros in a vector that of dimension up to 10^7

Wow!, I didn’t realize they were that sparse. To address your question, I’ll explain two different ways of handling sparse arrays.

1. Normally for medium sparse data ~0.01%-10% occupant (which your case doesn’t fit), the trick is to have two side-by-side sorted arrays—one for values and one for indices. Depending on the algorithm you’re performing on the sparse data, having the indices sorted and SIMD processing them in from lowest to largest indice reveals all kinds of optimization opportunities that can (again depending on what you’re doing) approach the speed of dense data processing. Then, instead of scattering, we redesign the algorithm to alternate intersect+accumulate (where we process the sparse data start-to-end, intersecting and substituting new updated sparse indices from the previous round) and sort+dedupe (where we take the gathered sparse output of the last round and put it into sorted order for the next.)
2. In your maximally sparse use-case, however, the far better approach is to specialize two versions of the code. One version handles the typical case of (say) <=16 occupant doubles, whereas the other case handles >16 doubles. The fast case with <=16 can be many times faster as we can eliminate all looping constructs and always perform 16 start-to-end SIMD calculations. When the actual occupancy is <16, we simply discard those unused junk calculations at the end.

Sorry if I’m not doing a good job of explaining this. It’s very difficult to describe as this is a microoptimization of how the C++ is written, not a fundamental change in the algorithm being performed.

I'm not suggesting that everything you suggest lacks insight into the computational reality of linear optimization,

I skimmed HiGHs very quickly and wrote my comment in a few strokes. I’m surprised you didn’t find more errors, ha ha!

I'll share your observations with our developers

I’d be happy to lend a hand in any way I can! (At their direction of course-I’ll work on whatever part they want me to work on however they want it to be.) I enjoy contributing to FOSS projects.

[jajhall]

Wow!, I didn’t realize they were that sparse.

The linear system solves in simplex are exceptionally sparse. I once tried seeing whether a parallel system solver (MUMPS - which is used commercially) might give better performance. It didn't, and its author was astonished at how sparse our computations were.

With $m$ constraints, the cost of a single simplex iteration with dense arithmetic is $O(m^2)$. When exploiting sparsity, that normally drops to at most $O(m)$, and for many "hyper-sparse" problems - with 10s of nonzeros in the linear system solutions, regardless of $m$ - it's $O(10)$.

To address your question, I’ll explain two different ways of handling sparse arrays.

1. Normally for medium sparse data ~0.01%-10% occupant (which your case doesn’t fit), the trick is to have two side-by-side sorted arrays—one for values and one for indices. Depending on the algorithm you’re performing on the sparse data, having the indices sorted and SIMD processing them in from lowest to largest indice reveals all kinds of optimization opportunities that can (again depending on what you’re doing) approach the speed of dense data processing. Then, instead of scattering, we redesign the algorithm to alternate intersect+accumulate (where we process the sparse data start-to-end, intersecting and substituting new updated sparse indices from the previous round) and sort+dedupe (where we take the gathered sparse output of the last round and put it into sorted order for the next.)

This is essentially what we do at the moment. Below 10% sparsity for the solution vectors, we just hold the values in a full-length array, as the overhead of tracking the indices becomes costly.

For the constraint matrix, the occasional dense (or even full) columns are held as packed vectors. Yes, we might gain a bit by having specialised data structure for them, but at the cost of code complexity. This is the sort of optimization that the commercial solvers may employ. However, the gains aren't going to big.

Essentially, if what you assert applied to linear optimization broadly, the multi-million dollar commercial companies like Gurobi and COPT (whose solvers have seen person-decades of expert, highly-paid development in a very competitive sphere) would use it, and their solvers would be 100x-1000x faster than HiGHS. They aren't.

Still, you may have suggested some local optimizations that might be good to apply "in principle". I'm far from fluent in C++, so I'll wait to see what our developers think.

[anonyco]

I hear everything you’re saying, and I’m going out on a limb to challenge you to consider how much branching was involved. You can do all the smart tricks and data structures out the wazoo, and they’ll only net negligible speedups if you don’t take care of branching and carefully design your code to dance around a fast path with few-to-no branches. This is the two-fold reason why simd gives such ridiculous speedups: (1.) simd processes many items simultaneously and (2.) writing simd means minimizing branching, which speeds up all computation. You’d be shocked at the numbers of large projects I’ve seen with person-decades of investment that don’t understand basic optimization like this. E.g. I sped up OpenJDK11 by 1.2x(? might have been 1.3x?) faster a few years ago with just a perl script to reorganize the worst offending branches across its codebase.

The only other thing I’ll say is that you don’t need ridiculous code complexity specializing every condition to be ridiculously fast. That’s a demonstrably false fallacy. E.g. in rare cases I do have to practically specialize a bazillion custom variants of my code to handle all logic efficiently, I write one version of the code with everything, then I leverage compiler optimizations to automatically rewrite subsections of it optimally for specific cases.

Let’s wait to hear back what your devs say. 👍

[filikat]

Thanks for your comments. I agree that many of the things that you mention can be done better in HiGHS. However, I would also agree with Julian that the expected speed-up is smaller than what you suggest.
Some of the things that you mentioned caught my attention. Do you have any advice/guideline/further reading regarding:

• how to obtain SIMD vectorization, without using intrinsics specific to a given architecture.
• how to restructure the code to exploit tailcall optimization.
• what are "noinline-coerced tailcall chaining" and "recycling hoisted argument packing classes optimization".
• what you mean by "collectivized memory allocation discipline".

Regarding the bug in HighsSplitDeque.h, I would need to ask Leona, who wrote that part of the code. My understanding of the scheduler is not good enough to comment on this. However, at the moment all loads and stores to splitRequest are relaxed. I suppose you mean that some of the stores to it are wrong as well, and should be marked release.