Lemke's algorithm and pivoting

[Numsgil]

I've implemented a naive version of Lemke's algorithm for solving LCPs based on this paper: Linear Complementarity and Mathematical (Non-linear) Programming.

To handle degeneracy issues (cases where we have a choice of pivot), which might cause issues with cycling, it suggests "perturbing" each equation in the system by a power of an epsilon. eg: ε, ε^2, ε^3, etc. I need to add this to my solver still, but I'm not quite clear how to do that.

Let's say that it just so happens that in the course of solving the system, you arrive at one equation with a constant term of (ε) * ε^2 and another with a constant term ε^3. They'll be identical and we'll have a choice of pivot again, wouldn't we? Which would break things again since we'd have a choice of pivot.

So I guess you have to keep track of the ε terms symbolically (maybe by adding extra columns to the "dictionary" for each power of ε). But then when it comes time to actually chose a pivot, you still have to actually choose a value for ε so you can evaluate the terms. So how do you chose a value for ε when you evaluate it?

Another source dealing with using Lemke to solve for Nash equilibrium suggested that you could just randomly select a pivot. The idea being that even if you ended up cycling once or twice, eventually the cycle would chose the other pivot and escape. Which is an interesting/clever idea. But it would make debugging sort of a pain since it would break determinacy (maybe there would also be numerical considerations? Like each time through the cycle you lose a bit of precision. So depending on the route taken, you could end up with numerically different answers for the same problem on different run throughs). Are there any other schemes to avoid degeneracy that I could use?

[bone]

Another source dealing with using Lemke to solve for Nash equilibrium suggested that you could just randomly select a pivot.

Well one could store a counter for the system, always initialized to a specific value (like zero). Then every time you need to choose between pivots, increment the counter and take the modulus of the number of choices. For a given system starting at a given initial state, this keeps the determinism for debugging purposes.

EDIT: this is the same as seeding the random number generator with the same seed each time (and making sure nobody else is sharing the generator)

[ngbinh]

I think this paper http://www.roboticsproceedings.org/rss04/p12.html could answer most of your questions. Also, you can find the implementation online ( google for it ).

[Numsgil]

I think this paper http://www.roboticsproceedings.org/rss04/p12.html could answer most of your question

Thanks, I'll take a look.

Also, you can find the implementation online ( google for it ).

I found Hecker's implementation, but it doesn't seem to handle pivoting at all from what I can tell.

I found a MATLAB implementation online, which seems to do what I want and more, but it's using backslash every cycle, which I think does a matrix solve, which is too expensive to be in the center of an iteration. Although as far as I can tell the matrix is only updated by a column every cycle, so I might be able to efficiently handle this by updating the matrix decomposition every frame... Anyway, as far as pivoting in this version, it chooses either the artificial variable every time if it's a candidate, else it chooses randomly from among the candidates.

I implemented both the powers of epsilon (I had to do this symbolically and just say that whatever epsilon actually is, n * ε^2 < ε, for all n) and randomly choosing pivots, and both seem to work as far as a pivoting strategy. Randomly choosing is obviously easier to code and maintain, and requires far less memory, so I'll probably stick with that.

However, I'm running into a failing test and for the life of me I can't figure out what might be wrong in code. It has to do with a degeneracy right at the start, I'm sure. The matrix is a 2x2 identity and the 'q' vector is [-1, -1]. During solving, the code I have right now reaches a state where the dictionary becomes infeasible (specifically entries in the constant vector become negative). I literally rewrote the entire lemke solver from scratch and it still fails on this case, so I'm pretty sure it's something I'm misunderstanding about the algorithm and not a code flub. This also happens exactly the same way for both pivoting strategies. The matlab implementation I linked above can handle this case, though, so I'm rather perplexed since I more or less copied it's pivoting code for randomly choosing from among candidate pivots. Being able to step through the matlab implementation in a debugger would answer a lot of questions, but I don't know how to do that in MATLAB or Octave (free MATLAB clone).

So I dunno, if anyone has a working implementation that they can step through each iteration with the example I gave above and give the basic variables and tableau/dictionary at each step that would be really helpful to compare with. Or even work it through by hand (that's what I'm trying to do now... chart out the different branches the pivot might chose).

[ngbinh]

I meant the implementation of the work on the paper.
You can find it here: http://www.openrtp.jp/openhrp3/en/index.html

Last time I checked it contains C++ implementation of their Lemke. But keep it mind that the method is quite expensive and not suitable for game physics.

[Numsgil]

I meant the implementation of the work on the paper.
You can find it here: http://www.openrtp.jp/openhrp3/en/index.html

Last time I checked it contains C++ implementation of their Lemke.

Hmm, I'll download the packages and see if I can find an actual implementation in there.

But keep it mind that the method is quite expensive and not suitable for game physics.

It's for a simulation of self organizing cells. Something like Tower of Goo but sturdy instead of goo-y and instead of manipulating the cells with a mouse you're programming them like little robots. So it doesn't have to run in real time (though fast is obviously nice), but it does have to be stable, which traditionally is the weak point of "projected gauss seidel" or whatever else you want to call the one-constraint-at-a-time iterative solvers.

My plan is to build up the matrix (which most physics engines don't do obviously) and try out different solvers (Lemke, the more common PGS, some interior point methods, etc.) so I can compare speed and behavior. My thinking right now is that a hybrid solver that starts out with an interior point method can be finished off with a pivoting algorithm like Lemke, and I should be able to get some reasonable speeds (by which I mean single digit frame rates). I might also play with something like PGS, but instead of doing one constraint at a time you solve an entire contact manifold between two bodies, or several bodies, and iterate like that.

And actually after that my goal is to try and build sparse versions of my dense solvers, which should run in more like linear time instead of cubic time and give me something that might work in a real time 60 hz game type situation. But that's a bit in the distant future.

[ngbinh]

It seems like you want to use PATH: http://pages.cs.wisc.edu/~ferris/path/

It should be much more robust and faster than Lemke.

[Numsgil]

It seems like you want to use PATH: http://pages.cs.wisc.edu/~ferris/path/

It should be much more robust and faster than Lemke.

Licensing issues make that unworkable, along with a bunch of other available library options, too. (LGPL and less restrictive would work).