Is it time for a native Julia LP(Simplex) solver?

For the last few months, I have been coming up to speed on both Linear Programming with Jump.jl using Clp.jl, Cbc.jl, and GLPK.jl. I have made some minor PRs to a few of the projects in order to improve the capability to solve my own problem.

But as a dig into each of those problems, I keep getting the feeling that Julia needs a native LP solver. Initially, this would most likely take the form of some flavor simplex solver. But with that foundation, I think a lot of progress could be made. From there a MILP solver could be built, and after that, there already exists higher level julia solvers. But to be performant it is clear that you need very low-level interoperability to start building those blocks. Just like Cbc has its own interfaces inside Clp.

Now a few years ago it appears that @miles.lubin started on an initial pass with jlSimplex. Though Julia has progressed a ton since then and I feel like with v1.0 out the capabilities are there finally there to make a native state of the art solver.

At the end of the day, I ask this because it is far outside my capabilities to start a project like this, but I would love to contribute. I am sure many of you have had a similar thought and it would be great to hear your opinions.

11 Likes

TL;DR: GitHub - JuliaLinearOptimizers/Simplex.jl

From a few conversations about this subject, the main problem is if it’s really needed. There are decent open source alternatives and professional softwares, all wrapped for JuMP, and dedicating time to developing a LP solver in Julia is not easily justifiable.

That being said, I decided to start a solver with my student @RenanOD, so that I could supervise a few students on slightly more advanced topics on LP. Hopefully, it becomes stable, efficient and robust enough to be used in research. Any help is welcome, short- or long-term.

The ominous organization name JuliaLinearOptimizers is to follow JuliaSmoothOptimizers naming, which is where most of my focus is right now.

13 Likes

Some already exist: https://github.com/oxfordcontrol/COSMO.jl, GitHub - ds4dm/Tulip.jl: Interior-point solver in pure Julia

5 Likes

I’ve toyed with the idea a little bit. I agree with @abelsiqueira in that it’s a little hard to justify considering how many alternatives there are, so if it happens, it’ll probably start out at least as being purely for the fun of it.

Writing performant simplex and interior point solvers for linear problems shouldn’t be all that difficult (indeed, it’s nice to see how few lines of code are in Simplex.jl). The problem is, usually what you want is something a little more complicated. In my job I mostly deal with mixed integer problems, and writing really good branch and bound routines seems like an aboslute nightmare, my naive guess is that it would be about 100 times the level of effort of writing the relaxation solvers to get something competitive with existing implementations (though hopefully the wonder of Julia would reduce that somewhat). There may also be lots of really sophisticated pre-solving routines that I’m not aware of, so I have no idea the amount of effort that would be involved in writing routines to simplify the initial problem, which in many cases can make a huge difference.

I do find myself extremely frustrated at how the open source solvers compare to the commercial ones (they tend not even to be very well documented), so I’ve spent a fair amount of time lately fantasizing about a really good open source solver written in pure Julia.

I think the a good starting point would be having really solid pure Julia implementations of simplex and interior point methods for purely linear problems that also have many of the important features that would be needed for them to ultimately be used in branch and bound such as warm starts. Looks like we are slowly getting there.

(By the way, I am most certainly not an expert in this subject, just someone who has at least a basic idea of what writing a solver would entail.)

4 Likes

It is funny I had the exact opposite view. When I dug into how to make a performant simplex solver I quickly felt over my head, but looking at branch and cut routines that made clear sense.

Beyond, that I get that I am effectively proposing that we should recreate/reinvent what is already in Clp, which is about 170k lines of c++. But If there is a good nugget to start around things can grow quickly.

@abelsiqueira makes a great point about open source alternative, but the same could have been said Julia when octave, matlab, and python were already around. If there was a vision to have a state of the art solver, the academic efforts could be contributing to that. Because I am sure you have noticed that the ease of hacking on Julia makes the learning curve for other domains come down as well.

2 Likes

Beyond, that I get that I am effectively proposing that we should recreate/reinvent what is already in Clp, which is about 170k lines of c++

Writing a competitive Clp solver is very difficult. Here is a (C++) solver that aims to be a Clp replacement GitHub - ERGO-Code/HiGHS: Linear optimization software

2 Likes

As you can imagine, I’m fan of simplex implementations in Julia. However, you haven’t really made clear the motivation for asking someone else to start a PhD-thesis worth of work on a native Julia simplex solver. What would a native Julia version do that existing solvers can’t? If there are small tweaks needed to existing solvers to make them useful to you (e.g., patches to expose functionality you need from Julia), why not do those instead?

2 Likes

I think being able to easily change the arithmetic precision of the algorithm would be a big advantage for a pure Julia implementation. From time to time it would be very helpful to have an LP solver with exact rational arithmetic.

7 Likes

Just to be clear my motivation here is essential to be proactive, not in the troll sense, but in the sense that I wanted to plant this in everyones minds. Because I know many of you are headed to JuMPcon, I will not, but I figure that will be a great opportunity for a provocative idea like this to be discussed.

That said alot of what drove me to finally ask this question is @ jeff.bezanson, @ StefanKarpinski, & @ viralbshah award presentation: SIAM CSE19: Solving the Two Language Problem in Scientific Computing and Machine Learning with Julia - YouTube

In there they make a few great points for example, they mention tensorflow vs flux:
tensorflow( ~1.5million LOC) ----> Flux.jl( 4300 LOC)

To your point when I can I am making those improvements to the existing solvers.

I agree somewhat with the branch-and-bound, but cuts are also difficult and technical to implement.
Also, don’t underestimate the sophistication of presolving and bound tightening routines, when implemented efficiently.

I feel that the most obvious MIP-solver ingredient for Julia implementation are primal heuristics, which can often be added to existing solvers via callbacks (so, no separate branch-and-bound implementation is required).

1 Like

My understanding is that the simplest possible implementation of a branch and bound solver should not be that hard, but for performance those usually rely on a wide variety of different heuristics, all of which must be implemented (very well) for good performance in the general case.

It’s possible that my perception that writing a continuous linear simplex or interior point solver wouldn’t be all that bad is ignoring major performance issues. I feel like I could probably figure one out in a few weeks, but perhaps whatever I’d come up with would have terrible performance (I also find interior point easier to understand than simplex).

Simplex in exact arithmetic is only a third of a PhD thesis :slight_smile: .

4 Likes

There is https://github.com/IainNZ/RationalSimplex.jl

2 Likes

Regarding the implementation of a competitive dual simplex method, Bob Bixby has held a 3-part lecture some years ago (video). It does not explain all critical ingredients, but references additional papers.

I also considered a Julia implementation, mostly as a learning exercise, but I now think that it would be more valuable to implement more MOI bridges or JuMP extensions, instead.
Or else, maybe support the maintenance of existing open-source solvers (e.g. Cbc).

Interesting thread…

It strikes me that given the broad interest in an open source performant MIP solver… why is this not manifesting as improvements to, say, Cbc for example? If all this interest was focused here, surely significant improvements could be made.

The fundamental here is that in the last decade, enormous strides have been made in open source tools… now a decent MIP solver is, perhaps, the last barrier to full open source optimisation for real world problems… now is definitely the time to focus efforts into the right place to get the progress we need.

2 Likes

If someone wants to give it a try, can I suggest e.g. Crash in small IP instance using the C interface and Cbc_setInitialSolution · Issue #202 · coin-or/Cbc · GitHub, which has a quite minimal reproducer and direct impact on the Julia bindings?

I did make a stab at it myself but lost heart when I had tracked the crash into CbcMain1 and found out that the function was 9000 lines of C++ code. At that point I rather wished that Cbc had been written in Julia instead.

7 Likes

Reading this thread at a very high level, it seems that there’s two perspectives:

  1. Really amazing solvers already exist that are not written in Julia, they’re usually written in C++. Why should someone go to the effort of reimplementing all of that in Julia rather than just making improvements to the existing projects?

  2. C++ is a a great zero-overhead language but so can Julia be and Julia tends to be much more concise and much closer to mathematical pseudocode, so a pure Julia implementation of a fast solver would potentially be much easier to maintain and contribute to.

Both perspectives are valid and there are going to be people who sympathize with one more than the other. It seems to me that like so many things, especially in open source, those who feel that having a really great pure Julia solver would be worth the effort are precisely the people who can and should make that happen. But of course, they should also be forewarned that it’s a huge undertaking and will be a lot of work—probably a PhD thesis worth, possibly even more to make it really production quality software.

Another question to ask is: Is it possible to make a solver that’s better than existing ones by a sufficient amount to motivate people to switch? In some areas pure Julia implementations of things have been sufficiently much faster or more generic or whatever to really be compelling. Is that the case in the area of LP solvers? What features of Julia make it a great language for making a better solution? Julia is great at code generation. It’s also great at higher order programming—at a level that C++ can only accomplish with templates but which Julia’s JIT can do with arbitrary user-defined functions. Can JIT be leveraged in some other powerful ways that would be too hard to accomplish in other languages?

5 Likes

Frankly some of the commercial solvers are really good, I don’t think it’s realistic to have something competitive with them any time soon unless someone were to actually throw money at this. On the other hand, the open source options currently available are disappointing, both in terms of performance and because it’s really annoying to mess around with their internals (though I imagine it’s probably not exactly easy to beat their performance).

I think if there were a nice Julia solver that were unambitious enough to be achievable with a realistic amount of effort people who are very focused on whatever specific use case it is developed for and who don’t have access to commercial solvers may well switch just from the code and interface being much cleaner (though perhaps that’s a very small audience). Perhaps some of the projects mentioned in this thread already fit that bill. The nice thing about Julia is that those initial efforts would likely be far easier to expand upon than if they had been written in C++, so perhaps in a few years there will be enough of a basis for a larger effort.

1 Like

I think there’s a big hole in the market for a simplex method that works with big floats: I once had a constrained L1 minimisation problem where the matrix was ill-conditioned which broke every solver I tried (CVX was the only one to give an answer), even though the problem it self was well-conditioned. (Think 0.0000001*abs(x). 0 is the unique minimum even with round off even though there are many almost-minimums).

It’s kind of like the LAPACK/generic_matmul! situation. It’s ok that generic_matmul! is order of magnitudes slower than LAPACK because it’s usecase is different.

3 Likes