[ODE] LCP solution methods

[Sergiy Migdalskiy]

Hello:

I'm working on an LCP solver, using JM^-1J^T method in Baraff's terminology 
(the LCP w=Mx+q, w,x>=0, wx=0 with matrix M being non singular and strictly 
positive definite). I have some pretty good progress, but I'm still far from 
a high-quality physics package. I wanted to ask for any suggestions how to 
implement iterative O(n) expected running time solver.

I'm looking at running up to 10000 contact problem in real time on the 
modern high-end PC with low memory footprint (don't ask me why :)). I'm 
pretty successful with this goal so far for frictionless balls :) but I hope 
to make it work for everything. The collision detection is not my focus now, 
as with this number of constraints the biggest problem is collision response 
(that is, when I'm using something more complicated than balls)

The LCP matrix is pretty sparse, it's relatively easy to represent it. But 
still, I have about 64000 float entries in the matrix itself, plus control 
structures. In essence, it means I'm limited to 50-100 sweeps through the 
matrix (matrix-vector multiplications). I tried constructing no-fill 
Cholesky factors for the matrix, doesn't seem to work well. Obviously I 
can't use any significant fill-in due to the prohibitive size of the 
problem.

For some simpler configurations (a pile of 1000 balls, tightly packed), the 
matrix is "almost diagonally dominant" in some sense, which means simple 
gauss-seidel iteration works just fine and I only need 15-30 matrix-vector 
multiplications, and it works beautifully fast.

For more complicated situations however, there are big off-diagonal 
elements, which brings gauss-seidel relaxation scheme I'm using to solve the 
LCP to its knees, just as the theory suggests. I tried block Gauss-Seidel, 
but it has essentially the same problems with the non-diagonally-dominant 
matrix.

Can anyone suggest a way to precondition the JM'J'-type matrix so that it 
becomes diagonal dominant? I really want to stick to a very simple iterative 
method, since it works so fast and well for simple configurations. Does 
anyone know any projected iteration methods, besides projected 
SOR/Gauss-Seidel, which work for positive-definite LCPs? I'd kill for a 
projective analog of Conjugate Gradient for LCP! CG solves the analogous 
linear problem (just find x:Mx+q=0) in 30 steps even without 
preconditioning - and my budget is 50-100 steps.

Has anyone tried to use multigrid/multilevel to solve contact forces? The 
literature on multigrid for LCP is sparse to say the least, any references 
would be very much appreciated.

Another way is to use some simplex method. I'm not a specialist in this, but 
to my best knowledge simplex methods (Lemke, criss-cross) don't work on 
sparse matrices, because the matrix gets more and more filled in during 
pivots. And memory is an issue for me. Am I missing something here? Has 
anyone tried simplex methods to solve contact forces in big systems (i.e. 
where it's impossible to fit the whole dense LCP in memory)?

Dantzig method requires keeping and updating a factor of variable size 
(depending on the active basis). I think I heard ODE uses or used this? 
Anyway, there's no way for me to keep a dense factor of a submatrix of a 
10000x10000 JM'J' matrix with a budget of 50-100 matrix-vector 
multiplications. Again, am I missing something? I didn't ever programmed 
updating of a factor by adding/deleting a row/column, but I did program a 
few sparse (incomplete) factors for this problem and even that was a huge 
pain, since the quality of the preconditioner falls dramatically as I try to 
shrink memory usage by raising drop tolerances or implementing some moderate 
fill-in heuristics. In any case, I suspect Dantzig method requires too many 
steps for a big (10000 vars) LCP. Am I wrong?

Another way to solve LCP is to use interior point methods. I experimented 
with them using alternative formulations (QP with such objective as energy 
minimization). It works pretty well, but again, I need at least 10 
outer-loop steps in bad cases, each step consisting of forming a 
preconditioner and solving at least one or two linear problems (e.g. for 
predictor-corrector type algorithm). I can only use cheap preconditioners, 
and it means at least 20-30 steps (very optimistically!) for each linear 
problem (I used MinRes, LSQR, QMR, SQMR, CG, CGS to solve the linear 
problem, all of them with a few different preconditioners - all just to make 
sure they work as expected by theory; I used both the Augmented System and 
Normal Equations approaches; AS worked better for me). So, in total I'm 
looking at least at 400 matrix-vector computations, which is out of budget. 
Besides, in reality, it'll be closer to 1200 due to ill-conditioned systems, 
especially closer to smaller complementarity gaps when the interior point 
converges to the solution. Plus the time to construct each preconditioner..

Oh yeah, to top it all, the nature of the engine is quite dynamic, and usage 
of warm-start (using the solution from the previous frame) is hardly 
possible to the extent that can make a big difference...

I'm starting to run out of sane ideas... Any suggestions would be greatly 
appreciated.

Sincerely,
Sergiy