In the comments to my last post, pozorvlak pointed out that what I was describing was awfully close to the branch and bound algorithm that an integer linear solver would be using anyway and asked if I had tried just adding integer constraints to the LP.

I did actually already know that what I was doing was a form of branch and bound. I’d come up with it independently before realising this, but alas there are no prizes for reinventing a 50 year old algorithm. The question was whether the way I was applying it was a good one or if I’d be better off just using the standard one in the ILP solver. So how do they compare?

Well, uh, there’s a funny story on that one which explains why I hadn’t tried just adding integer constraints to the LP. I actually tried it ages back, found it intolerably slow, and gave up on it as a bad idea.

Basically I only realised how effective the LP was even in large circumstances quite late because I spent the first chunk of working on this problem thinking that the LP couldn’t possibly scale to more than about 10 elements even without the integer constraints.

The reason for this is that I was using PuLP as an interface to the LP solver with its GLPK bindings. This gave me huge performance problems – it’s probably fine if you’re setting up an LP once or twice, maybe, but its performance at building the LP is terrible. It’s several orders of magnitude slower than actually running the LP.

So, I ripped out PuLP and wrote my own interface to lp_solve because it was easy enough to just generate its text format and do that (thanks to Bill Mill for the suggestion). Even that was quite slow at large n though – because the size of the LP is \(O(n^3)\) at large \(n\) you’re generating and parsing several MB to do a decent job of it. So even then the size of the LP I could hand off to it was strictly bounded.

Anyway, to cut a long story short, I’m now using the C API to lp_solve directly so setting up the model has become a lot cheaper (it’s still surprisingly slow for n bigger than around 30, but not intolerably so) and yeah it turns out that lp_solve’s integer constraints scale a lot better than my  custom tailored branch and bound algorithm does.

Moreover, after a while tweaking things and trying various clever stuff and huffily insisting that of course my approach was better and I just needed to optimise it more, I think I’ve convinced myself that the ILP branch and bound algorithm is intrinsically better here, for a reason that is both straightforward and kinda interesting.

The bit about it that is somewhat interesting is that it seems empirically that the LP that results from FAST in some sense “wants” to be 0/1 valued. The vast majority of the time the solution you get out of the relaxed LP is just integer valued automatically and you can call it a day as soon as you get the answer back from it.

The interesting thing is what happens when you don’t. You then have to branch, compute a new LP for each branch, and see what happens there. If each of those sub-problems is integral (or costs at least as much as your current best integer solution) then you’re done, else you’ll have to recurse some more.

And when this happens, my solution branches to n sub problems and the normal ILP branch and bound branches to two. If each is equally likely to produce integer valued sub-problems then the latter is obviously better.

My solution could be better if it chooses sub-problems that are substantially more likely to be integer valued. However this doesn’t seem to be the case – empirically it seems to almost always need to go at least two levels deep if it doesn’t find an integer solution immediately. In fact, it seems plausible to me that the greater freedom in choosing how to branch available to the ILP solution is more likely to let it achieve integer valued sub-problems. I don’t know if it does this, but it could just choose n/2 pairs to branch on, see if any of them produce both sides as integers and then be done if they do. At that point it would have calculated the same number of LPs as my solution but have had a substantially higher chance of not having to recurse further.

So, in conclusion, the better algorithm has won. It was a fun try in which I learned a few things but ultimately as a research avenue not very productive. The only thing from it that might be interesting to pursue further is something I actually didn’t mention in the last post, which is how to use the relaxed LP to get a good heuristic ordering, but I suspect that the quality of that ordering is directly correlated with how easy it is to just get the right answer from the ILP solver so even that is probably a dead end.